feat: add optional parallel feature for population-evaluation parallelism

Adds a `parallel` Cargo feature that pulls in rayon and parallelizes
the only step that's actually expensive in practice — calls to
`Problem::evaluate` — across the population. RNG-driven steps (parent
and donor selection, variation, replacement decisions) stay serial, so
seeded runs remain deterministic regardless of feature state, and the
default and `--features parallel` builds produce bit-identical
results.

Wiring:
- New `algorithms::parallel_eval::evaluate_batch` helper with two
  cfg-gated implementations (rayon's `into_par_iter` when the feature
  is on, plain `into_iter` otherwise). Both preserve input order, so
  pareto_front and crowding-distance decisions remain reproducible.
- `RandomSearch`, `Nsga2`, and `DifferentialEvolution` now route
  population/offspring evaluation through the helper. NSGA-II's main
  loop is restructured into a serial selection-and-variation phase
  followed by a parallel-friendly batch evaluation phase.
- DE's per-target loop is restructured into three phases (serial trial
  construction → batch evaluation → serial replacement). Side effect
  of the restructuring: DE is now the canonical synchronous DE/rand/1/bin
  rather than the asynchronous variant where target `i+1` sees `i`'s
  in-flight update. Synchronous is the textbook formulation, so this
  is a small correctness improvement on top of the parallelism enable.
- PAES stays serial — its main loop has a sequential dependency on the
  current candidate and would gain nothing from rayon.

Cost: algorithm impls now require `P: Sync` and `P::Decision: Send`
unconditionally so a single impl serves both feature modes. This is a
small bound tightening that any plain-data Problem already satisfies; in
return the public `Problem` trait itself stays unchanged and the
default build picks up no new dependencies.

Verified:
- `cargo test` and `cargo test --features parallel` both pass; the
  Nsga2 `deterministic_with_same_seed` test confirms reproducibility.
- `cargo run --release --example benchmarks` and the same with
  `--features parallel` produce bit-identical ZDT1 / Rastrigin
  results.

Cargo.toml

@@ -9,8 +9,10 @@ readme = "README.md"
 [features]
 default = []
 serde = ["dep:serde"]
+parallel = ["dep:rayon"]
 
 [dependencies]
 rand = "0.9"
 rand_distr = "0.5"
+rayon = { version = "1", optional = true }
 serde = { version = "1", features = ["derive"], optional = true }

src/algorithms/differential_evolution.rs

@@ -2,6 +2,7 @@
 
 use rand::Rng as _;
 
+use crate::algorithms::parallel_eval::evaluate_batch;
 use crate::core::candidate::Candidate;
 use crate::core::objective::Direction;
 use crate::core::population::Population;
