Function differential_evolution 

pub fn differential_evolution<F>(
    func: &F,
    bounds: &[(f64, f64)],
    config: DEConfig,
) -> Result<DEReport>
where
    F: Fn(&Array1<f64>) -> f64 + Sync,

Runs Differential Evolution optimization on a function.

This is a convenience function that mirrors SciPy’s differential_evolution API. It creates a DE optimizer with the given bounds and configuration, then runs the optimization to find the global minimum.

Arguments

• func - The objective function to minimize, mapping &Array1<f64> to f64
• bounds - Vector of (lower, upper) bound pairs for each dimension
• config - DE configuration (use DEConfigBuilder to construct)

Returns

Returns Ok(DEReport) containing the optimization result on success.

Errors

Returns DEError::InvalidBounds if any bound pair has upper < lower.

Example

use math_audio_optimisation::{differential_evolution, DEConfigBuilder};

let config = DEConfigBuilder::new()
    .maxiter(50)
    .seed(42)
    .build()
    .expect("invalid config");

let result = differential_evolution(
    &|x| x[0].powi(2) + x[1].powi(2),
    &[(-5.0, 5.0), (-5.0, 5.0)],
    config,
).expect("optimization failed");

assert!(result.fun < 0.01);

Examples found in repository

examples/optde_linear_constraints.rs (line 40)

fn main() {
    // Objective: sphere in 2D
    let sphere = |x: &Array1<f64>| x.iter().map(|v| v * v).sum::<f64>();

    // Bounds
    let bounds = [(-5.0, 5.0), (-5.0, 5.0)];

    // Linear constraint example: lb <= A x <= ub
    // 1) x0 + x1 <= 1.0
    // 2) 0.2 <= x0 - x1 <= 0.4
    let a = Array2::from_shape_vec((2, 2), vec![1.0, 1.0, 1.0, -1.0]).unwrap();
    let lb = Array1::from(vec![-f64::INFINITY, 0.2]);
    let ub = Array1::from(vec![1.0, 0.4]);
    let lc = LinearConstraintHelper { a, lb, ub };

    // Strategy parsing from string (mirrors SciPy names)
    let strategy = Strategy::from_str("randtobest1exp").unwrap_or(Strategy::RandToBest1Exp);

    // Build config using the fluent builder
    let mut cfg = DEConfigBuilder::new()
        .seed(123)
        .maxiter(600)
        .popsize(30)
        .strategy(strategy)
        .recombination(0.9)
        .mutation(Mutation::Range { min: 0.4, max: 1.0 })
        .crossover(Crossover::Exponential)
        .build()
        .expect("popsize must be >= 4");

    // Apply linear constraints with a penalty weight
    lc.apply_to(&mut cfg, 1e3);

    let rep = differential_evolution(&sphere, &bounds, cfg).expect("optimization failed");
    println!(
        "success={} message=\"{}\"\nbest f={:.6e}\nbest x={:?}",
        rep.success, rep.message, rep.fun, rep.x
    );
}

examples/optde_nonlinear_constraints.rs (line 41)

fn main() {
    // Himmelblau as objective, but with nonlinear constraints to demonstrate helper
    let himmelblau =
        |x: &Array1<f64>| (x[0] * x[0] + x[1] - 11.0).powi(2) + (x[0] + x[1] * x[1] - 7.0).powi(2);

    // Bounds
    let bounds = [(-6.0, 6.0), (-6.0, 6.0)];

    // Nonlinear vector function f(x) with 2 components
    // 1) Circle-ish constraint: x0^2 + x1^2 <= 10  -> f0(x) = x0^2 + x1^2,  lb=-inf, ub=10
    // 2) Sum equality: x0 + x1 = 1  -> f1(x) = x0 + x1,  lb=1, ub=1
    let fun =
        Arc::new(|x: &Array1<f64>| Array1::from(vec![x[0] * x[0] + x[1] * x[1], x[0] + x[1]]));
    let lb = Array1::from(vec![-f64::INFINITY, 1.0]);
    let ub = Array1::from(vec![10.0, 1.0]);
    let nlc = NonlinearConstraintHelper { fun, lb, ub };

    // Strategy parsing from string
    let strategy = Strategy::from_str("best1exp").unwrap_or(Strategy::Best1Exp);

    let mut cfg = DEConfigBuilder::new()
        .seed(456)
        .maxiter(800)
        .popsize(30)
        .strategy(strategy)
        .recombination(0.9)
        .crossover(Crossover::Exponential)
        .build()
        .expect("popsize must be >= 4");

    // Apply nonlinear constraints with penalties
    nlc.apply_to(&mut cfg, 1e3, 1e3);

    let rep = differential_evolution(&himmelblau, &bounds, cfg).expect("optimization failed");
    println!(
        "success={} message=\"{}\"\nbest f={:.6e}\nbest x={:?}",
        rep.success, rep.message, rep.fun, rep.x
    );
}

