scirs2-optimize 0.4.2

//! Tests for global optimization algorithms

use crate::global::{
    basinhopping, differential_evolution, dual_annealing, multi_start, particle_swarm,
    simulated_annealing, BasinHoppingOptions, DifferentialEvolutionOptions, DualAnnealingOptions,
    MultiStartOptions, ParticleSwarmOptions, SimulatedAnnealingOptions, StartingPointStrategy,
};
use crate::parallel::ParallelOptions;
use scirs2_core::ndarray::{array, ArrayView1};

#[test]
#[allow(dead_code)]
fn test_differential_evolution_rosenbrock() {
    // Rosenbrock function
    let rosenbrock = |x: &ArrayView1<f64>| {
        let a = 1.0;
        let b = 100.0;
        let x0 = x[0];
        let x1 = x[1];
        (a - x0).powi(2) + b * (x1 - x0.powi(2)).powi(2)
    };

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

    let options = DifferentialEvolutionOptions {
        maxiter: 100,
        popsize: 15,
        seed: Some(42),
        ..Default::default()
    };

    let result =
        differential_evolution(rosenbrock, bounds, Some(options), None).expect("Operation failed");

    assert!(result.success);
    assert!((result.x[0] - 1.0).abs() < 0.1);
    assert!((result.x[1] - 1.0).abs() < 0.1);
    assert!(result.fun < 0.01);
}

#[test]
#[allow(dead_code)]
fn test_basinhopping_sphere() {
    // Sphere function with minimum at origin
    let sphere = |x: &ArrayView1<f64>| x.iter().map(|&xi| xi.powi(2)).sum::<f64>();

    let x0 = array![1.0, 1.0];

    let options = BasinHoppingOptions {
        niter: 50,
        seed: Some(42),
        ..Default::default()
    };

    let result = basinhopping(sphere, x0, Some(options), None, None).expect("Operation failed");

    assert!(result.success);
    assert!(result.x.iter().all(|&xi| xi.abs() < 0.2));
    assert!(result.fun < 0.1);
}

#[test]
#[allow(dead_code)]
fn test_dual_annealing_rastrigin() {
    // Rastrigin function with many local minima
    let rastrigin = |x: &ArrayView1<f64>| {
        let n = x.len() as f64;
        10.0 * n
            + x.iter()
                .map(|&xi| xi.powi(2) - 10.0 * (2.0 * std::f64::consts::PI * xi).cos())
                .sum::<f64>()
    };

    let x0 = array![2.0, 2.0];
    let bounds = vec![(-5.12, 5.12), (-5.12, 5.12)];

    let options = DualAnnealingOptions {
        maxiter: 100,
        seed: Some(42),
        ..Default::default()
    };

    let result = dual_annealing(rastrigin, x0, bounds, Some(options)).expect("Operation failed");

    // Rastrigin global minimum is at origin with value 0
    assert!(result.x.iter().all(|&xi| xi.abs() < 0.5));
    assert!(result.fun < 1.0);
}

#[test]
#[allow(dead_code)]
fn test_differential_evolution_bounds() {
    // Function with minimum at (-1, -1), but bounded to positive quadrant
    let func = |x: &ArrayView1<f64>| (x[0] + 1.0).powi(2) + (x[1] + 1.0).powi(2);

    let bounds = vec![(0.0, 2.0), (0.0, 2.0)];

    let options = DifferentialEvolutionOptions {
        maxiter: 100,
        seed: Some(42),
        polish: true, // Re-enable polishing with bounds fix
        ..Default::default()
    };

    let result =
        differential_evolution(func, bounds, Some(options), None).expect("Operation failed");

    // Constrained minimum should be at (0, 0)
    assert!(result.success);
    // The bounds handling is working correctly - minimum should be near (0, 0)
    assert!(
        (result.fun - 2.0).abs() < 0.1,
        "Expected fun ≈ 2.0, got {}",
        result.fun
    );
}

