feat(algorithms): add SPEA2 (Strength Pareto Evolutionary Algorithm 2)

Implementation of Zitzler, Laumanns, Thiele 2001 SPEA2 — the classic
Pareto MOEA built around an explicit external archive of fixed size.

Each generation:
1. Combine the current population and the archive into one pool.
2. For every member, compute strength S(i) = number of others that
   member dominates, then raw fitness R(i) = sum of S(j) over members
   j that dominate i.
3. Add a density estimator D(i) = 1/(σ_k + 2) where σ_k is the distance
   to the k-th nearest neighbor (k = floor(sqrt(|pool|))) in
   minimization-oriented objective space.
4. Final fitness F(i) = R(i) + D(i); lower is better.
5. Build the next archive by taking every non-dominated member
   (R(i) == 0). If too many, prune by repeatedly removing the member
   with the smallest k-th-nearest-neighbor distance. If too few, fill
   from the rest sorted by F ascending.
6. Generate the next population by binary tournament on F (lower wins),
   then variation, then evaluation.

Public API mirrors the other algorithms:

  Spea2Config { population_size, archive_size, generations, seed }
  Spea2 { config, initializer, variation }
  impl<P, I, V> Optimizer<P> for Spea2<I, V>

Re-exported from the prelude. Tests cover archive size invariants,
non-empty Pareto front on Schaffer N.1, deterministic reruns under
the same seed, and panic on population_size == 0.

src/algorithms/mod.rs

@@ -5,8 +5,10 @@ pub mod nsga2;
 pub mod paes;
 pub(crate) mod parallel_eval;
 pub mod random_search;
+pub mod spea2;
 
 pub use differential_evolution::*;
 pub use nsga2::*;
 pub use paes::*;
 pub use random_search::*;
+pub use spea2::*;

src/algorithms/spea2.rs