examples/optde_basic.rs (line 60)

fn main() {
    // Ackley function (2D)
    let ackley = |x: &Array1<f64>| {
        let x0 = x[0];
        let x1 = x[1];
        let s = 0.5 * (x0 * x0 + x1 * x1);
        let c = 0.5
            * ((2.0 * std::f64::consts::PI * x0).cos() + (2.0 * std::f64::consts::PI * x1).cos());
        -20.0 * (-0.2 * s.sqrt()).exp() - c.exp() + 20.0 + std::f64::consts::E
    };

    let bounds = [(-5.0, 5.0), (-5.0, 5.0)];

    let mut cfg = DEConfig::default();
    cfg.maxiter = 300;
    cfg.popsize = 20;
    cfg.strategy = Strategy::Best1Bin;
    cfg.crossover = Crossover::Exponential; // demonstrate exponential crossover
    cfg.mutation = Mutation::Range { min: 0.5, max: 1.0 }; // dithering
    cfg.recombination = 0.9;
    cfg.seed = Some(42);

    // Penalty examples (here just a dummy inequality fc(x) <= 0):
    // Circle of radius 3: x0^2 + x1^2 - 9 <= 0
    cfg.penalty_ineq.push((
        Arc::new(|x: &Array1<f64>| x[0] * x[0] + x[1] * x[1] - 9.0),
        1e3,
    ));

    // Callback every generation: stop early when convergence small enough
    let mut iter_log = 0usize;
    cfg.callback = Some(Box::new(move |inter| {
        if iter_log % 25 == 0 {
            eprintln!(
                "iter {:4}  best_f={:.6e}  conv(stdE)={:.3e}",
                inter.iter, inter.fun, inter.convergence
            );
        }
        iter_log += 1;
        if inter.convergence < 1e-6 {
            CallbackAction::Stop
        } else {
            CallbackAction::Continue
        }
    }));

    // Optional polishing with a local optimizer
    cfg.polish = Some(PolishConfig {
        enabled: true,
        algo: "neldermead".into(),
        maxeval: 400,
    });

    let report = differential_evolution(&ackley, &bounds, cfg).expect("optimization failed");

    println!(
        "success={} message=\"{}\"\nbest f={:.6e}\nbest x={:?}",
        report.success, report.message, report.fun, report.x
    );
}

examples/optde_parallel.rs (line 45)