#[test]
#[allow(dead_code)]
fn test_basinhopping_with_bounds() {
    let func = |x: &ArrayView1<f64>| x[0].powi(2) + x[1].powi(2);
    let x0 = array![1.0, 1.0];

    let options = BasinHoppingOptions {
        niter: 20,
        bounds: Some(vec![(0.5, 2.0), (0.5, 2.0)]),
        seed: Some(42),
        ..Default::default()
    };

    let result = basinhopping(func, x0, Some(options), None, None).expect("Operation failed");

    // With bounds, minimum should be at (0.5, 0.5)
    assert!(result.success);
    // The bounds handling is working correctly - minimum should be near (0.5, 0.5)
    assert!(result.fun < 0.6); // Should be 0.5 at the constrained minimum
}

#[test]
#[allow(dead_code)]
fn test_differential_evolution_different_strategies() {
    let sphere = |x: &ArrayView1<f64>| x.iter().map(|&xi| xi.powi(2)).sum::<f64>();
    let bounds = vec![(-10.0, 10.0), (-10.0, 10.0)];

    let strategies = ["best1bin", "rand1bin", "best2bin", "currenttobest1bin"];

    for strategy in &strategies {
        let options = DifferentialEvolutionOptions {
            maxiter: 50,
            popsize: 10,
            seed: Some(42),
            ..Default::default()
        };

        let result = differential_evolution(sphere, bounds.clone(), Some(options), Some(strategy))
            .expect("Operation failed");
        assert!(result.fun < 0.1, "Strategy {} failed to converge", strategy);
    }
}

#[test]
#[allow(dead_code)]
fn test_particle_swarm_sphere() {
    // Sphere function with minimum at origin
    let sphere = |x: &ArrayView1<f64>| x.iter().map(|&xi| xi.powi(2)).sum::<f64>();
    let bounds = vec![(-5.0, 5.0), (-5.0, 5.0)];

    let options = ParticleSwarmOptions {
        swarm_size: 30,
        maxiter: 100,
        seed: Some(42),
        ..Default::default()
    };

    let result = particle_swarm(sphere, bounds, Some(options)).expect("Operation failed");

    assert!(result.x.iter().all(|&xi| xi.abs() < 0.1));
    assert!(result.fun < 0.01);
}

#[test]
#[allow(dead_code)]
fn test_simulated_annealing_rosenbrock() {
    // Rosenbrock function
    let rosenbrock = |x: &ArrayView1<f64>| {
        let a = 1.0;
        let b = 100.0;
        let x0 = x[0];
        let x1 = x[1];
        (a - x0).powi(2) + b * (x1 - x0.powi(2)).powi(2)
    };

    let x0 = array![0.0, 0.0];
    let bounds = Some(vec![(-5.0, 5.0), (-5.0, 5.0)]);

    let options = SimulatedAnnealingOptions {
        maxiter: 1000,
        initial_temp: 10.0,
        final_temp: 0.001,
        seed: Some(42),
        ..Default::default()
    };

    let result =
        simulated_annealing(rosenbrock, x0, bounds, Some(options)).expect("Operation failed");

    // Simulated annealing might not always find exact solution
    assert!((result.x[0] - 1.0).abs() < 0.5);
    assert!((result.x[1] - 1.0).abs() < 0.5);
    assert!(result.fun < 5.0);
}

#[test]
#[allow(dead_code)]
fn test_particle_swarm_with_bounds() {
    // Function with minimum at (-1, -1), but bounded to positive quadrant
    let func = |x: &ArrayView1<f64>| (x[0] + 1.0).powi(2) + (x[1] + 1.0).powi(2);
    let bounds = vec![(0.0, 2.0), (0.0, 2.0)];

    let options = ParticleSwarmOptions {
        swarm_size: 30,
        maxiter: 200,
        seed: Some(42),
        tol: 1e-4,
        ..Default::default()
    };

    let result = particle_swarm(func, bounds, Some(options)).expect("Operation failed");

    // Constrained minimum should be at (0, 0) with f(0,0) = 2.0
    // Even if convergence criterion not met, verify solution quality
    assert!(result.x[0].abs() < 0.5);
    assert!(result.x[1].abs() < 0.5);
    assert!((result.fun - 2.0).abs() < 0.5);
}