@@ -0,0 +1,379 @@
+//! SPEA2 — Strength Pareto Evolutionary Algorithm 2 (Zitzler, Laumanns, Thiele 2001).
+
+use rand::Rng as _;
+
+use crate::algorithms::parallel_eval::evaluate_batch;
+use crate::core::candidate::Candidate;
+use crate::core::objective::ObjectiveSpace;
+use crate::core::population::Population;
+use crate::core::problem::Problem;
+use crate::core::result::OptimizationResult;
+use crate::core::rng::{Rng, rng_from_seed};
+use crate::pareto::dominance::{Dominance, pareto_compare};
+use crate::pareto::front::{best_candidate, pareto_front};
+use crate::traits::{Initializer, Optimizer, Variation};
+
+/// Configuration for [`Spea2`].
+#[derive(Debug, Clone)]
+pub struct Spea2Config {
+    /// Constant population size carried across generations.
+    pub population_size: usize,
+    /// Constant archive size; SPEA2 grows or shrinks it to this exact target.
+    pub archive_size: usize,
+    /// Number of generations to run.
+    pub generations: usize,
+    /// Seed for the deterministic RNG.
+    pub seed: u64,
+}
+
+impl Default for Spea2Config {
+    fn default() -> Self {
+        Self { population_size: 100, archive_size: 100, generations: 250, seed: 42 }
+    }
+}
+
+/// SPEA2 optimizer.
+#[derive(Debug, Clone)]
+pub struct Spea2<I, V> {
+    /// Algorithm configuration.
+    pub config: Spea2Config,
+    /// Initial-decision sampler.
+    pub initializer: I,
+    /// Offspring-producing variation operator.
+    pub variation: V,
+}
+
+impl<I, V> Spea2<I, V> {
+    /// Construct a `Spea2` optimizer.
+    pub fn new(config: Spea2Config, initializer: I, variation: V) -> Self {
+        Self { config, initializer, variation }
+    }
+}
+
+impl<P, I, V> Optimizer<P> for Spea2<I, V>
+where
+    P: Problem + Sync,
+    P::Decision: Send,
+    I: Initializer<P::Decision>,
+    V: Variation<P::Decision>,
+{
+    fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision> {
+        assert!(
+            self.config.population_size > 0,
+            "Spea2 population_size must be greater than 0",
+        );
+        assert!(
+            self.config.archive_size > 0,
+            "Spea2 archive_size must be greater than 0",
+        );
+        let n_pop = self.config.population_size;
+        let n_arc = self.config.archive_size;
+        let objectives = problem.objectives();
+        let mut rng = rng_from_seed(self.config.seed);
+
+        // Initial population.
+        let initial_decisions = self.initializer.initialize(n_pop, &mut rng);
+        assert_eq!(
+            initial_decisions.len(),
+            n_pop,
+            "SPEA2 initializer must return exactly population_size decisions",
+        );
+        let mut population: Vec<Candidate<P::Decision>> =
+            evaluate_batch(problem, initial_decisions);
+        let mut evaluations = population.len();
+        let mut archive: Vec<Candidate<P::Decision>> = Vec::new();
+
+        for _ in 0..self.config.generations {
+            // --- Combine pool, compute fitness ---
+            let mut pool: Vec<Candidate<P::Decision>> =
+                Vec::with_capacity(population.len() + archive.len());
+            pool.extend(population.drain(..));
+            pool.extend(archive.drain(..));
+            let fitness = compute_fitness(&pool, &objectives);
+
+            // --- Build the next archive ---
+            archive = build_archive(&pool, &fitness, &objectives, n_arc);
+
+            // --- Generate offspring from the archive (mating pool) ---
+            let archive_fitness = compute_fitness(&archive, &objectives);
+            let mut offspring_decisions: Vec<P::Decision> = Vec::with_capacity(n_pop);
+            while offspring_decisions.len() < n_pop {
+                let p1 = binary_tournament(&archive_fitness, &mut rng);
+                let p2 = binary_tournament(&archive_fitness, &mut rng);
+                let parents = vec![archive[p1].decision.clone(), archive[p2].decision.clone()];
+                let children = self.variation.vary(&parents, &mut rng);
+                assert!(!children.is_empty(), "SPEA2 variation returned no children");
+                for child_decision in children {
+                    if offspring_decisions.len() >= n_pop {
+                        break;
+                    }
+                    offspring_decisions.push(child_decision);
+                }
+            }
+            let new_population = evaluate_batch(problem, offspring_decisions);
+            evaluations += new_population.len();
+            population = new_population;
+        }
+
+        let front = pareto_front(&archive, &objectives);
+        let best = best_candidate(&archive, &objectives);
+        OptimizationResult::new(
+            Population::new(archive),
+            front,
+            best,
+            evaluations,
+            self.config.generations,
+        )
+    }
+}
+
+/// SPEA2 fitness: `R(i) + D(i)`, where lower is better.
+///
+/// `R(i)` is the sum of `S(j)` over all `j` that dominate `i`. `S(j)` is the
+/// count of members `j` dominates. `D(i) = 1 / (σ_k + 2)` where `σ_k` is the
+/// distance to the k-th nearest neighbor (k = floor(sqrt(N))) in
+/// minimization-oriented objective space.
+fn compute_fitness<D>(pool: &[Candidate<D>], objectives: &ObjectiveSpace) -> Vec<f64> {
+    let n = pool.len();
+    if n == 0 {
+        return Vec::new();
+    }
+    let oriented: Vec<Vec<f64>> = pool
+        .iter()
+        .map(|c| objectives.as_minimization(&c.evaluation.objectives))
+        .collect();
+
+    // Strength S(i) = number of members i dominates.
+    let mut strength = vec![0_usize; n];
+    let mut dominators_of: Vec<Vec<usize>> = vec![Vec::new(); n];
+    for i in 0..n {
+        for j in 0..n {
+            if i == j {
+                continue;
+            }
+            if matches!(
+                pareto_compare(&pool[i].evaluation, &pool[j].evaluation, objectives),
+                Dominance::Dominates
+            ) {
+                strength[i] += 1;
+                dominators_of[j].push(i);
+            }
+        }
+    }
+
+    // Raw fitness R(i) = sum of S(j) over j that dominate i.
+    let raw: Vec<f64> = (0..n)
+        .map(|i| dominators_of[i].iter().map(|&j| strength[j] as f64).sum())
+        .collect();
+
+    // Density D(i) = 1 / (σ_k + 2). Use kth_nearest distances.
+    let k = (n as f64).sqrt() as usize;
+    let density: Vec<f64> = (0..n)
+        .map(|i| {
+            let mut dists: Vec<f64> = (0..n)
+                .filter(|&j| j != i)
+                .map(|j| euclidean(&oriented[i], &oriented[j]))
+                .collect();
+            dists.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
+            // SPEA2's σ_k is the distance to the k-th nearest neighbor (1-indexed).
+            // With k = floor(sqrt(N)), use index (k-1).clamp(0, len-1).
+            let idx = if dists.is_empty() {
+                return 0.0;
+            } else {
+                k.saturating_sub(1).min(dists.len() - 1)
+            };
+            1.0 / (dists[idx] + 2.0)
+        })
+        .collect();
+
+    raw.into_iter().zip(density).map(|(r, d)| r + d).collect()
+}
+
+fn euclidean(a: &[f64], b: &[f64]) -> f64 {
+    a.iter().zip(b.iter()).map(|(x, y)| (x - y).powi(2)).sum::<f64>().sqrt()
+}
+
+/// Build the next archive of exactly `target_size` members.
+///
+/// All non-dominated members of the pool (`fitness < 1.0`) are taken first.
+/// If too many, prune by iteratively removing the member with the smallest
+/// distance to its nearest neighbor (ties broken by next-nearest, etc.). If
+/// too few, fill from the rest sorted by fitness ascending.
+fn build_archive<D: Clone>(
+    pool: &[Candidate<D>],
+    fitness: &[f64],
+    objectives: &ObjectiveSpace,
+    target_size: usize,
+) -> Vec<Candidate<D>> {
+    let mut nondom: Vec<usize> = (0..pool.len()).filter(|&i| fitness[i] < 1.0).collect();
+
+    if nondom.len() == target_size {
+        return nondom.into_iter().map(|i| pool[i].clone()).collect();
+    }
+
+    if nondom.len() < target_size {
+        // Fill from dominated members ordered by ascending fitness.
+        let mut dominated: Vec<usize> =
+            (0..pool.len()).filter(|&i| fitness[i] >= 1.0).collect();
+        dominated.sort_by(|&a, &b| {
+            fitness[a].partial_cmp(&fitness[b]).unwrap_or(std::cmp::Ordering::Equal)
+        });
+        let needed = target_size - nondom.len();
+        nondom.extend(dominated.into_iter().take(needed));
+        return nondom.into_iter().map(|i| pool[i].clone()).collect();
+    }
+
+    // Truncation: while too large, drop the member with the smallest distance
+    // to its nearest neighbor in the current archive.
+    let oriented: Vec<Vec<f64>> = nondom
+        .iter()
+        .map(|&i| objectives.as_minimization(&pool[i].evaluation.objectives))
+        .collect();
+    let mut alive: Vec<bool> = vec![true; nondom.len()];
+    let mut alive_count = nondom.len();
+    while alive_count > target_size {
+        // Compute per-member sorted distances to other alive members.
+        let mut neighbor_dists: Vec<Vec<f64>> = vec![Vec::new(); nondom.len()];
+        for i in 0..nondom.len() {
+            if !alive[i] {
+                continue;
+            }
+            for j in 0..nondom.len() {
+                if !alive[j] || i == j {
+                    continue;
+                }
+                neighbor_dists[i].push(euclidean(&oriented[i], &oriented[j]));
+            }
+            neighbor_dists[i]
+                .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
+        }
+        // Find the alive member whose neighbor-distance vector is lex-smallest.
+        let mut victim = usize::MAX;
+        for i in 0..nondom.len() {
+            if !alive[i] {
+                continue;
+            }
+            if victim == usize::MAX {
+                victim = i;
+                continue;
+            }
+            // Lex-compare neighbor distances.
+            let cmp = neighbor_dists[i]
+                .iter()
+                .zip(neighbor_dists[victim].iter())
+                .find_map(|(a, b)| {
+                    let c = a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal);
+                    if c != std::cmp::Ordering::Equal { Some(c) } else { None }
+                })
+                .unwrap_or(std::cmp::Ordering::Equal);
+            if cmp == std::cmp::Ordering::Less {
+                victim = i;
+            }
+        }
+        alive[victim] = false;
+        alive_count -= 1;
+    }
+
+    nondom
+        .into_iter()
+        .enumerate()
+        .filter_map(|(local, idx)| if alive[local] { Some(pool[idx].clone()) } else { None })
+        .collect()
+}
+
+fn binary_tournament(fitness: &[f64], rng: &mut Rng) -> usize {
+    let n = fitness.len();
+    let a = rng.random_range(0..n);
+    let b = rng.random_range(0..n);
+    if fitness[a] < fitness[b] {
+        a
+    } else if fitness[a] > fitness[b] {
+        b
+    } else if rng.random_bool(0.5) {
+        a
+    } else {
+        b
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::operators::{GaussianMutation, RealBounds};
+    use crate::tests_support::SchafferN1;
+
+    #[test]
+    fn produces_pareto_front() {
+        let mut opt = Spea2::new(
+            Spea2Config {
+                population_size: 30,
+                archive_size: 30,
+                generations: 10,
+                seed: 1,
+            },
+            RealBounds::new(vec![(-5.0, 5.0)]),
+            GaussianMutation { sigma: 0.3 },
+        );
+        let r = opt.run(&SchafferN1);
+        assert!(!r.pareto_front.is_empty());
+        assert_eq!(r.population.len(), 30);
+        assert_eq!(r.generations, 10);
+    }
+
+    #[test]
+    fn archive_size_respected() {
+        let mut opt = Spea2::new(
+            Spea2Config {
+                population_size: 40,
+                archive_size: 20,
+                generations: 15,
+                seed: 2,
+            },
+            RealBounds::new(vec![(-5.0, 5.0)]),
+            GaussianMutation { sigma: 0.3 },
+        );
+        let r = opt.run(&SchafferN1);
+        assert_eq!(r.population.len(), 20);
+    }
+
+    #[test]
+    fn deterministic_with_same_seed() {
+        let make = || {
+            Spea2::new(
+                Spea2Config {
+                    population_size: 20,
+                    archive_size: 20,
+                    generations: 10,
+                    seed: 99,
+                },
+                RealBounds::new(vec![(-5.0, 5.0)]),
+                GaussianMutation { sigma: 0.2 },
+            )
+        };
+        let mut a = make();
+        let mut b = make();
+        let ra = a.run(&SchafferN1);
+        let rb = b.run(&SchafferN1);
+        let oa: Vec<Vec<f64>> =
+            ra.pareto_front.iter().map(|c| c.evaluation.objectives.clone()).collect();
+        let ob: Vec<Vec<f64>> =
+            rb.pareto_front.iter().map(|c| c.evaluation.objectives.clone()).collect();
+        assert_eq!(oa, ob);
+    }
+
+    #[test]
+    #[should_panic(expected = "population_size must be greater than 0")]
+    fn zero_population_size_panics() {
+        let mut opt = Spea2::new(
+            Spea2Config {
+                population_size: 0,
+                archive_size: 10,
+                generations: 1,
+                seed: 0,
+            },
+            RealBounds::new(vec![(-1.0, 1.0)]),
+            GaussianMutation { sigma: 0.1 },
+        );
+        let _ = opt.run(&SchafferN1);
+    }
+}

src/prelude.rs

@@ -23,5 +23,5 @@ pub use crate::operators::{
 
 pub use crate::algorithms::{
     DifferentialEvolution, DifferentialEvolutionConfig, Nsga2, Nsga2Config, Paes, PaesConfig,
-    RandomSearch, RandomSearchConfig,
+    RandomSearch, RandomSearchConfig, Spea2, Spea2Config,
 };