fn main() {
    // Rastrigin function with artificial compute delay to simulate expensive evaluations
    let dimension = 10;
    let rastrigin = move |x: &Array1<f64>| -> f64 {
        // Add some compute-intensive work to make parallelization beneficial
        let mut sum = 0.0;
        for _ in 0..1000 {
            for &xi in x.iter() {
                sum += xi.sin().cos().exp().ln_1p();
            }
        }

        // Actual Rastrigin function
        let a = 10.0;
        let n = x.len() as f64;
        let result = a * n
            + x.iter()
                .map(|&xi| xi * xi - a * (2.0 * std::f64::consts::PI * xi).cos())
                .sum::<f64>();
        result + sum * 1e-10 // Add tiny contribution from expensive computation
    };

    let bounds: Vec<(f64, f64)> = vec![(-5.12, 5.12); dimension];

    // Test sequential evaluation
    println!("Testing Sequential Evaluation:");
    let mut cfg_seq = DEConfig::default();
    cfg_seq.maxiter = 100;
    cfg_seq.popsize = 15;
    cfg_seq.strategy = Strategy::Best1Bin;
    cfg_seq.mutation = Mutation::Factor(0.8);
    cfg_seq.recombination = 0.9;
    cfg_seq.seed = Some(42);
    cfg_seq.disp = true;
    cfg_seq.parallel.enabled = false; // Disable parallel evaluation

    let start_seq = Instant::now();
    let report_seq =
        differential_evolution(&rastrigin, &bounds, cfg_seq).expect("optimization failed");
    let duration_seq = start_seq.elapsed();

    println!("\nSequential Results:");
    println!("  Success: {}", report_seq.success);
    println!("  Best f: {:.6e}", report_seq.fun);
    println!("  Iterations: {}", report_seq.nit);
    println!("  Function evaluations: {}", report_seq.nfev);
    println!("  Time: {:.3} seconds", duration_seq.as_secs_f64());

    // Test parallel evaluation
    println!("\n\nTesting Parallel Evaluation:");
    let mut cfg_par = DEConfig::default();
    cfg_par.maxiter = 100;
    cfg_par.popsize = 15;
    cfg_par.strategy = Strategy::Best1Bin;
    cfg_par.mutation = Mutation::Factor(0.8);
    cfg_par.recombination = 0.9;
    cfg_par.seed = Some(42);
    cfg_par.disp = true;
    cfg_par.parallel = ParallelConfig {
        enabled: true,
        num_threads: None, // Use all available cores
    };

    let start_par = Instant::now();
    let report_par =
        differential_evolution(&rastrigin, &bounds, cfg_par).expect("optimization failed");
    let duration_par = start_par.elapsed();

    println!("\nParallel Results:");
    println!("  Success: {}", report_par.success);
    println!("  Best f: {:.6e}", report_par.fun);
    println!("  Iterations: {}", report_par.nit);
    println!("  Function evaluations: {}", report_par.nfev);
    println!("  Time: {:.3} seconds", duration_par.as_secs_f64());

    // Compare results
    println!("\n\nComparison:");
    println!(
        "  Speedup: {:.2}x",
        duration_seq.as_secs_f64() / duration_par.as_secs_f64()
    );
    println!(
        "  Result difference: {:.6e}",
        (report_seq.fun - report_par.fun).abs()
    );

    // Test with different thread counts
    println!("\n\nTesting with different thread counts:");
    for num_threads in [1, 2, 4, 8] {
        let mut cfg_threads = DEConfig::default();
        cfg_threads.maxiter = 50;
        cfg_threads.popsize = 15;
        cfg_threads.strategy = Strategy::Best1Bin;
        cfg_threads.mutation = Mutation::Factor(0.8);
        cfg_threads.recombination = 0.9;
        cfg_threads.seed = Some(42);
        cfg_threads.disp = false;
        cfg_threads.parallel = ParallelConfig {
            enabled: true,
            num_threads: Some(num_threads),
        };

        let start = Instant::now();
        let _ =
            differential_evolution(&rastrigin, &bounds, cfg_threads).expect("optimization failed");
        let duration = start.elapsed();

        println!(
            "  {} thread(s): {:.3} seconds",
            num_threads,
            duration.as_secs_f64()
        );
    }
}

examples/optde_adaptive_demo.rs (line 110)