#[test]
#[allow(dead_code)]
fn test_multi_start_basic() {
    // Simple sphere function
    let sphere = |x: &ArrayView1<f64>| x.iter().map(|&xi| xi.powi(2)).sum::<f64>();
    let bounds = vec![(-5.0, 5.0), (-5.0, 5.0)];

    let options = MultiStartOptions {
        n_starts: 5,
        parallel: false,
        seed: Some(42),
        strategy: StartingPointStrategy::Random,
        ..Default::default()
    };

    let result = multi_start(sphere, bounds, Some(options)).expect("Operation failed");

    assert!(result.success);
    assert!(result.x.iter().all(|&xi| xi.abs() < 0.1));
    assert!(result.fun < 0.01);
}

#[test]
#[allow(dead_code)]
fn test_multi_start_latin_hypercube() {
    // Multimodal function
    let func = |x: &ArrayView1<f64>| {
        x[0].powi(2)
            + x[1].powi(2)
            + 10.0 * (1.0 - (2.0 * x[0]).cos())
            + 10.0 * (1.0 - (2.0 * x[1]).cos())
    };
    let bounds = vec![(-5.0, 5.0), (-5.0, 5.0)];

    let options = MultiStartOptions {
        n_starts: 10,
        strategy: StartingPointStrategy::LatinHypercube,
        seed: Some(42),
        ..Default::default()
    };

    let result = multi_start(func, bounds, Some(options)).expect("Operation failed");

    // Should find a good minimum
    assert!(result.success);
    assert!(result.fun < 5.0);
}

#[test]
#[allow(dead_code)]
fn test_multi_start_grid() {
    // Function with multiple minima
    let himmelblau = |x: &ArrayView1<f64>| {
        let x0 = x[0];
        let x1 = x[1];
        (x0.powi(2) + x1 - 11.0).powi(2) + (x0 + x1.powi(2) - 7.0).powi(2)
    };

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

    let options = MultiStartOptions {
        n_starts: 16, // 4x4 grid
        strategy: StartingPointStrategy::Grid,
        seed: Some(42),
        ..Default::default()
    };

    let result = multi_start(himmelblau, bounds, Some(options)).expect("Operation failed");

    // Himmelblau's function has 4 global minima, all with value 0
    assert!(result.success);
    assert!(result.fun < 0.1);
}

#[test]
#[allow(dead_code)]
fn test_differential_evolution_parallel_basic() {
    // Simple test that parallel execution works without error
    let bounds = vec![(-5.0, 5.0), (-5.0, 5.0)];

    // Simple quadratic function
    let sphere = |x: &ArrayView1<f64>| x[0].powi(2) + x[1].powi(2);

    let mut options = DifferentialEvolutionOptions::default();
    options.popsize = 10;
    options.maxiter = 10;
    options.seed = Some(42);
    options.parallel = Some(ParallelOptions {
        num_workers: Some(2), // Use 2 workers
        min_parallel_size: 4,
        chunk_size: 1,
        parallel_evaluations: true,
        parallel_gradient: true,
    });

    let result =
        differential_evolution(sphere, bounds, Some(options), None).expect("Operation failed");

    // Should find minimum near origin
    assert!(result.x[0].abs() < 1.0);
    assert!(result.x[1].abs() < 1.0);
    assert!(result.fun < 2.0);
}

#[test]
#[allow(dead_code)]
fn test_differential_evolution_parallel_correctness() {
    // Test that parallel execution produces correct results
    let bounds = vec![(-10.0, 10.0), (-10.0, 10.0)];

    // Rosenbrock function with minimum at (1, 1)
    let rosenbrock = |x: &ArrayView1<f64>| {
        let x0 = x[0];
        let x1 = x[1];
        (1.0 - x0).powi(2) + 100.0 * (x1 - x0.powi(2)).powi(2)
    };

    let mut options = DifferentialEvolutionOptions::default();
    options.popsize = 20;
    options.maxiter = 50;
    options.seed = Some(123);
    options.parallel = Some(ParallelOptions {
        num_workers: None, // Use all available cores
        min_parallel_size: 5,
        chunk_size: 1,
        parallel_evaluations: true,
        parallel_gradient: true,
    });

    let result =
        differential_evolution(rosenbrock, bounds, Some(options), None).expect("Operation failed");

    // Check that we found the minimum
    assert!(result.fun < 0.1); // Should be very close to 0
    assert!((result.x[0] - 1.0).abs() < 0.2);
    assert!((result.x[1] - 1.0).abs() < 0.2);
}