math_audio_optimisation/
mod.rs

1//! Differential Evolution optimization library.
2//!
3//! This crate provides a Rust implementation of the Differential Evolution (DE)
4//! algorithm, a population-based stochastic optimizer for continuous optimization
5//! problems. It also includes CMA-ES for coupled continuous landscapes and
6//! NSGA-II/NSGA-III for real Pareto multi-objective search.
7//!
8//! # Features
9//!
10//! - Multiple mutation strategies (Best1, Rand1, CurrentToBest1, etc.)
11//! - Binomial and exponential crossover
12//! - Adaptive parameter control (SHADE-style)
13//! - **L-SHADE**: Linear population size reduction for faster convergence
14//! - **External Archive**: Maintains diversity by storing discarded solutions
15//! - **Current-to-pbest/1**: SHADE mutation strategy for balanced exploration
16//! - Parallel population evaluation
17//! - Constraint handling via penalty methods
18//! - Latin Hypercube initialization
19//! - CMA-ES with full covariance adaptation
20//! - NSGA-II/III Pareto-front optimizers
21//!
22//! # Example
23//!
24//! ```rust
25//! use math_audio_optimisation::{differential_evolution, DEConfigBuilder};
26//!
27//! // Minimize the sphere function: f(x) = sum(x_i^2)
28//! let bounds = vec![(-5.0, 5.0), (-5.0, 5.0)];
29//! let config = DEConfigBuilder::new()
30//!     .maxiter(100)
31//!     .seed(42)
32//!     .build()
33//!     .expect("invalid config");
34//!
35//! let result = differential_evolution(
36//!     &|x| x.iter().map(|&xi| xi * xi).sum(),
37//!     &bounds,
38//!     config,
39//! ).expect("optimization should succeed");
40//!
41//! assert!(result.fun < 1e-6);
42//! ```
43//!
44//! # L-SHADE Example
45//!
46//! ```rust
47//! use math_audio_optimisation::{differential_evolution, DEConfigBuilder, Strategy, LShadeConfig};
48//!
49//! let bounds = vec![(-5.0, 5.0), (-5.0, 5.0)];
50//! let lshade_config = LShadeConfig {
51//!     np_init: 18,
52//!     np_final: 4,
53//!     p: 0.11,
54//!     arc_rate: 2.1,
55//!     ..Default::default()
56//! };
57//!
58//! let config = DEConfigBuilder::new()
59//!     .maxiter(100)
60//!     .strategy(Strategy::LShadeBin)
61//!     .lshade(lshade_config)
62//!     .seed(42)
63//!     .build()
64//!     .expect("invalid config");
65//!
66//! let result = differential_evolution(
67//!     &|x| x.iter().map(|&xi| xi * xi).sum(),
68//!     &bounds,
69//!     config,
70//! ).expect("optimization should succeed");
71//!
72//! assert!(result.fun < 1e-6);
73//! ```
74#![doc = include_str!("../README.md")]
75#![doc = include_str!("../REFERENCES.md")]
76#![warn(missing_docs)]
77
78pub mod error;
79pub mod levenberg_marquardt;
80pub use error::{DEError, Result};
81pub use levenberg_marquardt::{
82    LMCallbackAction, LMConfig, LMConfigBuilder, LMError, LMIntermediate, LMReport, LMResult,
83    levenberg_marquardt,
84};
85
86/// Pure-Rust COBYLA local optimizer.
87///
88/// Hand-translated from NLopt 2.7.1's `cobyla.c` (Powell 1994) — see
89/// [`cobyla_native`] for the implementation. The public API lives in
90/// [`cobyla`].
91pub mod cobyla;
92mod cobyla_native;
93pub use cobyla::{CobylaConfig, CobylaConstraint, CobylaReport, CobylaStopTols};
94
95/// Pure-Rust ISRES global constrained optimizer.
96pub mod isres;
97pub use isres::{IsresConfig, IsresConstraint, IsresReport, isres};
98
99/// Pure-Rust CMA-ES global optimizer.
100pub mod cmaes;
101pub use cmaes::{CmaEsConfig, CmaEsIntermediate, CmaEsReport, cma_es};
102
103/// Pure-Rust NSGA-II/III multi-objective optimizers.
104pub mod nsga;
105pub use nsga::{NsgaConfig, NsgaReport, NsgaVariant, ParetoSolution, nsga, nsga2, nsga3};
106
107/// Gaussian-process Bayesian optimizers for expensive objectives.
108pub mod bayesian;
109
110/// Continuous-prior loss integration for spatial / area-based optimization
111/// (expected value, worst-case, CVaR over a [`continuous_area::Prior`]).
112pub mod continuous_area;
113pub use bayesian::{
114    BayesAcquisition, BayesOptCallback, BayesOptConfig, BayesOptIntermediate, BayesOptParetoReport,
115    BayesOptReport, BayesParetoSolution, bayesian_multi_objective, bayesian_optimization,
116};
117pub use continuous_area::{
118    AreaError, AreaScalarisation, Prior, Quadrature, build_quadrature_points, evaluate_area_loss,
119    try_evaluate_area_loss,
120};
121
122use std::fmt;
123use std::str::FromStr;
124use std::sync::Arc;
125use std::sync::RwLock;
126
127use ndarray::{Array1, Array2, Zip};
128use oxiblas_ndarray::blas::matvec;
129use rand::rngs::StdRng;
130use rand::{Rng, SeedableRng};
131use std::time::Instant;
132
133/// Linear penalty stacking utilities for combining multiple constraints.
134pub mod stack_linear_penalty;
135
136/// Integer variable handling for mixed-integer optimization.
137pub mod apply_integrality;
138/// Wrapper Local Search (WLS) application for local refinement.
139pub mod apply_wls;
140
141/// Utilities for selecting distinct random indices from a population.
142pub mod distinct_indices;
143/// External archive for L-SHADE algorithm.
144pub mod external_archive;
145/// Latin Hypercube Sampling initialization strategy.
146pub mod init_latin_hypercube;
147/// Random uniform initialization strategy.
148pub mod init_random;
149/// Sobol quasi-random initialization strategy.
150pub mod init_sobol;
151/// L-SHADE configuration.
152pub mod lshade;
153
154/// Adaptive mutation strategy with dynamic parameter control.
155pub mod mutant_adaptive;
156/// Best/1 mutation strategy: uses best individual plus one difference vector.
157pub mod mutant_best1;
158/// Best/2 mutation strategy: uses best individual plus two difference vectors.
159pub mod mutant_best2;
160/// Current-to-best/1 mutation: blends current with best individual.
161pub mod mutant_current_to_best1;
162/// Current-to-pbest/1 mutation: blends current with p-best individual (SHADE).
163pub mod mutant_current_to_pbest1;
164/// Rand/1 mutation strategy: uses random individual plus one difference vector.
165pub mod mutant_rand1;
166/// Rand/2 mutation strategy: uses random individual plus two difference vectors.
167pub mod mutant_rand2;
168/// Rand-to-best/1 mutation: blends random with best individual.
169pub mod mutant_rand_to_best1;
170
171/// Binomial (uniform) crossover implementation.
172pub mod crossover_binomial;
173/// Exponential crossover implementation.
174pub mod crossover_exponential;
175
176/// Comprehensive tests for DE strategies and features.
177#[cfg(test)]
178mod de_tests;
179/// Main differential evolution algorithm implementation.
180pub mod differential_evolution;
181/// Registry of standard test functions for benchmarking.
182pub mod function_registry;
183/// Internal helper functions for DE implementation.
184pub mod impl_helpers;
185/// Metadata-driven optimization examples and tests.
186pub mod metadata;
187/// Parallel population evaluation support.
188pub mod parallel_eval;
189/// Optimization recording for analysis and debugging.
190pub mod recorder;
191/// Recorded optimization wrapper for testing.
192pub mod run_recorded;
193pub use differential_evolution::differential_evolution;
194pub use external_archive::ExternalArchive;
195pub use lshade::LShadeConfig;
196pub use parallel_eval::ParallelConfig;
197pub use recorder::{OptimizationRecord, OptimizationRecorder};
198pub use run_recorded::run_recorded_differential_evolution;
199
200// Type aliases to reduce complexity
201/// Scalar constraint function type
202pub type ScalarConstraintFn = Arc<dyn Fn(&Array1<f64>) -> f64 + Send + Sync>;
203/// Vector constraint function type
204pub type VectorConstraintFn = Arc<dyn Fn(&Array1<f64>) -> Array1<f64> + Send + Sync>;
205/// Penalty tuple type (function, weight)
206pub type PenaltyTuple = (ScalarConstraintFn, f64);
207/// Callback function type
208pub type CallbackFn = Box<dyn FnMut(&DEIntermediate) -> CallbackAction>;
209
210mod argmin;
211pub(crate) use argmin::argmin;
212
213/// Differential Evolution mutation/crossover strategy.
214///
215/// The strategy name follows the pattern `{mutation}{n}{crossover}` where:
216/// - `mutation`: Base vector selection (Best, Rand, CurrentToBest, RandToBest, Adaptive)
217/// - `n`: Number of difference vectors (1 or 2)
218/// - `crossover`: Crossover type (Bin = binomial, Exp = exponential)
219#[derive(Debug, Clone, Copy)]
220pub enum Strategy {
221    /// Best/1/Bin: Best individual + 1 difference vector, binomial crossover
222    Best1Bin,
223    /// Best/1/Exp: Best individual + 1 difference vector, exponential crossover
224    Best1Exp,
225    /// Rand/1/Bin: Random individual + 1 difference vector, binomial crossover
226    Rand1Bin,
227    /// Rand/1/Exp: Random individual + 1 difference vector, exponential crossover
228    Rand1Exp,
229    /// Rand/2/Bin: Random individual + 2 difference vectors, binomial crossover
230    Rand2Bin,
231    /// Rand/2/Exp: Random individual + 2 difference vectors, exponential crossover
232    Rand2Exp,
233    /// Current-to-best/1/Bin: Blend of current and best + 1 diff, binomial crossover
234    CurrentToBest1Bin,
235    /// Current-to-best/1/Exp: Blend of current and best + 1 diff, exponential crossover
236    CurrentToBest1Exp,
237    /// Best/2/Bin: Best individual + 2 difference vectors, binomial crossover
238    Best2Bin,
239    /// Best/2/Exp: Best individual + 2 difference vectors, exponential crossover
240    Best2Exp,
241    /// Rand-to-best/1/Bin: Blend of random and best + 1 diff, binomial crossover
242    RandToBest1Bin,
243    /// Rand-to-best/1/Exp: Blend of random and best + 1 diff, exponential crossover
244    RandToBest1Exp,
245    /// Adaptive mutation with binomial crossover (SAM approach)
246    AdaptiveBin,
247    /// Adaptive mutation with exponential crossover
248    AdaptiveExp,
249    /// L-SHADE with binomial crossover: current-to-pbest/1 + linear pop reduction + archive
250    LShadeBin,
251    /// L-SHADE with exponential crossover: current-to-pbest/1 + linear pop reduction + archive
252    LShadeExp,
253}
254
255impl FromStr for Strategy {
256    type Err = String;
257    fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
258        let t = s.to_lowercase();
259        match t.as_str() {
260            "best1bin" | "best1" => Ok(Strategy::Best1Bin),
261            "best1exp" => Ok(Strategy::Best1Exp),
262            "rand1bin" | "rand1" => Ok(Strategy::Rand1Bin),
263            "rand1exp" => Ok(Strategy::Rand1Exp),
264            "rand2bin" | "rand2" => Ok(Strategy::Rand2Bin),
265            "rand2exp" => Ok(Strategy::Rand2Exp),
266            "currenttobest1bin" | "current-to-best1bin" | "current_to_best1bin" => {
267                Ok(Strategy::CurrentToBest1Bin)
268            }
269            "currenttobest1exp" | "current-to-best1exp" | "current_to_best1exp" => {
270                Ok(Strategy::CurrentToBest1Exp)
271            }
272            "best2bin" | "best2" => Ok(Strategy::Best2Bin),
273            "best2exp" => Ok(Strategy::Best2Exp),
274            "randtobest1bin" | "rand-to-best1bin" | "rand_to_best1bin" => {
275                Ok(Strategy::RandToBest1Bin)
276            }
277            "randtobest1exp" | "rand-to-best1exp" | "rand_to_best1exp" => {
278                Ok(Strategy::RandToBest1Exp)
279            }
280            "adaptivebin" | "adaptive-bin" | "adaptive_bin" | "adaptive" => {
281                Ok(Strategy::AdaptiveBin)
282            }
283            "adaptiveexp" | "adaptive-exp" | "adaptive_exp" => Ok(Strategy::AdaptiveExp),
284            "lshadebin" | "lshade-bin" | "lshade_bin" | "lshade" => Ok(Strategy::LShadeBin),
285            "lshadeexp" | "lshade-exp" | "lshade_exp" => Ok(Strategy::LShadeExp),
286            _ => Err(format!("unknown strategy: {}", s)),
287        }
288    }
289}
290
291/// Crossover type
292#[derive(Debug, Clone, Copy, Default)]
293pub enum Crossover {
294    /// Binomial (uniform) crossover
295    #[default]
296    Binomial,
297    /// Exponential crossover
298    Exponential,
299}
300
301/// Mutation setting: either a fixed factor, a uniform range (dithering), or adaptive.
302#[derive(Debug, Clone, Copy)]
303pub enum Mutation {
304    /// Fixed mutation factor F in [0, 2).
305    Factor(f64),
306    /// Dithering range [min, max) with 0 <= min < max <= 2.
307    Range {
308        /// Minimum mutation factor.
309        min: f64,
310        /// Maximum mutation factor.
311        max: f64,
312    },
313    /// Adaptive mutation factor using Cauchy distribution.
314    Adaptive {
315        /// Initial mutation factor before adaptation.
316        initial_f: f64,
317    },
318}
319
320impl Default for Mutation {
321    fn default() -> Self {
322        let _ = Mutation::Factor(0.8);
323        Mutation::Range { min: 0.0, max: 2.0 }
324    }
325}
326
327impl Mutation {
328    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
329        match *self {
330            Mutation::Factor(f) => f,
331            Mutation::Range { min, max } => rng.random_range(min..max),
332            Mutation::Adaptive { initial_f } => initial_f, // Will be overridden by adaptive logic
333        }
334    }
335
336    /// Sample from Cauchy distribution for adaptive mutation (F parameter)
337    #[allow(dead_code)]
338    fn sample_cauchy<R: Rng + ?Sized>(&self, f_m: f64, _scale: f64, rng: &mut R) -> f64 {
339        // Simplified version using normal random for now
340        let perturbation = (rng.random::<f64>() - 0.5) * 0.2; // Small perturbation
341        (f_m + perturbation).clamp(0.0, 2.0) // Clamp to valid range
342    }
343}
344
345/// Initialization scheme for the population.
346#[derive(Debug, Clone, Copy, Default)]
347pub enum Init {
348    /// Latin Hypercube Sampling for better space coverage.
349    #[default]
350    LatinHypercube,
351    /// Uniform random initialization.
352    Random,
353}
354
355/// Whether best updates during a generation (we use Deferred only).
356#[derive(Debug, Clone, Copy, Default)]
357pub enum Updating {
358    /// Deferred update: best is updated after all trials are evaluated.
359    #[default]
360    Deferred,
361}
362
363/// Linear penalty specification: lb <= A x <= ub (component-wise).
364#[derive(Debug, Clone)]
365pub struct LinearPenalty {
366    /// Constraint matrix A (m x n).
367    pub a: Array2<f64>,
368    /// Lower bounds vector (m elements).
369    pub lb: Array1<f64>,
370    /// Upper bounds vector (m elements).
371    pub ub: Array1<f64>,
372    /// Penalty weight for constraint violations.
373    pub weight: f64,
374}
375
376/// SciPy-like linear constraint helper: lb <= A x <= ub.
377#[derive(Debug, Clone)]
378pub struct LinearConstraintHelper {
379    /// Constraint matrix A (m x n).
380    pub a: Array2<f64>,
381    /// Lower bounds vector (m elements).
382    pub lb: Array1<f64>,
383    /// Upper bounds vector (m elements).
384    pub ub: Array1<f64>,
385}
386
387impl LinearConstraintHelper {
388    /// Apply helper by merging into DEConfig.linear_penalty (stacking rows if already present)
389    pub fn apply_to(&self, cfg: &mut DEConfig, weight: f64) {
390        use stack_linear_penalty::stack_linear_penalty;
391
392        let new_lp = LinearPenalty {
393            a: self.a.clone(),
394            lb: self.lb.clone(),
395            ub: self.ub.clone(),
396            weight,
397        };
398        match &mut cfg.linear_penalty {
399            Some(existing) => stack_linear_penalty(existing, &new_lp),
400            None => cfg.linear_penalty = Some(new_lp),
401        }
402    }
403}
404
405/// SciPy-like nonlinear constraint helper: vector-valued fun(x) with lb <= fun(x) <= ub.
406#[derive(Clone)]
407pub struct NonlinearConstraintHelper {
408    /// Vector-valued constraint function.
409    pub fun: VectorConstraintFn,
410    /// Lower bounds for each constraint component.
411    pub lb: Array1<f64>,
412    /// Upper bounds for each constraint component.
413    pub ub: Array1<f64>,
414}
415
416impl NonlinearConstraintHelper {
417    /// Apply helper by emitting penalty closures per component.
418    /// lb <= f_i(x) <= ub becomes two inequalities: f_i(x)-ub <= 0 and lb - f_i(x) <= 0.
419    /// If lb==ub, emit an equality penalty for f_i(x)-lb.
420    pub fn apply_to(&self, cfg: &mut DEConfig, weight_ineq: f64, weight_eq: f64) {
421        let f = self.fun.clone();
422        let lb = self.lb.clone();
423        let ub = self.ub.clone();
424        let m = lb.len().min(ub.len());
425        for i in 0..m {
426            let l = lb[i];
427            let u = ub[i];
428            let tol = 1e-12 * (l.abs() + u.abs()).max(1.0);
429            if (u - l).abs() <= tol {
430                let fi = f.clone();
431                cfg.penalty_eq.push((
432                    Arc::new(move |x: &Array1<f64>| {
433                        let y = (fi)(x);
434                        y[i] - l
435                    }),
436                    weight_eq,
437                ));
438            } else {
439                let fi_u = f.clone();
440                cfg.penalty_ineq.push((
441                    Arc::new(move |x: &Array1<f64>| {
442                        let y = (fi_u)(x);
443                        y[i] - u
444                    }),
445                    weight_ineq,
446                ));
447                let fi_l = f.clone();
448                cfg.penalty_ineq.push((
449                    Arc::new(move |x: &Array1<f64>| {
450                        let y = (fi_l)(x);
451                        l - y[i]
452                    }),
453                    weight_ineq,
454                ));
455            }
456        }
457    }
458}
459
460/// Structures for tracking adaptive parameters
461#[derive(Debug, Clone)]
462struct AdaptiveState {
463    /// Current F_m parameter for Cauchy distribution (mutation)
464    f_m: f64,
465    /// Current CR_m parameter for Gaussian distribution (crossover)
466    cr_m: f64,
467    /// Successful F values from this generation
468    successful_f: Vec<f64>,
469    /// Successful CR values from this generation
470    successful_cr: Vec<f64>,
471    /// Current linearly decreasing weight for adaptive mutation
472    current_w: f64,
473}
474
475impl AdaptiveState {
476    fn new(config: &AdaptiveConfig) -> Self {
477        Self {
478            f_m: config.f_m,
479            cr_m: config.cr_m,
480            successful_f: Vec::new(),
481            successful_cr: Vec::new(),
482            current_w: config.w_max, // Start with maximum weight
483        }
484    }
485
486    /// Update adaptive parameters based on successful trials
487    fn update(&mut self, config: &AdaptiveConfig, iter: usize, max_iter: usize) {
488        // Update linearly decreasing weight (Equation 19 from the paper)
489        let iter_ratio = iter as f64 / max_iter as f64;
490        self.current_w = config.w_max - (config.w_max - config.w_min) * iter_ratio;
491
492        // Update F_m using Lehmer mean of successful F values
493        if !self.successful_f.is_empty() {
494            let mean_f = self.compute_lehmer_mean(&self.successful_f);
495            self.f_m = (1.0 - config.w_f) * self.f_m + config.w_f * mean_f;
496            self.f_m = self.f_m.clamp(0.2, 1.2);
497        }
498
499        // Update CR_m using arithmetic mean of successful CR values
500        if !self.successful_cr.is_empty() {
501            let mean_cr = self.compute_arithmetic_mean(&self.successful_cr);
502            self.cr_m = (1.0 - config.w_cr) * self.cr_m + config.w_cr * mean_cr;
503            self.cr_m = self.cr_m.clamp(0.1, 0.9);
504        }
505
506        // Clear successful values for next generation
507        self.successful_f.clear();
508        self.successful_cr.clear();
509    }
510
511    /// Compute Lehmer mean for successful F values (p=2)
512    fn compute_lehmer_mean(&self, values: &[f64]) -> f64 {
513        if values.is_empty() {
514            return 0.5; // Default fallback
515        }
516
517        let sum_sq: f64 = values.iter().map(|&x| x * x).sum();
518        let sum: f64 = values.iter().sum();
519
520        if sum > 0.0 {
521            sum_sq / sum
522        } else {
523            0.5 // Fallback if sum is zero
524        }
525    }
526
527    /// Compute arithmetic mean for successful CR values
528    fn compute_arithmetic_mean(&self, values: &[f64]) -> f64 {
529        if values.is_empty() {
530            return 0.5; // Default fallback
531        }
532        values.iter().sum::<f64>() / values.len() as f64
533    }
534
535    /// Record successful parameter values
536    fn record_success(&mut self, f_val: f64, cr_val: f64) {
537        self.successful_f.push(f_val);
538        self.successful_cr.push(cr_val);
539    }
540
541    /// Sample adaptive F parameter using Cauchy distribution (heavy tails).
542    /// JADE/SHADE literature uses Cauchy for F and Normal for CR.
543    fn sample_f<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
544        let u = rng.random::<f64>().clamp(1e-10, 1.0 - 1e-10);
545        let cauchy = (std::f64::consts::PI * (u - 0.5)).tan();
546        let sample = self.f_m + 0.1 * cauchy;
547
548        sample.clamp(0.1, 1.5)
549    }
550
551    /// Sample adaptive CR parameter using conservative Gaussian distribution
552    fn sample_cr<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
553        let u1: f64 = rng.random::<f64>().max(1e-15);
554        let u2: f64 = rng.random::<f64>();
555
556        let normal = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
557        let sample = self.cr_m + 0.05 * normal;
558
559        sample.clamp(0.1, 0.9)
560    }
561}
562
563/// Adaptive differential evolution configuration
564#[derive(Debug, Clone)]
565pub struct AdaptiveConfig {
566    /// Enable adaptive mutation strategy
567    pub adaptive_mutation: bool,
568    /// Enable Wrapper Local Search (WLS)
569    pub wls_enabled: bool,
570    /// Maximum weight for adaptive mutation (w_max)
571    pub w_max: f64,
572    /// Minimum weight for adaptive mutation (w_min)
573    pub w_min: f64,
574    /// Weight factor for F parameter adaptation (between 0.8 and 1.0)
575    pub w_f: f64,
576    /// Weight factor for CR parameter adaptation (between 0.9 and 1.0)
577    pub w_cr: f64,
578    /// Initial location parameter for Cauchy distribution (F_m)
579    pub f_m: f64,
580    /// Initial location parameter for Gaussian distribution (CR_m)
581    pub cr_m: f64,
582    /// WLS probability (what fraction of population to apply WLS to)
583    pub wls_prob: f64,
584    /// WLS Cauchy scale parameter
585    pub wls_scale: f64,
586}
587
588impl Default for AdaptiveConfig {
589    fn default() -> Self {
590        Self {
591            adaptive_mutation: false,
592            wls_enabled: false,
593            w_max: 0.9,
594            w_min: 0.1,
595            w_f: 0.9,
596            w_cr: 0.9,
597            f_m: 0.5,
598            cr_m: 0.6,
599            wls_prob: 0.1,
600            wls_scale: 0.1,
601        }
602    }
603}
604
605/// Polishing configuration for coordinate-wise local search within bounds.
606#[derive(Debug, Clone)]
607pub struct PolishConfig {
608    /// Whether polishing is enabled.
609    pub enabled: bool,
610    /// Maximum function evaluations for polishing (e.g., 200*n).
611    pub maxeval: usize,
612}
613
614/// Configuration for the Differential Evolution optimizer.
615///
616/// This struct holds all parameters controlling the DE algorithm behavior,
617/// including population size, mutation/crossover settings, constraints, and
618/// convergence criteria.
619pub struct DEConfig {
620    /// Maximum number of generations (iterations).
621    pub maxiter: usize,
622    /// Population size multiplier (total NP = popsize * n_params_free).
623    pub popsize: usize,
624    /// Relative tolerance for convergence (population energy std dev).
625    pub tol: f64,
626    /// Absolute tolerance for convergence on best fitness.
627    pub atol: f64,
628    /// Mutation factor setting.
629    pub mutation: Mutation,
630    /// Crossover probability CR in [0, 1].
631    pub recombination: f64,
632    /// Mutation/crossover strategy.
633    pub strategy: Strategy,
634    /// Crossover type (binomial or exponential).
635    pub crossover: Crossover,
636    /// Population initialization scheme.
637    pub init: Init,
638    /// Update timing (deferred).
639    pub updating: Updating,
640    /// Optional random seed for reproducibility.
641    pub seed: Option<u64>,
642    /// Optional integrality mask; true => variable is integer-constrained.
643    pub integrality: Option<Vec<bool>>,
644    /// Optional initial guess used to replace the best member after init.
645    pub x0: Option<Array1<f64>>,
646    /// Print objective best at each iteration.
647    pub disp: bool,
648    /// Optional per-iteration callback (may stop early).
649    pub callback: Option<CallbackFn>,
650    /// Penalty-based inequality constraints: fc(x) <= 0.
651    pub penalty_ineq: Vec<PenaltyTuple>,
652    /// Penalty-based equality constraints: h(x) = 0.
653    pub penalty_eq: Vec<PenaltyTuple>,
654    /// Optional linear constraints treated by penalty: lb <= A x <= ub.
655    pub linear_penalty: Option<LinearPenalty>,
656    /// Polishing configuration (optional).
657    pub polish: Option<PolishConfig>,
658    /// Adaptive differential evolution configuration.
659    pub adaptive: AdaptiveConfig,
660    /// Parallel evaluation configuration.
661    pub parallel: parallel_eval::ParallelConfig,
662    /// L-SHADE configuration for population reduction and pbest selection.
663    pub lshade: lshade::LShadeConfig,
664}
665
666impl Default for DEConfig {
667    fn default() -> Self {
668        Self {
669            maxiter: 1000,
670            popsize: 15,
671            tol: 1e-2,
672            atol: 0.0,
673            mutation: Mutation::default(),
674            recombination: 0.7,
675            strategy: Strategy::Best1Bin,
676            crossover: Crossover::default(),
677            init: Init::default(),
678            updating: Updating::default(),
679            seed: None,
680            integrality: None,
681            x0: None,
682            disp: false,
683            callback: None,
684            penalty_ineq: Vec::new(),
685            penalty_eq: Vec::new(),
686            linear_penalty: None,
687            polish: None,
688            adaptive: AdaptiveConfig::default(),
689            parallel: parallel_eval::ParallelConfig::default(),
690            lshade: lshade::LShadeConfig::default(),
691        }
692    }
693}
694
695/// Fluent builder for `DEConfig` for ergonomic configuration.
696///
697/// # Example
698///
699/// ```rust
700/// use math_audio_optimisation::{DEConfigBuilder, Strategy, Mutation};
701///
702/// let config = DEConfigBuilder::new()
703///     .maxiter(500)
704///     .popsize(20)
705///     .strategy(Strategy::Best1Bin)
706///     .mutation(Mutation::Factor(0.8))
707///     .recombination(0.9)
708///     .seed(42)
709///     .build();
710/// ```
711pub struct DEConfigBuilder {
712    cfg: DEConfig,
713}
714impl Default for DEConfigBuilder {
715    fn default() -> Self {
716        Self::new()
717    }
718}
719
720impl DEConfigBuilder {
721    /// Creates a new builder with default configuration.
722    pub fn new() -> Self {
723        Self {
724            cfg: DEConfig::default(),
725        }
726    }
727    /// Sets the maximum number of iterations.
728    pub fn maxiter(mut self, v: usize) -> Self {
729        self.cfg.maxiter = v;
730        self
731    }
732    /// Sets the population size multiplier.
733    ///
734    /// # Panics
735    ///
736    /// Panics if `v < 4` since DE requires at least 4 individuals for
737    /// rand/1 and rand/2 mutation strategies.
738    pub fn popsize(mut self, v: usize) -> Self {
739        self.cfg.popsize = v;
740        self
741    }
742    /// Sets the relative convergence tolerance.
743    pub fn tol(mut self, v: f64) -> Self {
744        self.cfg.tol = v;
745        self
746    }
747    /// Sets the absolute convergence tolerance.
748    pub fn atol(mut self, v: f64) -> Self {
749        self.cfg.atol = v;
750        self
751    }
752    /// Sets the mutation factor configuration.
753    pub fn mutation(mut self, v: Mutation) -> Self {
754        self.cfg.mutation = v;
755        self
756    }
757    /// Sets the crossover probability (CR).
758    pub fn recombination(mut self, v: f64) -> Self {
759        self.cfg.recombination = v;
760        self
761    }
762    /// Sets the mutation/crossover strategy.
763    pub fn strategy(mut self, v: Strategy) -> Self {
764        self.cfg.strategy = v;
765        self
766    }
767    /// Sets the crossover type.
768    pub fn crossover(mut self, v: Crossover) -> Self {
769        self.cfg.crossover = v;
770        self
771    }
772    /// Sets the population initialization scheme.
773    pub fn init(mut self, v: Init) -> Self {
774        self.cfg.init = v;
775        self
776    }
777    /// Sets the random seed for reproducibility.
778    pub fn seed(mut self, v: u64) -> Self {
779        self.cfg.seed = Some(v);
780        self
781    }
782    /// Sets the integrality mask for mixed-integer optimization.
783    pub fn integrality(mut self, v: Vec<bool>) -> Self {
784        self.cfg.integrality = Some(v);
785        self
786    }
787    /// Sets an initial guess to seed the population.
788    pub fn x0(mut self, v: Array1<f64>) -> Self {
789        self.cfg.x0 = Some(v);
790        self
791    }
792    /// Enables/disables progress display.
793    pub fn disp(mut self, v: bool) -> Self {
794        self.cfg.disp = v;
795        self
796    }
797    /// Sets a per-iteration callback function.
798    pub fn callback(mut self, cb: Box<dyn FnMut(&DEIntermediate) -> CallbackAction>) -> Self {
799        self.cfg.callback = Some(cb);
800        self
801    }
802    /// Adds an inequality constraint penalty function.
803    pub fn add_penalty_ineq<FN>(mut self, f: FN, w: f64) -> Self
804    where
805        FN: Fn(&Array1<f64>) -> f64 + Send + Sync + 'static,
806    {
807        self.cfg.penalty_ineq.push((Arc::new(f), w));
808        self
809    }
810    /// Adds an equality constraint penalty function.
811    pub fn add_penalty_eq<FN>(mut self, f: FN, w: f64) -> Self
812    where
813        FN: Fn(&Array1<f64>) -> f64 + Send + Sync + 'static,
814    {
815        self.cfg.penalty_eq.push((Arc::new(f), w));
816        self
817    }
818    /// Sets a linear constraint penalty.
819    pub fn linear_penalty(mut self, lp: LinearPenalty) -> Self {
820        self.cfg.linear_penalty = Some(lp);
821        self
822    }
823    /// Sets the polishing configuration.
824    pub fn polish(mut self, pol: PolishConfig) -> Self {
825        self.cfg.polish = Some(pol);
826        self
827    }
828    /// Sets the adaptive DE configuration.
829    pub fn adaptive(mut self, adaptive: AdaptiveConfig) -> Self {
830        self.cfg.adaptive = adaptive;
831        self
832    }
833    /// Enables/disables adaptive mutation.
834    pub fn enable_adaptive_mutation(mut self, enable: bool) -> Self {
835        self.cfg.adaptive.adaptive_mutation = enable;
836        self
837    }
838    /// Enables/disables Wrapper Local Search.
839    pub fn enable_wls(mut self, enable: bool) -> Self {
840        self.cfg.adaptive.wls_enabled = enable;
841        self
842    }
843    /// Sets the adaptive weight bounds.
844    pub fn adaptive_weights(mut self, w_max: f64, w_min: f64) -> Self {
845        self.cfg.adaptive.w_max = w_max;
846        self.cfg.adaptive.w_min = w_min;
847        self
848    }
849    /// Sets the parallel evaluation configuration.
850    pub fn parallel(mut self, parallel: parallel_eval::ParallelConfig) -> Self {
851        self.cfg.parallel = parallel;
852        self
853    }
854    /// Enables/disables parallel evaluation.
855    pub fn enable_parallel(mut self, enable: bool) -> Self {
856        self.cfg.parallel.enabled = enable;
857        self
858    }
859    /// Sets the number of parallel threads.
860    pub fn parallel_threads(mut self, num_threads: usize) -> Self {
861        self.cfg.parallel.num_threads = Some(num_threads);
862        self
863    }
864    /// Sets the L-SHADE configuration.
865    pub fn lshade(mut self, lshade: lshade::LShadeConfig) -> Self {
866        self.cfg.lshade = lshade;
867        self
868    }
869    /// Builds and returns the configuration.
870    ///
871    /// # Errors
872    ///
873    /// Returns `DEError::PopulationTooSmall` if `popsize < 4`.
874    pub fn build(self) -> error::Result<DEConfig> {
875        if self.cfg.popsize < 4 {
876            return Err(DEError::PopulationTooSmall {
877                pop_size: self.cfg.popsize,
878            });
879        }
880        Ok(self.cfg)
881    }
882}
883
884/// Result/report of a DE optimization run.
885///
886/// Contains the optimal solution, convergence status, and statistics.
887#[derive(Clone)]
888pub struct DEReport {
889    /// The optimal solution vector.
890    pub x: Array1<f64>,
891    /// The objective function value at the optimal solution.
892    pub fun: f64,
893    /// Whether the optimization converged successfully.
894    pub success: bool,
895    /// Human-readable status message.
896    pub message: String,
897    /// Number of iterations (generations) performed.
898    pub nit: usize,
899    /// Number of function evaluations performed.
900    pub nfev: usize,
901    /// Final population matrix (NP x n).
902    pub population: Array2<f64>,
903    /// Fitness values for each population member.
904    pub population_energies: Array1<f64>,
905}
906
907impl fmt::Debug for DEReport {
908    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
909        f.debug_struct("DEReport")
910            .field("x", &format!("len={}", self.x.len()))
911            .field("fun", &self.fun)
912            .field("success", &self.success)
913            .field("message", &self.message)
914            .field("nit", &self.nit)
915            .field("nfev", &self.nfev)
916            .field(
917                "population",
918                &format!("{}x{}", self.population.nrows(), self.population.ncols()),
919            )
920            .field(
921                "population_energies",
922                &format!("len={}", self.population_energies.len()),
923            )
924            .finish()
925    }
926}
927
928/// Information passed to callback after each generation.
929pub struct DEIntermediate {
930    /// Current best solution vector.
931    pub x: Array1<f64>,
932    /// Current best objective value.
933    pub fun: f64,
934    /// Convergence measure (population fitness std dev).
935    pub convergence: f64,
936    /// Current iteration number.
937    pub iter: usize,
938}
939
940/// Action returned by callback to control optimization flow.
941pub enum CallbackAction {
942    /// Continue optimization.
943    Continue,
944    /// Stop optimization early.
945    Stop,
946}
947
948/// Differential Evolution optimizer.
949///
950/// A population-based stochastic optimizer for continuous functions.
951/// Use [`DifferentialEvolution::new`] to create an instance, configure
952/// with [`config_mut`](Self::config_mut), then call [`solve`](Self::solve).
953pub struct DifferentialEvolution<'a, F>
954where
955    F: Fn(&Array1<f64>) -> f64 + Sync,
956{
957    func: &'a F,
958    lower: Array1<f64>,
959    upper: Array1<f64>,
960    config: DEConfig,
961}
962
963impl<'a, F> DifferentialEvolution<'a, F>
964where
965    F: Fn(&Array1<f64>) -> f64 + Sync,
966{
967    /// Creates a new DE optimizer with objective `func` and bounds [lower, upper].
968    ///
969    /// # Errors
970    ///
971    /// Returns `DEError::BoundsMismatch` if `lower` and `upper` have different lengths.
972    /// Returns `DEError::InvalidBounds` if any lower bound exceeds its corresponding upper bound.
973    pub fn new(func: &'a F, lower: Array1<f64>, upper: Array1<f64>) -> Result<Self> {
974        if lower.len() != upper.len() {
975            return Err(DEError::BoundsMismatch {
976                lower_len: lower.len(),
977                upper_len: upper.len(),
978            });
979        }
980
981        // Validate that lower <= upper for all dimensions
982        for i in 0..lower.len() {
983            if lower[i] > upper[i] {
984                return Err(DEError::InvalidBounds {
985                    index: i,
986                    lower: lower[i],
987                    upper: upper[i],
988                });
989            }
990        }
991
992        Ok(Self {
993            func,
994            lower,
995            upper,
996            config: DEConfig::default(),
997        })
998    }
999
1000    /// Mutable access to configuration
1001    pub fn config_mut(&mut self) -> &mut DEConfig {
1002        &mut self.config
1003    }
1004
1005    /// Run the optimization and return a report
1006    pub fn solve(&mut self) -> DEReport {
1007        use apply_integrality::apply_integrality;
1008        use apply_wls::apply_wls;
1009        use crossover_binomial::binomial_crossover;
1010        use crossover_exponential::exponential_crossover;
1011        use init_latin_hypercube::init_latin_hypercube;
1012        use init_random::init_random;
1013        use mutant_adaptive::mutant_adaptive;
1014        use mutant_best1::mutant_best1;
1015        use mutant_best2::mutant_best2;
1016        use mutant_current_to_best1::mutant_current_to_best1;
1017        use mutant_current_to_pbest1;
1018        use mutant_rand_to_best1::mutant_rand_to_best1;
1019        use mutant_rand1::mutant_rand1;
1020        use mutant_rand2::mutant_rand2;
1021        use parallel_eval::evaluate_trials_parallel;
1022        use std::sync::Arc;
1023
1024        let n = self.lower.len();
1025
1026        // Identify fixed (equal-bounds) and free variables
1027        let mut is_free: Vec<bool> = Vec::with_capacity(n);
1028        for i in 0..n {
1029            is_free.push((self.upper[i] - self.lower[i]).abs() > 0.0);
1030        }
1031        let n_free = is_free.iter().filter(|&&b| b).count();
1032        let _n_equal = n - n_free;
1033        if n_free == 0 {
1034            // All fixed; just evaluate x = lower
1035            let x_fixed = self.lower.clone();
1036            let mut x_eval = x_fixed.clone();
1037            if let Some(mask) = &self.config.integrality {
1038                apply_integrality(&mut x_eval, mask, &self.lower, &self.upper);
1039            }
1040            let f = (self.func)(&x_eval);
1041            return DEReport {
1042                x: x_eval,
1043                fun: f,
1044                success: true,
1045                message: "All variables fixed by bounds".into(),
1046                nit: 0,
1047                nfev: 1,
1048                population: Array2::zeros((1, n)),
1049                population_energies: Array1::from(vec![f]),
1050            };
1051        }
1052
1053        let is_lshade = matches!(
1054            self.config.strategy,
1055            Strategy::LShadeBin | Strategy::LShadeExp
1056        );
1057        let mut npop = if is_lshade {
1058            self.config.lshade.initial_population_size(n_free)
1059        } else {
1060            self.config.popsize * n_free
1061        };
1062        let max_nfev = if is_lshade {
1063            self.config.maxiter * npop
1064        } else {
1065            0
1066        };
1067        let _bounds_span = &self.upper - &self.lower;
1068
1069        if self.config.disp {
1070            eprintln!(
1071                "DE Init: {} dimensions ({} free), population={}, maxiter={}",
1072                n, n_free, npop, self.config.maxiter
1073            );
1074            eprintln!(
1075                "  Strategy: {:?}, Mutation: {:?}, Crossover: CR={:.3}",
1076                self.config.strategy, self.config.mutation, self.config.recombination
1077            );
1078            eprintln!(
1079                "  Tolerances: tol={:.2e}, atol={:.2e}",
1080                self.config.tol, self.config.atol
1081            );
1082        }
1083
1084        // Timing toggle via env var
1085        let timing_enabled = std::env::var("AUTOEQ_DE_TIMING")
1086            .map(|v| v != "0")
1087            .unwrap_or(false);
1088
1089        // Configure global rayon thread pool once if requested
1090        if let Some(n) = self.config.parallel.num_threads {
1091            // Ignore error if global pool already set
1092            let _ = rayon::ThreadPoolBuilder::new()
1093                .num_threads(n)
1094                .build_global();
1095        }
1096
1097        // RNG
1098        let mut rng: StdRng = match self.config.seed {
1099            Some(s) => StdRng::seed_from_u64(s),
1100            None => {
1101                let mut thread_rng = rand::rng();
1102                StdRng::from_rng(&mut thread_rng)
1103            }
1104        };
1105
1106        // Initialize population in [lower, upper]
1107        let mut pop = match self.config.init {
1108            Init::LatinHypercube => {
1109                if self.config.disp {
1110                    eprintln!("  Using Latin Hypercube initialization");
1111                }
1112                init_latin_hypercube(n, npop, &self.lower, &self.upper, &is_free, &mut rng)
1113            }
1114            Init::Random => {
1115                if self.config.disp {
1116                    eprintln!("  Using Random initialization");
1117                }
1118                init_random(n, npop, &self.lower, &self.upper, &is_free, &mut rng)
1119            }
1120        };
1121
1122        // Evaluate energies (objective + penalties)
1123        let mut nfev: usize = 0;
1124        if self.config.disp {
1125            eprintln!("  Evaluating initial population of {} individuals...", npop);
1126        }
1127
1128        // Prepare population for evaluation (apply integrality constraints)
1129        let mut eval_pop = pop.clone();
1130        let t_integrality0 = Instant::now();
1131        if let Some(mask) = &self.config.integrality {
1132            for i in 0..npop {
1133                let mut row = eval_pop.row_mut(i);
1134                let mut x_eval = row.to_owned();
1135                apply_integrality(&mut x_eval, mask, &self.lower, &self.upper);
1136                row.assign(&x_eval);
1137            }
1138        }
1139        let t_integrality = t_integrality0.elapsed();
1140
1141        // Build thread-safe energy function that includes penalties
1142        let func_ref = self.func;
1143        let penalty_ineq_vec: Vec<PenaltyTuple> = self
1144            .config
1145            .penalty_ineq
1146            .iter()
1147            .map(|(f, w)| (f.clone(), *w))
1148            .collect();
1149        let penalty_eq_vec: Vec<PenaltyTuple> = self
1150            .config
1151            .penalty_eq
1152            .iter()
1153            .map(|(f, w)| (f.clone(), *w))
1154            .collect();
1155        let linear_penalty = self.config.linear_penalty.clone();
1156
1157        let energy_fn = Arc::new(move |x: &Array1<f64>| -> f64 {
1158            let base = (func_ref)(x);
1159            let mut p = 0.0;
1160            for (f, w) in &penalty_ineq_vec {
1161                let v = f(x);
1162                let viol = v.max(0.0);
1163                p += w * viol * viol;
1164            }
1165            for (h, w) in &penalty_eq_vec {
1166                let v = h(x);
1167                p += w * v * v;
1168            }
1169            if let Some(ref lp) = linear_penalty {
1170                let ax = matvec(&lp.a, &x.to_owned());
1171                Zip::from(&ax)
1172                    .and(&lp.lb)
1173                    .and(&lp.ub)
1174                    .for_each(|&v, &lo, &hi| {
1175                        if v < lo {
1176                            let d = lo - v;
1177                            p += lp.weight * d * d;
1178                        } else if v > hi {
1179                            let d = v - hi;
1180                            p += lp.weight * d * d;
1181                        }
1182                    });
1183            }
1184            let energy = base + p;
1185            if energy.is_finite() {
1186                energy
1187            } else {
1188                f64::INFINITY
1189            }
1190        });
1191
1192        let t_eval0 = Instant::now();
1193        let mut energies = parallel_eval::evaluate_population_parallel(
1194            &eval_pop,
1195            energy_fn.clone(),
1196            &self.config.parallel,
1197        );
1198        let t_eval_init = t_eval0.elapsed();
1199        nfev += npop;
1200        if timing_enabled {
1201            eprintln!(
1202                "TIMING init: integrality={:.3} ms, eval={:.3} ms",
1203                t_integrality.as_secs_f64() * 1e3,
1204                t_eval_init.as_secs_f64() * 1e3
1205            );
1206        }
1207
1208        // Report initial population statistics
1209        let pop_mean = energies.mean().unwrap_or(0.0);
1210        let pop_std = energies.std(0.0);
1211        if self.config.disp {
1212            eprintln!(
1213                "  Initial population: mean={:.6e}, std={:.6e}",
1214                pop_mean, pop_std
1215            );
1216        }
1217
1218        // If x0 provided, override the best member
1219        if let Some(x0) = &self.config.x0 {
1220            let mut x0c = x0.clone();
1221            // Clip to bounds using ndarray
1222            for i in 0..x0c.len() {
1223                x0c[i] = x0c[i].clamp(self.lower[i], self.upper[i]);
1224            }
1225            if let Some(mask) = &self.config.integrality {
1226                apply_integrality(&mut x0c, mask, &self.lower, &self.upper);
1227            }
1228            let f0 = self.energy(&x0c);
1229            nfev += 1;
1230            // find current best
1231            let (best_idx, _best_f) = argmin(&energies);
1232            pop.row_mut(best_idx).assign(&x0c.view());
1233            energies[best_idx] = f0;
1234        }
1235
1236        let (mut best_idx, mut best_f) = argmin(&energies);
1237        let mut best_x = pop.row(best_idx).to_owned();
1238
1239        if self.config.disp {
1240            eprintln!(
1241                "  Initial best: fitness={:.6e} at index {}",
1242                best_f, best_idx
1243            );
1244            let param_summary: Vec<String> = (0..best_x.len() / 3)
1245                .map(|i| {
1246                    let freq = 10f64.powf(best_x[i * 3]);
1247                    let q = best_x[i * 3 + 1];
1248                    let gain = best_x[i * 3 + 2];
1249                    format!("f{:.0}Hz/Q{:.2}/G{:.2}dB", freq, q, gain)
1250                })
1251                .collect();
1252            eprintln!("  Initial best params: [{}]", param_summary.join(", "));
1253        }
1254
1255        if self.config.disp {
1256            eprintln!("DE iter {:4}  best_f={:.6e}", 0, best_f);
1257        }
1258
1259        // Initialize adaptive state if adaptive strategies are enabled
1260        let mut adaptive_state = if matches!(
1261            self.config.strategy,
1262            Strategy::AdaptiveBin | Strategy::AdaptiveExp
1263        ) || self.config.adaptive.adaptive_mutation
1264        {
1265            Some(AdaptiveState::new(&self.config.adaptive))
1266        } else {
1267            None
1268        };
1269
1270        // Initialize external archive for L-SHADE strategies
1271        let external_archive: Option<Arc<RwLock<external_archive::ExternalArchive>>> = if matches!(
1272            self.config.strategy,
1273            Strategy::LShadeBin | Strategy::LShadeExp
1274        ) {
1275            Some(Arc::new(RwLock::new(
1276                external_archive::ExternalArchive::with_population_size(
1277                    npop,
1278                    self.config.lshade.arc_rate,
1279                ),
1280            )))
1281        } else {
1282            None
1283        };
1284
1285        // Main loop
1286        let mut success = false;
1287        let mut message = String::new();
1288        let mut nit = 0;
1289        let mut accepted_trials;
1290        let mut improvement_count;
1291
1292        let mut t_build_tot = std::time::Duration::ZERO;
1293        let mut t_eval_tot = std::time::Duration::ZERO;
1294        let mut t_select_tot = std::time::Duration::ZERO;
1295        let mut t_iter_tot = std::time::Duration::ZERO;
1296
1297        for iter in 1..=self.config.maxiter {
1298            nit = iter;
1299            accepted_trials = 0;
1300            improvement_count = 0;
1301
1302            let iter_start = Instant::now();
1303
1304            // Pre-sort indices for adaptive/L-SHADE strategies to avoid re-sorting in the loop
1305            let sorted_indices = if matches!(
1306                self.config.strategy,
1307                Strategy::AdaptiveBin
1308                    | Strategy::AdaptiveExp
1309                    | Strategy::LShadeBin
1310                    | Strategy::LShadeExp
1311            ) {
1312                let mut indices: Vec<usize> = (0..npop).collect();
1313                indices.sort_by(|&a, &b| {
1314                    energies[a]
1315                        .partial_cmp(&energies[b])
1316                        .unwrap_or(std::cmp::Ordering::Equal)
1317                });
1318                indices
1319            } else {
1320                Vec::new() // Not needed for other strategies
1321            };
1322
1323            // Generate all trials first, then evaluate in parallel
1324            let t_build0 = Instant::now();
1325
1326            // Parallelize trial generation using rayon
1327            use rayon::prelude::*;
1328            let (trials, trial_params): (Vec<_>, Vec<_>) = (0..npop)
1329                .into_par_iter()
1330                .map(|i| {
1331                    // Create thread-local RNG from base seed + iteration + individual index
1332                    let mut local_rng: StdRng = if let Some(base_seed) = self.config.seed {
1333                        StdRng::seed_from_u64(
1334                            base_seed
1335                                .wrapping_add((iter as u64) << 32)
1336                                .wrapping_add(i as u64),
1337                        )
1338                    } else {
1339                        // Use thread_rng for unseeded runs
1340                        let mut thread_rng = rand::rng();
1341                        StdRng::from_rng(&mut thread_rng)
1342                    };
1343
1344                    // Sample mutation factor and crossover rate (adaptive or fixed)
1345                    let (f, cr) = if let Some(ref adaptive) = adaptive_state {
1346                        // Use adaptive parameter sampling
1347                        let adaptive_f = adaptive.sample_f(&mut local_rng);
1348                        let adaptive_cr = adaptive.sample_cr(&mut local_rng);
1349                        (adaptive_f, adaptive_cr)
1350                    } else {
1351                        // Use fixed or dithered parameters
1352                        (
1353                            self.config.mutation.sample(&mut local_rng),
1354                            self.config.recombination,
1355                        )
1356                    };
1357
1358                    // Generate mutant and apply crossover based on strategy
1359                    let (mutant, cross) = match self.config.strategy {
1360                        Strategy::Best1Bin => (
1361                            mutant_best1(i, &pop, best_idx, f, &mut local_rng),
1362                            Crossover::Binomial,
1363                        ),
1364                        Strategy::Best1Exp => (
1365                            mutant_best1(i, &pop, best_idx, f, &mut local_rng),
1366                            Crossover::Exponential,
1367                        ),
1368                        Strategy::Rand1Bin => (
1369                            mutant_rand1(i, &pop, f, &mut local_rng),
1370                            Crossover::Binomial,
1371                        ),
1372                        Strategy::Rand1Exp => (
1373                            mutant_rand1(i, &pop, f, &mut local_rng),
1374                            Crossover::Exponential,
1375                        ),
1376                        Strategy::Rand2Bin => (
1377                            mutant_rand2(i, &pop, f, &mut local_rng),
1378                            Crossover::Binomial,
1379                        ),
1380                        Strategy::Rand2Exp => (
1381                            mutant_rand2(i, &pop, f, &mut local_rng),
1382                            Crossover::Exponential,
1383                        ),
1384                        Strategy::CurrentToBest1Bin => (
1385                            mutant_current_to_best1(i, &pop, best_idx, f, &mut local_rng),
1386                            Crossover::Binomial,
1387                        ),
1388                        Strategy::CurrentToBest1Exp => (
1389                            mutant_current_to_best1(i, &pop, best_idx, f, &mut local_rng),
1390                            Crossover::Exponential,
1391                        ),
1392                        Strategy::Best2Bin => (
1393                            mutant_best2(i, &pop, best_idx, f, &mut local_rng),
1394                            Crossover::Binomial,
1395                        ),
1396                        Strategy::Best2Exp => (
1397                            mutant_best2(i, &pop, best_idx, f, &mut local_rng),
1398                            Crossover::Exponential,
1399                        ),
1400                        Strategy::RandToBest1Bin => (
1401                            mutant_rand_to_best1(i, &pop, best_idx, f, &mut local_rng),
1402                            Crossover::Binomial,
1403                        ),
1404                        Strategy::RandToBest1Exp => (
1405                            mutant_rand_to_best1(i, &pop, best_idx, f, &mut local_rng),
1406                            Crossover::Exponential,
1407                        ),
1408                        Strategy::AdaptiveBin => {
1409                            if let Some(ref adaptive) = adaptive_state {
1410                                (
1411                                    mutant_adaptive(
1412                                        i,
1413                                        &pop,
1414                                        &sorted_indices,
1415                                        adaptive.current_w,
1416                                        f,
1417                                        &mut local_rng,
1418                                    ),
1419                                    Crossover::Binomial,
1420                                )
1421                            } else {
1422                                // Fallback to rand1 if adaptive state not available
1423                                (
1424                                    mutant_rand1(i, &pop, f, &mut local_rng),
1425                                    Crossover::Binomial,
1426                                )
1427                            }
1428                        }
1429                        Strategy::AdaptiveExp => {
1430                            if let Some(ref adaptive) = adaptive_state {
1431                                (
1432                                    mutant_adaptive(
1433                                        i,
1434                                        &pop,
1435                                        &sorted_indices,
1436                                        adaptive.current_w,
1437                                        f,
1438                                        &mut local_rng,
1439                                    ),
1440                                    Crossover::Exponential,
1441                                )
1442                            } else {
1443                                // Fallback to rand1 if adaptive state not available
1444                                (
1445                                    mutant_rand1(i, &pop, f, &mut local_rng),
1446                                    Crossover::Exponential,
1447                                )
1448                            }
1449                        }
1450                        Strategy::LShadeBin => {
1451                            let pbest_size = mutant_current_to_pbest1::compute_pbest_size(
1452                                self.config.lshade.p,
1453                                pop.nrows(),
1454                            );
1455                            let archive_ref = external_archive.as_ref().and_then(|a| a.read().ok());
1456                            (
1457                                mutant_current_to_pbest1::mutant_current_to_pbest1(
1458                                    i,
1459                                    &pop,
1460                                    &sorted_indices,
1461                                    pbest_size,
1462                                    archive_ref.as_deref(),
1463                                    f,
1464                                    &mut local_rng,
1465                                ),
1466                                Crossover::Binomial,
1467                            )
1468                        }
1469                        Strategy::LShadeExp => {
1470                            let pbest_size = mutant_current_to_pbest1::compute_pbest_size(
1471                                self.config.lshade.p,
1472                                pop.nrows(),
1473                            );
1474                            let archive_ref = external_archive.as_ref().and_then(|a| a.read().ok());
1475                            (
1476                                mutant_current_to_pbest1::mutant_current_to_pbest1(
1477                                    i,
1478                                    &pop,
1479                                    &sorted_indices,
1480                                    pbest_size,
1481                                    archive_ref.as_deref(),
1482                                    f,
1483                                    &mut local_rng,
1484                                ),
1485                                Crossover::Exponential,
1486                            )
1487                        }
1488                    };
1489
1490                    // If strategy didn't dictate crossover, fallback to config
1491                    let crossover = cross;
1492                    let trial = match crossover {
1493                        Crossover::Binomial => {
1494                            binomial_crossover(&pop.row(i).to_owned(), &mutant, cr, &mut local_rng)
1495                        }
1496                        Crossover::Exponential => exponential_crossover(
1497                            &pop.row(i).to_owned(),
1498                            &mutant,
1499                            cr,
1500                            &mut local_rng,
1501                        ),
1502                    };
1503
1504                    // Apply WLS if enabled
1505                    let wls_trial = if self.config.adaptive.wls_enabled
1506                        && local_rng.random::<f64>() < self.config.adaptive.wls_prob
1507                    {
1508                        apply_wls(
1509                            &trial,
1510                            &self.lower,
1511                            &self.upper,
1512                            self.config.adaptive.wls_scale,
1513                            &mut local_rng,
1514                        )
1515                    } else {
1516                        trial.clone()
1517                    };
1518
1519                    // Clip to bounds using vectorized operation
1520                    let mut trial_clipped = wls_trial;
1521                    Zip::from(&mut trial_clipped)
1522                        .and(&self.lower)
1523                        .and(&self.upper)
1524                        .for_each(|x, lo, hi| *x = x.clamp(*lo, *hi));
1525
1526                    // Apply integrality if provided
1527                    if let Some(mask) = &self.config.integrality {
1528                        apply_integrality(&mut trial_clipped, mask, &self.lower, &self.upper);
1529                    }
1530
1531                    // Return trial and parameters
1532                    (trial_clipped, (f, cr))
1533                })
1534                .unzip();
1535
1536            let t_build = t_build0.elapsed();
1537            let t_eval0 = Instant::now();
1538            let trial_energies =
1539                evaluate_trials_parallel(&trials, energy_fn.clone(), &self.config.parallel);
1540            let t_eval = t_eval0.elapsed();
1541            nfev += npop;
1542
1543            let t_select0 = Instant::now();
1544            // Selection phase: update population based on trial results
1545            for (i, (trial, trial_energy)) in
1546                trials.into_iter().zip(trial_energies.iter()).enumerate()
1547            {
1548                let (f, cr) = trial_params[i];
1549
1550                // Selection: replace if better
1551                if *trial_energy <= energies[i] {
1552                    // For L-SHADE: add replaced solution to archive
1553                    if let Some(ref archive) = external_archive
1554                        && *trial_energy < energies[i]
1555                        && let Ok(mut arch) = archive.write()
1556                    {
1557                        arch.add(pop.row(i).to_owned());
1558                    }
1559                    pop.row_mut(i).assign(&trial.view());
1560                    energies[i] = *trial_energy;
1561                    accepted_trials += 1;
1562
1563                    // Update adaptive parameters if improvement
1564                    if let Some(ref mut adaptive) = adaptive_state {
1565                        adaptive.record_success(f, cr);
1566                    }
1567
1568                    // Track if this is an improvement over the current best
1569                    if *trial_energy < best_f {
1570                        improvement_count += 1;
1571                    }
1572                }
1573            }
1574            let t_select = t_select0.elapsed();
1575
1576            // L-SHADE: linear population size reduction
1577            if is_lshade {
1578                let current_npop = self
1579                    .config
1580                    .lshade
1581                    .current_population_size(n_free, nfev, max_nfev);
1582                if current_npop < pop.nrows() {
1583                    let mut indices: Vec<usize> = (0..pop.nrows()).collect();
1584                    indices.sort_by(|&a, &b| {
1585                        energies[a]
1586                            .partial_cmp(&energies[b])
1587                            .unwrap_or(std::cmp::Ordering::Equal)
1588                    });
1589                    let keep: Vec<usize> = indices.into_iter().take(current_npop).collect();
1590
1591                    let mut new_pop = Array2::zeros((current_npop, n));
1592                    let mut new_energies = Array1::zeros(current_npop);
1593                    for (new_i, &old_i) in keep.iter().enumerate() {
1594                        new_pop.row_mut(new_i).assign(&pop.row(old_i));
1595                        new_energies[new_i] = energies[old_i];
1596                    }
1597                    pop = new_pop;
1598                    energies = new_energies;
1599                    npop = current_npop;
1600
1601                    if let Some(ref archive) = external_archive
1602                        && let Ok(mut arch) = archive.write()
1603                    {
1604                        arch.resize(self.config.lshade.current_archive_size(current_npop));
1605                    }
1606                }
1607            }
1608
1609            t_build_tot += t_build;
1610            t_eval_tot += t_eval;
1611            t_select_tot += t_select;
1612            let iter_dur = iter_start.elapsed();
1613            t_iter_tot += iter_dur;
1614
1615            if timing_enabled && (iter <= 5 || iter % 10 == 0) {
1616                eprintln!(
1617                    "TIMING iter {:4}: build={:.3} ms, eval={:.3} ms, select={:.3} ms, total={:.3} ms",
1618                    iter,
1619                    t_build.as_secs_f64() * 1e3,
1620                    t_eval.as_secs_f64() * 1e3,
1621                    t_select.as_secs_f64() * 1e3,
1622                    iter_dur.as_secs_f64() * 1e3,
1623                );
1624            }
1625
1626            // Update adaptive parameters after each generation
1627            if let Some(ref mut adaptive) = adaptive_state {
1628                adaptive.update(&self.config.adaptive, iter, self.config.maxiter);
1629            }
1630
1631            // Update best solution after generation
1632            let (new_best_idx, new_best_f) = argmin(&energies);
1633            if new_best_f < best_f {
1634                best_idx = new_best_idx;
1635                best_f = new_best_f;
1636                best_x = pop.row(best_idx).to_owned();
1637            }
1638
1639            // Convergence check
1640            let pop_mean = energies.mean().unwrap_or(0.0);
1641            let pop_std = energies.std(0.0);
1642            let convergence_threshold = self.config.atol + self.config.tol * pop_mean.abs();
1643
1644            if self.config.disp {
1645                eprintln!(
1646                    "DE iter {:4}  best_f={:.6e}  std={:.3e}  accepted={}/{}, improved={}",
1647                    iter, best_f, pop_std, accepted_trials, npop, improvement_count
1648                );
1649            }
1650
1651            // Callback
1652            if let Some(ref mut cb) = self.config.callback {
1653                let intermediate = DEIntermediate {
1654                    x: best_x.clone(),
1655                    fun: best_f,
1656                    convergence: pop_std,
1657                    iter,
1658                };
1659                match cb(&intermediate) {
1660                    CallbackAction::Stop => {
1661                        success = true;
1662                        message = "Optimization stopped by callback".to_string();
1663                        break;
1664                    }
1665                    CallbackAction::Continue => {}
1666                }
1667            }
1668
1669            if pop_std <= convergence_threshold {
1670                success = true;
1671                message = format!(
1672                    "Converged: std(pop_f)={:.3e} <= threshold={:.3e}",
1673                    pop_std, convergence_threshold
1674                );
1675                break;
1676            }
1677        }
1678
1679        if !success {
1680            message = format!("Maximum iterations reached: {}", self.config.maxiter);
1681        }
1682
1683        if self.config.disp {
1684            eprintln!("DE finished: {}", message);
1685        }
1686
1687        // Polish if configured
1688        let (final_x, final_f, polish_nfev) = if let Some(ref polish_cfg) = self.config.polish {
1689            if polish_cfg.enabled {
1690                self.polish(&best_x)
1691            } else {
1692                (best_x.clone(), best_f, 0)
1693            }
1694        } else {
1695            (best_x.clone(), best_f, 0)
1696        };
1697
1698        if timing_enabled {
1699            eprintln!(
1700                "TIMING total: build={:.3} s, eval={:.3} s, select={:.3} s, iter_total={:.3} s",
1701                t_build_tot.as_secs_f64(),
1702                t_eval_tot.as_secs_f64(),
1703                t_select_tot.as_secs_f64(),
1704                t_iter_tot.as_secs_f64()
1705            );
1706        }
1707
1708        self.finish_report(
1709            pop,
1710            energies,
1711            final_x,
1712            final_f,
1713            success,
1714            message,
1715            nit,
1716            nfev + polish_nfev,
1717        )
1718    }
1719}
1720
1721#[cfg(test)]
1722mod strategy_tests {
1723    use super::*;
1724    use rand::SeedableRng;
1725
1726    #[test]
1727    fn test_parse_strategy_variants() {
1728        assert!(matches!(
1729            "best1exp".parse::<Strategy>().unwrap(),
1730            Strategy::Best1Exp
1731        ));
1732        assert!(matches!(
1733            "rand1bin".parse::<Strategy>().unwrap(),
1734            Strategy::Rand1Bin
1735        ));
1736        assert!(matches!(
1737            "randtobest1exp".parse::<Strategy>().unwrap(),
1738            Strategy::RandToBest1Exp
1739        ));
1740    }
1741
1742    #[test]
1743    fn sample_f_uses_cauchy_heavy_tails() {
1744        // JADE/SHADE literature requires Cauchy sampling for F so that
1745        // occasional large mutation factors can escape local minima.
1746        // A Normal(0.5, 0.05) has negligible probability outside [0.1, 0.9]
1747        // (≈8σ). A Cauchy(0.5, 0.1) has ~15% mass outside that range.
1748        let config = AdaptiveConfig::default();
1749        let state = AdaptiveState::new(&config);
1750        let mut rng = StdRng::seed_from_u64(42);
1751        let mut extreme = 0;
1752        for _ in 0..200 {
1753            let f = state.sample_f(&mut rng);
1754            if !(0.1..=0.9).contains(&f) {
1755                extreme += 1;
1756            }
1757        }
1758        assert!(
1759            extreme > 5,
1760            "Cauchy sampling should produce heavy tails, got {} extreme values",
1761            extreme
1762        );
1763    }
1764}
math_audio_optimisation/mod.rs

math_audio_optimisation/
mod.rs