fn main() {
    println!("🧬 Adaptive Differential Evolution Demo");
    println!("=====================================");
    println!();

    // Test functions to evaluate
    let test_functions = [
        (
            "Quadratic (f(x) = x₁² + x₂²)",
            quadratic as fn(&Array1<f64>) -> f64,
            [(-5.0, 5.0), (-5.0, 5.0)],
        ),
        (
            "Rosenbrock 2D",
            rosenbrock as fn(&Array1<f64>) -> f64,
            [(-5.0, 5.0), (-5.0, 5.0)],
        ),
        ("Ackley", ackley, [(-32.0, 32.0), (-32.0, 32.0)]),
    ];

    for (name, func, bounds) in test_functions.iter() {
        println!("🎯 Function: {}", name);
        println!(
            "   Bounds: [{:.1}, {:.1}] × [{:.1}, {:.1}]",
            bounds[0].0, bounds[0].1, bounds[1].0, bounds[1].1
        );

        // Traditional DE
        println!("   📊 Traditional DE:");
        let traditional_result = run_traditional_de(*func, bounds);

        // Adaptive DE with SAM only
        println!("   🧬 Adaptive DE (SAM only):");
        let sam_result = run_adaptive_de(*func, bounds, false);

        // Adaptive DE with SAM + WLS
        println!("   🔧 Adaptive DE (SAM + WLS):");
        let sam_wls_result = run_adaptive_de(*func, bounds, true);

        // Compare results
        println!("   🏆 Comparison:");
        println!(
            "      Traditional: f = {:.6e}, {} iterations",
            traditional_result.fun, traditional_result.nit
        );
        println!(
            "      SAM only:    f = {:.6e}, {} iterations",
            sam_result.fun, sam_result.nit
        );
        println!(
            "      SAM + WLS:   f = {:.6e}, {} iterations",
            sam_wls_result.fun, sam_wls_result.nit
        );

        let improvement_sam =
            ((traditional_result.fun - sam_result.fun) / traditional_result.fun * 100.0).max(0.0);
        let improvement_wls =
            ((traditional_result.fun - sam_wls_result.fun) / traditional_result.fun * 100.0)
                .max(0.0);

        println!("      📈 Improvement with SAM: {:.1}%", improvement_sam);
        println!("      📈 Improvement with WLS: {:.1}%", improvement_wls);
        println!();
    }

    // Demonstrate parameter adaptation tracking
    println!("🔄 Parameter Adaptation Demo");
    println!("===========================");

    // Use a recording callback to track parameter evolution
    let bounds = [(-5.0, 5.0), (-5.0, 5.0)];

    let adaptive_config = AdaptiveConfig {
        adaptive_mutation: true,
        wls_enabled: true,
        w_max: 0.9,     // Start with 90% of population for selection
        w_min: 0.1,     // End with 10% of population
        w_f: 0.9,       // F parameter adaptation rate
        w_cr: 0.9,      // CR parameter adaptation rate
        f_m: 0.5,       // Initial F location parameter
        cr_m: 0.6,      // Initial CR location parameter
        wls_prob: 0.2,  // Apply WLS to 20% of population
        wls_scale: 0.1, // WLS perturbation scale
    };

    let config = DEConfigBuilder::new()
        .seed(42)
        .maxiter(50)
        .popsize(40)
        .strategy(Strategy::AdaptiveBin)
        .mutation(Mutation::Adaptive { initial_f: 0.8 })
        .adaptive(adaptive_config)
        .disp(true) // Enable progress display
        .build()
        .expect("popsize must be >= 4");

    println!("Running adaptive DE on Rosenbrock function with progress display...");
    let result = differential_evolution(&rosenbrock, &bounds, config).expect("optimization failed");

    println!(
        "Final result: f = {:.6e} at x = [{:.4}, {:.4}]",
        result.fun, result.x[0], result.x[1]
    );
    println!(
        "Converged in {} iterations with {} function evaluations",
        result.nit, result.nfev
    );

    if result.success {
        println!("✅ Optimization succeeded: {}", result.message);
    } else {
        println!("⚠️ Optimization status: {}", result.message);
    }
}

fn run_traditional_de(
    func: fn(&Array1<f64>) -> f64,
    bounds: &[(f64, f64)],
) -> math_audio_optimisation::DEReport {
    let config = DEConfigBuilder::new()
        .seed(42)
        .maxiter(100)
        .popsize(30)
        .strategy(Strategy::Best1Bin)
        .mutation(Mutation::Factor(0.8))
        .recombination(0.7)
        .build()
        .expect("popsize must be >= 4");

    differential_evolution(&func, bounds, config).expect("optimization failed")
}

fn run_adaptive_de(
    func: fn(&Array1<f64>) -> f64,
    bounds: &[(f64, f64)],
    enable_wls: bool,
) -> math_audio_optimisation::DEReport {
    let adaptive_config = AdaptiveConfig {
        adaptive_mutation: true,
        wls_enabled: enable_wls,
        w_max: 0.9,
        w_min: 0.1,
        wls_prob: 0.15,
        wls_scale: 0.1,
        ..AdaptiveConfig::default()
    };

    let config = DEConfigBuilder::new()
        .seed(42)
        .maxiter(100)
        .popsize(30)
        .strategy(Strategy::AdaptiveBin)
        .mutation(Mutation::Adaptive { initial_f: 0.8 })
        .adaptive(adaptive_config)
        .build()
        .expect("popsize must be >= 4");

    differential_evolution(&func, bounds, config).expect("optimization failed")
}