@@ -60,7 +61,7 @@ impl DifferentialEvolution {
 
 impl<P> Optimizer<P> for DifferentialEvolution
 where
-    P: Problem<Decision = Vec<f64>>,
+    P: Problem<Decision = Vec<f64>> + Sync,
 {
     fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision> {
         assert!(
@@ -88,57 +89,56 @@ where
             use crate::traits::Initializer as _;
             self.bounds.initialize(n, &mut rng)
         };
-        let mut evaluations = 0usize;
+        let initial_pop = evaluate_batch(problem, decisions.clone());
-        let mut evals: Vec<f64> = decisions
+        let mut evaluations = initial_pop.len();
-            .iter()
+        let mut evals: Vec<f64> =
-            .map(|d| {
+            initial_pop.iter().map(|c| c.evaluation.objectives[0]).collect();
-                let e = problem.evaluate(d);
-                evaluations += 1;
-                e.objectives[0]
-            })
-            .collect();
 
         for _gen in 0..self.config.generations {
-            for i in 0..n {
+            // Phase 1 (serial): construct one trial per target. RNG state is
-                let (r1, r2, r3) = pick_three_distinct(n, i, &mut rng);
+            // consumed in deterministic order so seeded runs reproduce
-                let j_rand = rng.random_range(0..dim);
+            // exactly regardless of the `parallel` feature.
-                let mut trial = decisions[i].clone();
+            let trials: Vec<Vec<f64>> = (0..n)
-                for j in 0..dim {
+                .map(|i| {
-                    let take_donor =
+                    let (r1, r2, r3) = pick_three_distinct(n, i, &mut rng);
-                        rng.random_bool(self.config.crossover_probability) || j == j_rand;
+                    let j_rand = rng.random_range(0..dim);
-                    if take_donor {
+                    let mut trial = decisions[i].clone();
-                        let mutant = decisions[r1][j]
+                    for j in 0..dim {
-                            + self.config.differential_weight
+                        let take_donor =
-                                * (decisions[r2][j] - decisions[r3][j]);
+                            rng.random_bool(self.config.crossover_probability) || j == j_rand;
-                        let (lo, hi) = self.bounds.bounds[j];
+                        if take_donor {
-                        trial[j] = mutant.clamp(lo, hi);
+                            let mutant = decisions[r1][j]
+                                + self.config.differential_weight
+                                    * (decisions[r2][j] - decisions[r3][j]);
+                            let (lo, hi) = self.bounds.bounds[j];
+                            trial[j] = mutant.clamp(lo, hi);
+                        }
                     }
-                }
+                    trial
-                let trial_obj = {
+                })
-                    let e = problem.evaluate(&trial);
+                .collect();
-                    evaluations += 1;
+
-                    e.objectives[0]
+            // Phase 2 (parallel-friendly): evaluate every trial.
-                };
+            let trial_cands = evaluate_batch(problem, trials);
+            evaluations += trial_cands.len();
+
+            // Phase 3 (serial): greedy replacement.
+            for (i, trial_cand) in trial_cands.into_iter().enumerate() {
+                let trial_obj = trial_cand.evaluation.objectives[0];
                 let target_obj = evals[i];
                 let trial_better = match direction {
                     Direction::Minimize => trial_obj <= target_obj,
                     Direction::Maximize => trial_obj >= target_obj,
                 };
                 if trial_better {
-                    decisions[i] = trial;
+                    decisions[i] = trial_cand.decision;
                     evals[i] = trial_obj;
                 }
             }
         }
 
-        let final_pop: Vec<Candidate<Vec<f64>>> = decisions
+        let final_pop: Vec<Candidate<Vec<f64>>> = evaluate_batch(problem, decisions);
-            .into_iter()
+        evaluations += final_pop.len();
-            .map(|d| {
-                let e = problem.evaluate(&d);
-                evaluations += 1;
-                Candidate::new(d, e)
-            })
-            .collect();
         let front = pareto_front(&final_pop, &objectives);
         let best = best_candidate(&final_pop, &objectives);
         OptimizationResult::new(

src/algorithms/mod.rs

@@ -3,6 +3,7 @@
 pub mod differential_evolution;
 pub mod nsga2;
 pub mod paes;
+pub(crate) mod parallel_eval;
 pub mod random_search;
 
 pub use differential_evolution::*;

src/algorithms/nsga2.rs

@@ -2,6 +2,7 @@
 
 use rand::Rng as _;
 
+use crate::algorithms::parallel_eval::evaluate_batch;
 use crate::core::candidate::Candidate;
 use crate::core::population::Population;
 use crate::core::problem::Problem;
@@ -56,7 +57,8 @@ struct Nsga2Entry<D> {
 
 impl<P, I, V> Optimizer<P> for Nsga2<I, V>
 where
-    P: Problem,
+    P: Problem + Sync,
+    P::Decision: Send,
     I: Initializer<P::Decision>,
     V: Variation<P::Decision>,
 {
@@ -76,13 +78,8 @@ where
             n,
             "NSGA-II initializer must return exactly population_size decisions",
         );
-        let mut population: Vec<Candidate<P::Decision>> = initial_decisions
+        let population: Vec<Candidate<P::Decision>> =
-            .into_iter()
+            evaluate_batch(problem, initial_decisions);
-            .map(|d| {
-                let e = problem.evaluate(&d);
-                Candidate::new(d, e)
-            })
-            .collect();
         let mut evaluations = population.len();
 
         // Annotate the starting population with rank and crowding so the first
@@ -90,9 +87,9 @@ where
         let mut annotated = annotate(population, &objectives);
 
         for _ in 0..self.config.generations {
-            // --- Parent selection + offspring generation ---
+            // --- Phase 1: serial parent selection + variation ---
-            let mut offspring: Vec<Candidate<P::Decision>> = Vec::with_capacity(n);
+            let mut offspring_decisions: Vec<P::Decision> = Vec::with_capacity(n);
-            while offspring.len() < n {
+            while offspring_decisions.len() < n {
                 let p1 = binary_tournament(&annotated, &mut rng);
                 let p2 = binary_tournament(&annotated, &mut rng);
                 let parents = vec![
@@ -105,14 +102,16 @@ where
                     "NSGA-II variation returned no children",
                 );
                 for child_decision in children {
-                    if offspring.len() >= n {
+                    if offspring_decisions.len() >= n {
                         break;
                     }
-                    let eval = problem.evaluate(&child_decision);
+                    offspring_decisions.push(child_decision);
-                    evaluations += 1;
-                    offspring.push(Candidate::new(child_decision, eval));
                 }
             }
+            // --- Phase 2: parallel-friendly batch evaluation ---
+            let offspring: Vec<Candidate<P::Decision>> =
+                evaluate_batch(problem, offspring_decisions);
+            evaluations += offspring.len();
 
             // --- Combine + survival selection ---
             let mut combined: Vec<Candidate<P::Decision>> = Vec::with_capacity(2 * n);
@@ -143,8 +142,7 @@ where
                     break;
                 }
             }
-            population = next;
+            annotated = annotate(next, &objectives);
-            annotated = annotate(population, &objectives);
         }
 
         // Return final state.

src/algorithms/parallel_eval.rs

@@ -0,0 +1,55 @@
+//! Population-wide evaluation helper used by the population-based algorithms.
+//!
+//! Two cfg-gated implementations:
+//!
+//! - With the `parallel` feature: rayon's `into_par_iter` evaluates decisions
+//!   on the global thread pool. The result preserves input order so any
+//!   downstream sort/dominance/crowding decisions remain reproducible.
+//! - Without the feature: plain serial `into_iter`.
+//!
+//! Both implementations require `P: Problem + Sync` and `P::Decision: Send`,
+//! so each algorithm's `Optimizer<P>` impl carries the same bounds regardless
+//! of feature state. This keeps the public `Problem` trait itself unchanged.
+
+use crate::core::candidate::Candidate;
+use crate::core::problem::Problem;
+
+/// Evaluate every decision in `decisions` against `problem` and return the
+/// resulting candidates in the same order.
+#[cfg(feature = "parallel")]
+pub(crate) fn evaluate_batch<P>(
+    problem: &P,
+    decisions: Vec<P::Decision>,
+) -> Vec<Candidate<P::Decision>>
+where
+    P: Problem + Sync,
+    P::Decision: Send,
+{
+    use rayon::prelude::*;
+    decisions
+        .into_par_iter()
+        .map(|d| {
+            let e = problem.evaluate(&d);
+            Candidate::new(d, e)
+        })
+        .collect()
+}
+
+/// Serial fallback used when the `parallel` feature is disabled.
+#[cfg(not(feature = "parallel"))]
+pub(crate) fn evaluate_batch<P>(
+    problem: &P,
+    decisions: Vec<P::Decision>,
+) -> Vec<Candidate<P::Decision>>
+where
+    P: Problem + Sync,
+    P::Decision: Send,
+{
+    decisions
+        .into_iter()
+        .map(|d| {
+            let e = problem.evaluate(&d);
+            Candidate::new(d, e)
+        })
+        .collect()
+}

src/algorithms/random_search.rs

@@ -3,6 +3,7 @@
 //! This is the reference example for spec §2.4 / §12.1: read this file before
 //! writing your own optimizer.
 
+use crate::algorithms::parallel_eval::evaluate_batch;
 use crate::core::candidate::Candidate;
 use crate::core::population::Population;
 use crate::core::problem::Problem;
@@ -50,7 +51,8 @@ impl<I> RandomSearch<I> {
 
 impl<P, I> Optimizer<P> for RandomSearch<I>
 where
-    P: Problem,
+    P: Problem + Sync,
+    P::Decision: Send,
     I: Initializer<P::Decision>,
 {
     fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision> {
@@ -61,11 +63,8 @@ where
 
         for _ in 0..self.config.iterations {
             let decisions = self.initializer.initialize(self.config.batch_size, &mut rng);
-            for decision in decisions {
+            evaluations += decisions.len();
-                let eval = problem.evaluate(&decision);
+            all.extend(evaluate_batch(problem, decisions));
-                evaluations += 1;
-                all.push(Candidate::new(decision, eval));
-            }
         }
 
         let front = pareto_front(&all, &objectives);