chore: release v0.9.1

This reverts commit ed80c59795.

# Conflicts:
#	src/lib.rs
#	src/storage/sqlite.rs

.github/workflows/ci.yml

@@ -42,11 +42,6 @@ jobs:
   test:
     name: Test
     runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        features:
-          - ""
-          - "async"
     steps:
       - uses: actions/checkout@v6
       - name: Install Rust
@@ -55,12 +50,7 @@ jobs:
           rustup update stable
       - uses: Swatinem/rust-cache@v2
       - name: Run tests
-        run: |
-          if [ -z "${{ matrix.features }}" ]; then
-            cargo test --verbose
-          else
-            cargo test --verbose --features "${{ matrix.features }}"
-          fi
+        run: cargo test --verbose --all-features --all-targets
 
   examples:
     name: Examples
@@ -81,7 +71,7 @@ jobs:
       - name: Run async example (async_api_optimization)
         run: cargo run --example async_api_optimization --features async
       - name: Run visualization example (visualization_demo)
-        run: cargo run --example visualization_demo --features visualization
+        run: cargo run --example visualization_demo
 
   docs:
     name: Docs
@@ -98,8 +88,12 @@ jobs:
           RUSTDOCFLAGS: -D warnings
 
   feature-check:
-    name: Feature Combinations
+    name: Feature Combinations (${{ matrix.partition }}/4)
     runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        partition: [1, 2, 3, 4]
     steps:
       - uses: actions/checkout@v6
       - name: Install Rust
@@ -110,7 +104,7 @@ jobs:
         uses: taiki-e/install-action@cargo-hack
       - uses: Swatinem/rust-cache@v2
       - name: Check feature powerset
-        run: cargo hack check --feature-powerset --no-dev-deps
+        run: cargo hack check --feature-powerset --no-dev-deps --partition ${{ matrix.partition }}/4
 
   coverage:
     name: Coverage
@@ -135,14 +129,14 @@ jobs:
           fail_ci_if_error: true
 
   msrv:
-    name: MSRV (1.88)
+    name: MSRV (1.89)
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v6
-      - name: Install Rust 1.88
+      - name: Install Rust 1.89
         run: |
-          rustup override set 1.88
-          rustup update 1.88
+          rustup override set 1.89
+          rustup update 1.89
       - uses: Swatinem/rust-cache@v2
       - name: Check compilation
         run: cargo check --all-features
@@ -370,6 +364,7 @@ jobs:
       - cross-platform
       - cross-targets
       - typos
+      - bench-check
     steps:
       - uses: actions/checkout@v6
       - name: Install Rust

Cargo.toml

@@ -3,9 +3,9 @@ members = ["optimizer-derive"]
 
 [package]
 name = "optimizer"
-version = "0.8.1"
+version = "0.9.1"
 edition = "2024"
-rust-version = "1.88"
+rust-version = "1.89"
 license = "MIT"
 authors = ["Manuel Raimann <raimannma@outlook.de"]
 description = "A Rust library for optimization algorithms."
@@ -16,7 +16,7 @@ categories = ["algorithm", "science", "data-structures"]
 readme = "README.md"
 
 [dependencies]
-rand = "0.10"
+fastrand = "2.3"
 thiserror = "2"
 parking_lot = "0.12"
 tokio = { version = "1", features = ["sync", "rt-multi-thread"], optional = true }
@@ -26,17 +26,18 @@ serde_json = { version = "1", optional = true }
 tracing = { version = "0.1.29", optional = true }
 sobol_burley = { version = "0.5", optional = true }
 nalgebra = { version = "0.34", optional = true }
+fs2 = { version = "0.4", optional = true }
 
 [features]
 default = []
 async = ["dep:tokio"]
 derive = ["dep:optimizer-derive"]
 serde = ["dep:serde", "dep:serde_json"]
+journal = ["dep:fs2", "serde"]
 tracing = ["dep:tracing"]
 sobol = ["dep:sobol_burley"]
 cma-es = ["dep:nalgebra"]
-visualization = []
-fanova = []
+gp = ["dep:nalgebra"]
 
 [dev-dependencies]
 tokio = { version = "1", features = ["rt-multi-thread", "macros", "time"] }
@@ -69,4 +70,7 @@ required-features = ["derive"]
 [[example]]
 name = "visualization_demo"
 path = "examples/visualization_demo.rs"
-required-features = ["visualization"]
+
+[[test]]
+name = "journal_tests"
+required-features = ["journal"]

deny.toml

@@ -14,6 +14,7 @@ allow = [
     "Apache-2.0",
     "Apache-2.0 WITH LLVM-exception",
     "Unicode-3.0",
+    "Zlib",
 ]
 confidence-threshold = 0.8
 

optimizer-derive/Cargo.toml

@@ -2,7 +2,7 @@
 name = "optimizer-derive"
 version = "0.1.0"
 edition = "2024"
-rust-version = "1.88"
+rust-version = "1.89"
 license = "MIT"
 description = "Derive macros for the optimizer crate"
 repository = "https://github.com/raimannma/rust-optimizer"

src/error.rs

@@ -93,6 +93,11 @@ pub enum Error {
     #[cfg(feature = "async")]
     #[error("async task error: {0}")]
     TaskError(String),
+
+    /// Returned when a storage operation fails.
+    #[cfg(feature = "journal")]
+    #[error("storage error: {0}")]
+    Storage(String),
 }
 
 pub type Result<T> = core::result::Result<T, Error>;

src/fanova.rs

@@ -9,9 +9,6 @@
 //! 3. Computes main effects (single-parameter importance)
 //! 4. Computes interaction effects (pairwise parameter importance)
 
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
-
 /// Result of fANOVA analysis.
 #[derive(Debug, Clone)]
 pub struct FanovaResult {
@@ -81,7 +78,7 @@ impl DecisionTree {
         targets: &[f64],
         indices: &[usize],
         config: &FanovaConfig,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> Self {
         let mut tree = Self { nodes: Vec::new() };
         tree.build_node(data, targets, indices, 0, config, rng);
@@ -96,7 +93,7 @@ impl DecisionTree {
         indices: &[usize],
         depth: usize,
         config: &FanovaConfig,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> usize {
         let n = indices.len();
         let mean = indices.iter().map(|&i| targets[i]).sum::<f64>() / n as f64;
@@ -264,11 +261,11 @@ impl DecisionTree {
 // --- Helper Functions ---
 
 /// Select `k` random indices from `0..n` using partial Fisher-Yates shuffle.
-fn partial_shuffle(n: usize, k: usize, rng: &mut StdRng) -> Vec<usize> {
+fn partial_shuffle(n: usize, k: usize, rng: &mut fastrand::Rng) -> Vec<usize> {
     let mut indices: Vec<usize> = (0..n).collect();
     let k = k.min(n);
     for i in 0..k {
-        let j = rng.random_range(i..n);
+        let j = rng.usize(i..n);
         indices.swap(i, j);
     }
     indices.truncate(k);
@@ -331,16 +328,14 @@ pub(crate) fn compute_fanova(
     let n_samples = data.len();
     let n_features = data[0].len();
 
-    let mut rng: StdRng = config
+    let mut rng: fastrand::Rng = config
         .seed
-        .map_or_else(rand::make_rng, StdRng::seed_from_u64);
+        .map_or_else(fastrand::Rng::new, fastrand::Rng::with_seed);
 
     // Build random forest with bootstrap sampling
     let trees: Vec<DecisionTree> = (0..config.n_trees)
         .map(|_| {
-            let bootstrap: Vec<usize> = (0..n_samples)
-                .map(|_| rng.random_range(0..n_samples))
-                .collect();
+            let bootstrap: Vec<usize> = (0..n_samples).map(|_| rng.usize(0..n_samples)).collect();
             DecisionTree::build(data, targets, &bootstrap, config, &mut rng)
         })
         .collect();
@@ -414,14 +409,20 @@ pub(crate) fn compute_fanova(
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::rng_util;
 
     #[test]
     fn single_dominant_parameter() {
         // f(x, y) = x — only x matters
-        let mut rng = StdRng::seed_from_u64(0);
+        let mut rng = fastrand::Rng::with_seed(0);
         let n = 100;
         let data: Vec<Vec<f64>> = (0..n)
-            .map(|_| vec![rng.random_range(0.0..10.0), rng.random_range(0.0..10.0)])
+            .map(|_| {
+                vec![
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                ]
+            })
             .collect();
         let targets: Vec<f64> = data.iter().map(|row| row[0]).collect();
 
@@ -443,10 +444,15 @@ mod tests {
     #[test]
     fn interaction_detection() {
         // f(x, y) = x * y — both matter and interact
-        let mut rng = StdRng::seed_from_u64(0);
+        let mut rng = fastrand::Rng::with_seed(42);
         let n = 200;
         let data: Vec<Vec<f64>> = (0..n)
-            .map(|_| vec![rng.random_range(0.0..10.0), rng.random_range(0.0..10.0)])
+            .map(|_| {
+                vec![
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                ]
+            })
             .collect();
         let targets: Vec<f64> = data.iter().map(|row| row[0] * row[1]).collect();
 
@@ -477,14 +483,14 @@ mod tests {
     #[test]
     fn three_params_one_dominant() {
         // f(x, y, z) = 3*x + 0.1*y + 0*z
-        let mut rng = StdRng::seed_from_u64(7);
+        let mut rng = fastrand::Rng::with_seed(7);
         let n = 150;
         let data: Vec<Vec<f64>> = (0..n)
             .map(|_| {
                 vec![
-                    rng.random_range(0.0..10.0),
-                    rng.random_range(0.0..10.0),
-                    rng.random_range(0.0..10.0),
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
                 ]
             })
             .collect();
@@ -512,10 +518,15 @@ mod tests {
 
     #[test]
     fn importances_sum_to_one() {
-        let mut rng = StdRng::seed_from_u64(3);
+        let mut rng = fastrand::Rng::with_seed(3);
         let n = 100;
         let data: Vec<Vec<f64>> = (0..n)
-            .map(|_| vec![rng.random_range(0.0..10.0), rng.random_range(0.0..10.0)])
+            .map(|_| {
+                vec![
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                    rng_util::f64_range(&mut rng, 0.0, 10.0),
+                ]
+            })
             .collect();
         let targets: Vec<f64> = data.iter().map(|r| r[0] + r[1]).collect();
 

src/kde/multivariate.rs

@@ -5,8 +5,6 @@
 //! parameter independently, the multivariate KDE models the joint distribution
 //! to better capture correlations between parameters.
 
-use rand::{Rng, RngExt};
-
 use crate::error::{Error, Result};
 
 /// A multivariate Gaussian kernel density estimator for joint distributions.
@@ -308,9 +306,9 @@ impl MultivariateKDE {
     /// # Returns
     ///
     /// A `Vec<f64>` of length `n_dims` representing a sample from the KDE.
-    pub(crate) fn sample<R: Rng>(&self, rng: &mut R) -> Vec<f64> {
+    pub(crate) fn sample(&self, rng: &mut fastrand::Rng) -> Vec<f64> {
         // Select a random sample to center the kernel on
-        let idx = rng.random_range(0..self.samples.len());
+        let idx = rng.usize(0..self.samples.len());
         let center = &self.samples[idx];
 
         // Add independent Gaussian noise to each dimension
@@ -319,8 +317,8 @@ impl MultivariateKDE {
             .iter()
             .zip(self.bandwidths.iter())
             .map(|(&center_j, &bandwidth_j)| {
-                let u1: f64 = rng.random();
-                let u2: f64 = rng.random();
+                let u1: f64 = rng.f64();
+                let u2: f64 = rng.f64();
 
                 // Box-Muller transform: generates standard normal variate
                 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * core::f64::consts::PI * u2).cos();
@@ -724,7 +722,7 @@ mod tests {
     fn test_multivariate_kde_sample_basic() {
         let samples = vec![vec![0.0, 0.0], vec![1.0, 1.0], vec![2.0, 2.0]];
         let kde = MultivariateKDE::new(samples).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Sample should have correct dimensionality
         let sample = kde.sample(&mut rng);
@@ -741,7 +739,7 @@ mod tests {
             vec![4.0, 4.0],
         ];
         let kde = MultivariateKDE::new(samples).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Samples should generally be in a reasonable range around the data
         for _ in 0..100 {
@@ -766,7 +764,7 @@ mod tests {
         // When KDE has only one sample, all samples should be centered around it
         let samples = vec![vec![5.0, 10.0]];
         let kde = MultivariateKDE::new(samples).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Generate many samples and check they cluster around (5.0, 10.0)
         let n_samples = 100;
@@ -803,7 +801,7 @@ mod tests {
             })
             .collect();
         let kde = MultivariateKDE::new(samples).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Sample should have correct dimensionality
         for _ in 0..50 {
@@ -823,7 +821,7 @@ mod tests {
         let data = vec![vec![0.0, 0.0], vec![0.0, 0.0], vec![0.0, 0.0]];
         let bandwidths = vec![0.1, 10.0]; // Small bandwidth in x, large in y
         let kde = MultivariateKDE::with_bandwidths(data, bandwidths).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Generate samples and check variance in each dimension
         let n_samples = 1000;
@@ -868,7 +866,7 @@ mod tests {
             vec![4.0, 4.0],
         ];
         let kde = MultivariateKDE::new(data).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Sample many points and verify the mean is near the center
         let n_samples = 500;
@@ -896,14 +894,12 @@ mod tests {
 
     #[test]
     fn test_multivariate_kde_sample_deterministic_with_seeded_rng() {
-        use rand::SeedableRng;
-
         let data = vec![vec![1.0, 2.0], vec![3.0, 4.0]];
         let kde = MultivariateKDE::new(data).unwrap();
 
         // Use a seeded RNG for reproducibility
-        let mut rng1 = rand::rngs::StdRng::seed_from_u64(42);
-        let mut rng2 = rand::rngs::StdRng::seed_from_u64(42);
+        let mut rng1 = fastrand::Rng::with_seed(42);
+        let mut rng2 = fastrand::Rng::with_seed(42);
 
         // Same seed should produce same samples
         let result1 = kde.sample(&mut rng1);

src/kde/univariate.rs

@@ -3,8 +3,6 @@
 //! This module provides a Gaussian kernel density estimator used by the TPE
 //! sampler to model probability distributions over good and bad trial regions.
 
-use rand::{Rng, RngExt};
-
 use crate::error::{Error, Result};
 
 /// A Gaussian kernel density estimator for continuous distributions.
@@ -26,7 +24,7 @@ use crate::error::{Error, Result};
 /// assert!(density > 0.0);
 ///
 /// // Sample from the estimated distribution
-/// let mut rng = rand::rng();
+/// let mut rng = fastrand::Rng::new();
 /// let sample = kde.sample(&mut rng);
 /// ```
 #[derive(Clone, Debug)]
@@ -132,15 +130,15 @@ impl KernelDensityEstimator {
     /// Sampling works by:
     /// 1. Uniformly selecting one of the kernel centers (samples)
     /// 2. Adding Gaussian noise with the bandwidth as standard deviation
-    pub(crate) fn sample<R: Rng>(&self, rng: &mut R) -> f64 {
+    pub(crate) fn sample(&self, rng: &mut fastrand::Rng) -> f64 {
         // Select a random sample to center the kernel on
-        let idx = rng.random_range(0..self.samples.len());
+        let idx = rng.usize(0..self.samples.len());
         let center = self.samples[idx];
 
         // Add Gaussian noise with bandwidth as standard deviation
         // Using Box-Muller transform for Gaussian sampling
-        let u1: f64 = rng.random();
-        let u2: f64 = rng.random();
+        let u1: f64 = rng.f64();
+        let u2: f64 = rng.f64();
 
         let z = (-2.0 * u1.ln()).sqrt() * (2.0 * core::f64::consts::PI * u2).cos();
         center + z * self.bandwidth
@@ -211,7 +209,7 @@ mod tests {
     fn test_kde_sample_in_reasonable_range() {
         let samples = vec![0.0, 1.0, 2.0, 3.0, 4.0];
         let kde = KernelDensityEstimator::new(samples).unwrap();
-        let mut rng = rand::rng();
+        let mut rng = fastrand::Rng::new();
 
         // Samples should generally be in a reasonable range around the data
         for _ in 0..100 {

src/lib.rs

@@ -19,8 +19,12 @@
 //! - **Grid Search** - Exhaustive search over a specified parameter grid
 //! - **Sobol (QMC)** - Quasi-random sampling for better space coverage (requires `sobol` feature)
 //! - **CMA-ES** - Covariance Matrix Adaptation Evolution Strategy for continuous optimization (requires `cma-es` feature)
+//! - **DE** - Differential Evolution for population-based global optimization
+//! - **GP** - Gaussian Process Bayesian optimization with Expected Improvement (requires `gp` feature)
 //! - **BOHB** - Bayesian Optimization + `HyperBand` for budget-aware TPE sampling
 //! - **NSGA-II** - Non-dominated Sorting Genetic Algorithm II for multi-objective optimization
+//! - **NSGA-III** - Reference-point-based NSGA for many-objective (3+) optimization
+//! - **MOEA/D** - Decomposition-based multi-objective with Tchebycheff, Weighted Sum, or PBI
 //! - **MOTPE** - Multi-Objective Tree-Parzen Estimator for Bayesian multi-objective optimization
 //!
 //! Additional features include:
@@ -190,6 +194,7 @@
 //! - `serde`: Enable `Serialize`/`Deserialize` on public types and `Study::save()`/`Study::load()`
 //! - `sobol`: Enable the Sobol quasi-random sampler for better space coverage
 //! - `cma-es`: Enable the CMA-ES sampler for continuous optimization
+//! - `gp`: Enable the Gaussian Process sampler for Bayesian optimization
 //! - `visualization`: Generate self-contained HTML reports with interactive Plotly.js charts
 //! - `tracing`: Emit structured log events via the [`tracing`](https://docs.rs/tracing) crate at key optimization points
 
@@ -219,7 +224,6 @@ macro_rules! trace_debug {
 
 mod distribution;
 mod error;
-#[cfg(feature = "fanova")]
 mod fanova;
 mod importance;
 mod kde;
@@ -228,15 +232,15 @@ mod param;
 pub mod parameter;
 pub mod pareto;
 pub mod pruner;
+mod rng_util;
 pub mod sampler;
+pub mod storage;
 mod study;
 mod trial;
 mod types;
-#[cfg(feature = "visualization")]
 mod visualization;
 
 pub use error::{Error, Result, TrialPruned};
-#[cfg(feature = "fanova")]
 pub use fanova::{FanovaConfig, FanovaResult};
 pub use multi_objective::{MultiObjectiveSampler, MultiObjectiveStudy, MultiObjectiveTrial};
 #[cfg(feature = "derive")]
@@ -253,19 +257,28 @@ pub use sampler::CompletedTrial;
 pub use sampler::bohb::BohbSampler;
 #[cfg(feature = "cma-es")]
 pub use sampler::cma_es::CmaEsSampler;
+pub use sampler::differential_evolution::{
+    DifferentialEvolutionSampler, DifferentialEvolutionStrategy,
+};
+#[cfg(feature = "gp")]
+pub use sampler::gp::GpSampler;
 pub use sampler::grid::GridSearchSampler;
+pub use sampler::moead::{Decomposition, MoeadSampler};
 pub use sampler::motpe::MotpeSampler;
 pub use sampler::nsga2::Nsga2Sampler;
+pub use sampler::nsga3::Nsga3Sampler;
 pub use sampler::random::RandomSampler;
 #[cfg(feature = "sobol")]
 pub use sampler::sobol::SobolSampler;
 pub use sampler::tpe::TpeSampler;
-pub use study::Study;
+#[cfg(feature = "journal")]
+pub use storage::JournalStorage;
+pub use storage::{MemoryStorage, Storage};
 #[cfg(feature = "serde")]
 pub use study::StudySnapshot;
+pub use study::{Study, StudyBuilder};
 pub use trial::{AttrValue, Trial};
 pub use types::{Direction, TrialState};
-#[cfg(feature = "visualization")]
 pub use visualization::generate_html_report;
 
 /// Convenient wildcard import for the most common types.
@@ -278,7 +291,6 @@ pub mod prelude {
     pub use optimizer_derive::Categorical as DeriveCategory;
 
     pub use crate::error::{Error, Result, TrialPruned};
-    #[cfg(feature = "fanova")]
     pub use crate::fanova::{FanovaConfig, FanovaResult};
     pub use crate::multi_objective::{MultiObjectiveStudy, MultiObjectiveTrial};
     pub use crate::param::ParamValue;
@@ -293,18 +305,27 @@ pub mod prelude {
     pub use crate::sampler::bohb::BohbSampler;
     #[cfg(feature = "cma-es")]
     pub use crate::sampler::cma_es::CmaEsSampler;
+    pub use crate::sampler::differential_evolution::{
+        DifferentialEvolutionSampler, DifferentialEvolutionStrategy,
+    };
+    #[cfg(feature = "gp")]
+    pub use crate::sampler::gp::GpSampler;
     pub use crate::sampler::grid::GridSearchSampler;
+    pub use crate::sampler::moead::{Decomposition, MoeadSampler};
     pub use crate::sampler::motpe::MotpeSampler;
     pub use crate::sampler::nsga2::Nsga2Sampler;
+    pub use crate::sampler::nsga3::Nsga3Sampler;
     pub use crate::sampler::random::RandomSampler;
     #[cfg(feature = "sobol")]
     pub use crate::sampler::sobol::SobolSampler;
     pub use crate::sampler::tpe::TpeSampler;
-    pub use crate::study::Study;
+    #[cfg(feature = "journal")]
+    pub use crate::storage::JournalStorage;
+    pub use crate::storage::{MemoryStorage, Storage};
     #[cfg(feature = "serde")]
     pub use crate::study::StudySnapshot;
+    pub use crate::study::{Study, StudyBuilder};
     pub use crate::trial::{AttrValue, Trial};
     pub use crate::types::Direction;
-    #[cfg(feature = "visualization")]
     pub use crate::visualization::generate_html_report;
 }

src/rng_util.rs

@@ -0,0 +1,5 @@
+/// Generate a random `f64` in the range `[low, high)`.
+#[inline]
+pub(crate) fn f64_range(rng: &mut fastrand::Rng, low: f64, high: f64) -> f64 {
+    low + rng.f64() * (high - low)
+}

src/sampler/cma_es.rs

@@ -25,11 +25,10 @@ use std::collections::HashMap;
 
 use nalgebra::{DMatrix, DVector};
 use parking_lot::Mutex;
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
 
 use crate::distribution::Distribution;
 use crate::param::ParamValue;
+use crate::rng_util;
 use crate::sampler::{CompletedTrial, Sampler};
 
 /// CMA-ES sampler for continuous optimization.
@@ -267,7 +266,7 @@ enum Phase {
 /// Top-level mutable state behind the `Mutex`.
 struct CmaEsState {
     /// The RNG used for sampling.
-    rng: StdRng,
+    rng: fastrand::Rng,
     /// User-provided initial sigma (None = auto).
     sigma0: Option<f64>,
     /// User-provided population size (None = auto).
@@ -290,7 +289,7 @@ struct CmaEsState {
 
 impl CmaEsState {
     fn new(sigma0: Option<f64>, user_lambda: Option<usize>, seed: Option<u64>) -> Self {
-        let rng = seed.map_or_else(rand::make_rng, StdRng::seed_from_u64);
+        let rng = seed.map_or_else(fastrand::Rng::new, fastrand::Rng::with_seed);
         Self {
             rng,
             sigma0,
@@ -421,7 +420,7 @@ impl CmaEsAlgorithm {
     /// Generate `lambda` candidate vectors from the current distribution.
     fn generate_candidates(
         &self,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
         dimensions: &[DimensionInfo],
     ) -> Vec<Candidate> {
         let n = self.constants.n;
@@ -439,7 +438,7 @@ impl CmaEsAlgorithm {
     /// Generate a single candidate from the current distribution.
     fn generate_single_candidate(
         &self,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
         dimensions: &[DimensionInfo],
         n: usize,
     ) -> Candidate {
@@ -452,7 +451,7 @@ impl CmaEsAlgorithm {
             if !dim.is_continuous
                 && let Distribution::Categorical(cat) = &dim.distribution
             {
-                categorical_values.insert(i, rng.random_range(0..cat.n_choices));
+                categorical_values.insert(i, rng.usize(0..cat.n_choices));
             }
         }
 
@@ -465,7 +464,7 @@ impl CmaEsAlgorithm {
     /// Sample a candidate vector with rejection sampling for bounds.
     fn sample_with_rejection(
         &self,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
         dimensions: &[DimensionInfo],
         n: usize,
     ) -> DVector<f64> {
@@ -679,36 +678,36 @@ fn internal_bounds(distribution: &Distribution) -> Option<(f64, f64)> {
 }
 
 /// Sample a value from the standard normal distribution using Box-Muller transform.
-fn sample_standard_normal(rng: &mut StdRng) -> f64 {
+fn sample_standard_normal(rng: &mut fastrand::Rng) -> f64 {
     // Box-Muller transform
-    let u1: f64 = rng.random_range(f64::EPSILON..=1.0);
-    let u2: f64 = rng.random_range(0.0_f64..=core::f64::consts::TAU);
+    let u1: f64 = rng_util::f64_range(rng, f64::EPSILON, 1.0);
+    let u2: f64 = rng_util::f64_range(rng, 0.0_f64, core::f64::consts::TAU);
     (-2.0 * u1.ln()).sqrt() * u2.cos()
 }
 
 /// Sample a categorical value randomly.
-fn sample_random_categorical(rng: &mut StdRng, distribution: &Distribution) -> ParamValue {
+fn sample_random_categorical(rng: &mut fastrand::Rng, distribution: &Distribution) -> ParamValue {
     match distribution {
-        Distribution::Categorical(d) => ParamValue::Categorical(rng.random_range(0..d.n_choices)),
+        Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
         _ => unreachable!("sample_random_categorical called with non-categorical distribution"),
     }
 }
 
 /// Sample a random value for any distribution (used during discovery phase).
 #[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
-fn sample_random(rng: &mut StdRng, distribution: &Distribution) -> ParamValue {
+fn sample_random(rng: &mut fastrand::Rng, distribution: &Distribution) -> ParamValue {
     match distribution {
         Distribution::Float(d) => {
             let value = if d.log_scale {
                 let log_low = d.low.ln();
                 let log_high = d.high.ln();
-                rng.random_range(log_low..=log_high).exp()
+                rng_util::f64_range(rng, log_low, log_high).exp()
             } else if let Some(step) = d.step {
                 let n_steps = ((d.high - d.low) / step).floor() as i64;
-                let k = rng.random_range(0..=n_steps);
+                let k = rng.i64(0..=n_steps);
                 d.low + (k as f64) * step
             } else {
-                rng.random_range(d.low..=d.high)
+                rng_util::f64_range(rng, d.low, d.high)
             };
             ParamValue::Float(value)
         }
@@ -716,18 +715,18 @@ fn sample_random(rng: &mut StdRng, distribution: &Distribution) -> ParamValue {
             let value = if d.log_scale {
                 let log_low = (d.low as f64).ln();
                 let log_high = (d.high as f64).ln();
-                let raw = rng.random_range(log_low..=log_high).exp().round() as i64;
+                let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
                 raw.clamp(d.low, d.high)
             } else if let Some(step) = d.step {
                 let n_steps = (d.high - d.low) / step;
-                let k = rng.random_range(0..=n_steps);
+                let k = rng.i64(0..=n_steps);
                 d.low + k * step
             } else {
-                rng.random_range(d.low..=d.high)
+                rng.i64(d.low..=d.high)
             };
             ParamValue::Int(value)
         }
-        Distribution::Categorical(d) => ParamValue::Categorical(rng.random_range(0..d.n_choices)),
+        Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
     }
 }
 
@@ -811,7 +810,7 @@ fn finalize_discovery(state: &mut CmaEsState) {
 
 /// Generate candidates that are purely categorical (no continuous dims).
 fn generate_pure_categorical_candidates(
-    rng: &mut StdRng,
+    rng: &mut fastrand::Rng,
     dimensions: &[DimensionInfo],
     lambda: usize,
 ) -> Vec<Candidate> {
@@ -820,7 +819,7 @@ fn generate_pure_categorical_candidates(
             let mut categorical_values = HashMap::new();
             for (i, dim) in dimensions.iter().enumerate() {
                 if let Distribution::Categorical(cat) = &dim.distribution {
-                    categorical_values.insert(i, rng.random_range(0..cat.n_choices));
+                    categorical_values.insert(i, rng.usize(0..cat.n_choices));
                 }
             }
             Candidate {
@@ -913,7 +912,7 @@ fn generate_overflow_candidate(state: &mut CmaEsState) -> Candidate {
         let mut categorical_values = HashMap::new();
         for (i, dim) in state.dimensions.iter().enumerate() {
             if let Distribution::Categorical(cat) = &dim.distribution {
-                categorical_values.insert(i, state.rng.random_range(0..cat.n_choices));
+                categorical_values.insert(i, state.rng.usize(0..cat.n_choices));
             }
         }
         return Candidate {

src/sampler/differential_evolution.rs

@@ -0,0 +1,995 @@
+//! Differential Evolution (DE) sampler.
+//!
+//! DE is a population-based metaheuristic that maintains a population of
+//! candidate solutions and creates new candidates by combining (mutating +
+//! crossing over) existing ones. It is competitive with CMA-ES on many
+//! problems and simpler to implement.
+//!
+//! Categorical parameters are sampled uniformly at random (not part of the
+//! DE vector). If all parameters are categorical, the sampler falls back to
+//! pure random sampling.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::sampler::differential_evolution::DifferentialEvolutionSampler;
+//! use optimizer::{Direction, Study};
+//!
+//! let sampler = DifferentialEvolutionSampler::with_seed(42);
+//! let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+//! ```
+
+use std::collections::HashMap;
+
+use parking_lot::Mutex;
+
+use crate::distribution::Distribution;
+use crate::param::ParamValue;
+use crate::rng_util;
+use crate::sampler::{CompletedTrial, Sampler};
+
+/// Differential Evolution mutation strategy.
+///
+/// Controls how mutant vectors are created from the current population.
+#[derive(Clone, Copy, Debug, Default)]
+pub enum DifferentialEvolutionStrategy {
+    /// DE/rand/1: `v = x_r1 + F * (x_r2 - x_r3)`
+    ///
+    /// The most robust strategy. Uses three random population members.
+    #[default]
+    Rand1,
+    /// DE/best/1: `v = x_best + F * (x_r1 - x_r2)`
+    ///
+    /// Greedier strategy that biases toward the current best solution.
+    Best1,
+    /// DE/current-to-best/1: `v = x_i + F * (x_best - x_i) + F * (x_r1 - x_r2)`
+    ///
+    /// Balances exploration and exploitation by blending the current
+    /// individual with the best.
+    CurrentToBest1,
+}
+
+/// Differential Evolution sampler for continuous global optimization.
+///
+/// Maintains a population of candidate solutions. New candidates are
+/// created by combining (mutating + crossing over) existing members.
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::differential_evolution::DifferentialEvolutionSampler;
+/// use optimizer::{Direction, Study};
+///
+/// // Default configuration
+/// let study: Study<f64> =
+///     Study::with_sampler(Direction::Minimize, DifferentialEvolutionSampler::new());
+///
+/// // With seed for reproducibility
+/// let study: Study<f64> = Study::with_sampler(
+///     Direction::Minimize,
+///     DifferentialEvolutionSampler::with_seed(42),
+/// );
+///
+/// // Custom configuration via builder
+/// use optimizer::sampler::differential_evolution::DifferentialEvolutionStrategy;
+/// let sampler = DifferentialEvolutionSampler::builder()
+///     .mutation_factor(0.8)
+///     .crossover_rate(0.9)
+///     .strategy(DifferentialEvolutionStrategy::Best1)
+///     .population_size(30)
+///     .seed(42)
+///     .build();
+/// let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+/// ```
+pub struct DifferentialEvolutionSampler {
+    state: Mutex<State>,
+}
+
+impl DifferentialEvolutionSampler {
+    /// Creates a new DE sampler with default settings and a random seed.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            state: Mutex::new(State::new(
+                None,
+                0.8,
+                0.9,
+                DifferentialEvolutionStrategy::Rand1,
+                None,
+            )),
+        }
+    }
+
+    /// Creates a new DE sampler with a fixed seed for reproducibility.
+    #[must_use]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            state: Mutex::new(State::new(
+                None,
+                0.8,
+                0.9,
+                DifferentialEvolutionStrategy::Rand1,
+                Some(seed),
+            )),
+        }
+    }
+
+    /// Creates a builder for configuring a `DifferentialEvolutionSampler`.
+    #[must_use]
+    pub fn builder() -> DifferentialEvolutionSamplerBuilder {
+        DifferentialEvolutionSamplerBuilder::new()
+    }
+}
+
+impl Default for DifferentialEvolutionSampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+/// Builder for configuring a [`DifferentialEvolutionSampler`].
+///
+/// All options have sensible defaults:
+/// - `population_size`: `max(10 * n_dims, 15)` (auto-computed from parameter count)
+/// - `mutation_factor` (F): 0.8
+/// - `crossover_rate` (CR): 0.9
+/// - `strategy`: `Rand1`
+/// - `seed`: random
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::differential_evolution::{
+///     DifferentialEvolutionSamplerBuilder, DifferentialEvolutionStrategy,
+/// };
+///
+/// let sampler = DifferentialEvolutionSamplerBuilder::new()
+///     .mutation_factor(0.5)
+///     .crossover_rate(0.7)
+///     .strategy(DifferentialEvolutionStrategy::CurrentToBest1)
+///     .population_size(20)
+///     .seed(42)
+///     .build();
+/// ```
+#[derive(Debug, Clone)]
+pub struct DifferentialEvolutionSamplerBuilder {
+    population_size: Option<usize>,
+    mutation_factor: f64,
+    crossover_rate: f64,
+    strategy: DifferentialEvolutionStrategy,
+    seed: Option<u64>,
+}
+
+impl Default for DifferentialEvolutionSamplerBuilder {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl DifferentialEvolutionSamplerBuilder {
+    /// Creates a new builder with default settings.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            population_size: None,
+            mutation_factor: 0.8,
+            crossover_rate: 0.9,
+            strategy: DifferentialEvolutionStrategy::Rand1,
+            seed: None,
+        }
+    }
+
+    /// Sets the population size.
+    ///
+    /// Number of candidate solutions maintained across generations.
+    /// Larger populations improve robustness but require more evaluations
+    /// per generation.
+    ///
+    /// Default: `max(10 * n_continuous_dims, 15)`.
+    #[must_use]
+    pub fn population_size(mut self, size: usize) -> Self {
+        self.population_size = Some(size);
+        self
+    }
+
+    /// Sets the mutation factor (F).
+    ///
+    /// Controls the amplification of differential variation.
+    /// Typical values are in `[0.5, 1.0]`. Higher values increase
+    /// exploration; lower values favor exploitation.
+    ///
+    /// Default: 0.8.
+    #[must_use]
+    pub fn mutation_factor(mut self, f: f64) -> Self {
+        self.mutation_factor = f;
+        self
+    }
+
+    /// Sets the crossover rate (CR).
+    ///
+    /// Probability of each dimension being taken from the mutant vector
+    /// rather than the parent. Typical values are in `[0.7, 1.0]`.
+    ///
+    /// Default: 0.9.
+    #[must_use]
+    pub fn crossover_rate(mut self, cr: f64) -> Self {
+        self.crossover_rate = cr;
+        self
+    }
+
+    /// Sets the mutation strategy.
+    ///
+    /// Default: [`DifferentialEvolutionStrategy::Rand1`].
+    #[must_use]
+    pub fn strategy(mut self, strategy: DifferentialEvolutionStrategy) -> Self {
+        self.strategy = strategy;
+        self
+    }
+
+    /// Sets the random seed for reproducibility.
+    #[must_use]
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.seed = Some(seed);
+        self
+    }
+
+    /// Builds the configured [`DifferentialEvolutionSampler`].
+    #[must_use]
+    pub fn build(self) -> DifferentialEvolutionSampler {
+        DifferentialEvolutionSampler {
+            state: Mutex::new(State::new(
+                self.population_size,
+                self.mutation_factor,
+                self.crossover_rate,
+                self.strategy,
+                self.seed,
+            )),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Internal types
+// ---------------------------------------------------------------------------
+
+/// Describes how a parameter dimension maps into the DE internal vector.
+#[derive(Clone, Debug)]
+struct DimensionInfo {
+    /// The distribution for this dimension (stored for decoding).
+    distribution: Distribution,
+    /// Whether this dimension participates in DE (Float/Int = true, Categorical = false).
+    is_continuous: bool,
+    /// Internal-space bounds for continuous dimensions: `(low, high)`.
+    /// For log-scale parameters these are in log-space.
+    bounds: Option<(f64, f64)>,
+}
+
+/// A candidate solution produced by mutation + crossover.
+#[derive(Clone, Debug)]
+struct Candidate {
+    /// Internal-space vector (only continuous dimensions).
+    x: Vec<f64>,
+    /// Values for categorical dimensions (index in `dimensions` -> categorical index).
+    categorical_values: HashMap<usize, usize>,
+    /// Index of the population member this candidate competes against.
+    target_idx: usize,
+}
+
+/// Tracks per-trial sampling progress.
+#[derive(Clone, Debug)]
+struct TrialProgress {
+    /// Index of the candidate assigned to this trial.
+    candidate_idx: usize,
+    /// Next dimension to return for this trial.
+    next_dim: usize,
+}
+
+/// Phase of the DE state machine.
+enum Phase {
+    /// Discovering the search space structure (first trial).
+    Discovery,
+    /// Active sampling and evolving.
+    Active,
+}
+
+/// Top-level mutable state behind the `Mutex`.
+struct State {
+    /// The RNG used for sampling.
+    rng: fastrand::Rng,
+    /// User-provided population size (None = auto).
+    user_population_size: Option<usize>,
+    /// Mutation factor (F).
+    mutation_factor: f64,
+    /// Crossover rate (CR).
+    crossover_rate: f64,
+    /// Mutation strategy.
+    strategy: DifferentialEvolutionStrategy,
+    /// Current phase.
+    phase: Phase,
+    /// Discovered dimension info (populated during discovery).
+    dimensions: Vec<DimensionInfo>,
+    /// Last `trial_id` seen during discovery.
+    discovery_trial_id: Option<u64>,
+
+    // --- Population state ---
+    /// Current population (internal-space vectors, continuous dims only).
+    population: Vec<Vec<f64>>,
+    /// Categorical values for each population member.
+    population_categorical: Vec<HashMap<usize, usize>>,
+    /// Objective values for the current population.
+    population_values: Vec<f64>,
+    /// Index of the best population member.
+    best_idx: usize,
+    /// Whether the initial population has been evaluated.
+    initialized: bool,
+    /// Effective population size (resolved after discovery).
+    population_size: usize,
+
+    // --- Current generation ---
+    /// Current generation's candidates.
+    candidates: Vec<Candidate>,
+    /// Mapping from `trial_id` to its progress.
+    trial_progress: HashMap<u64, TrialProgress>,
+    /// Number of candidates assigned so far in the current generation.
+    assigned_count: usize,
+    /// Trial IDs assigned in the current generation.
+    generation_trial_ids: Vec<u64>,
+}
+
+impl State {
+    fn new(
+        user_population_size: Option<usize>,
+        mutation_factor: f64,
+        crossover_rate: f64,
+        strategy: DifferentialEvolutionStrategy,
+        seed: Option<u64>,
+    ) -> Self {
+        let rng = seed.map_or_else(fastrand::Rng::new, fastrand::Rng::with_seed);
+        Self {
+            rng,
+            user_population_size,
+            mutation_factor,
+            crossover_rate,
+            strategy,
+            phase: Phase::Discovery,
+            dimensions: Vec::new(),
+            discovery_trial_id: None,
+            population: Vec::new(),
+            population_categorical: Vec::new(),
+            population_values: Vec::new(),
+            best_idx: 0,
+            initialized: false,
+            population_size: 0,
+            candidates: Vec::new(),
+            trial_progress: HashMap::new(),
+            assigned_count: 0,
+            generation_trial_ids: Vec::new(),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Helpers
+// ---------------------------------------------------------------------------
+
+/// Compute internal-space bounds for a distribution.
+#[allow(clippy::cast_precision_loss)]
+fn internal_bounds(distribution: &Distribution) -> Option<(f64, f64)> {
+    match distribution {
+        Distribution::Float(d) => {
+            if d.log_scale {
+                Some((d.low.ln(), d.high.ln()))
+            } else {
+                Some((d.low, d.high))
+            }
+        }
+        Distribution::Int(d) => {
+            if d.log_scale {
+                Some(((d.low as f64).ln(), (d.high as f64).ln()))
+            } else {
+                Some((d.low as f64, d.high as f64))
+            }
+        }
+        Distribution::Categorical(_) => None,
+    }
+}
+
+/// Convert an internal-space value back to a `ParamValue`.
+#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+fn from_internal(value: f64, distribution: &Distribution) -> ParamValue {
+    match distribution {
+        Distribution::Float(d) => {
+            let v = if d.log_scale { value.exp() } else { value };
+            let v = if let Some(step) = d.step {
+                let k = ((v - d.low) / step).round();
+                d.low + k * step
+            } else {
+                v
+            };
+            ParamValue::Float(v.clamp(d.low, d.high))
+        }
+        Distribution::Int(d) => {
+            let v = if d.log_scale { value.exp() } else { value };
+            let v = if let Some(step) = d.step {
+                let k = ((v - d.low as f64) / step as f64).round() as i64;
+                d.low + k * step
+            } else {
+                v.round() as i64
+            };
+            ParamValue::Int(v.clamp(d.low, d.high))
+        }
+        Distribution::Categorical(_) => {
+            unreachable!("from_internal should not be called for categorical distributions")
+        }
+    }
+}
+
+/// Sample a random value for any distribution.
+#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+fn sample_random(rng: &mut fastrand::Rng, distribution: &Distribution) -> ParamValue {
+    match distribution {
+        Distribution::Float(d) => {
+            let value = if d.log_scale {
+                let log_low = d.low.ln();
+                let log_high = d.high.ln();
+                rng_util::f64_range(rng, log_low, log_high).exp()
+            } else if let Some(step) = d.step {
+                let n_steps = ((d.high - d.low) / step).floor() as i64;
+                let k = rng.i64(0..=n_steps);
+                d.low + (k as f64) * step
+            } else {
+                rng_util::f64_range(rng, d.low, d.high)
+            };
+            ParamValue::Float(value)
+        }
+        Distribution::Int(d) => {
+            let value = if d.log_scale {
+                let log_low = (d.low as f64).ln();
+                let log_high = (d.high as f64).ln();
+                let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
+                raw.clamp(d.low, d.high)
+            } else if let Some(step) = d.step {
+                let n_steps = (d.high - d.low) / step;
+                let k = rng.i64(0..=n_steps);
+                d.low + k * step
+            } else {
+                rng.i64(d.low..=d.high)
+            };
+            ParamValue::Int(value)
+        }
+        Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
+    }
+}
+
+/// Sample a random value in internal space for a continuous dimension.
+fn sample_random_internal(rng: &mut fastrand::Rng, bounds: (f64, f64)) -> f64 {
+    rng_util::f64_range(rng, bounds.0, bounds.1)
+}
+
+/// Clamp a value to the given bounds.
+fn clamp_to_bounds(value: f64, bounds: Option<(f64, f64)>) -> f64 {
+    if let Some((lo, hi)) = bounds {
+        value.clamp(lo, hi)
+    } else {
+        value
+    }
+}
+
+// ---------------------------------------------------------------------------
+// DE algorithm
+// ---------------------------------------------------------------------------
+
+/// Select `count` distinct random indices from `0..n`, all different from `exclude`.
+fn select_random_indices(
+    rng: &mut fastrand::Rng,
+    n: usize,
+    count: usize,
+    exclude: &[usize],
+) -> Vec<usize> {
+    let mut selected = Vec::with_capacity(count);
+    while selected.len() < count {
+        let idx = rng.usize(0..n);
+        if !exclude.contains(&idx) && !selected.contains(&idx) {
+            selected.push(idx);
+        }
+    }
+    selected
+}
+
+/// Generate trial vectors (mutation + crossover) for the current population.
+fn generate_trial_vectors(state: &mut State) -> Vec<Candidate> {
+    let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+    let pop_size = state.population_size;
+
+    let mut candidates = Vec::with_capacity(pop_size);
+
+    for i in 0..pop_size {
+        // Mutation
+        let mutant = create_mutant_with_rng(state, i, n_continuous);
+
+        // Crossover (binomial)
+        let j_rand = state.rng.usize(0..n_continuous.max(1));
+        let trial_x: Vec<f64> = if n_continuous > 0 {
+            (0..n_continuous)
+                .map(|j| {
+                    let use_mutant = j == j_rand || state.rng.f64() < state.crossover_rate;
+                    let val = if use_mutant {
+                        mutant[j]
+                    } else {
+                        state.population[i][j]
+                    };
+                    // Clamp to bounds
+                    let dim_bounds = continuous_dim_bounds(&state.dimensions, j);
+                    clamp_to_bounds(val, dim_bounds)
+                })
+                .collect()
+        } else {
+            Vec::new()
+        };
+
+        // Categorical: randomly sample (DE doesn't optimize categoricals)
+        let mut categorical_values = HashMap::new();
+        for (dim_idx, dim) in state.dimensions.iter().enumerate() {
+            if !dim.is_continuous
+                && let Distribution::Categorical(cat) = &dim.distribution
+            {
+                categorical_values.insert(dim_idx, state.rng.usize(0..cat.n_choices));
+            }
+        }
+
+        candidates.push(Candidate {
+            x: trial_x,
+            categorical_values,
+            target_idx: i,
+        });
+    }
+
+    candidates
+}
+
+/// Create a mutant vector, consuming RNG from state.
+fn create_mutant_with_rng(state: &mut State, target_idx: usize, n_continuous: usize) -> Vec<f64> {
+    if n_continuous == 0 {
+        return Vec::new();
+    }
+
+    let pop = &state.population;
+    let best_idx = state.best_idx;
+    let f = state.mutation_factor;
+    let pop_size = state.population_size;
+
+    match state.strategy {
+        DifferentialEvolutionStrategy::Rand1 => {
+            let indices = select_random_indices(&mut state.rng, pop_size, 3, &[target_idx]);
+            let (r1, r2, r3) = (indices[0], indices[1], indices[2]);
+            (0..n_continuous)
+                .map(|j| pop[r1][j] + f * (pop[r2][j] - pop[r3][j]))
+                .collect()
+        }
+        DifferentialEvolutionStrategy::Best1 => {
+            let indices = select_random_indices(&mut state.rng, pop_size, 2, &[target_idx]);
+            let (r1, r2) = (indices[0], indices[1]);
+            (0..n_continuous)
+                .map(|j| pop[best_idx][j] + f * (pop[r1][j] - pop[r2][j]))
+                .collect()
+        }
+        DifferentialEvolutionStrategy::CurrentToBest1 => {
+            let indices = select_random_indices(&mut state.rng, pop_size, 2, &[target_idx]);
+            let (r1, r2) = (indices[0], indices[1]);
+            (0..n_continuous)
+                .map(|j| {
+                    pop[target_idx][j]
+                        + f * (pop[best_idx][j] - pop[target_idx][j])
+                        + f * (pop[r1][j] - pop[r2][j])
+                })
+                .collect()
+        }
+    }
+}
+
+/// Get the bounds for the j-th continuous dimension.
+fn continuous_dim_bounds(
+    dimensions: &[DimensionInfo],
+    continuous_idx: usize,
+) -> Option<(f64, f64)> {
+    let mut ci = 0;
+    for dim in dimensions {
+        if dim.is_continuous {
+            if ci == continuous_idx {
+                return dim.bounds;
+            }
+            ci += 1;
+        }
+    }
+    None
+}
+
+/// Generate the initial random population.
+fn generate_initial_population(state: &mut State) -> Vec<Candidate> {
+    let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+
+    let mut candidates = Vec::with_capacity(state.population_size);
+
+    for i in 0..state.population_size {
+        let x: Vec<f64> = if n_continuous > 0 {
+            let mut v = Vec::with_capacity(n_continuous);
+            for dim in &state.dimensions {
+                if dim.is_continuous {
+                    let val = if let Some(bounds) = dim.bounds {
+                        sample_random_internal(&mut state.rng, bounds)
+                    } else {
+                        0.0
+                    };
+                    v.push(val);
+                }
+            }
+            v
+        } else {
+            Vec::new()
+        };
+
+        let mut categorical_values = HashMap::new();
+        for (dim_idx, dim) in state.dimensions.iter().enumerate() {
+            if !dim.is_continuous
+                && let Distribution::Categorical(cat) = &dim.distribution
+            {
+                categorical_values.insert(dim_idx, state.rng.usize(0..cat.n_choices));
+            }
+        }
+
+        candidates.push(Candidate {
+            x,
+            categorical_values,
+            target_idx: i,
+        });
+    }
+
+    candidates
+}
+
+// ---------------------------------------------------------------------------
+// Sampler trait implementation
+// ---------------------------------------------------------------------------
+
+impl Sampler for DifferentialEvolutionSampler {
+    #[allow(clippy::cast_precision_loss)]
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[CompletedTrial],
+    ) -> ParamValue {
+        let mut state = self.state.lock();
+
+        match &state.phase {
+            Phase::Discovery => sample_discovery(&mut state, distribution, trial_id),
+            Phase::Active => sample_active(&mut state, distribution, trial_id, history),
+        }
+    }
+}
+
+/// Handle sampling during the discovery phase.
+fn sample_discovery(state: &mut State, distribution: &Distribution, trial_id: u64) -> ParamValue {
+    // Check if this is a new trial (discovery phase ended for previous trial)
+    if let Some(prev_id) = state.discovery_trial_id
+        && trial_id != prev_id
+    {
+        // First trial is done; we know the search space. Initialize DE.
+        finalize_discovery(state);
+        return sample_active(state, distribution, trial_id, &[]);
+    }
+
+    // Record this trial_id
+    state.discovery_trial_id = Some(trial_id);
+
+    // Record this dimension
+    let is_continuous = !matches!(distribution, Distribution::Categorical(_));
+    let bounds = internal_bounds(distribution);
+    state.dimensions.push(DimensionInfo {
+        distribution: distribution.clone(),
+        is_continuous,
+        bounds,
+    });
+
+    // Sample randomly for the discovery trial
+    sample_random(&mut state.rng, distribution)
+}
+
+/// Finalize discovery and transition to the active phase.
+#[allow(clippy::cast_precision_loss)]
+fn finalize_discovery(state: &mut State) {
+    let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+
+    // Resolve population size
+    state.population_size = state
+        .user_population_size
+        .unwrap_or_else(|| (10 * n_continuous).max(15));
+
+    // Ensure population size is at least 4 (DE needs distinct random indices)
+    state.population_size = state.population_size.max(4);
+
+    // Generate initial random population
+    state.candidates = generate_initial_population(state);
+    state.assigned_count = 0;
+    state.generation_trial_ids.clear();
+    state.trial_progress.clear();
+    state.phase = Phase::Active;
+}
+
+/// Handle sampling during the active phase.
+fn sample_active(
+    state: &mut State,
+    distribution: &Distribution,
+    trial_id: u64,
+    history: &[CompletedTrial],
+) -> ParamValue {
+    // Check if we need to process completed trials and start a new generation
+    maybe_update_generation(state, history);
+
+    // Assign a candidate to this trial if not yet done
+    if !state.trial_progress.contains_key(&trial_id) {
+        assign_candidate(state, trial_id);
+    }
+
+    let progress = state.trial_progress.get_mut(&trial_id).unwrap();
+    let dim_idx = progress.next_dim;
+    progress.next_dim += 1;
+
+    // Safety check
+    if dim_idx >= state.dimensions.len() {
+        return sample_random(&mut state.rng, distribution);
+    }
+
+    let candidate = &state.candidates[progress.candidate_idx];
+    let dim_info = &state.dimensions[dim_idx];
+
+    if dim_info.is_continuous {
+        // Map from overall dimension index to continuous index
+        let ci = state.dimensions[..dim_idx]
+            .iter()
+            .filter(|d| d.is_continuous)
+            .count();
+        if ci < candidate.x.len() {
+            from_internal(candidate.x[ci], &dim_info.distribution)
+        } else {
+            sample_random(&mut state.rng, distribution)
+        }
+    } else {
+        // Categorical: use pre-sampled value
+        if let Some(&cat_idx) = candidate.categorical_values.get(&dim_idx) {
+            ParamValue::Categorical(cat_idx)
+        } else {
+            sample_random(&mut state.rng, distribution)
+        }
+    }
+}
+
+/// Assign a candidate to a trial.
+fn assign_candidate(state: &mut State, trial_id: u64) {
+    let candidate_idx = if state.assigned_count < state.candidates.len() {
+        let idx = state.assigned_count;
+        state.assigned_count += 1;
+        idx
+    } else {
+        // Overflow: generate an extra random candidate
+        let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+        let x: Vec<f64> = (0..n_continuous)
+            .map(|j| {
+                let bounds = continuous_dim_bounds(&state.dimensions, j);
+                if let Some(b) = bounds {
+                    sample_random_internal(&mut state.rng, b)
+                } else {
+                    0.0
+                }
+            })
+            .collect();
+        let mut categorical_values = HashMap::new();
+        for (dim_idx, dim) in state.dimensions.iter().enumerate() {
+            if !dim.is_continuous
+                && let Distribution::Categorical(cat) = &dim.distribution
+            {
+                categorical_values.insert(dim_idx, state.rng.usize(0..cat.n_choices));
+            }
+        }
+        state.candidates.push(Candidate {
+            x,
+            categorical_values,
+            target_idx: 0, // overflow candidates don't compete
+        });
+        let idx = state.candidates.len() - 1;
+        state.assigned_count = state.candidates.len();
+        idx
+    };
+
+    state.trial_progress.insert(
+        trial_id,
+        TrialProgress {
+            candidate_idx,
+            next_dim: 0,
+        },
+    );
+    state.generation_trial_ids.push(trial_id);
+}
+
+/// Check if we should process completed trials and start a new generation.
+fn maybe_update_generation(state: &mut State, history: &[CompletedTrial]) {
+    let pop_size = state.population_size;
+
+    // Only update when at least pop_size candidates have been assigned
+    if state.generation_trial_ids.len() < pop_size {
+        return;
+    }
+
+    // Check if the first pop_size trial IDs are all completed
+    let trial_ids: Vec<u64> = state
+        .generation_trial_ids
+        .iter()
+        .take(pop_size)
+        .copied()
+        .collect();
+    let history_map: HashMap<u64, f64> = history.iter().map(|t| (t.id, t.value)).collect();
+
+    let all_completed = trial_ids.iter().all(|id| history_map.contains_key(id));
+    if !all_completed {
+        return;
+    }
+
+    let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+
+    if state.initialized {
+        // Subsequent generations: selection
+        perform_selection(state, &trial_ids, &history_map);
+    } else {
+        // First generation: store as initial population
+        initialize_population(state, &trial_ids, &history_map, n_continuous);
+    }
+
+    // Generate next generation's trial vectors
+    state.candidates = if state.initialized && n_continuous > 0 {
+        generate_trial_vectors(state)
+    } else {
+        generate_initial_population(state)
+    };
+    state.assigned_count = 0;
+    state.generation_trial_ids.clear();
+    state.trial_progress.clear();
+}
+
+/// Initialize the population from the first generation's results.
+fn initialize_population(
+    state: &mut State,
+    trial_ids: &[u64],
+    history_map: &HashMap<u64, f64>,
+    _n_continuous: usize,
+) {
+    state.population.clear();
+    state.population_categorical.clear();
+    state.population_values.clear();
+
+    let mut best_value = f64::INFINITY;
+    let mut best_idx = 0;
+
+    for (i, &trial_id) in trial_ids.iter().enumerate() {
+        let progress = &state.trial_progress[&trial_id];
+        let candidate = &state.candidates[progress.candidate_idx];
+        let value = history_map[&trial_id];
+
+        state.population.push(candidate.x.clone());
+        state
+            .population_categorical
+            .push(candidate.categorical_values.clone());
+        state.population_values.push(value);
+
+        if value < best_value {
+            best_value = value;
+            best_idx = i;
+        }
+    }
+
+    state.best_idx = best_idx;
+    state.initialized = true;
+}
+
+/// Perform DE selection: replace parent if trial vector is better.
+fn perform_selection(state: &mut State, trial_ids: &[u64], history_map: &HashMap<u64, f64>) {
+    for &trial_id in trial_ids {
+        let progress = &state.trial_progress[&trial_id];
+        let candidate = &state.candidates[progress.candidate_idx];
+        let trial_value = history_map[&trial_id];
+        let target_idx = candidate.target_idx;
+
+        if target_idx < state.population_size && trial_value <= state.population_values[target_idx]
+        {
+            state.population[target_idx] = candidate.x.clone();
+            state.population_categorical[target_idx] = candidate.categorical_values.clone();
+            state.population_values[target_idx] = trial_value;
+        }
+    }
+
+    // Update best index
+    let mut best_value = f64::INFINITY;
+    let mut best_idx = 0;
+    for (i, &val) in state.population_values.iter().enumerate() {
+        if val < best_value {
+            best_value = val;
+            best_idx = i;
+        }
+    }
+    state.best_idx = best_idx;
+}
+
+#[cfg(test)]
+#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+mod tests {
+    use super::*;
+    use crate::distribution::FloatDistribution;
+
+    #[test]
+    fn test_de_sampler_basic_float() {
+        let sampler = DifferentialEvolutionSampler::with_seed(42);
+        let dist = Distribution::Float(FloatDistribution {
+            low: -5.0,
+            high: 5.0,
+            log_scale: false,
+            step: None,
+        });
+
+        // Sample many values and check bounds
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Float(v) = value {
+                assert!(
+                    (-5.0..=5.0).contains(&v),
+                    "value {v} out of bounds at trial {i}"
+                );
+            } else {
+                panic!("Expected Float value");
+            }
+        }
+    }
+
+    #[test]
+    fn test_de_sampler_reproducibility() {
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        let sample_values = |seed: u64| {
+            let sampler = DifferentialEvolutionSampler::with_seed(seed);
+            (0..20)
+                .map(|i| sampler.sample(&dist, i, &[]))
+                .collect::<Vec<_>>()
+        };
+
+        let v1 = sample_values(42);
+        let v2 = sample_values(42);
+        assert_eq!(v1, v2, "same seed should produce same results");
+
+        let v3 = sample_values(99);
+        assert_ne!(v1, v3, "different seeds should produce different results");
+    }
+
+    #[test]
+    fn test_de_strategy_default() {
+        assert!(matches!(
+            DifferentialEvolutionStrategy::default(),
+            DifferentialEvolutionStrategy::Rand1
+        ));
+    }
+
+    #[test]
+    fn test_builder_defaults() {
+        let builder = DifferentialEvolutionSamplerBuilder::new();
+        assert!(builder.population_size.is_none());
+        assert!((builder.mutation_factor - 0.8).abs() < f64::EPSILON);
+        assert!((builder.crossover_rate - 0.9).abs() < f64::EPSILON);
+        assert!(matches!(
+            builder.strategy,
+            DifferentialEvolutionStrategy::Rand1
+        ));
+        assert!(builder.seed.is_none());
+    }
+}

src/sampler/genetic.rs

@@ -0,0 +1,572 @@
+//! Shared types and genetic operators for evolutionary multi-objective samplers.
+//!
+//! This module extracts common functionality used by NSGA-II, NSGA-III, and MOEA/D:
+//! candidate management, discovery/active phase logic, SBX crossover,
+//! polynomial mutation, and Das-Dennis reference point generation.
+
+use std::collections::HashMap;
+
+use crate::distribution::Distribution;
+use crate::multi_objective::MultiObjectiveTrial;
+use crate::param::ParamValue;
+use crate::rng_util;
+
+/// Describes a parameter dimension discovered during the first trial.
+#[derive(Clone, Debug)]
+pub(crate) struct DimensionInfo {
+    pub distribution: Distribution,
+}
+
+/// A candidate solution: one value per dimension.
+#[derive(Clone, Debug)]
+pub(crate) struct Candidate {
+    pub params: Vec<ParamValue>,
+}
+
+/// Tracks per-trial sampling progress (which candidate, which dimension next).
+#[derive(Clone, Debug)]
+pub(crate) struct TrialProgress {
+    pub candidate_idx: usize,
+    pub next_dim: usize,
+}
+
+/// Phase of an evolutionary sampler.
+pub(crate) enum Phase {
+    /// First trial reveals parameter dimensions.
+    Discovery,
+    /// Evolutionary optimisation.
+    Active,
+}
+
+/// Common state shared by all evolutionary multi-objective samplers.
+pub(crate) struct EvolutionaryState {
+    pub rng: fastrand::Rng,
+    pub phase: Phase,
+    pub dimensions: Vec<DimensionInfo>,
+    pub population_size: usize,
+    pub candidates: Vec<Candidate>,
+    pub trial_progress: HashMap<u64, TrialProgress>,
+    pub assigned_count: usize,
+    pub generation_trial_ids: Vec<u64>,
+    pub discovery_trial_id: Option<u64>,
+    pub generation: usize,
+}
+
+impl EvolutionaryState {
+    pub(crate) fn new(seed: Option<u64>) -> Self {
+        let rng = seed.map_or_else(fastrand::Rng::new, fastrand::Rng::with_seed);
+        Self {
+            rng,
+            phase: Phase::Discovery,
+            dimensions: Vec::new(),
+            population_size: 4,
+            candidates: Vec::new(),
+            trial_progress: HashMap::new(),
+            assigned_count: 0,
+            generation_trial_ids: Vec::new(),
+            discovery_trial_id: None,
+            generation: 0,
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Discovery phase helpers
+// ---------------------------------------------------------------------------
+
+/// Handle sampling during the discovery phase.
+///
+/// Returns `Some(value)` if the discovery phase handled the sample,
+/// or `None` if it transitioned to active phase and the caller should
+/// generate candidates and sample from them.
+pub(crate) fn sample_discovery(
+    evo: &mut EvolutionaryState,
+    distribution: &Distribution,
+    trial_id: u64,
+) -> Option<ParamValue> {
+    if let Some(prev_id) = evo.discovery_trial_id
+        && trial_id != prev_id
+    {
+        // A new trial arrived — transition to active phase
+        return None;
+    }
+
+    evo.discovery_trial_id = Some(trial_id);
+    evo.dimensions.push(DimensionInfo {
+        distribution: distribution.clone(),
+    });
+
+    Some(sample_random(&mut evo.rng, distribution))
+}
+
+/// Compute population size from dimensions and optional user override.
+#[allow(
+    clippy::cast_precision_loss,
+    clippy::cast_possible_truncation,
+    clippy::cast_sign_loss
+)]
+pub(crate) fn compute_population_size(
+    n_dims: usize,
+    user_pop_size: Option<usize>,
+    minimum: usize,
+) -> usize {
+    user_pop_size
+        .unwrap_or_else(|| (4.0 + 3.0 * (n_dims as f64).ln().max(0.0)).floor() as usize)
+        .max(minimum)
+}
+
+/// Transition from discovery to active phase.
+pub(crate) fn finalize_discovery(evo: &mut EvolutionaryState, user_pop_size: Option<usize>) {
+    evo.population_size = compute_population_size(evo.dimensions.len(), user_pop_size, 4);
+    evo.phase = Phase::Active;
+}
+
+/// Generate `population_size` random candidates.
+pub(crate) fn generate_random_candidates(evo: &mut EvolutionaryState) {
+    let pop = evo.population_size;
+    evo.candidates = (0..pop)
+        .map(|_| {
+            let params: Vec<ParamValue> = evo
+                .dimensions
+                .iter()
+                .map(|d| sample_random(&mut evo.rng, &d.distribution))
+                .collect();
+            Candidate { params }
+        })
+        .collect();
+    evo.assigned_count = 0;
+    evo.generation_trial_ids.clear();
+    evo.trial_progress.clear();
+}
+
+/// Assign a candidate to a trial and return the next dimension value.
+pub(crate) fn sample_from_candidate(evo: &mut EvolutionaryState, trial_id: u64) -> ParamValue {
+    if !evo.trial_progress.contains_key(&trial_id) {
+        let candidate_idx = if evo.assigned_count < evo.candidates.len() {
+            let idx = evo.assigned_count;
+            evo.assigned_count += 1;
+            idx
+        } else {
+            // Overflow: generate a random candidate
+            let params: Vec<ParamValue> = evo
+                .dimensions
+                .iter()
+                .map(|d| sample_random(&mut evo.rng, &d.distribution))
+                .collect();
+            evo.candidates.push(Candidate { params });
+            let idx = evo.candidates.len() - 1;
+            evo.assigned_count = evo.candidates.len();
+            idx
+        };
+
+        evo.trial_progress.insert(
+            trial_id,
+            TrialProgress {
+                candidate_idx,
+                next_dim: 0,
+            },
+        );
+        evo.generation_trial_ids.push(trial_id);
+    }
+
+    let progress = evo.trial_progress.get_mut(&trial_id).unwrap();
+    let dim_idx = progress.next_dim;
+    progress.next_dim += 1;
+
+    if dim_idx >= evo.dimensions.len() {
+        return sample_random(&mut evo.rng, &evo.dimensions.last().unwrap().distribution);
+    }
+
+    evo.candidates[progress.candidate_idx].params[dim_idx].clone()
+}
+
+/// Extract parameter values from a trial, ordered by dimension index.
+pub(crate) fn extract_trial_params(
+    trial: &MultiObjectiveTrial,
+    dimensions: &[DimensionInfo],
+    rng: &mut fastrand::Rng,
+) -> Vec<ParamValue> {
+    let mut param_pairs: Vec<_> = trial.params.iter().collect();
+    param_pairs.sort_by_key(|(id, _)| *id);
+
+    dimensions
+        .iter()
+        .enumerate()
+        .map(|(dim_idx, dim_info)| {
+            if dim_idx < param_pairs.len() {
+                param_pairs[dim_idx].1.clone()
+            } else {
+                sample_random(rng, &dim_info.distribution)
+            }
+        })
+        .collect()
+}
+
+/// Install new offspring as the next generation's candidates.
+pub(crate) fn advance_generation(evo: &mut EvolutionaryState, offspring: Vec<Candidate>) {
+    evo.candidates = offspring;
+    evo.assigned_count = 0;
+    evo.generation_trial_ids.clear();
+    evo.trial_progress.clear();
+    evo.generation += 1;
+}
+
+/// Check if the current generation is fully evaluated and return the
+/// evaluated trials if so.
+pub(crate) fn collect_evaluated_generation<'a>(
+    evo: &EvolutionaryState,
+    history: &'a [MultiObjectiveTrial],
+) -> Option<Vec<&'a MultiObjectiveTrial>> {
+    let pop_size = evo.population_size;
+
+    if evo.generation_trial_ids.len() < pop_size {
+        return None;
+    }
+
+    let gen_ids: Vec<u64> = evo
+        .generation_trial_ids
+        .iter()
+        .take(pop_size)
+        .copied()
+        .collect();
+    let history_map: HashMap<u64, &MultiObjectiveTrial> =
+        history.iter().map(|t| (t.id, t)).collect();
+
+    if !gen_ids.iter().all(|id| history_map.contains_key(id)) {
+        return None;
+    }
+
+    Some(
+        gen_ids
+            .iter()
+            .filter_map(|id| history_map.get(id).copied())
+            .collect(),
+    )
+}
+
+// ---------------------------------------------------------------------------
+// Genetic operators
+// ---------------------------------------------------------------------------
+
+/// SBX crossover for continuous params, uniform crossover for categorical.
+pub(crate) fn crossover(
+    rng: &mut fastrand::Rng,
+    parent1: &[ParamValue],
+    parent2: &[ParamValue],
+    dimensions: &[DimensionInfo],
+    crossover_prob: f64,
+    eta: f64,
+) -> (Vec<ParamValue>, Vec<ParamValue>) {
+    let n = parent1.len();
+    let mut child1 = parent1.to_vec();
+    let mut child2 = parent2.to_vec();
+
+    let u: f64 = rng_util::f64_range(rng, 0.0, 1.0);
+    if u > crossover_prob {
+        return (child1, child2);
+    }
+
+    for i in 0..n {
+        match (&parent1[i], &parent2[i], &dimensions[i].distribution) {
+            (ParamValue::Float(p1), ParamValue::Float(p2), Distribution::Float(d)) => {
+                if (p1 - p2).abs() < 1e-14 {
+                    continue;
+                }
+                let (c1, c2) = sbx_crossover_f64(rng, *p1, *p2, d.low, d.high, eta);
+                child1[i] = ParamValue::Float(c1);
+                child2[i] = ParamValue::Float(c2);
+            }
+            (ParamValue::Int(p1), ParamValue::Int(p2), Distribution::Int(d)) => {
+                if p1 == p2 {
+                    continue;
+                }
+                #[allow(clippy::cast_precision_loss)]
+                let (c1, c2) = sbx_crossover_f64(
+                    rng,
+                    *p1 as f64,
+                    *p2 as f64,
+                    d.low as f64,
+                    d.high as f64,
+                    eta,
+                );
+                #[allow(clippy::cast_possible_truncation)]
+                {
+                    child1[i] = ParamValue::Int((c1.round() as i64).clamp(d.low, d.high));
+                    child2[i] = ParamValue::Int((c2.round() as i64).clamp(d.low, d.high));
+                }
+            }
+            (ParamValue::Categorical(_), ParamValue::Categorical(_), _) => {
+                if rng_util::f64_range(rng, 0.0, 1.0) < 0.5 {
+                    core::mem::swap(&mut child1[i], &mut child2[i]);
+                }
+            }
+            _ => {}
+        }
+    }
+
+    (child1, child2)
+}
+
+/// SBX crossover for a single float dimension.
+pub(crate) fn sbx_crossover_f64(
+    rng: &mut fastrand::Rng,
+    p1: f64,
+    p2: f64,
+    low: f64,
+    high: f64,
+    eta: f64,
+) -> (f64, f64) {
+    let u: f64 = rng_util::f64_range(rng, 0.0, 1.0);
+
+    let beta = if u <= 0.5 {
+        (2.0 * u).powf(1.0 / (eta + 1.0))
+    } else {
+        (1.0 / (2.0 * (1.0 - u))).powf(1.0 / (eta + 1.0))
+    };
+
+    let c1 = 0.5 * ((1.0 + beta) * p1 + (1.0 - beta) * p2);
+    let c2 = 0.5 * ((1.0 - beta) * p1 + (1.0 + beta) * p2);
+
+    (c1.clamp(low, high), c2.clamp(low, high))
+}
+
+/// Polynomial mutation for each dimension.
+#[allow(clippy::cast_precision_loss)]
+pub(crate) fn mutate(
+    rng: &mut fastrand::Rng,
+    individual: &mut [ParamValue],
+    dimensions: &[DimensionInfo],
+    eta: f64,
+) {
+    let n = individual.len();
+    if n == 0 {
+        return;
+    }
+    let mutation_prob = 1.0 / n as f64;
+
+    for (i, value) in individual.iter_mut().enumerate() {
+        if rng_util::f64_range(rng, 0.0, 1.0) >= mutation_prob {
+            continue;
+        }
+
+        match (value, &dimensions[i].distribution) {
+            (v @ ParamValue::Float(_), Distribution::Float(d)) => {
+                let ParamValue::Float(x) = *v else {
+                    unreachable!();
+                };
+                let mutated = polynomial_mutation_f64(rng, x, d.low, d.high, eta);
+                *v = ParamValue::Float(mutated);
+            }
+            (v @ ParamValue::Int(_), Distribution::Int(d)) => {
+                let ParamValue::Int(x) = *v else {
+                    unreachable!();
+                };
+                #[allow(clippy::cast_possible_truncation)]
+                {
+                    let mutated =
+                        polynomial_mutation_f64(rng, x as f64, d.low as f64, d.high as f64, eta);
+                    *v = ParamValue::Int((mutated.round() as i64).clamp(d.low, d.high));
+                }
+            }
+            (v @ ParamValue::Categorical(_), Distribution::Categorical(d)) => {
+                *v = ParamValue::Categorical(rng.usize(0..d.n_choices));
+            }
+            _ => {}
+        }
+    }
+}
+
+/// Polynomial mutation for a single float value.
+pub(crate) fn polynomial_mutation_f64(
+    rng: &mut fastrand::Rng,
+    x: f64,
+    low: f64,
+    high: f64,
+    eta: f64,
+) -> f64 {
+    let u: f64 = rng_util::f64_range(rng, 0.0, 1.0);
+    let range = high - low;
+    if range <= 0.0 {
+        return x;
+    }
+
+    let delta1 = (x - low) / range;
+    let delta2 = (high - x) / range;
+
+    let delta_q = if u < 0.5 {
+        let xy = 1.0 - delta1;
+        let val = 2.0 * u + (1.0 - 2.0 * u) * xy.powf(eta + 1.0);
+        val.powf(1.0 / (eta + 1.0)) - 1.0
+    } else {
+        let xy = 1.0 - delta2;
+        let val = 2.0 * (1.0 - u) + 2.0 * (u - 0.5) * xy.powf(eta + 1.0);
+        1.0 - val.powf(1.0 / (eta + 1.0))
+    };
+
+    (x + delta_q * range).clamp(low, high)
+}
+
+/// Random sampling for a single distribution.
+#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+pub(crate) fn sample_random(rng: &mut fastrand::Rng, distribution: &Distribution) -> ParamValue {
+    match distribution {
+        Distribution::Float(d) => {
+            let value = if d.log_scale {
+                let log_low = d.low.ln();
+                let log_high = d.high.ln();
+                rng_util::f64_range(rng, log_low, log_high).exp()
+            } else if let Some(step) = d.step {
+                let n_steps = ((d.high - d.low) / step).floor() as i64;
+                let k = rng.i64(0..=n_steps);
+                d.low + (k as f64) * step
+            } else {
+                rng_util::f64_range(rng, d.low, d.high)
+            };
+            ParamValue::Float(value)
+        }
+        Distribution::Int(d) => {
+            let value = if d.log_scale {
+                let log_low = (d.low as f64).ln();
+                let log_high = (d.high as f64).ln();
+                let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
+                raw.clamp(d.low, d.high)
+            } else if let Some(step) = d.step {
+                let n_steps = (d.high - d.low) / step;
+                let k = rng.i64(0..=n_steps);
+                d.low + k * step
+            } else {
+                rng.i64(d.low..=d.high)
+            };
+            ParamValue::Int(value)
+        }
+        Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Das-Dennis reference point generation
+// ---------------------------------------------------------------------------
+
+/// Generate Das-Dennis (simplex-lattice) reference points.
+///
+/// Returns `C(H + M - 1, M - 1)` uniformly spaced points on the
+/// `M`-dimensional unit simplex, where `M = n_objectives` and
+/// `H = divisions`.
+pub(crate) fn das_dennis(n_objectives: usize, divisions: usize) -> Vec<Vec<f64>> {
+    let mut points = Vec::new();
+    let mut point = vec![0.0_f64; n_objectives];
+    das_dennis_recursive(
+        n_objectives,
+        divisions,
+        0,
+        divisions,
+        &mut point,
+        &mut points,
+    );
+    points
+}
+
+#[allow(clippy::cast_precision_loss)]
+fn das_dennis_recursive(
+    n_objectives: usize,
+    divisions: usize,
+    depth: usize,
+    remaining: usize,
+    current: &mut Vec<f64>,
+    result: &mut Vec<Vec<f64>>,
+) {
+    if depth == n_objectives - 1 {
+        current[depth] = remaining as f64 / divisions as f64;
+        result.push(current.clone());
+        return;
+    }
+
+    for i in 0..=remaining {
+        current[depth] = i as f64 / divisions as f64;
+        das_dennis_recursive(
+            n_objectives,
+            divisions,
+            depth + 1,
+            remaining - i,
+            current,
+            result,
+        );
+    }
+}
+
+/// Choose the number of divisions for Das-Dennis to get close to a target
+/// population size.
+///
+/// The number of reference points is `C(H + M - 1, M - 1)`. This function
+/// finds the smallest `H` such that the number of points >= `target_pop`.
+pub(crate) fn auto_divisions(n_objectives: usize, target_pop: usize) -> usize {
+    let m = n_objectives;
+    for h in 1..200 {
+        let n_points = n_combinations(h + m - 1, m - 1);
+        if n_points >= target_pop {
+            return h;
+        }
+    }
+    12
+}
+
+/// Compute `C(n, k)` = n! / (k! * (n-k)!).
+fn n_combinations(n: usize, k: usize) -> usize {
+    if k > n {
+        return 0;
+    }
+    let k = k.min(n - k);
+    let mut result: usize = 1;
+    for i in 0..k {
+        result = result.saturating_mul(n - i) / (i + 1);
+    }
+    result
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_das_dennis_2d() {
+        let points = das_dennis(2, 4);
+        // C(4+1, 1) = 5 points
+        assert_eq!(points.len(), 5);
+        for p in &points {
+            let sum: f64 = p.iter().sum();
+            assert!((sum - 1.0).abs() < 1e-10, "point {p:?} doesn't sum to 1");
+        }
+    }
+
+    #[test]
+    fn test_das_dennis_3d() {
+        let points = das_dennis(3, 4);
+        // C(4+2, 2) = 15 points
+        assert_eq!(points.len(), 15);
+        for p in &points {
+            let sum: f64 = p.iter().sum();
+            assert!((sum - 1.0).abs() < 1e-10);
+        }
+    }
+
+    #[test]
+    fn test_auto_divisions() {
+        // For 2 objectives targeting 10 points: H=9 gives C(10,1)=10
+        let h = auto_divisions(2, 10);
+        let n = n_combinations(h + 1, 1);
+        assert!(n >= 10);
+
+        // For 3 objectives targeting ~91 points: H=12 gives C(14,2)=91
+        let h3 = auto_divisions(3, 91);
+        let n3 = n_combinations(h3 + 2, 2);
+        assert!(n3 >= 91);
+    }
+
+    #[test]
+    fn test_n_combinations() {
+        assert_eq!(n_combinations(5, 2), 10);
+        assert_eq!(n_combinations(4, 0), 1);
+        assert_eq!(n_combinations(4, 4), 1);
+        assert_eq!(n_combinations(6, 3), 20);
+    }
+}

src/sampler/gp.rs

@@ -0,0 +1,848 @@
+//! Gaussian Process (GP) sampler with Expected Improvement acquisition.
+//!
+//! A classical Bayesian optimization sampler that uses a Gaussian Process
+//! surrogate model with a Matérn 5/2 kernel and Expected Improvement (EI)
+//! acquisition function. Best suited for small, expensive evaluations in
+//! low-dimensional continuous spaces (d ≤ 20).
+//!
+//! Categorical parameters are sampled uniformly at random (not part of
+//! the GP model). If all parameters are categorical, the sampler falls
+//! back to pure random sampling.
+//!
+//! Requires the `gp` feature flag.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::sampler::gp::GpSampler;
+//! use optimizer::{Direction, Study};
+//!
+//! let sampler = GpSampler::with_seed(42);
+//! let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+//! ```
+
+use std::collections::HashMap;
+
+use nalgebra::DMatrix;
+use parking_lot::Mutex;
+
+use crate::distribution::Distribution;
+use crate::param::ParamValue;
+use crate::rng_util;
+use crate::sampler::{CompletedTrial, Sampler};
+
+// ---------------------------------------------------------------------------
+// Public API
+// ---------------------------------------------------------------------------
+
+/// Gaussian Process sampler for Bayesian optimization.
+///
+/// Uses a GP surrogate with Matérn 5/2 kernel and Expected Improvement
+/// acquisition to guide sampling toward promising regions of the search
+/// space. Best suited for continuous (float/int) parameters in low
+/// dimensions (up to ~20).
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::gp::GpSampler;
+/// use optimizer::{Direction, Study};
+///
+/// // Default configuration
+/// let study: Study<f64> = Study::with_sampler(Direction::Minimize, GpSampler::new());
+///
+/// // With seed for reproducibility
+/// let study: Study<f64> = Study::with_sampler(Direction::Minimize, GpSampler::with_seed(42));
+///
+/// // Custom configuration via builder
+/// let sampler = GpSampler::builder()
+///     .n_startup_trials(15)
+///     .n_candidates(2000)
+///     .noise_variance(1e-4)
+///     .seed(42)
+///     .build();
+/// let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+/// ```
+pub struct GpSampler {
+    state: Mutex<GpState>,
+}
+
+impl GpSampler {
+    /// Creates a new GP sampler with a random seed.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            state: Mutex::new(GpState::new(None, None, None, None)),
+        }
+    }
+
+    /// Creates a new GP sampler with a fixed seed for reproducibility.
+    #[must_use]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            state: Mutex::new(GpState::new(None, None, None, Some(seed))),
+        }
+    }
+
+    /// Creates a builder for configuring a `GpSampler`.
+    #[must_use]
+    pub fn builder() -> GpSamplerBuilder {
+        GpSamplerBuilder::new()
+    }
+}
+
+impl Default for GpSampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+/// Builder for configuring a [`GpSampler`].
+///
+/// All options have sensible defaults:
+/// - `n_startup_trials`: 10
+/// - `n_candidates`: 1000
+/// - `noise_variance`: 1e-6
+/// - `seed`: random
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::gp::GpSamplerBuilder;
+///
+/// let sampler = GpSamplerBuilder::new()
+///     .n_startup_trials(15)
+///     .n_candidates(2000)
+///     .noise_variance(1e-4)
+///     .seed(42)
+///     .build();
+/// ```
+#[derive(Debug, Clone, Default)]
+pub struct GpSamplerBuilder {
+    n_startup_trials: Option<usize>,
+    n_candidates: Option<usize>,
+    noise_variance: Option<f64>,
+    seed: Option<u64>,
+}
+
+impl GpSamplerBuilder {
+    /// Creates a new builder with default settings.
+    #[must_use]
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    /// Sets the number of random trials before GP-guided sampling begins.
+    ///
+    /// Default: 10.
+    #[must_use]
+    pub fn n_startup_trials(mut self, n: usize) -> Self {
+        self.n_startup_trials = Some(n);
+        self
+    }
+
+    /// Sets the number of random candidate points for acquisition optimization.
+    ///
+    /// More candidates improve the quality of the acquisition maximum
+    /// at the cost of more GP predictions per trial.
+    ///
+    /// Default: 1000.
+    #[must_use]
+    pub fn n_candidates(mut self, n: usize) -> Self {
+        self.n_candidates = Some(n);
+        self
+    }
+
+    /// Sets the observation noise variance added to the kernel diagonal.
+    ///
+    /// Controls the assumed noise level. Larger values make the GP smoother.
+    ///
+    /// Default: 1e-6 (near-noiseless).
+    #[must_use]
+    pub fn noise_variance(mut self, v: f64) -> Self {
+        self.noise_variance = Some(v);
+        self
+    }
+
+    /// Sets the random seed for reproducibility.
+    #[must_use]
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.seed = Some(seed);
+        self
+    }
+
+    /// Builds the configured [`GpSampler`].
+    #[must_use]
+    pub fn build(self) -> GpSampler {
+        GpSampler {
+            state: Mutex::new(GpState::new(
+                self.n_startup_trials,
+                self.n_candidates,
+                self.noise_variance,
+                self.seed,
+            )),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Internal types
+// ---------------------------------------------------------------------------
+
+/// Default number of random startup trials before GP kicks in.
+const DEFAULT_N_STARTUP: usize = 10;
+/// Default number of candidate points for EI optimization.
+const DEFAULT_N_CANDIDATES: usize = 1000;
+/// Default observation noise variance.
+const DEFAULT_NOISE_VAR: f64 = 1e-6;
+
+/// Describes how a parameter dimension maps into the GP internal vector.
+#[derive(Clone, Debug)]
+struct DimensionInfo {
+    distribution: Distribution,
+    is_continuous: bool,
+    bounds: Option<(f64, f64)>,
+}
+
+/// Tracks per-trial sampling progress.
+#[derive(Clone, Debug)]
+struct TrialProgress {
+    /// The candidate values for each dimension.
+    values: Vec<ParamValue>,
+    /// Next dimension to return.
+    next_dim: usize,
+}
+
+/// Phase of the GP state machine.
+enum GpPhase {
+    /// Discovering the search space (first trial).
+    Discovery,
+    /// Steady-state sampling.
+    Active,
+}
+
+/// A fitted GP model ready for predictions.
+struct GpModel {
+    /// Cholesky factor L of K + σ²I.
+    cholesky: nalgebra::linalg::Cholesky<f64, nalgebra::Dyn>,
+    /// α = (K + σ²I)^{-1} y.
+    alpha: nalgebra::DVector<f64>,
+    /// Training inputs (each row is a data point, normalized to [0, 1]).
+    x_train: Vec<Vec<f64>>,
+    /// ARD lengthscales per dimension.
+    lengthscales: Vec<f64>,
+    /// Signal variance.
+    signal_var: f64,
+    /// Mean of original y values (for un-standardization, unused but kept for diagnostics).
+    _y_mean: f64,
+    /// Std dev of original y values (unused but kept for diagnostics).
+    _y_std: f64,
+    /// Best observed (standardized) y.
+    f_best: f64,
+}
+
+/// Top-level mutable state behind the `Mutex`.
+struct GpState {
+    rng: fastrand::Rng,
+    n_startup_trials: usize,
+    n_candidates: usize,
+    noise_variance: f64,
+    phase: GpPhase,
+    dimensions: Vec<DimensionInfo>,
+    trial_progress: HashMap<u64, TrialProgress>,
+    discovery_trial_id: Option<u64>,
+}
+
+impl GpState {
+    fn new(
+        n_startup: Option<usize>,
+        n_candidates: Option<usize>,
+        noise_var: Option<f64>,
+        seed: Option<u64>,
+    ) -> Self {
+        let rng = seed.map_or_else(fastrand::Rng::new, fastrand::Rng::with_seed);
+        Self {
+            rng,
+            n_startup_trials: n_startup.unwrap_or(DEFAULT_N_STARTUP),
+            n_candidates: n_candidates.unwrap_or(DEFAULT_N_CANDIDATES),
+            noise_variance: noise_var.unwrap_or(DEFAULT_NOISE_VAR),
+            phase: GpPhase::Discovery,
+            dimensions: Vec::new(),
+            trial_progress: HashMap::new(),
+            discovery_trial_id: None,
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Matérn 5/2 kernel
+// ---------------------------------------------------------------------------
+
+/// Matérn 5/2 kernel with ARD lengthscales.
+///
+/// `k(x1, x2) = σ² (1 + √5 r + 5/3 r²) exp(-√5 r)`
+/// where `r = sqrt(Σ ((x1_i - x2_i) / l_i)²)`
+fn matern52(x1: &[f64], x2: &[f64], lengthscales: &[f64], signal_var: f64) -> f64 {
+    let mut r_sq = 0.0;
+    for i in 0..x1.len() {
+        let diff = (x1[i] - x2[i]) / lengthscales[i];
+        r_sq += diff * diff;
+    }
+    let r = r_sq.sqrt();
+    let sqrt5_r = SQRT_5 * r;
+    signal_var * (1.0 + sqrt5_r + 5.0 / 3.0 * r_sq) * (-sqrt5_r).exp()
+}
+
+/// Build the kernel matrix `K + σ²I`.
+fn kernel_matrix(
+    x: &[Vec<f64>],
+    lengthscales: &[f64],
+    signal_var: f64,
+    noise_var: f64,
+) -> DMatrix<f64> {
+    let n = x.len();
+    DMatrix::from_fn(n, n, |i, j| {
+        let k = matern52(&x[i], &x[j], lengthscales, signal_var);
+        if i == j { k + noise_var } else { k }
+    })
+}
+
+/// Compute the kernel vector k(x*, X) for a test point.
+fn kernel_vector(
+    x_star: &[f64],
+    x_train: &[Vec<f64>],
+    lengthscales: &[f64],
+    signal_var: f64,
+) -> nalgebra::DVector<f64> {
+    nalgebra::DVector::from_fn(x_train.len(), |i, _| {
+        matern52(x_star, &x_train[i], lengthscales, signal_var)
+    })
+}
+
+/// Precomputed √5 constant.
+const SQRT_5: f64 = 2.236_213_562_373_095;
+
+// ---------------------------------------------------------------------------
+// GP fitting and prediction
+// ---------------------------------------------------------------------------
+
+/// Fit a GP model to the training data.
+///
+/// Returns `None` if fitting fails (e.g. Cholesky decomposition failure).
+#[allow(clippy::cast_precision_loss)]
+fn fit_gp(x_train: &[Vec<f64>], y_train: &[f64], noise_var: f64) -> Option<GpModel> {
+    let n = y_train.len();
+    if n == 0 {
+        return None;
+    }
+
+    // Standardize y
+    let y_mean = y_train.iter().sum::<f64>() / n as f64;
+    let y_var = if n > 1 {
+        y_train.iter().map(|&y| (y - y_mean).powi(2)).sum::<f64>() / (n - 1) as f64
+    } else {
+        1.0
+    };
+    let y_std = y_var.sqrt().max(1e-10);
+    let y_standardized: Vec<f64> = y_train.iter().map(|&y| (y - y_mean) / y_std).collect();
+
+    let f_best = y_standardized.iter().copied().fold(f64::INFINITY, f64::min);
+
+    // ARD lengthscales: per-dimension std dev of training X, clamped
+    let d = if x_train.is_empty() {
+        0
+    } else {
+        x_train[0].len()
+    };
+    let lengthscales: Vec<f64> = (0..d)
+        .map(|j| {
+            let vals: Vec<f64> = x_train.iter().map(|x| x[j]).collect();
+            let mean_j = vals.iter().sum::<f64>() / n as f64;
+            let var_j = vals.iter().map(|&v| (v - mean_j).powi(2)).sum::<f64>() / n as f64;
+            var_j.sqrt().max(0.01)
+        })
+        .collect();
+
+    // Signal variance = 1.0 (data is standardized)
+    let signal_var = 1.0;
+
+    let k = kernel_matrix(x_train, &lengthscales, signal_var, noise_var);
+    let cholesky = nalgebra::linalg::Cholesky::new(k)?;
+
+    // α = (K + σ²I)^{-1} y
+    let y_vec = nalgebra::DVector::from_column_slice(&y_standardized);
+    let alpha = cholesky.solve(&y_vec);
+
+    Some(GpModel {
+        cholesky,
+        alpha,
+        x_train: x_train.to_vec(),
+        lengthscales,
+        signal_var,
+        _y_mean: y_mean,
+        _y_std: y_std,
+        f_best,
+    })
+}
+
+/// Predict mean and standard deviation at a test point.
+fn predict(model: &GpModel, x: &[f64]) -> (f64, f64) {
+    let k_star = kernel_vector(x, &model.x_train, &model.lengthscales, model.signal_var);
+
+    // Mean: k*^T α
+    let mean = k_star.dot(&model.alpha);
+
+    // Variance: k(x*, x*) - k*^T (K + σ²I)^{-1} k*
+    let k_self = model.signal_var;
+    let v = model.cholesky.solve(&k_star);
+    let var = (k_self - k_star.dot(&v)).max(0.0);
+
+    (mean, var.sqrt())
+}
+
+// ---------------------------------------------------------------------------
+// Normal distribution helpers (Abramowitz-Stegun approximation)
+// ---------------------------------------------------------------------------
+
+/// Standard normal PDF.
+fn norm_pdf(x: f64) -> f64 {
+    const INV_SQRT_2PI: f64 = 0.398_942_280_401_432_7;
+    INV_SQRT_2PI * (-0.5 * x * x).exp()
+}
+
+/// Standard normal CDF (Abramowitz-Stegun rational approximation).
+fn norm_cdf(x: f64) -> f64 {
+    // Hart approximation (higher precision than basic A&S)
+    if x < -8.0 {
+        return 0.0;
+    }
+    if x > 8.0 {
+        return 1.0;
+    }
+
+    let abs_x = x.abs();
+    let t = 1.0 / (1.0 + 0.231_641_9 * abs_x);
+    let t2 = t * t;
+    let t3 = t2 * t;
+    let t4 = t3 * t;
+    let t5 = t4 * t;
+
+    let poly = 0.319_381_530 * t - 0.356_563_782 * t2 + 1.781_477_937 * t3 - 1.821_255_978 * t4
+        + 1.330_274_429 * t5;
+    let pdf = norm_pdf(abs_x);
+    let cdf = 1.0 - pdf * poly;
+
+    if x >= 0.0 { cdf } else { 1.0 - cdf }
+}
+
+// ---------------------------------------------------------------------------
+// Expected Improvement
+// ---------------------------------------------------------------------------
+
+/// Compute Expected Improvement at a point.
+///
+/// `EI(x) = (f_best - mean) Φ(z) + std φ(z)`
+/// where `z = (f_best - mean) / std`
+fn expected_improvement(mean: f64, std: f64, f_best: f64) -> f64 {
+    if std < 1e-12 {
+        return (f_best - mean).max(0.0);
+    }
+    let z = (f_best - mean) / std;
+    let improvement = (f_best - mean) * norm_cdf(z) + std * norm_pdf(z);
+    improvement.max(0.0)
+}
+
+// ---------------------------------------------------------------------------
+// Acquisition optimization
+// ---------------------------------------------------------------------------
+
+/// Find the point in [0, 1]^d that maximizes EI via multi-start random search.
+fn optimize_acquisition(
+    model: &GpModel,
+    n_dims: usize,
+    n_candidates: usize,
+    rng: &mut fastrand::Rng,
+) -> Vec<f64> {
+    let mut best_ei = f64::NEG_INFINITY;
+    let mut best_x = vec![0.5; n_dims];
+
+    for _ in 0..n_candidates {
+        let x: Vec<f64> = (0..n_dims)
+            .map(|_| rng_util::f64_range(rng, 0.0, 1.0))
+            .collect();
+        let (mean, std) = predict(model, &x);
+        let ei = expected_improvement(mean, std, model.f_best);
+        if ei > best_ei {
+            best_ei = ei;
+            best_x = x;
+        }
+    }
+
+    best_x
+}
+
+// ---------------------------------------------------------------------------
+// Data preprocessing helpers
+// ---------------------------------------------------------------------------
+
+/// Compute internal-space bounds for a distribution.
+#[allow(clippy::cast_precision_loss)]
+fn internal_bounds(distribution: &Distribution) -> Option<(f64, f64)> {
+    match distribution {
+        Distribution::Float(d) => {
+            if d.log_scale {
+                Some((d.low.ln(), d.high.ln()))
+            } else {
+                Some((d.low, d.high))
+            }
+        }
+        Distribution::Int(d) => {
+            if d.log_scale {
+                Some(((d.low as f64).ln(), (d.high as f64).ln()))
+            } else {
+                Some((d.low as f64, d.high as f64))
+            }
+        }
+        Distribution::Categorical(_) => None,
+    }
+}
+
+/// Convert a value from internal space to a `ParamValue` in original space.
+#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+fn from_internal(value: f64, distribution: &Distribution) -> ParamValue {
+    match distribution {
+        Distribution::Float(d) => {
+            let v = if d.log_scale { value.exp() } else { value };
+            let v = if let Some(step) = d.step {
+                let k = ((v - d.low) / step).round();
+                d.low + k * step
+            } else {
+                v
+            };
+            ParamValue::Float(v.clamp(d.low, d.high))
+        }
+        Distribution::Int(d) => {
+            let v = if d.log_scale { value.exp() } else { value };
+            let v = if let Some(step) = d.step {
+                let k = ((v - d.low as f64) / step as f64).round() as i64;
+                d.low + k * step
+            } else {
+                v.round() as i64
+            };
+            ParamValue::Int(v.clamp(d.low, d.high))
+        }
+        Distribution::Categorical(_) => {
+            unreachable!("from_internal should not be called for categorical distributions")
+        }
+    }
+}
+
+/// Convert an internal-space value to normalized [0, 1] using bounds.
+fn to_normalized(value: f64, lo: f64, hi: f64) -> f64 {
+    if (hi - lo).abs() < 1e-15 {
+        0.5
+    } else {
+        (value - lo) / (hi - lo)
+    }
+}
+
+/// Convert a normalized [0, 1] value back to internal space.
+fn from_normalized(value: f64, lo: f64, hi: f64) -> f64 {
+    lo + value * (hi - lo)
+}
+
+/// Convert a `ParamValue` to its internal-space representation.
+#[allow(clippy::cast_precision_loss)]
+fn to_internal(value: &ParamValue, distribution: &Distribution) -> f64 {
+    match (value, distribution) {
+        (ParamValue::Float(v), Distribution::Float(d)) => {
+            if d.log_scale {
+                v.ln()
+            } else {
+                *v
+            }
+        }
+        (ParamValue::Int(v), Distribution::Int(d)) => {
+            if d.log_scale {
+                (*v as f64).ln()
+            } else {
+                *v as f64
+            }
+        }
+        _ => 0.0,
+    }
+}
+
+/// Sample a random value for any distribution.
+#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+fn sample_random(rng: &mut fastrand::Rng, distribution: &Distribution) -> ParamValue {
+    match distribution {
+        Distribution::Float(d) => {
+            let value = if d.log_scale {
+                let log_low = d.low.ln();
+                let log_high = d.high.ln();
+                rng_util::f64_range(rng, log_low, log_high).exp()
+            } else if let Some(step) = d.step {
+                let n_steps = ((d.high - d.low) / step).floor() as i64;
+                let k = rng.i64(0..=n_steps);
+                d.low + (k as f64) * step
+            } else {
+                rng_util::f64_range(rng, d.low, d.high)
+            };
+            ParamValue::Float(value)
+        }
+        Distribution::Int(d) => {
+            let value = if d.log_scale {
+                let log_low = (d.low as f64).ln();
+                let log_high = (d.high as f64).ln();
+                let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
+                raw.clamp(d.low, d.high)
+            } else if let Some(step) = d.step {
+                let n_steps = (d.high - d.low) / step;
+                let k = rng.i64(0..=n_steps);
+                d.low + k * step
+            } else {
+                rng.i64(d.low..=d.high)
+            };
+            ParamValue::Int(value)
+        }
+        Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Extract training data from history
+// ---------------------------------------------------------------------------
+
+/// Maximum number of training points to use for the GP.
+/// Caps computational cost at O(`MAX_TRAIN_POINTS`^3) per trial.
+const MAX_TRAIN_POINTS: usize = 100;
+
+/// Establish a deterministic mapping from dimension index to `ParamId`
+/// using the first trial in history.
+///
+/// Matches dimensions to params by distribution equality, consuming
+/// matched params to correctly handle duplicate distributions.
+fn establish_param_mapping(
+    trial: &CompletedTrial,
+    dimensions: &[DimensionInfo],
+) -> Vec<Option<crate::parameter::ParamId>> {
+    use crate::parameter::ParamId;
+
+    let mut available: Vec<(ParamId, &Distribution)> =
+        trial.distributions.iter().map(|(id, d)| (*id, d)).collect();
+    // Sort for deterministic matching order
+    available.sort_by_key(|(id, _)| *id);
+
+    let mut mapping = Vec::with_capacity(dimensions.len());
+    for dim in dimensions {
+        let pos = available.iter().position(|(_, d)| **d == dim.distribution);
+        if let Some(pos) = pos {
+            mapping.push(Some(available.remove(pos).0));
+        } else {
+            mapping.push(None);
+        }
+    }
+    mapping
+}
+
+/// Build normalized training data from completed trials.
+///
+/// Returns `(x_train, y_train)` where x values are normalized to [0, 1]
+/// per dimension using the bounds from `dimensions`. Only continuous
+/// dimensions are included. Uses at most [`MAX_TRAIN_POINTS`] most recent
+/// trials.
+#[allow(clippy::cast_precision_loss)]
+fn build_training_data(
+    history: &[CompletedTrial],
+    dimensions: &[DimensionInfo],
+) -> (Vec<Vec<f64>>, Vec<f64>) {
+    if history.is_empty() {
+        return (Vec::new(), Vec::new());
+    }
+
+    // Use only the most recent trials to cap GP fitting cost
+    let start = history.len().saturating_sub(MAX_TRAIN_POINTS);
+    let recent = &history[start..];
+
+    // Establish dimension → ParamId mapping from the first trial
+    let param_mapping = establish_param_mapping(&recent[0], dimensions);
+
+    let continuous_indices: Vec<usize> = dimensions
+        .iter()
+        .enumerate()
+        .filter(|(_, d)| d.is_continuous)
+        .map(|(i, _)| i)
+        .collect();
+
+    let mut x_train = Vec::with_capacity(recent.len());
+    let mut y_train = Vec::with_capacity(recent.len());
+
+    for trial in recent {
+        let mut x_row = Vec::with_capacity(continuous_indices.len());
+        let mut valid = true;
+
+        for &dim_idx in &continuous_indices {
+            let dim_info = &dimensions[dim_idx];
+            if let Some(param_id) = param_mapping[dim_idx] {
+                if let Some(param_val) = trial.params.get(&param_id) {
+                    let internal = to_internal(param_val, &dim_info.distribution);
+                    let (lo, hi) = dim_info.bounds.unwrap_or((0.0, 1.0));
+                    x_row.push(to_normalized(internal, lo, hi));
+                } else {
+                    valid = false;
+                    break;
+                }
+            } else {
+                valid = false;
+                break;
+            }
+        }
+
+        if valid && x_row.len() == continuous_indices.len() {
+            x_train.push(x_row);
+            y_train.push(trial.value);
+        }
+    }
+
+    (x_train, y_train)
+}
+
+// ---------------------------------------------------------------------------
+// Sampler trait implementation
+// ---------------------------------------------------------------------------
+
+impl Sampler for GpSampler {
+    #[allow(clippy::cast_precision_loss, clippy::cast_possible_truncation)]
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[CompletedTrial],
+    ) -> ParamValue {
+        let mut state = self.state.lock();
+
+        match &state.phase {
+            GpPhase::Discovery => sample_discovery(&mut state, distribution, trial_id),
+            GpPhase::Active => sample_active(&mut state, distribution, trial_id, history),
+        }
+    }
+}
+
+/// Handle sampling during the discovery phase.
+fn sample_discovery(state: &mut GpState, distribution: &Distribution, trial_id: u64) -> ParamValue {
+    // A new trial_id means discovery is done
+    if let Some(prev_id) = state.discovery_trial_id
+        && trial_id != prev_id
+    {
+        finalize_discovery(state);
+        return sample_active(state, distribution, trial_id, &[]);
+    }
+
+    state.discovery_trial_id = Some(trial_id);
+
+    let is_continuous = !matches!(distribution, Distribution::Categorical(_));
+    let bounds = internal_bounds(distribution);
+    state.dimensions.push(DimensionInfo {
+        distribution: distribution.clone(),
+        is_continuous,
+        bounds,
+    });
+
+    sample_random(&mut state.rng, distribution)
+}
+
+/// Finalize discovery and transition to the active phase.
+fn finalize_discovery(state: &mut GpState) {
+    state.phase = GpPhase::Active;
+    state.trial_progress.clear();
+}
+
+/// Handle sampling during the active phase.
+fn sample_active(
+    state: &mut GpState,
+    distribution: &Distribution,
+    trial_id: u64,
+    history: &[CompletedTrial],
+) -> ParamValue {
+    // If this trial already has progress, return the next pre-computed value
+    if let Some(progress) = state.trial_progress.get_mut(&trial_id) {
+        let dim_idx = progress.next_dim;
+        progress.next_dim += 1;
+        if dim_idx < progress.values.len() {
+            return progress.values[dim_idx].clone();
+        }
+        // Extra dimension not seen during discovery
+        return sample_random(&mut state.rng, distribution);
+    }
+
+    // New trial: compute all dimension values at once
+    let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+    let use_gp = n_continuous > 0 && history.len() >= state.n_startup_trials;
+
+    let values = if use_gp {
+        compute_gp_candidate(state, history)
+    } else {
+        compute_random_candidate(state)
+    };
+
+    let first_value = values
+        .first()
+        .cloned()
+        .unwrap_or_else(|| sample_random(&mut state.rng, distribution));
+
+    state.trial_progress.insert(
+        trial_id,
+        TrialProgress {
+            values,
+            next_dim: 1,
+        },
+    );
+
+    first_value
+}
+
+/// Compute a candidate using the GP model.
+fn compute_gp_candidate(state: &mut GpState, history: &[CompletedTrial]) -> Vec<ParamValue> {
+    let (x_train, y_train) = build_training_data(history, &state.dimensions);
+
+    // Try to fit GP; fall back to random if it fails
+    let model = fit_gp(&x_train, &y_train, state.noise_variance);
+
+    let n_continuous = state.dimensions.iter().filter(|d| d.is_continuous).count();
+
+    let normalized_candidate = if let Some(ref model) = model {
+        optimize_acquisition(model, n_continuous, state.n_candidates, &mut state.rng)
+    } else {
+        // GP fitting failed; use random
+        (0..n_continuous)
+            .map(|_| rng_util::f64_range(&mut state.rng, 0.0, 1.0))
+            .collect()
+    };
+
+    // Convert normalized candidate back to parameter values
+    let mut values = Vec::with_capacity(state.dimensions.len());
+    let mut ci = 0; // continuous dimension index
+
+    for dim in &state.dimensions {
+        if dim.is_continuous {
+            let (lo, hi) = dim.bounds.unwrap_or((0.0, 1.0));
+            let internal_val = from_normalized(normalized_candidate[ci], lo, hi);
+            values.push(from_internal(internal_val, &dim.distribution));
+            ci += 1;
+        } else {
+            values.push(sample_random(&mut state.rng, &dim.distribution));
+        }
+    }
+
+    values
+}
+
+/// Compute a random candidate for all dimensions.
+fn compute_random_candidate(state: &mut GpState) -> Vec<ParamValue> {
+    state
+        .dimensions
+        .iter()
+        .map(|dim| sample_random(&mut state.rng, &dim.distribution))
+        .collect()
+}

src/sampler/mod.rs

@@ -3,9 +3,15 @@
 pub mod bohb;
 #[cfg(feature = "cma-es")]
 pub mod cma_es;
+pub mod differential_evolution;
+pub(crate) mod genetic;
+#[cfg(feature = "gp")]
+pub mod gp;
 pub mod grid;
+pub mod moead;
 pub mod motpe;
 pub mod nsga2;
+pub mod nsga3;
 pub mod random;
 #[cfg(feature = "sobol")]
 pub mod sobol;

src/sampler/moead.rs

@@ -0,0 +1,603 @@
+//! MOEA/D (Multi-Objective Evolutionary Algorithm based on Decomposition) sampler.
+//!
+//! Decomposes a multi-objective problem into scalar subproblems using
+//! weight vectors and solves them collaboratively. Supports Weighted Sum,
+//! Tchebycheff, and Penalty-based Boundary Intersection (PBI) scalarization.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::Direction;
+//! use optimizer::multi_objective::MultiObjectiveStudy;
+//! use optimizer::parameter::{FloatParam, Parameter};
+//! use optimizer::sampler::moead::MoeadSampler;
+//!
+//! let sampler = MoeadSampler::with_seed(42);
+//! let study =
+//!     MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+//!
+//! let x = FloatParam::new(0.0, 1.0);
+//! study
+//!     .optimize(100, |trial| {
+//!         let xv = x.suggest(trial)?;
+//!         Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+//!     })
+//!     .unwrap();
+//! ```
+
+use parking_lot::Mutex;
+
+use super::genetic::{
+    self, Candidate, EvolutionaryState, Phase, advance_generation, auto_divisions,
+    collect_evaluated_generation, crossover, das_dennis, extract_trial_params,
+    generate_random_candidates, mutate, sample_from_candidate, sample_random,
+};
+use crate::distribution::Distribution;
+use crate::multi_objective::MultiObjectiveTrial;
+use crate::param::ParamValue;
+use crate::types::Direction;
+
+/// Decomposition (scalarization) method for MOEA/D.
+#[derive(Debug, Clone, Default)]
+pub enum Decomposition {
+    /// Weighted sum: `sum(w_i * f_i)`.
+    WeightedSum,
+    /// Tchebycheff: `max(w_i * |f_i - z_i*|)`.
+    #[default]
+    Tchebycheff,
+    /// Penalty-based Boundary Intersection with parameter theta.
+    Pbi {
+        /// Penalty parameter controlling the balance between convergence
+        /// and diversity. Default: 5.0.
+        theta: f64,
+    },
+}
+
+/// MOEA/D sampler for multi-objective optimization.
+///
+/// Decomposes the multi-objective problem into scalar subproblems
+/// using weight vectors, solving them collaboratively via
+/// neighborhood-based mating and replacement.
+pub struct MoeadSampler {
+    state: Mutex<MoeadState>,
+}
+
+impl MoeadSampler {
+    /// Creates a new MOEA/D sampler with a random seed.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            state: Mutex::new(MoeadState::new(MoeadConfig::default(), None)),
+        }
+    }
+
+    /// Creates a new MOEA/D sampler with a fixed seed.
+    #[must_use]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            state: Mutex::new(MoeadState::new(MoeadConfig::default(), Some(seed))),
+        }
+    }
+
+    /// Creates a builder for configuring a `MoeadSampler`.
+    #[must_use]
+    pub fn builder() -> MoeadSamplerBuilder {
+        MoeadSamplerBuilder::default()
+    }
+}
+
+impl Default for MoeadSampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+/// Builder for [`MoeadSampler`].
+#[derive(Debug, Clone, Default)]
+pub struct MoeadSamplerBuilder {
+    population_size: Option<usize>,
+    neighborhood_size: Option<usize>,
+    decomposition: Decomposition,
+    crossover_prob: Option<f64>,
+    crossover_eta: Option<f64>,
+    mutation_eta: Option<f64>,
+    seed: Option<u64>,
+}
+
+impl MoeadSamplerBuilder {
+    /// Sets the population size. If unset, equals the number of
+    /// Das-Dennis weight vectors.
+    #[must_use]
+    pub fn population_size(mut self, size: usize) -> Self {
+        self.population_size = Some(size);
+        self
+    }
+
+    /// Sets the neighborhood size (T). Default: `min(20, pop_size)`.
+    #[must_use]
+    pub fn neighborhood_size(mut self, size: usize) -> Self {
+        self.neighborhood_size = Some(size);
+        self
+    }
+
+    /// Sets the decomposition method. Default: Tchebycheff.
+    #[must_use]
+    pub fn decomposition(mut self, decomp: Decomposition) -> Self {
+        self.decomposition = decomp;
+        self
+    }
+
+    /// Sets the crossover probability. Default: 1.0.
+    #[must_use]
+    pub fn crossover_prob(mut self, prob: f64) -> Self {
+        self.crossover_prob = Some(prob);
+        self
+    }
+
+    /// Sets the SBX distribution index. Default: 20.0.
+    #[must_use]
+    pub fn crossover_eta(mut self, eta: f64) -> Self {
+        self.crossover_eta = Some(eta);
+        self
+    }
+
+    /// Sets the polynomial mutation distribution index. Default: 20.0.
+    #[must_use]
+    pub fn mutation_eta(mut self, eta: f64) -> Self {
+        self.mutation_eta = Some(eta);
+        self
+    }
+
+    /// Sets the random seed for reproducibility.
+    #[must_use]
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.seed = Some(seed);
+        self
+    }
+
+    /// Builds the configured [`MoeadSampler`].
+    #[must_use]
+    pub fn build(self) -> MoeadSampler {
+        let config = MoeadConfig {
+            user_population_size: self.population_size,
+            neighborhood_size: self.neighborhood_size,
+            decomposition: self.decomposition,
+            crossover_prob: self.crossover_prob.unwrap_or(1.0),
+            crossover_eta: self.crossover_eta.unwrap_or(20.0),
+            mutation_eta: self.mutation_eta.unwrap_or(20.0),
+        };
+        MoeadSampler {
+            state: Mutex::new(MoeadState::new(config, self.seed)),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Internal types
+// ---------------------------------------------------------------------------
+
+#[derive(Debug, Clone)]
+struct MoeadConfig {
+    user_population_size: Option<usize>,
+    neighborhood_size: Option<usize>,
+    decomposition: Decomposition,
+    crossover_prob: f64,
+    crossover_eta: f64,
+    mutation_eta: f64,
+}
+
+impl Default for MoeadConfig {
+    fn default() -> Self {
+        Self {
+            user_population_size: None,
+            neighborhood_size: None,
+            decomposition: Decomposition::default(),
+            crossover_prob: 1.0,
+            crossover_eta: 20.0,
+            mutation_eta: 20.0,
+        }
+    }
+}
+
+struct MoeadState {
+    evo: EvolutionaryState,
+    config: MoeadConfig,
+    /// Weight vectors (Das-Dennis), one per subproblem.
+    weight_vectors: Vec<Vec<f64>>,
+    /// Neighborhoods: for each subproblem, indices of T nearest weight vectors.
+    neighborhoods: Vec<Vec<usize>>,
+    /// Ideal point z* (best per-objective in minimize-space).
+    ideal_point: Vec<f64>,
+    /// Current population's objective values in minimize-space (one per subproblem).
+    population_values: Vec<Vec<f64>>,
+    /// Current population's parameter vectors (one per subproblem).
+    population_params: Vec<Vec<ParamValue>>,
+    /// Whether the MOEA/D state has been initialized.
+    initialized: bool,
+}
+
+impl MoeadState {
+    fn new(config: MoeadConfig, seed: Option<u64>) -> Self {
+        Self {
+            evo: EvolutionaryState::new(seed),
+            config,
+            weight_vectors: Vec::new(),
+            neighborhoods: Vec::new(),
+            ideal_point: Vec::new(),
+            population_values: Vec::new(),
+            population_params: Vec::new(),
+            initialized: false,
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MultiObjectiveSampler implementation
+// ---------------------------------------------------------------------------
+
+impl crate::multi_objective::MultiObjectiveSampler for MoeadSampler {
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[MultiObjectiveTrial],
+        directions: &[Direction],
+    ) -> ParamValue {
+        let mut state = self.state.lock();
+
+        match &state.evo.phase {
+            Phase::Discovery => {
+                if let Some(value) =
+                    genetic::sample_discovery(&mut state.evo, distribution, trial_id)
+                {
+                    return value;
+                }
+                // Transitioned to active phase
+                initialize_moead(&mut state, directions);
+                generate_random_candidates(&mut state.evo);
+                sample_from_candidate(&mut state.evo, trial_id)
+            }
+            Phase::Active => {
+                maybe_generate_new_generation(&mut state, history, directions);
+                sample_from_candidate(&mut state.evo, trial_id)
+            }
+        }
+    }
+}
+
+/// Initialize MOEA/D: weight vectors, neighborhoods, ideal point.
+fn initialize_moead(state: &mut MoeadState, directions: &[Direction]) {
+    let n_obj = directions.len();
+
+    // Generate weight vectors
+    let divisions = auto_divisions(n_obj, state.config.user_population_size.unwrap_or(100));
+    state.weight_vectors = das_dennis(n_obj, divisions);
+
+    let pop_size = state
+        .config
+        .user_population_size
+        .unwrap_or(state.weight_vectors.len())
+        .max(4);
+
+    // Trim or pad weight vectors to match population size
+    state.weight_vectors.truncate(pop_size);
+    while state.weight_vectors.len() < pop_size {
+        // Duplicate random existing weight vectors
+        let idx = state.evo.rng.usize(0..state.weight_vectors.len());
+        let w = state.weight_vectors[idx].clone();
+        state.weight_vectors.push(w);
+    }
+
+    // Compute neighborhoods
+    let t = state
+        .config
+        .neighborhood_size
+        .unwrap_or_else(|| 20.min(pop_size));
+    let t = t.min(pop_size);
+    state.neighborhoods = compute_neighborhoods(&state.weight_vectors, t);
+
+    state.evo.population_size = pop_size;
+    state.evo.phase = Phase::Active;
+    state.ideal_point = vec![f64::INFINITY; n_obj];
+    state.initialized = true;
+}
+
+/// Compute T-nearest neighborhoods by Euclidean distance between weight vectors.
+fn compute_neighborhoods(weights: &[Vec<f64>], t: usize) -> Vec<Vec<usize>> {
+    let n = weights.len();
+    weights
+        .iter()
+        .map(|wi| {
+            let mut distances: Vec<(usize, f64)> = (0..n)
+                .map(|j| {
+                    let d: f64 = wi
+                        .iter()
+                        .zip(&weights[j])
+                        .map(|(&a, &b)| (a - b).powi(2))
+                        .sum::<f64>()
+                        .sqrt();
+                    (j, d)
+                })
+                .collect();
+            distances.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(core::cmp::Ordering::Equal));
+            distances.into_iter().take(t).map(|(idx, _)| idx).collect()
+        })
+        .collect()
+}
+
+/// Convert values to minimize-space.
+fn to_minimize_space(values: &[f64], directions: &[Direction]) -> Vec<f64> {
+    values
+        .iter()
+        .zip(directions)
+        .map(|(&v, d)| match d {
+            Direction::Minimize => v,
+            Direction::Maximize => -v,
+        })
+        .collect()
+}
+
+fn maybe_generate_new_generation(
+    state: &mut MoeadState,
+    history: &[MultiObjectiveTrial],
+    directions: &[Direction],
+) {
+    if state.evo.candidates.is_empty() {
+        generate_random_candidates(&mut state.evo);
+        return;
+    }
+
+    if let Some(evaluated) = collect_evaluated_generation(&state.evo, history) {
+        let offspring = moead_generate_offspring(state, &evaluated, directions);
+        advance_generation(&mut state.evo, offspring);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Scalarization functions
+// ---------------------------------------------------------------------------
+
+/// Weighted sum scalarization: `sum(w_i * f_i)`.
+fn scalarize_weighted_sum(values: &[f64], weight: &[f64]) -> f64 {
+    values.iter().zip(weight).map(|(&v, &w)| w * v).sum()
+}
+
+/// Tchebycheff scalarization: `max(w_i * |f_i - z_i*|)`.
+fn scalarize_tchebycheff(values: &[f64], weight: &[f64], ideal: &[f64]) -> f64 {
+    values
+        .iter()
+        .zip(weight)
+        .zip(ideal)
+        .map(|((&v, &w), &z)| {
+            let w = if w < 1e-6 { 1e-6 } else { w };
+            w * (v - z).abs()
+        })
+        .fold(f64::NEG_INFINITY, f64::max)
+}
+
+/// PBI scalarization: `d1 + theta * d2`.
+///
+/// d1 = projection onto weight direction, d2 = perpendicular distance.
+fn scalarize_pbi(values: &[f64], weight: &[f64], ideal: &[f64], theta: f64) -> f64 {
+    let n = values.len();
+
+    // Direction from ideal to the point
+    let diff: Vec<f64> = values.iter().zip(ideal).map(|(&v, &z)| v - z).collect();
+
+    // Normalize weight vector
+    let w_norm: f64 = weight.iter().map(|&w| w * w).sum::<f64>().sqrt();
+    if w_norm < 1e-30 {
+        return f64::INFINITY;
+    }
+    let w_unit: Vec<f64> = weight.iter().map(|&w| w / w_norm).collect();
+
+    // d1 = projection of diff onto weight direction
+    let d1: f64 = diff.iter().zip(&w_unit).map(|(&d, &w)| d * w).sum();
+
+    // d2 = perpendicular distance
+    let d2_sq: f64 = (0..n)
+        .map(|i| {
+            let proj = d1 * w_unit[i];
+            (diff[i] - proj).powi(2)
+        })
+        .sum::<f64>();
+
+    d1 + theta * d2_sq.sqrt()
+}
+
+/// Evaluate scalarization for a given decomposition method.
+fn scalarize(values: &[f64], weight: &[f64], ideal: &[f64], decomposition: &Decomposition) -> f64 {
+    match decomposition {
+        Decomposition::WeightedSum => scalarize_weighted_sum(values, weight),
+        Decomposition::Tchebycheff => scalarize_tchebycheff(values, weight, ideal),
+        Decomposition::Pbi { theta } => scalarize_pbi(values, weight, ideal, *theta),
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MOEA/D generation algorithm
+// ---------------------------------------------------------------------------
+
+fn moead_generate_offspring(
+    state: &mut MoeadState,
+    population: &[&MultiObjectiveTrial],
+    directions: &[Direction],
+) -> Vec<Candidate> {
+    let pop_size = state.evo.population_size;
+
+    if population.len() < 2 {
+        return (0..pop_size)
+            .map(|_| {
+                let params = state
+                    .evo
+                    .dimensions
+                    .iter()
+                    .map(|d| sample_random(&mut state.evo.rng, &d.distribution))
+                    .collect();
+                Candidate { params }
+            })
+            .collect();
+    }
+
+    // Extract current population parameters and objective values
+    let current_params: Vec<Vec<ParamValue>> = population
+        .iter()
+        .map(|t| extract_trial_params(t, &state.evo.dimensions, &mut state.evo.rng))
+        .collect();
+
+    let current_values: Vec<Vec<f64>> = population
+        .iter()
+        .map(|t| to_minimize_space(&t.values, directions))
+        .collect();
+
+    // Update ideal point
+    for vals in &current_values {
+        for (i, &v) in vals.iter().enumerate() {
+            if i < state.ideal_point.len() && v < state.ideal_point[i] {
+                state.ideal_point[i] = v;
+            }
+        }
+    }
+
+    // Assign each solution to its best subproblem via scalarization
+    // and select the best solution for each subproblem as its representative
+    let n_weights = state.weight_vectors.len();
+    let mut best_for_subproblem: Vec<usize> = Vec::with_capacity(n_weights);
+
+    for j in 0..n_weights {
+        let mut best_idx = 0;
+        let mut best_val = f64::INFINITY;
+        for (k, vals) in current_values.iter().enumerate() {
+            let s = scalarize(
+                vals,
+                &state.weight_vectors[j],
+                &state.ideal_point,
+                &state.config.decomposition,
+            );
+            if s < best_val {
+                best_val = s;
+                best_idx = k;
+            }
+        }
+        best_for_subproblem.push(best_idx);
+    }
+
+    // Store current population state
+    state.population_values = current_values;
+    state.population_params = current_params;
+
+    // Generate offspring: for each subproblem, mate from neighborhood
+    let mut offspring = Vec::with_capacity(pop_size);
+
+    for i in 0..pop_size.min(state.neighborhoods.len()) {
+        let neighborhood = &state.neighborhoods[i];
+
+        // Pick two parents from the neighborhood using subproblem assignments
+        let n1 = neighborhood[state.evo.rng.usize(0..neighborhood.len())];
+        let n2 = neighborhood[state.evo.rng.usize(0..neighborhood.len())];
+
+        let p1_idx = best_for_subproblem[n1 % best_for_subproblem.len()];
+        let p2_idx = best_for_subproblem[n2 % best_for_subproblem.len()];
+
+        let p1 = &state.population_params[p1_idx];
+        let p2 = &state.population_params[p2_idx];
+
+        let (mut child1, _child2) = crossover(
+            &mut state.evo.rng,
+            p1,
+            p2,
+            &state.evo.dimensions,
+            state.config.crossover_prob,
+            state.config.crossover_eta,
+        );
+
+        mutate(
+            &mut state.evo.rng,
+            &mut child1,
+            &state.evo.dimensions,
+            state.config.mutation_eta,
+        );
+
+        offspring.push(Candidate { params: child1 });
+    }
+
+    // If pop_size > neighborhoods, fill remaining with random neighborhood crossover
+    while offspring.len() < pop_size {
+        let i = state.evo.rng.usize(0..state.neighborhoods.len());
+        let neighborhood = &state.neighborhoods[i];
+        let n1 = neighborhood[state.evo.rng.usize(0..neighborhood.len())];
+        let n2 = neighborhood[state.evo.rng.usize(0..neighborhood.len())];
+
+        let p1_idx = best_for_subproblem[n1 % best_for_subproblem.len()];
+        let p2_idx = best_for_subproblem[n2 % best_for_subproblem.len()];
+
+        let (mut child1, _) = crossover(
+            &mut state.evo.rng,
+            &state.population_params[p1_idx],
+            &state.population_params[p2_idx],
+            &state.evo.dimensions,
+            state.config.crossover_prob,
+            state.config.crossover_eta,
+        );
+
+        mutate(
+            &mut state.evo.rng,
+            &mut child1,
+            &state.evo.dimensions,
+            state.config.mutation_eta,
+        );
+
+        offspring.push(Candidate { params: child1 });
+    }
+
+    offspring
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_scalarize_weighted_sum() {
+        let values = [1.0, 2.0, 3.0];
+        let weight = [0.5, 0.3, 0.2];
+        let result = scalarize_weighted_sum(&values, &weight);
+        assert!((result - (0.5 + 0.6 + 0.6)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_scalarize_tchebycheff() {
+        let values = [3.0, 2.0];
+        let weight = [0.5, 0.5];
+        let ideal = [1.0, 1.0];
+        let result = scalarize_tchebycheff(&values, &weight, &ideal);
+        // max(0.5 * |3-1|, 0.5 * |2-1|) = max(1.0, 0.5) = 1.0
+        assert!((result - 1.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_scalarize_pbi() {
+        let values = [2.0, 2.0];
+        let weight = [1.0, 1.0];
+        let ideal = [0.0, 0.0];
+        let result = scalarize_pbi(&values, &weight, &ideal, 5.0);
+        // d1 = projection of (2,2) onto (1/√2, 1/√2) = 2*√2
+        // d2 = 0 (point is on the weight direction)
+        let expected_d1 = 2.0 * (2.0_f64).sqrt();
+        assert!((result - expected_d1).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_compute_neighborhoods() {
+        let weights = vec![vec![1.0, 0.0], vec![0.5, 0.5], vec![0.0, 1.0]];
+        let neighborhoods = compute_neighborhoods(&weights, 2);
+        assert_eq!(neighborhoods.len(), 3);
+        // Each neighborhood should have 2 entries
+        for n in &neighborhoods {
+            assert_eq!(n.len(), 2);
+        }
+        // First weight [1,0] should be closest to itself and [0.5,0.5]
+        assert_eq!(neighborhoods[0][0], 0); // itself
+        assert_eq!(neighborhoods[0][1], 1); // nearest neighbor
+    }
+}

src/sampler/motpe.rs

@@ -40,15 +40,13 @@
 //! ```
 
 use parking_lot::Mutex;
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
 
 use crate::distribution::Distribution;
 use crate::kde::KernelDensityEstimator;
 use crate::multi_objective::{MultiObjectiveSampler, MultiObjectiveTrial};
 use crate::param::ParamValue;
-use crate::pareto;
 use crate::types::{Direction, TrialState};
+use crate::{pareto, rng_util};
 
 /// Multi-Objective TPE (MOTPE) sampler for multi-objective Bayesian optimization.
 ///
@@ -85,7 +83,7 @@ pub struct MotpeSampler {
     /// Optional fixed bandwidth for KDE. If None, uses Scott's rule.
     kde_bandwidth: Option<f64>,
     /// Thread-safe RNG for sampling.
-    rng: Mutex<StdRng>,
+    rng: Mutex<fastrand::Rng>,
 }
 
 impl MotpeSampler {
@@ -101,7 +99,7 @@ impl MotpeSampler {
             n_startup_trials: 11,
             n_ei_candidates: 24,
             kde_bandwidth: None,
-            rng: Mutex::new(rand::make_rng()),
+            rng: Mutex::new(fastrand::Rng::new()),
         }
     }
 
@@ -112,7 +110,7 @@ impl MotpeSampler {
             n_startup_trials: 11,
             n_ei_candidates: 24,
             kde_bandwidth: None,
-            rng: Mutex::new(StdRng::seed_from_u64(seed)),
+            rng: Mutex::new(fastrand::Rng::with_seed(seed)),
         }
     }
 
@@ -173,19 +171,19 @@ impl MotpeSampler {
         clippy::cast_precision_loss,
         clippy::unused_self
     )]
-    fn sample_uniform(distribution: &Distribution, rng: &mut StdRng) -> ParamValue {
+    fn sample_uniform(distribution: &Distribution, rng: &mut fastrand::Rng) -> ParamValue {
         match distribution {
             Distribution::Float(d) => {
                 let value = if d.log_scale {
                     let log_low = d.low.ln();
                     let log_high = d.high.ln();
-                    rng.random_range(log_low..=log_high).exp()
+                    rng_util::f64_range(rng, log_low, log_high).exp()
                 } else if let Some(step) = d.step {
                     let n_steps = ((d.high - d.low) / step).floor() as i64;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + (k as f64) * step
                 } else {
-                    rng.random_range(d.low..=d.high)
+                    rng_util::f64_range(rng, d.low, d.high)
                 };
                 ParamValue::Float(value)
             }
@@ -193,20 +191,18 @@ impl MotpeSampler {
                 let value = if d.log_scale {
                     let log_low = (d.low as f64).ln();
                     let log_high = (d.high as f64).ln();
-                    let raw = rng.random_range(log_low..=log_high).exp().round() as i64;
+                    let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
                     raw.clamp(d.low, d.high)
                 } else if let Some(step) = d.step {
                     let n_steps = (d.high - d.low) / step;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + k * step
                 } else {
-                    rng.random_range(d.low..=d.high)
+                    rng.i64(d.low..=d.high)
                 };
                 ParamValue::Int(value)
             }
-            Distribution::Categorical(d) => {
-                ParamValue::Categorical(rng.random_range(0..d.n_choices))
-            }
+            Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
         }
     }
 
@@ -220,7 +216,7 @@ impl MotpeSampler {
         step: Option<f64>,
         good_values: Vec<f64>,
         bad_values: Vec<f64>,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> f64 {
         // Transform to internal space (log space if needed)
         let (internal_low, internal_high, good_internal, bad_internal) = if log_scale {
@@ -245,7 +241,7 @@ impl MotpeSampler {
 
         // If KDE construction fails, fall back to uniform sampling
         let (Ok(l_kde), Ok(g_kde)) = (l_kde, g_kde) else {
-            return rng.random_range(low..=high);
+            return rng_util::f64_range(rng, low, high);
         };
 
         // Generate candidates from l(x) and select the one with best l(x)/g(x)
@@ -304,7 +300,7 @@ impl MotpeSampler {
         step: Option<i64>,
         good_values: &[i64],
         bad_values: &[i64],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> i64 {
         let good_floats: Vec<f64> = good_values.iter().map(|&v| v as f64).collect();
         let bad_floats: Vec<f64> = bad_values.iter().map(|&v| v as f64).collect();
@@ -336,7 +332,7 @@ impl MotpeSampler {
         n_choices: usize,
         good_indices: &[usize],
         bad_indices: &[usize],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> usize {
         let mut good_counts = vec![0usize; n_choices];
         let mut bad_counts = vec![0usize; n_choices];
@@ -365,7 +361,7 @@ impl MotpeSampler {
 
         // Sample proportionally to weights
         let total_weight: f64 = weights.iter().sum();
-        let threshold = rng.random::<f64>() * total_weight;
+        let threshold = rng.f64() * total_weight;
 
         let mut cumulative = 0.0;
         for (i, &w) in weights.iter().enumerate() {
@@ -589,8 +585,8 @@ impl MotpeSamplerBuilder {
     #[must_use]
     pub fn build(self) -> MotpeSampler {
         let rng = match self.seed {
-            Some(s) => StdRng::seed_from_u64(s),
-            None => rand::make_rng(),
+            Some(s) => fastrand::Rng::with_seed(s),
+            None => fastrand::Rng::new(),
         };
 
         MotpeSampler {

src/sampler/nsga2.rs

@@ -24,12 +24,13 @@
 //!     .unwrap();
 //! ```
 
-use std::collections::HashMap;
-
 use parking_lot::Mutex;
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
 
+use super::genetic::{
+    self, Candidate, EvolutionaryState, Phase, advance_generation, collect_evaluated_generation,
+    crossover, extract_trial_params, finalize_discovery, generate_random_candidates, mutate,
+    sample_from_candidate, sample_random,
+};
 use crate::distribution::Distribution;
 use crate::multi_objective::MultiObjectiveTrial;
 use crate::param::ParamValue;
@@ -158,62 +159,16 @@ impl Default for Nsga2Config {
     }
 }
 
-/// Describes a parameter dimension.
-#[derive(Clone, Debug)]
-struct DimensionInfo {
-    distribution: Distribution,
-}
-
-/// A candidate solution: one value per dimension.
-#[derive(Clone, Debug)]
-struct Candidate {
-    params: Vec<ParamValue>,
-}
-
-/// Tracks per-trial sampling progress.
-#[derive(Clone, Debug)]
-struct TrialProgress {
-    candidate_idx: usize,
-    next_dim: usize,
-}
-
-enum Phase {
-    /// First trial reveals parameter dimensions.
-    Discovery,
-    /// NSGA-II optimisation.
-    Active,
-}
-
 struct Nsga2State {
-    rng: StdRng,
+    evo: EvolutionaryState,
     config: Nsga2Config,
-    phase: Phase,
-    dimensions: Vec<DimensionInfo>,
-    population_size: usize,
-    candidates: Vec<Candidate>,
-    trial_progress: HashMap<u64, TrialProgress>,
-    assigned_count: usize,
-    generation_trial_ids: Vec<u64>,
-    discovery_trial_id: Option<u64>,
-    /// How many complete generations have been evaluated.
-    generation: usize,
 }
 
 impl Nsga2State {
     fn new(config: Nsga2Config, seed: Option<u64>) -> Self {
-        let rng = seed.map_or_else(rand::make_rng, StdRng::seed_from_u64);
         Self {
-            rng,
+            evo: EvolutionaryState::new(seed),
             config,
-            phase: Phase::Discovery,
-            dimensions: Vec::new(),
-            population_size: 4,
-            candidates: Vec::new(),
-            trial_progress: HashMap::new(),
-            assigned_count: 0,
-            generation_trial_ids: Vec::new(),
-            discovery_trial_id: None,
-            generation: 0,
         }
     }
 }
@@ -232,130 +187,27 @@ impl crate::multi_objective::MultiObjectiveSampler for Nsga2Sampler {
     ) -> ParamValue {
         let mut state = self.state.lock();
 
-        match &state.phase {
-            Phase::Discovery => sample_discovery(&mut state, distribution, trial_id),
-            Phase::Active => sample_active(&mut state, distribution, trial_id, history, directions),
+        match &state.evo.phase {
+            Phase::Discovery => {
+                if let Some(value) =
+                    genetic::sample_discovery(&mut state.evo, distribution, trial_id)
+                {
+                    return value;
+                }
+                // Transitioned to active phase
+                let user_pop = state.config.user_population_size;
+                finalize_discovery(&mut state.evo, user_pop);
+                generate_random_candidates(&mut state.evo);
+                sample_from_candidate(&mut state.evo, trial_id)
+            }
+            Phase::Active => {
+                maybe_generate_new_generation(&mut state, history, directions);
+                sample_from_candidate(&mut state.evo, trial_id)
+            }
         }
     }
 }
 
-/// Handle sampling during the discovery phase.
-fn sample_discovery(
-    state: &mut Nsga2State,
-    distribution: &Distribution,
-    trial_id: u64,
-) -> ParamValue {
-    if let Some(prev_id) = state.discovery_trial_id
-        && trial_id != prev_id
-    {
-        finalize_discovery(state);
-        // Assign this trial a random candidate (no history yet)
-        generate_random_candidates(state);
-        return sample_from_candidate(state, trial_id);
-    }
-
-    state.discovery_trial_id = Some(trial_id);
-    state.dimensions.push(DimensionInfo {
-        distribution: distribution.clone(),
-    });
-
-    sample_random(&mut state.rng, distribution)
-}
-
-/// Transition from discovery to active phase.
-#[allow(
-    clippy::cast_precision_loss,
-    clippy::cast_possible_truncation,
-    clippy::cast_sign_loss
-)]
-fn finalize_discovery(state: &mut Nsga2State) {
-    let n = state.dimensions.len();
-    state.population_size = state
-        .config
-        .user_population_size
-        .unwrap_or_else(|| (4.0 + 3.0 * (n as f64).ln().max(0.0)).floor() as usize)
-        .max(4);
-    state.phase = Phase::Active;
-}
-
-/// Generate `population_size` random candidates.
-fn generate_random_candidates(state: &mut Nsga2State) {
-    let pop = state.population_size;
-    state.candidates = (0..pop)
-        .map(|_| {
-            let params: Vec<ParamValue> = state
-                .dimensions
-                .iter()
-                .map(|d| sample_random(&mut state.rng, &d.distribution))
-                .collect();
-            Candidate { params }
-        })
-        .collect();
-    state.assigned_count = 0;
-    state.generation_trial_ids.clear();
-    state.trial_progress.clear();
-}
-
-/// Active-phase sampling.
-fn sample_active(
-    state: &mut Nsga2State,
-    _distribution: &Distribution,
-    trial_id: u64,
-    history: &[MultiObjectiveTrial],
-    directions: &[Direction],
-) -> ParamValue {
-    // Check if we need to generate a new generation
-    maybe_generate_new_generation(state, history, directions);
-
-    sample_from_candidate(state, trial_id)
-}
-
-/// Assign a candidate to a trial and return the next dimension value.
-fn sample_from_candidate(state: &mut Nsga2State, trial_id: u64) -> ParamValue {
-    // Assign candidate if not yet done
-    if !state.trial_progress.contains_key(&trial_id) {
-        let candidate_idx = if state.assigned_count < state.candidates.len() {
-            let idx = state.assigned_count;
-            state.assigned_count += 1;
-            idx
-        } else {
-            // Overflow: generate a random candidate
-            let params: Vec<ParamValue> = state
-                .dimensions
-                .iter()
-                .map(|d| sample_random(&mut state.rng, &d.distribution))
-                .collect();
-            state.candidates.push(Candidate { params });
-            let idx = state.candidates.len() - 1;
-            state.assigned_count = state.candidates.len();
-            idx
-        };
-
-        state.trial_progress.insert(
-            trial_id,
-            TrialProgress {
-                candidate_idx,
-                next_dim: 0,
-            },
-        );
-        state.generation_trial_ids.push(trial_id);
-    }
-
-    let progress = state.trial_progress.get_mut(&trial_id).unwrap();
-    let dim_idx = progress.next_dim;
-    progress.next_dim += 1;
-
-    if dim_idx >= state.dimensions.len() {
-        // Extra dimension: sample randomly
-        return sample_random(
-            &mut state.rng,
-            &state.dimensions.last().unwrap().distribution,
-        );
-    }
-
-    state.candidates[progress.candidate_idx].params[dim_idx].clone()
-}
-
 /// Check if all candidates in the current generation have been evaluated;
 /// if so, run NSGA-II selection and generate offspring.
 fn maybe_generate_new_generation(
@@ -363,45 +215,15 @@ fn maybe_generate_new_generation(
     history: &[MultiObjectiveTrial],
     directions: &[Direction],
 ) {
-    let pop_size = state.population_size;
-
-    // Need at least pop_size assigned trials
-    if state.generation_trial_ids.len() < pop_size {
-        // Not enough candidates assigned yet — check if we need initial candidates
-        if state.candidates.is_empty() {
-            generate_random_candidates(state);
-        }
+    if state.evo.candidates.is_empty() {
+        generate_random_candidates(&mut state.evo);
         return;
     }
 
-    // Check if the first pop_size trials are completed
-    let gen_ids: Vec<u64> = state
-        .generation_trial_ids
-        .iter()
-        .take(pop_size)
-        .copied()
-        .collect();
-    let history_map: HashMap<u64, &MultiObjectiveTrial> =
-        history.iter().map(|t| (t.id, t)).collect();
-
-    let all_completed = gen_ids.iter().all(|id| history_map.contains_key(id));
-    if !all_completed {
-        return;
+    if let Some(evaluated) = collect_evaluated_generation(&state.evo, history) {
+        let offspring = nsga2_generate_offspring(state, &evaluated, directions);
+        advance_generation(&mut state.evo, offspring);
     }
-
-    // Collect the evaluated population
-    let evaluated: Vec<&MultiObjectiveTrial> = gen_ids
-        .iter()
-        .filter_map(|id| history_map.get(id).copied())
-        .collect();
-
-    // Run NSGA-II to produce offspring
-    let offspring = nsga2_generate_offspring(state, &evaluated, directions);
-    state.candidates = offspring;
-    state.assigned_count = 0;
-    state.generation_trial_ids.clear();
-    state.trial_progress.clear();
-    state.generation += 1;
 }
 
 // ---------------------------------------------------------------------------
@@ -415,7 +237,7 @@ fn nsga2_select(
     population: &[&MultiObjectiveTrial],
     directions: &[Direction],
 ) -> (Vec<Vec<ParamValue>>, Vec<usize>, Vec<f64>) {
-    let pop_size = state.population_size;
+    let pop_size = state.evo.population_size;
 
     let values: Vec<Vec<f64>> = population.iter().map(|t| t.values.clone()).collect();
     let constraints: Vec<Vec<f64>> = population.iter().map(|t| t.constraints.clone()).collect();
@@ -457,13 +279,14 @@ fn nsga2_select(
     }
 
     while selected.len() < pop_size {
-        selected.push(state.rng.random_range(0..n));
+        selected.push(state.evo.rng.usize(0..n));
     }
 
-    // Extract parent parameter vectors ordered by dimension
     let parents: Vec<Vec<ParamValue>> = selected
         .iter()
-        .map(|&idx| extract_trial_params(population[idx], &state.dimensions, &mut state.rng))
+        .map(|&idx| {
+            extract_trial_params(population[idx], &state.evo.dimensions, &mut state.evo.rng)
+        })
         .collect();
 
     let sel_rank: Vec<usize> = selected.iter().map(|&i| rank[i]).collect();
@@ -472,43 +295,22 @@ fn nsga2_select(
     (parents, sel_rank, sel_crowding)
 }
 
-/// Extract parameter values from a trial, ordered by dimension index.
-fn extract_trial_params(
-    trial: &MultiObjectiveTrial,
-    dimensions: &[DimensionInfo],
-    rng: &mut StdRng,
-) -> Vec<ParamValue> {
-    let mut param_pairs: Vec<_> = trial.params.iter().collect();
-    param_pairs.sort_by_key(|(id, _)| *id);
-
-    dimensions
-        .iter()
-        .enumerate()
-        .map(|(dim_idx, dim_info)| {
-            if dim_idx < param_pairs.len() {
-                param_pairs[dim_idx].1.clone()
-            } else {
-                sample_random(rng, &dim_info.distribution)
-            }
-        })
-        .collect()
-}
-
 /// Runs NSGA-II selection and generates offspring candidates.
 fn nsga2_generate_offspring(
     state: &mut Nsga2State,
     population: &[&MultiObjectiveTrial],
     directions: &[Direction],
 ) -> Vec<Candidate> {
-    let pop_size = state.population_size;
+    let pop_size = state.evo.population_size;
 
     if population.len() < 2 {
         return (0..pop_size)
             .map(|_| {
                 let params = state
+                    .evo
                     .dimensions
                     .iter()
-                    .map(|d| sample_random(&mut state.rng, &d.distribution))
+                    .map(|d| sample_random(&mut state.evo.rng, &d.distribution))
                     .collect();
                 Candidate { params }
             })
@@ -519,28 +321,28 @@ fn nsga2_generate_offspring(
 
     let mut offspring = Vec::with_capacity(pop_size);
     while offspring.len() < pop_size {
-        let p1 = tournament_select(&mut state.rng, &sel_rank, &sel_crowding, parents.len());
-        let p2 = tournament_select(&mut state.rng, &sel_rank, &sel_crowding, parents.len());
+        let p1 = tournament_select(&mut state.evo.rng, &sel_rank, &sel_crowding, parents.len());
+        let p2 = tournament_select(&mut state.evo.rng, &sel_rank, &sel_crowding, parents.len());
 
         let (mut child1, mut child2) = crossover(
-            &mut state.rng,
+            &mut state.evo.rng,
             &parents[p1],
             &parents[p2],
-            &state.dimensions,
+            &state.evo.dimensions,
             state.config.crossover_prob,
             state.config.crossover_eta,
         );
 
         mutate(
-            &mut state.rng,
+            &mut state.evo.rng,
             &mut child1,
-            &state.dimensions,
+            &state.evo.dimensions,
             state.config.mutation_eta,
         );
         mutate(
-            &mut state.rng,
+            &mut state.evo.rng,
             &mut child2,
-            &state.dimensions,
+            &state.evo.dimensions,
             state.config.mutation_eta,
         );
 
@@ -553,15 +355,16 @@ fn nsga2_generate_offspring(
     offspring
 }
 
-// ---------------------------------------------------------------------------
-// Genetic operators
-// ---------------------------------------------------------------------------
-
 /// Tournament selection: pick 2 random individuals, return index of winner.
 /// Winner has lower rank; ties broken by higher crowding distance.
-fn tournament_select(rng: &mut StdRng, ranks: &[usize], crowding: &[f64], n: usize) -> usize {
-    let a = rng.random_range(0..n);
-    let b = rng.random_range(0..n);
+fn tournament_select(
+    rng: &mut fastrand::Rng,
+    ranks: &[usize],
+    crowding: &[f64],
+    n: usize,
+) -> usize {
+    let a = rng.usize(0..n);
+    let b = rng.usize(0..n);
 
     if ranks[a] < ranks[b] {
         a
@@ -573,191 +376,3 @@ fn tournament_select(rng: &mut StdRng, ranks: &[usize], crowding: &[f64], n: usi
         b
     }
 }
-
-/// SBX crossover for continuous params, uniform crossover for categorical.
-fn crossover(
-    rng: &mut StdRng,
-    parent1: &[ParamValue],
-    parent2: &[ParamValue],
-    dimensions: &[DimensionInfo],
-    crossover_prob: f64,
-    eta: f64,
-) -> (Vec<ParamValue>, Vec<ParamValue>) {
-    let n = parent1.len();
-    let mut child1 = parent1.to_vec();
-    let mut child2 = parent2.to_vec();
-
-    let u: f64 = rng.random_range(0.0..=1.0);
-    if u > crossover_prob {
-        return (child1, child2);
-    }
-
-    for i in 0..n {
-        match (&parent1[i], &parent2[i], &dimensions[i].distribution) {
-            (ParamValue::Float(p1), ParamValue::Float(p2), Distribution::Float(d)) => {
-                if (p1 - p2).abs() < 1e-14 {
-                    continue;
-                }
-                let (c1, c2) = sbx_crossover_f64(rng, *p1, *p2, d.low, d.high, eta);
-                child1[i] = ParamValue::Float(c1);
-                child2[i] = ParamValue::Float(c2);
-            }
-            (ParamValue::Int(p1), ParamValue::Int(p2), Distribution::Int(d)) => {
-                if p1 == p2 {
-                    continue;
-                }
-                #[allow(clippy::cast_precision_loss)]
-                let (c1, c2) = sbx_crossover_f64(
-                    rng,
-                    *p1 as f64,
-                    *p2 as f64,
-                    d.low as f64,
-                    d.high as f64,
-                    eta,
-                );
-                #[allow(clippy::cast_possible_truncation)]
-                {
-                    child1[i] = ParamValue::Int((c1.round() as i64).clamp(d.low, d.high));
-                    child2[i] = ParamValue::Int((c2.round() as i64).clamp(d.low, d.high));
-                }
-            }
-            (ParamValue::Categorical(_), ParamValue::Categorical(_), _) => {
-                // Uniform crossover: swap with 50% probability
-                if rng.random_range(0.0..=1.0) < 0.5 {
-                    core::mem::swap(&mut child1[i], &mut child2[i]);
-                }
-            }
-            _ => {}
-        }
-    }
-
-    (child1, child2)
-}
-
-/// SBX crossover for a single float dimension.
-fn sbx_crossover_f64(
-    rng: &mut StdRng,
-    p1: f64,
-    p2: f64,
-    low: f64,
-    high: f64,
-    eta: f64,
-) -> (f64, f64) {
-    let u: f64 = rng.random_range(0.0_f64..1.0_f64);
-
-    let beta = if u <= 0.5 {
-        (2.0 * u).powf(1.0 / (eta + 1.0))
-    } else {
-        (1.0 / (2.0 * (1.0 - u))).powf(1.0 / (eta + 1.0))
-    };
-
-    let c1 = 0.5 * ((1.0 + beta) * p1 + (1.0 - beta) * p2);
-    let c2 = 0.5 * ((1.0 - beta) * p1 + (1.0 + beta) * p2);
-
-    (c1.clamp(low, high), c2.clamp(low, high))
-}
-
-/// Polynomial mutation for each dimension.
-#[allow(clippy::cast_precision_loss)]
-fn mutate(rng: &mut StdRng, individual: &mut [ParamValue], dimensions: &[DimensionInfo], eta: f64) {
-    let n = individual.len();
-    if n == 0 {
-        return;
-    }
-    let mutation_prob = 1.0 / n as f64;
-
-    for (i, value) in individual.iter_mut().enumerate() {
-        if rng.random_range(0.0..=1.0) >= mutation_prob {
-            continue;
-        }
-
-        match (value, &dimensions[i].distribution) {
-            (v @ ParamValue::Float(_), Distribution::Float(d)) => {
-                let ParamValue::Float(x) = *v else {
-                    unreachable!();
-                };
-                let mutated = polynomial_mutation_f64(rng, x, d.low, d.high, eta);
-                *v = ParamValue::Float(mutated);
-            }
-            (v @ ParamValue::Int(_), Distribution::Int(d)) => {
-                let ParamValue::Int(x) = *v else {
-                    unreachable!();
-                };
-                #[allow(clippy::cast_possible_truncation)]
-                {
-                    let mutated =
-                        polynomial_mutation_f64(rng, x as f64, d.low as f64, d.high as f64, eta);
-                    *v = ParamValue::Int((mutated.round() as i64).clamp(d.low, d.high));
-                }
-            }
-            (v @ ParamValue::Categorical(_), Distribution::Categorical(d)) => {
-                *v = ParamValue::Categorical(rng.random_range(0..d.n_choices));
-            }
-            _ => {}
-        }
-    }
-}
-
-/// Polynomial mutation for a single float value.
-fn polynomial_mutation_f64(rng: &mut StdRng, x: f64, low: f64, high: f64, eta: f64) -> f64 {
-    let u: f64 = rng.random_range(0.0_f64..1.0_f64);
-    let range = high - low;
-    if range <= 0.0 {
-        return x;
-    }
-
-    let delta1 = (x - low) / range;
-    let delta2 = (high - x) / range;
-
-    let delta_q = if u < 0.5 {
-        let xy = 1.0 - delta1;
-        let val = 2.0 * u + (1.0 - 2.0 * u) * xy.powf(eta + 1.0);
-        val.powf(1.0 / (eta + 1.0)) - 1.0
-    } else {
-        let xy = 1.0 - delta2;
-        let val = 2.0 * (1.0 - u) + 2.0 * (u - 0.5) * xy.powf(eta + 1.0);
-        1.0 - val.powf(1.0 / (eta + 1.0))
-    };
-
-    (x + delta_q * range).clamp(low, high)
-}
-
-// ---------------------------------------------------------------------------
-// Random sampling helper (for discovery phase)
-// ---------------------------------------------------------------------------
-
-#[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
-fn sample_random(rng: &mut StdRng, distribution: &Distribution) -> ParamValue {
-    match distribution {
-        Distribution::Float(d) => {
-            let value = if d.log_scale {
-                let log_low = d.low.ln();
-                let log_high = d.high.ln();
-                rng.random_range(log_low..=log_high).exp()
-            } else if let Some(step) = d.step {
-                let n_steps = ((d.high - d.low) / step).floor() as i64;
-                let k = rng.random_range(0..=n_steps);
-                d.low + (k as f64) * step
-            } else {
-                rng.random_range(d.low..=d.high)
-            };
-            ParamValue::Float(value)
-        }
-        Distribution::Int(d) => {
-            let value = if d.log_scale {
-                let log_low = (d.low as f64).ln();
-                let log_high = (d.high as f64).ln();
-                let raw = rng.random_range(log_low..=log_high).exp().round() as i64;
-                raw.clamp(d.low, d.high)
-            } else if let Some(step) = d.step {
-                let n_steps = (d.high - d.low) / step;
-                let k = rng.random_range(0..=n_steps);
-                d.low + k * step
-            } else {
-                rng.random_range(d.low..=d.high)
-            };
-            ParamValue::Int(value)
-        }
-        Distribution::Categorical(d) => ParamValue::Categorical(rng.random_range(0..d.n_choices)),
-    }
-}

src/sampler/nsga3.rs

@@ -0,0 +1,720 @@
+//! NSGA-III (Non-dominated Sorting Genetic Algorithm III) sampler.
+//!
+//! Uses reference-point-based niching for better diversity in
+//! many-objective (3+) optimization problems. Das-Dennis structured
+//! reference points guide the search toward a well-distributed
+//! Pareto front.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::Direction;
+//! use optimizer::multi_objective::MultiObjectiveStudy;
+//! use optimizer::parameter::{FloatParam, Parameter};
+//! use optimizer::sampler::nsga3::Nsga3Sampler;
+//!
+//! let sampler = Nsga3Sampler::with_seed(42);
+//! let study = MultiObjectiveStudy::with_sampler(
+//!     vec![
+//!         Direction::Minimize,
+//!         Direction::Minimize,
+//!         Direction::Minimize,
+//!     ],
+//!     sampler,
+//! );
+//!
+//! let x = FloatParam::new(0.0, 1.0);
+//! let y = FloatParam::new(0.0, 1.0);
+//! study
+//!     .optimize(100, |trial| {
+//!         let xv = x.suggest(trial)?;
+//!         let yv = y.suggest(trial)?;
+//!         Ok::<_, optimizer::Error>(vec![xv, yv, (1.0 - xv - yv).abs()])
+//!     })
+//!     .unwrap();
+//! ```
+
+use parking_lot::Mutex;
+
+use super::genetic::{
+    self, Candidate, EvolutionaryState, Phase, advance_generation, auto_divisions,
+    collect_evaluated_generation, crossover, das_dennis, extract_trial_params,
+    generate_random_candidates, mutate, sample_from_candidate, sample_random,
+};
+use crate::distribution::Distribution;
+use crate::multi_objective::MultiObjectiveTrial;
+use crate::param::ParamValue;
+use crate::pareto;
+use crate::types::Direction;
+
+/// NSGA-III sampler for multi-objective optimization.
+///
+/// Uses reference-point-based niching to maintain diversity,
+/// especially effective for problems with 3 or more objectives.
+pub struct Nsga3Sampler {
+    state: Mutex<Nsga3State>,
+}
+
+impl Nsga3Sampler {
+    /// Creates a new NSGA-III sampler with a random seed.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            state: Mutex::new(Nsga3State::new(Nsga3Config::default(), None)),
+        }
+    }
+
+    /// Creates a new NSGA-III sampler with a fixed seed.
+    #[must_use]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            state: Mutex::new(Nsga3State::new(Nsga3Config::default(), Some(seed))),
+        }
+    }
+
+    /// Creates a builder for configuring an `Nsga3Sampler`.
+    #[must_use]
+    pub fn builder() -> Nsga3SamplerBuilder {
+        Nsga3SamplerBuilder::default()
+    }
+}
+
+impl Default for Nsga3Sampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+/// Builder for [`Nsga3Sampler`].
+#[derive(Debug, Clone, Default)]
+pub struct Nsga3SamplerBuilder {
+    population_size: Option<usize>,
+    n_divisions: Option<usize>,
+    crossover_prob: Option<f64>,
+    crossover_eta: Option<f64>,
+    mutation_eta: Option<f64>,
+    seed: Option<u64>,
+}
+
+impl Nsga3SamplerBuilder {
+    /// Sets the population size. If unset, equals the number of
+    /// Das-Dennis reference points.
+    #[must_use]
+    pub fn population_size(mut self, size: usize) -> Self {
+        self.population_size = Some(size);
+        self
+    }
+
+    /// Sets the number of divisions (H) for Das-Dennis reference points.
+    /// If unset, automatically chosen based on population size and number
+    /// of objectives.
+    #[must_use]
+    pub fn n_divisions(mut self, h: usize) -> Self {
+        self.n_divisions = Some(h);
+        self
+    }
+
+    /// Sets the crossover probability. Default: 1.0.
+    #[must_use]
+    pub fn crossover_prob(mut self, prob: f64) -> Self {
+        self.crossover_prob = Some(prob);
+        self
+    }
+
+    /// Sets the SBX distribution index. Default: 30.0.
+    #[must_use]
+    pub fn crossover_eta(mut self, eta: f64) -> Self {
+        self.crossover_eta = Some(eta);
+        self
+    }
+
+    /// Sets the polynomial mutation distribution index. Default: 20.0.
+    #[must_use]
+    pub fn mutation_eta(mut self, eta: f64) -> Self {
+        self.mutation_eta = Some(eta);
+        self
+    }
+
+    /// Sets the random seed for reproducibility.
+    #[must_use]
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.seed = Some(seed);
+        self
+    }
+
+    /// Builds the configured [`Nsga3Sampler`].
+    #[must_use]
+    pub fn build(self) -> Nsga3Sampler {
+        let config = Nsga3Config {
+            user_population_size: self.population_size,
+            n_divisions: self.n_divisions,
+            crossover_prob: self.crossover_prob.unwrap_or(1.0),
+            crossover_eta: self.crossover_eta.unwrap_or(30.0),
+            mutation_eta: self.mutation_eta.unwrap_or(20.0),
+        };
+        Nsga3Sampler {
+            state: Mutex::new(Nsga3State::new(config, self.seed)),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Internal types
+// ---------------------------------------------------------------------------
+
+#[derive(Clone, Debug)]
+struct Nsga3Config {
+    user_population_size: Option<usize>,
+    n_divisions: Option<usize>,
+    crossover_prob: f64,
+    crossover_eta: f64,
+    mutation_eta: f64,
+}
+
+impl Default for Nsga3Config {
+    fn default() -> Self {
+        Self {
+            user_population_size: None,
+            n_divisions: None,
+            crossover_prob: 1.0,
+            crossover_eta: 30.0,
+            mutation_eta: 20.0,
+        }
+    }
+}
+
+struct Nsga3State {
+    evo: EvolutionaryState,
+    config: Nsga3Config,
+    /// Das-Dennis reference points (lazily generated once objectives are known).
+    reference_points: Vec<Vec<f64>>,
+    /// Best value seen per objective (minimize-space).
+    ideal_point: Vec<f64>,
+    /// Whether reference points have been initialized.
+    initialized: bool,
+}
+
+impl Nsga3State {
+    fn new(config: Nsga3Config, seed: Option<u64>) -> Self {
+        Self {
+            evo: EvolutionaryState::new(seed),
+            config,
+            reference_points: Vec::new(),
+            ideal_point: Vec::new(),
+            initialized: false,
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MultiObjectiveSampler implementation
+// ---------------------------------------------------------------------------
+
+impl crate::multi_objective::MultiObjectiveSampler for Nsga3Sampler {
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[MultiObjectiveTrial],
+        directions: &[Direction],
+    ) -> ParamValue {
+        let mut state = self.state.lock();
+
+        match &state.evo.phase {
+            Phase::Discovery => {
+                if let Some(value) =
+                    genetic::sample_discovery(&mut state.evo, distribution, trial_id)
+                {
+                    return value;
+                }
+                // Transitioned to active phase
+                initialize_nsga3(&mut state, directions);
+                generate_random_candidates(&mut state.evo);
+                sample_from_candidate(&mut state.evo, trial_id)
+            }
+            Phase::Active => {
+                maybe_generate_new_generation(&mut state, history, directions);
+                sample_from_candidate(&mut state.evo, trial_id)
+            }
+        }
+    }
+}
+
+/// Initialize NSGA-III: generate reference points and set population size.
+fn initialize_nsga3(state: &mut Nsga3State, directions: &[Direction]) {
+    let n_obj = directions.len();
+
+    // Determine divisions
+    let divisions = state
+        .config
+        .n_divisions
+        .unwrap_or_else(|| auto_divisions(n_obj, state.config.user_population_size.unwrap_or(100)));
+
+    state.reference_points = das_dennis(n_obj, divisions);
+    let n_ref = state.reference_points.len();
+
+    // Population size = number of reference points (or user override, at least n_ref)
+    let pop_size = state.config.user_population_size.unwrap_or(n_ref).max(4);
+    state.evo.population_size = pop_size;
+    state.evo.phase = Phase::Active;
+    state.ideal_point = vec![f64::INFINITY; n_obj];
+    state.initialized = true;
+}
+
+fn maybe_generate_new_generation(
+    state: &mut Nsga3State,
+    history: &[MultiObjectiveTrial],
+    directions: &[Direction],
+) {
+    if state.evo.candidates.is_empty() {
+        generate_random_candidates(&mut state.evo);
+        return;
+    }
+
+    if let Some(evaluated) = collect_evaluated_generation(&state.evo, history) {
+        let offspring = nsga3_generate_offspring(state, &evaluated, directions);
+        advance_generation(&mut state.evo, offspring);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// NSGA-III selection algorithm
+// ---------------------------------------------------------------------------
+
+/// Normalize objectives to minimize-space.
+fn to_minimize_space(values: &[f64], directions: &[Direction]) -> Vec<f64> {
+    values
+        .iter()
+        .zip(directions)
+        .map(|(&v, d)| match d {
+            Direction::Minimize => v,
+            Direction::Maximize => -v,
+        })
+        .collect()
+}
+
+/// Update ideal point with new observations.
+fn update_ideal_point(ideal: &mut [f64], normalized_values: &[Vec<f64>]) {
+    for vals in normalized_values {
+        for (i, &v) in vals.iter().enumerate() {
+            if v < ideal[i] {
+                ideal[i] = v;
+            }
+        }
+    }
+}
+
+/// Compute Achievement Scalarizing Function (ASF) for extreme point finding.
+fn asf(point: &[f64], weight: &[f64], ideal: &[f64]) -> f64 {
+    point
+        .iter()
+        .zip(weight)
+        .zip(ideal)
+        .map(|((&p, &w), &z)| {
+            let w = if w < 1e-6 { 1e-6 } else { w };
+            (p - z) / w
+        })
+        .fold(f64::NEG_INFINITY, f64::max)
+}
+
+/// Find intercepts for normalization via extreme points.
+///
+/// For each objective, find the point with best ASF (using a weight vector
+/// that emphasizes that objective). The intercepts are where the hyperplane
+/// through the extreme points crosses each axis.
+fn find_intercepts(normalized_values: &[Vec<f64>], ideal: &[f64]) -> Vec<f64> {
+    let n_obj = ideal.len();
+    let n = normalized_values.len();
+
+    if n == 0 || n_obj == 0 {
+        return vec![1.0; n_obj];
+    }
+
+    // Find extreme points (one per objective)
+    let mut extreme_indices = Vec::with_capacity(n_obj);
+    for obj in 0..n_obj {
+        let mut weight = vec![1e-6; n_obj];
+        weight[obj] = 1.0;
+
+        let mut best_idx = 0;
+        let mut best_asf = f64::INFINITY;
+        for (i, vals) in normalized_values.iter().enumerate() {
+            let a = asf(vals, &weight, ideal);
+            if a < best_asf {
+                best_asf = a;
+                best_idx = i;
+            }
+        }
+        extreme_indices.push(best_idx);
+    }
+
+    // Try to compute hyperplane intercepts
+    // For stability, if the extreme points are degenerate, fall back to
+    // max - ideal per objective
+    let mut intercepts = Vec::with_capacity(n_obj);
+    for obj in 0..n_obj {
+        let max_val = normalized_values
+            .iter()
+            .map(|v| v[obj])
+            .fold(f64::NEG_INFINITY, f64::max);
+        let intercept = max_val - ideal[obj];
+        intercepts.push(if intercept > 1e-10 { intercept } else { 1.0 });
+    }
+
+    intercepts
+}
+
+/// Normalize objective values: subtract ideal, divide by intercepts.
+fn normalize_objectives(values: &[Vec<f64>], ideal: &[f64], intercepts: &[f64]) -> Vec<Vec<f64>> {
+    values
+        .iter()
+        .map(|v| {
+            v.iter()
+                .zip(ideal)
+                .zip(intercepts)
+                .map(|((&val, &z), &a)| {
+                    let norm = if a > 1e-10 { a } else { 1.0 };
+                    (val - z) / norm
+                })
+                .collect()
+        })
+        .collect()
+}
+
+/// Perpendicular distance from a point to a reference line (direction vector).
+fn perpendicular_distance(point: &[f64], reference: &[f64]) -> f64 {
+    let dot: f64 = point.iter().zip(reference).map(|(&p, &r)| p * r).sum();
+    let ref_norm_sq: f64 = reference.iter().map(|&r| r * r).sum();
+
+    if ref_norm_sq < 1e-30 {
+        return f64::INFINITY;
+    }
+
+    let proj_scalar = dot / ref_norm_sq;
+    let dist_sq: f64 = point
+        .iter()
+        .zip(reference)
+        .map(|(&p, &r)| {
+            let proj = proj_scalar * r;
+            (p - proj).powi(2)
+        })
+        .sum();
+
+    dist_sq.sqrt()
+}
+
+/// Associate each solution with its nearest reference point.
+/// Returns (`closest_ref_idx`, distance) for each solution.
+fn associate_to_reference_points(
+    normalized: &[Vec<f64>],
+    reference_points: &[Vec<f64>],
+) -> Vec<(usize, f64)> {
+    normalized
+        .iter()
+        .map(|point| {
+            let mut best_ref = 0;
+            let mut best_dist = f64::INFINITY;
+            for (j, rp) in reference_points.iter().enumerate() {
+                let d = perpendicular_distance(point, rp);
+                if d < best_dist {
+                    best_dist = d;
+                    best_ref = j;
+                }
+            }
+            (best_ref, best_dist)
+        })
+        .collect()
+}
+
+/// NSGA-III niching-based selection from the last front.
+///
+/// `already_selected` are indices into the combined population that are
+/// already accepted (from fronts 0..L-1). `last_front` contains indices
+/// from front L. We need to pick `remaining` more from `last_front`.
+fn niching_select(
+    rng: &mut fastrand::Rng,
+    associations: &[(usize, f64)],
+    already_selected: &[usize],
+    last_front: &[usize],
+    n_reference_points: usize,
+    remaining: usize,
+) -> Vec<usize> {
+    // Count niche per reference point for already selected
+    let mut niche_count = vec![0_usize; n_reference_points];
+    for &idx in already_selected {
+        niche_count[associations[idx].0] += 1;
+    }
+
+    // Build per-reference-point candidate lists from the last front
+    let mut ref_candidates: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n_reference_points];
+    for &idx in last_front {
+        let (ref_idx, dist) = associations[idx];
+        ref_candidates[ref_idx].push((idx, dist));
+    }
+
+    let mut selected = Vec::with_capacity(remaining);
+    let mut excluded = vec![false; associations.len()];
+
+    for _ in 0..remaining {
+        // Find minimum niche count among reference points that still have candidates
+        let min_count = (0..n_reference_points)
+            .filter(|&j| ref_candidates[j].iter().any(|&(idx, _)| !excluded[idx]))
+            .map(|j| niche_count[j])
+            .min();
+
+        let Some(min_count) = min_count else {
+            break;
+        };
+
+        // Collect reference points with this minimum count that have candidates
+        let min_refs: Vec<usize> = (0..n_reference_points)
+            .filter(|&j| {
+                niche_count[j] == min_count
+                    && ref_candidates[j].iter().any(|&(idx, _)| !excluded[idx])
+            })
+            .collect();
+
+        if min_refs.is_empty() {
+            break;
+        }
+
+        // Pick a random reference point from the minimum set
+        let chosen_ref = min_refs[rng.usize(0..min_refs.len())];
+
+        // Available candidates for this reference point
+        let available: Vec<(usize, f64)> = ref_candidates[chosen_ref]
+            .iter()
+            .filter(|&&(idx, _)| !excluded[idx])
+            .copied()
+            .collect();
+
+        if available.is_empty() {
+            continue;
+        }
+
+        let chosen_idx = if min_count == 0 {
+            // Pick closest to reference line
+            available
+                .iter()
+                .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(core::cmp::Ordering::Equal))
+                .unwrap()
+                .0
+        } else {
+            // Pick random
+            available[rng.usize(0..available.len())].0
+        };
+
+        selected.push(chosen_idx);
+        excluded[chosen_idx] = true;
+        niche_count[chosen_ref] += 1;
+    }
+
+    selected
+}
+
+/// Perform NSGA-III selection: non-dominated sort + reference-point niching.
+fn nsga3_select(
+    state: &mut Nsga3State,
+    population: &[&MultiObjectiveTrial],
+    directions: &[Direction],
+) -> Vec<Vec<ParamValue>> {
+    let pop_size = state.evo.population_size;
+    let n_obj = directions.len();
+
+    // Convert to minimize-space
+    let min_values: Vec<Vec<f64>> = population
+        .iter()
+        .map(|t| to_minimize_space(&t.values, directions))
+        .collect();
+
+    // Non-dominated sort
+    let constraints: Vec<Vec<f64>> = population.iter().map(|t| t.constraints.clone()).collect();
+    let has_constraints = constraints.iter().any(|c| !c.is_empty());
+    let fronts = if has_constraints {
+        pareto::fast_non_dominated_sort_constrained(
+            &min_values,
+            &vec![Direction::Minimize; n_obj],
+            &constraints,
+        )
+    } else {
+        pareto::fast_non_dominated_sort(&min_values, &vec![Direction::Minimize; n_obj])
+    };
+
+    // Fill front-by-front
+    let mut selected: Vec<usize> = Vec::with_capacity(pop_size);
+    let mut last_front_idx = None;
+
+    for (fi, front) in fronts.iter().enumerate() {
+        if selected.len() + front.len() <= pop_size {
+            selected.extend_from_slice(front);
+        } else {
+            last_front_idx = Some(fi);
+            break;
+        }
+    }
+
+    // If we filled exactly or all fronts fit, done
+    if selected.len() < pop_size
+        && let Some(lf_idx) = last_front_idx
+    {
+        // Need niching from the last partial front
+        let remaining = pop_size - selected.len();
+
+        // Update ideal point
+        update_ideal_point(&mut state.ideal_point, &min_values);
+
+        // Find intercepts and normalize
+        let intercepts = find_intercepts(&min_values, &state.ideal_point);
+        let normalized = normalize_objectives(&min_values, &state.ideal_point, &intercepts);
+
+        // Associate all solutions with reference points
+        let associations = associate_to_reference_points(&normalized, &state.reference_points);
+
+        // Select from last front using niching
+        let last_front = &fronts[lf_idx];
+        let additional = niching_select(
+            &mut state.evo.rng,
+            &associations,
+            &selected,
+            last_front,
+            state.reference_points.len(),
+            remaining,
+        );
+        selected.extend(additional);
+    }
+
+    // Pad if needed
+    let n = population.len();
+    while selected.len() < pop_size {
+        selected.push(state.evo.rng.usize(0..n));
+    }
+
+    selected
+        .iter()
+        .map(|&idx| {
+            extract_trial_params(population[idx], &state.evo.dimensions, &mut state.evo.rng)
+        })
+        .collect()
+}
+
+/// Tournament selection based on rank only (no crowding distance in NSGA-III).
+fn tournament_select_rank(rng: &mut fastrand::Rng, ranks: &[usize], n: usize) -> usize {
+    let a = rng.usize(0..n);
+    let b = rng.usize(0..n);
+
+    if ranks[a] <= ranks[b] { a } else { b }
+}
+
+fn nsga3_generate_offspring(
+    state: &mut Nsga3State,
+    population: &[&MultiObjectiveTrial],
+    directions: &[Direction],
+) -> Vec<Candidate> {
+    let pop_size = state.evo.population_size;
+
+    if population.len() < 2 {
+        return (0..pop_size)
+            .map(|_| {
+                let params = state
+                    .evo
+                    .dimensions
+                    .iter()
+                    .map(|d| sample_random(&mut state.evo.rng, &d.distribution))
+                    .collect();
+                Candidate { params }
+            })
+            .collect();
+    }
+
+    // Initialize reference points and ideal on first generation
+    if !state.initialized {
+        initialize_nsga3(state, directions);
+    }
+
+    let parents = nsga3_select(state, population, directions);
+
+    // Assign ranks for tournament selection
+    let n_obj = directions.len();
+    let min_values: Vec<Vec<f64>> = population
+        .iter()
+        .map(|t| to_minimize_space(&t.values, directions))
+        .collect();
+    let fronts = pareto::fast_non_dominated_sort(&min_values, &vec![Direction::Minimize; n_obj]);
+    let mut rank = vec![0_usize; parents.len()];
+    for (front_rank, front) in fronts.iter().enumerate() {
+        for &idx in front {
+            if idx < rank.len() {
+                rank[idx] = front_rank;
+            }
+        }
+    }
+    // Ranks for selected parents (simplified: use index order)
+    let parent_ranks: Vec<usize> = (0..parents.len())
+        .map(|i| i % (fronts.len().max(1)))
+        .collect();
+
+    let mut offspring = Vec::with_capacity(pop_size);
+    while offspring.len() < pop_size {
+        let p1 = tournament_select_rank(&mut state.evo.rng, &parent_ranks, parents.len());
+        let p2 = tournament_select_rank(&mut state.evo.rng, &parent_ranks, parents.len());
+
+        let (mut child1, mut child2) = crossover(
+            &mut state.evo.rng,
+            &parents[p1],
+            &parents[p2],
+            &state.evo.dimensions,
+            state.config.crossover_prob,
+            state.config.crossover_eta,
+        );
+
+        mutate(
+            &mut state.evo.rng,
+            &mut child1,
+            &state.evo.dimensions,
+            state.config.mutation_eta,
+        );
+        mutate(
+            &mut state.evo.rng,
+            &mut child2,
+            &state.evo.dimensions,
+            state.config.mutation_eta,
+        );
+
+        offspring.push(Candidate { params: child1 });
+        if offspring.len() < pop_size {
+            offspring.push(Candidate { params: child2 });
+        }
+    }
+
+    offspring
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_perpendicular_distance() {
+        // Point (1, 0) to reference line (1, 1) (45-degree line)
+        let d = perpendicular_distance(&[1.0, 0.0], &[1.0, 1.0]);
+        // Projection is (0.5, 0.5), distance = sqrt(0.25 + 0.25) = sqrt(0.5)
+        assert!((d - (0.5_f64).sqrt()).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_perpendicular_distance_on_line() {
+        // Point on the reference line
+        let d = perpendicular_distance(&[2.0, 2.0], &[1.0, 1.0]);
+        assert!(d < 1e-10);
+    }
+
+    #[test]
+    fn test_normalize_objectives() {
+        let values = vec![vec![2.0, 4.0], vec![4.0, 2.0]];
+        let ideal = vec![1.0, 1.0];
+        let intercepts = vec![3.0, 3.0];
+        let normalized = normalize_objectives(&values, &ideal, &intercepts);
+        assert!((normalized[0][0] - 1.0 / 3.0).abs() < 1e-10);
+        assert!((normalized[0][1] - 1.0).abs() < 1e-10);
+    }
+}

src/sampler/random.rs

@@ -1,11 +1,10 @@
 //! Random sampler implementation.
 
 use parking_lot::Mutex;
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
 
 use crate::distribution::Distribution;
 use crate::param::ParamValue;
+use crate::rng_util;
 use crate::sampler::{CompletedTrial, Sampler};
 
 /// A simple random sampler that samples uniformly from distributions.
@@ -26,7 +25,7 @@ use crate::sampler::{CompletedTrial, Sampler};
 /// let sampler = RandomSampler::with_seed(42);
 /// ```
 pub struct RandomSampler {
-    rng: Mutex<StdRng>,
+    rng: Mutex<fastrand::Rng>,
 }
 
 impl RandomSampler {
@@ -34,7 +33,7 @@ impl RandomSampler {
     #[must_use]
     pub fn new() -> Self {
         Self {
-            rng: Mutex::new(rand::make_rng()),
+            rng: Mutex::new(fastrand::Rng::new()),
         }
     }
 
@@ -44,7 +43,7 @@ impl RandomSampler {
     #[must_use]
     pub fn with_seed(seed: u64) -> Self {
         Self {
-            rng: Mutex::new(StdRng::seed_from_u64(seed)),
+            rng: Mutex::new(fastrand::Rng::with_seed(seed)),
         }
     }
 }
@@ -71,16 +70,16 @@ impl Sampler for RandomSampler {
                     // Sample uniformly in log space
                     let log_low = d.low.ln();
                     let log_high = d.high.ln();
-                    let log_value = rng.random_range(log_low..=log_high);
+                    let log_value = rng_util::f64_range(&mut rng, log_low, log_high);
                     log_value.exp()
                 } else if let Some(step) = d.step {
                     // Sample from step grid
                     let n_steps = ((d.high - d.low) / step).floor() as i64;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + (k as f64) * step
                 } else {
                     // Uniform sampling
-                    rng.random_range(d.low..=d.high)
+                    rng_util::f64_range(&mut rng, d.low, d.high)
                 };
                 ParamValue::Float(value)
             }
@@ -89,23 +88,23 @@ impl Sampler for RandomSampler {
                     // Sample uniformly in log space, then round
                     let log_low = (d.low as f64).ln();
                     let log_high = (d.high as f64).ln();
-                    let log_value = rng.random_range(log_low..=log_high);
+                    let log_value = rng_util::f64_range(&mut rng, log_low, log_high);
                     let raw = log_value.exp().round() as i64;
                     // Clamp to bounds since rounding might push outside
                     raw.clamp(d.low, d.high)
                 } else if let Some(step) = d.step {
                     // Sample from step grid
                     let n_steps = (d.high - d.low) / step;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + k * step
                 } else {
                     // Uniform sampling
-                    rng.random_range(d.low..=d.high)
+                    rng.i64(d.low..=d.high)
                 };
                 ParamValue::Int(value)
             }
             Distribution::Categorical(d) => {
-                let index = rng.random_range(0..d.n_choices);
+                let index = rng.usize(0..d.n_choices);
                 ParamValue::Categorical(index)
             }
         }

src/sampler/tpe/multivariate.rs

@@ -116,14 +116,13 @@ use std::collections::HashMap;
 use std::sync::Arc;
 
 use parking_lot::Mutex;
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
 
 use super::{FixedGamma, GammaStrategy};
 use crate::distribution::Distribution;
 use crate::error::Result;
 use crate::param::ParamValue;
 use crate::parameter::ParamId;
+use crate::rng_util;
 use crate::sampler::{CompletedTrial, PendingTrial, Sampler};
 
 /// Strategy for imputing objective values for pending/running trials during parallel optimization.
@@ -189,7 +188,7 @@ pub struct MultivariateTpeSampler {
     /// Strategy for imputing objective values for pending trials in parallel optimization.
     constant_liar: ConstantLiarStrategy,
     /// Thread-safe RNG for sampling.
-    rng: Mutex<StdRng>,
+    rng: Mutex<fastrand::Rng>,
     /// Cache for joint samples to maintain consistency across parameters within the same trial.
     /// The tuple contains (`trial_id`, cached joint sample).
     joint_sample_cache: Mutex<Option<(u64, HashMap<ParamId, ParamValue>)>>,
@@ -219,7 +218,7 @@ impl MultivariateTpeSampler {
             n_ei_candidates: 24,
             group: false,
             constant_liar: ConstantLiarStrategy::None,
-            rng: Mutex::new(rand::make_rng()),
+            rng: Mutex::new(fastrand::Rng::new()),
             joint_sample_cache: Mutex::new(None),
         }
     }
@@ -695,7 +694,7 @@ impl MultivariateTpeSampler {
 
         // Generate candidates from the good distribution
         let candidates: Vec<Vec<f64>> = (0..self.n_ei_candidates)
-            .map(|_| good_kde.sample(&mut *rng))
+            .map(|_| good_kde.sample(&mut rng))
             .collect();
 
         // Compute log(l(x)) - log(g(x)) for each candidate
@@ -731,7 +730,7 @@ impl MultivariateTpeSampler {
         &self,
         good_kde: &crate::kde::MultivariateKDE,
         bad_kde: &crate::kde::MultivariateKDE,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> Vec<f64> {
         // Generate candidates from the good distribution
         let candidates: Vec<Vec<f64>> = (0..self.n_ei_candidates)
@@ -769,7 +768,7 @@ impl MultivariateTpeSampler {
     fn sample_all_uniform(
         &self,
         search_space: &HashMap<ParamId, Distribution>,
-        rng: &mut rand::rngs::StdRng,
+        rng: &mut fastrand::Rng,
     ) -> HashMap<ParamId, ParamValue> {
         search_space
             .iter()
@@ -778,25 +777,22 @@ impl MultivariateTpeSampler {
     }
 
     /// Samples a single parameter uniformly at random from its distribution.
-    fn sample_uniform_single(
-        distribution: &Distribution,
-        rng: &mut rand::rngs::StdRng,
-    ) -> ParamValue {
+    fn sample_uniform_single(distribution: &Distribution, rng: &mut fastrand::Rng) -> ParamValue {
         match distribution {
             Distribution::Float(d) => {
                 let value = if d.log_scale {
                     let log_low = d.low.ln();
                     let log_high = d.high.ln();
-                    rng.random_range(log_low..=log_high).exp()
+                    rng_util::f64_range(rng, log_low, log_high).exp()
                 } else if let Some(step) = d.step {
                     #[allow(clippy::cast_possible_truncation)]
                     let n_steps = ((d.high - d.low) / step).floor() as i64;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     #[allow(clippy::cast_precision_loss)]
                     let result = d.low + (k as f64) * step;
                     result
                 } else {
-                    rng.random_range(d.low..=d.high)
+                    rng_util::f64_range(rng, d.low, d.high)
                 };
                 ParamValue::Float(value)
             }
@@ -806,20 +802,20 @@ impl MultivariateTpeSampler {
                     let log_low = (d.low as f64).ln();
                     let log_high = (d.high as f64).ln();
                     #[allow(clippy::cast_possible_truncation)]
-                    let raw = rng.random_range(log_low..=log_high).exp().round() as i64;
+                    let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
                     raw.clamp(d.low, d.high)
                 } else if let Some(step) = d.step {
                     #[allow(clippy::cast_possible_truncation)]
                     let n_steps = (d.high - d.low) / step;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + k * step
                 } else {
-                    rng.random_range(d.low..=d.high)
+                    rng.i64(d.low..=d.high)
                 };
                 ParamValue::Int(value)
             }
             Distribution::Categorical(d) => {
-                let index = rng.random_range(0..d.n_choices);
+                let index = rng.usize(0..d.n_choices);
                 ParamValue::Categorical(index)
             }
         }
@@ -1018,7 +1014,7 @@ impl MultivariateTpeSampler {
         &self,
         search_space: &HashMap<ParamId, Distribution>,
         history: &[CompletedTrial],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> HashMap<ParamId, ParamValue> {
         use super::IntersectionSearchSpace;
         use crate::kde::MultivariateKDE;
@@ -1220,7 +1216,7 @@ impl MultivariateTpeSampler {
         _intersection: &HashMap<ParamId, Distribution>,
         history: &[CompletedTrial],
         result: &mut HashMap<ParamId, ParamValue>,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) {
         // Identify parameters not in result (and not in intersection)
         let missing_params: Vec<(&ParamId, &Distribution)> = search_space
@@ -1253,7 +1249,7 @@ impl MultivariateTpeSampler {
         distribution: &Distribution,
         good_trials: &[&CompletedTrial],
         bad_trials: &[&CompletedTrial],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> ParamValue {
         match distribution {
             Distribution::Float(d) => {
@@ -1370,7 +1366,7 @@ impl MultivariateTpeSampler {
         step: Option<f64>,
         good_values: Vec<f64>,
         bad_values: Vec<f64>,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> f64 {
         use crate::kde::KernelDensityEstimator;
 
@@ -1391,7 +1387,7 @@ impl MultivariateTpeSampler {
 
         // If KDE construction fails, fall back to uniform sampling
         let (Ok(l_kde), Ok(g_kde)) = (l_kde, g_kde) else {
-            return rng.random_range(low..=high);
+            return rng_util::f64_range(rng, low, high);
         };
 
         // Generate candidates from l(x) and select the one with best l(x)/g(x) ratio
@@ -1455,7 +1451,7 @@ impl MultivariateTpeSampler {
         step: Option<i64>,
         good_values: &[i64],
         bad_values: &[i64],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> i64 {
         // Convert to floats for KDE
         let good_floats: Vec<f64> = good_values.iter().map(|&v| v as f64).collect();
@@ -1518,7 +1514,7 @@ impl MultivariateTpeSampler {
         &self,
         search_space: &HashMap<ParamId, Distribution>,
         history: &[CompletedTrial],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> HashMap<ParamId, ParamValue> {
         // Split trials for independent sampling
         let (good_trials, bad_trials) = self.split_trials(&history.iter().collect::<Vec<_>>());
@@ -1562,7 +1558,7 @@ impl MultivariateTpeSampler {
         n_choices: usize,
         good_indices: &[usize],
         bad_indices: &[usize],
-        rng: &mut rand::rngs::StdRng,
+        rng: &mut fastrand::Rng,
     ) -> usize {
         // Count occurrences in good and bad groups
         let mut good_counts = vec![0usize; n_choices];
@@ -1593,7 +1589,7 @@ impl MultivariateTpeSampler {
 
         // Sample proportionally to weights
         let total_weight: f64 = weights.iter().sum();
-        let threshold = rng.random::<f64>() * total_weight;
+        let threshold = rng.f64() * total_weight;
 
         let mut cumulative = 0.0;
         for (i, &w) in weights.iter().enumerate() {
@@ -2069,8 +2065,8 @@ impl MultivariateTpeSamplerBuilder {
         };
 
         let rng = match self.seed {
-            Some(s) => StdRng::seed_from_u64(s),
-            None => rand::make_rng(),
+            Some(s) => fastrand::Rng::with_seed(s),
+            None => fastrand::Rng::new(),
         };
 
         Ok(MultivariateTpeSampler {
@@ -4854,8 +4850,7 @@ mod tests {
 
         #[test]
         fn test_sample_tpe_categorical_basic() {
-            use rand::SeedableRng;
-            let mut rng = rand::rngs::StdRng::seed_from_u64(42);
+            let mut rng = fastrand::Rng::with_seed(42);
 
             // Category 0 is good (appears more in good trials)
             let good_indices = vec![0, 0, 0, 1];
@@ -4890,8 +4885,7 @@ mod tests {
 
         #[test]
         fn test_sample_tpe_categorical_laplace_smoothing() {
-            use rand::SeedableRng;
-            let mut rng = rand::rngs::StdRng::seed_from_u64(42);
+            let mut rng = fastrand::Rng::with_seed(42);
 
             // Category 2 never appears, but should still be sampled due to Laplace smoothing
             let good_indices = vec![0, 0, 1];
@@ -4919,8 +4913,7 @@ mod tests {
 
         #[test]
         fn test_sample_tpe_categorical_empty_good() {
-            use rand::SeedableRng;
-            let mut rng = rand::rngs::StdRng::seed_from_u64(42);
+            let mut rng = fastrand::Rng::with_seed(42);
 
             // Empty good group - all categories should have equal probability
             let good_indices: Vec<usize> = vec![];
@@ -4945,8 +4938,7 @@ mod tests {
 
         #[test]
         fn test_sample_tpe_categorical_all_indices_valid() {
-            use rand::SeedableRng;
-            let mut rng = rand::rngs::StdRng::seed_from_u64(42);
+            let mut rng = fastrand::Rng::with_seed(42);
 
             let n_choices = 4;
             let good_indices = vec![0, 1, 2, 3];

src/sampler/tpe/sampler.rs

@@ -59,13 +59,12 @@ use core::fmt::Debug;
 use std::sync::Arc;
 
 use parking_lot::Mutex;
-use rand::rngs::StdRng;
-use rand::{RngExt, SeedableRng};
 
 use crate::distribution::Distribution;
 use crate::error::{Error, Result};
 use crate::kde::KernelDensityEstimator;
 use crate::param::ParamValue;
+use crate::rng_util;
 use crate::sampler::tpe::gamma::{FixedGamma, GammaStrategy};
 use crate::sampler::{CompletedTrial, Sampler};
 
@@ -131,7 +130,7 @@ pub struct TpeSampler {
     /// Optional fixed bandwidth for KDE. If None, uses Scott's rule.
     kde_bandwidth: Option<f64>,
     /// Thread-safe RNG for sampling.
-    rng: Mutex<StdRng>,
+    rng: Mutex<fastrand::Rng>,
 }
 
 impl TpeSampler {
@@ -149,7 +148,7 @@ impl TpeSampler {
             n_startup_trials: 10,
             n_ei_candidates: 24,
             kde_bandwidth: None,
-            rng: Mutex::new(rand::make_rng()),
+            rng: Mutex::new(fastrand::Rng::new()),
         }
     }
 
@@ -250,8 +249,8 @@ impl TpeSampler {
         }
 
         let rng = match seed {
-            Some(s) => StdRng::seed_from_u64(s),
-            None => rand::make_rng(),
+            Some(s) => fastrand::Rng::with_seed(s),
+            None => fastrand::Rng::new(),
         };
 
         Ok(Self {
@@ -328,19 +327,19 @@ impl TpeSampler {
         clippy::cast_precision_loss,
         clippy::unused_self
     )]
-    fn sample_uniform(&self, distribution: &Distribution, rng: &mut StdRng) -> ParamValue {
+    fn sample_uniform(&self, distribution: &Distribution, rng: &mut fastrand::Rng) -> ParamValue {
         match distribution {
             Distribution::Float(d) => {
                 let value = if d.log_scale {
                     let log_low = d.low.ln();
                     let log_high = d.high.ln();
-                    rng.random_range(log_low..=log_high).exp()
+                    rng_util::f64_range(rng, log_low, log_high).exp()
                 } else if let Some(step) = d.step {
                     let n_steps = ((d.high - d.low) / step).floor() as i64;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + (k as f64) * step
                 } else {
-                    rng.random_range(d.low..=d.high)
+                    rng_util::f64_range(rng, d.low, d.high)
                 };
                 ParamValue::Float(value)
             }
@@ -348,20 +347,18 @@ impl TpeSampler {
                 let value = if d.log_scale {
                     let log_low = (d.low as f64).ln();
                     let log_high = (d.high as f64).ln();
-                    let raw = rng.random_range(log_low..=log_high).exp().round() as i64;
+                    let raw = rng_util::f64_range(rng, log_low, log_high).exp().round() as i64;
                     raw.clamp(d.low, d.high)
                 } else if let Some(step) = d.step {
                     let n_steps = (d.high - d.low) / step;
-                    let k = rng.random_range(0..=n_steps);
+                    let k = rng.i64(0..=n_steps);
                     d.low + k * step
                 } else {
-                    rng.random_range(d.low..=d.high)
+                    rng.i64(d.low..=d.high)
                 };
                 ParamValue::Int(value)
             }
-            Distribution::Categorical(d) => {
-                ParamValue::Categorical(rng.random_range(0..d.n_choices))
-            }
+            Distribution::Categorical(d) => ParamValue::Categorical(rng.usize(0..d.n_choices)),
         }
     }
 
@@ -375,7 +372,7 @@ impl TpeSampler {
         step: Option<f64>,
         good_values: Vec<f64>,
         bad_values: Vec<f64>,
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> f64 {
         // Transform to internal space (log space if needed)
         let (internal_low, internal_high, good_internal, bad_internal) = if log_scale {
@@ -400,7 +397,7 @@ impl TpeSampler {
 
         // If KDE construction fails, fall back to uniform sampling
         let (Ok(l_kde), Ok(g_kde)) = (l_kde, g_kde) else {
-            return rng.random_range(low..=high);
+            return rng_util::f64_range(rng, low, high);
         };
 
         // Generate candidates from l(x) and select the one with best l(x)/g(x) ratio
@@ -464,7 +461,7 @@ impl TpeSampler {
         step: Option<i64>,
         good_values: &[i64],
         bad_values: &[i64],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> i64 {
         // Convert to floats for KDE
         let good_floats: Vec<f64> = good_values.iter().map(|&v| v as f64).collect();
@@ -503,7 +500,7 @@ impl TpeSampler {
         n_choices: usize,
         good_indices: &[usize],
         bad_indices: &[usize],
-        rng: &mut StdRng,
+        rng: &mut fastrand::Rng,
     ) -> usize {
         // Count occurrences in good and bad groups
         let mut good_counts = vec![0usize; n_choices];
@@ -534,7 +531,7 @@ impl TpeSampler {
 
         // Sample proportionally to weights
         let total_weight: f64 = weights.iter().sum();
-        let threshold = rng.random::<f64>() * total_weight;
+        let threshold = rng.f64() * total_weight;
 
         let mut cumulative = 0.0;
         for (i, &w) in weights.iter().enumerate() {
@@ -856,8 +853,8 @@ impl TpeSamplerBuilder {
         }
 
         let rng = match self.seed {
-            Some(s) => StdRng::seed_from_u64(s),
-            None => rand::make_rng(),
+            Some(s) => fastrand::Rng::with_seed(s),
+            None => fastrand::Rng::new(),
         };
 
         Ok(TpeSampler {

src/storage/journal.rs

@@ -0,0 +1,176 @@
+//! JSONL-based journal storage backend.
+
+use core::marker::PhantomData;
+use std::fs::{File, OpenOptions};
+use std::io::{BufRead, BufReader, Seek, SeekFrom, Write};
+use std::path::{Path, PathBuf};
+use std::sync::Arc;
+
+use fs2::FileExt;
+use parking_lot::{Mutex, RwLock};
+use serde::Serialize;
+use serde::de::DeserializeOwned;
+
+use super::{MemoryStorage, Storage};
+use crate::sampler::CompletedTrial;
+
+/// A storage backend that appends completed trials as JSON lines to a file.
+///
+/// Trials are kept in memory for fast read access and simultaneously
+/// persisted to a JSONL file.  Multiple processes can safely share
+/// the same file: writes use an exclusive file lock, reads use a
+/// shared file lock.
+///
+/// The type parameter `V` is the objective value type (typically `f64`).
+/// It must be serializable so that trials can be written to disk.
+///
+/// # Examples
+///
+/// ```no_run
+/// use optimizer::storage::JournalStorage;
+///
+/// let storage: JournalStorage<f64> = JournalStorage::new("trials.jsonl");
+/// ```
+pub struct JournalStorage<V = f64> {
+    memory: MemoryStorage<V>,
+    path: PathBuf,
+    /// Serialise in-process writes so we only hold the file lock briefly.
+    write_lock: Mutex<()>,
+    _marker: PhantomData<V>,
+}
+
+impl<V: Serialize + DeserializeOwned + Send + Sync> JournalStorage<V> {
+    /// Creates a new journal storage that writes to the given path.
+    ///
+    /// The file does not need to exist yet — it will be created on the
+    /// first write.  Existing trials in the file are **not** loaded
+    /// until [`refresh`](Storage::refresh) is called (which happens
+    /// automatically at the start of each trial via the [`Study`](crate::Study)).
+    ///
+    /// To pre-load existing trials, use [`JournalStorage::open`].
+    #[must_use]
+    pub fn new(path: impl AsRef<Path>) -> Self {
+        Self {
+            memory: MemoryStorage::new(),
+            path: path.as_ref().to_path_buf(),
+            write_lock: Mutex::new(()),
+            _marker: PhantomData,
+        }
+    }
+
+    /// Opens an existing journal file and loads all stored trials.
+    ///
+    /// If the file does not exist, returns an empty storage (no error).
+    ///
+    /// # Errors
+    ///
+    /// Returns a [`Storage`](crate::Error::Storage) error if the file
+    /// exists but cannot be read or parsed.
+    pub fn open(path: impl AsRef<Path>) -> crate::Result<Self> {
+        let path = path.as_ref().to_path_buf();
+        let trials = load_trials_from_file(&path)?;
+        Ok(Self {
+            memory: MemoryStorage::with_trials(trials),
+            path,
+            write_lock: Mutex::new(()),
+            _marker: PhantomData,
+        })
+    }
+
+    /// Append a single trial to the JSONL file (best-effort).
+    fn write_to_file(&self, trial: &CompletedTrial<V>) -> crate::Result<()> {
+        let _guard = self.write_lock.lock();
+
+        let mut file = OpenOptions::new()
+            .create(true)
+            .truncate(false)
+            .read(true)
+            .write(true)
+            .open(&self.path)
+            .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+        file.lock_exclusive()
+            .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+        file.seek(SeekFrom::End(0))
+            .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+        let line =
+            serde_json::to_string(trial).map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+        writeln!(file, "{line}").map_err(|e| crate::Error::Storage(e.to_string()))?;
+        file.sync_data()
+            .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+        file.unlock()
+            .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+        Ok(())
+    }
+}
+
+impl<V: Serialize + DeserializeOwned + Send + Sync> Storage<V> for JournalStorage<V> {
+    fn push(&self, trial: CompletedTrial<V>) {
+        // Best-effort persist; the trial stays in memory regardless.
+        let _ = self.write_to_file(&trial);
+        self.memory.push(trial);
+    }
+
+    fn trials_arc(&self) -> &Arc<RwLock<Vec<CompletedTrial<V>>>> {
+        self.memory.trials_arc()
+    }
+
+    fn next_trial_id(&self) -> u64 {
+        self.memory.next_trial_id()
+    }
+
+    fn refresh(&self) -> bool {
+        let Ok(loaded) = load_trials_from_file::<V>(&self.path) else {
+            return false;
+        };
+        let mut guard = self.memory.trials_arc().write();
+        if loaded.len() > guard.len() {
+            if let Some(max_id) = loaded.iter().map(|t| t.id).max() {
+                self.memory.bump_next_id(max_id + 1);
+            }
+            *guard = loaded;
+            true
+        } else {
+            false
+        }
+    }
+}
+
+/// Read all trials from a JSONL file.  Returns an empty vec if the
+/// file does not exist.
+fn load_trials_from_file<V: DeserializeOwned>(
+    path: &Path,
+) -> crate::Result<Vec<CompletedTrial<V>>> {
+    let file = match File::open(path) {
+        Ok(f) => f,
+        Err(e) if e.kind() == std::io::ErrorKind::NotFound => return Ok(Vec::new()),
+        Err(e) => return Err(crate::Error::Storage(e.to_string())),
+    };
+
+    file.lock_shared()
+        .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+    let reader = BufReader::new(&file);
+    let mut trials = Vec::new();
+
+    for line in reader.lines() {
+        let line = line.map_err(|e| crate::Error::Storage(e.to_string()))?;
+        let line = line.trim();
+        if line.is_empty() {
+            continue;
+        }
+        let trial: CompletedTrial<V> =
+            serde_json::from_str(line).map_err(|e| crate::Error::Storage(e.to_string()))?;
+        trials.push(trial);
+    }
+
+    file.unlock()
+        .map_err(|e| crate::Error::Storage(e.to_string()))?;
+
+    Ok(trials)
+}

src/storage/memory.rs

@@ -0,0 +1,61 @@
+use core::sync::atomic::{AtomicU64, Ordering};
+use std::sync::Arc;
+
+use parking_lot::RwLock;
+
+use super::Storage;
+use crate::sampler::CompletedTrial;
+
+/// In-memory trial storage (the default).
+///
+/// This is a thin wrapper around `Arc<RwLock<Vec<CompletedTrial<V>>>>`.
+pub struct MemoryStorage<V> {
+    trials: Arc<RwLock<Vec<CompletedTrial<V>>>>,
+    next_id: AtomicU64,
+}
+
+impl<V> MemoryStorage<V> {
+    /// Creates a new, empty in-memory store.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            trials: Arc::new(RwLock::new(Vec::new())),
+            next_id: AtomicU64::new(0),
+        }
+    }
+
+    /// Creates an in-memory store pre-populated with `trials`.
+    #[must_use]
+    pub fn with_trials(trials: Vec<CompletedTrial<V>>) -> Self {
+        let next_id = trials.iter().map(|t| t.id).max().map_or(0, |id| id + 1);
+        Self {
+            trials: Arc::new(RwLock::new(trials)),
+            next_id: AtomicU64::new(next_id),
+        }
+    }
+
+    /// Ensures the ID counter is at least `min_value`.
+    pub(crate) fn bump_next_id(&self, min_value: u64) {
+        self.next_id.fetch_max(min_value, Ordering::SeqCst);
+    }
+}
+
+impl<V> Default for MemoryStorage<V> {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl<V: Send + Sync> Storage<V> for MemoryStorage<V> {
+    fn push(&self, trial: CompletedTrial<V>) {
+        self.trials.write().push(trial);
+    }
+
+    fn trials_arc(&self) -> &Arc<RwLock<Vec<CompletedTrial<V>>>> {
+        &self.trials
+    }
+
+    fn next_trial_id(&self) -> u64 {
+        self.next_id.fetch_add(1, Ordering::SeqCst)
+    }
+}

src/storage/mod.rs

@@ -0,0 +1,54 @@
+//! Trial storage backends.
+//!
+//! The [`Storage`] trait defines how completed trials are stored and
+//! accessed.  [`MemoryStorage`] keeps trials in memory (the default).
+//! With the `journal` feature enabled, [`JournalStorage`] appends
+//! trials to a JSONL file with file-level locking so multiple
+//! processes can safely share state.
+
+#[cfg(feature = "journal")]
+mod journal;
+
+use std::sync::Arc;
+
+#[cfg(feature = "journal")]
+pub use journal::JournalStorage;
+use parking_lot::RwLock;
+
+mod memory;
+pub use memory::MemoryStorage;
+
+use crate::sampler::CompletedTrial;
+
+/// Trait for storing and retrieving completed trials.
+///
+/// Every [`Study`](crate::Study) owns an `Arc<dyn Storage<V>>`.  The
+/// default implementation is [`MemoryStorage`], which keeps trials in
+/// a plain `Vec` behind a read-write lock.
+///
+/// Implementations must be safe to use from multiple threads.
+pub trait Storage<V>: Send + Sync {
+    /// Append a completed trial to the store.
+    fn push(&self, trial: CompletedTrial<V>);
+
+    /// Return a reference to the in-memory trial buffer.
+    ///
+    /// All implementations must maintain an `Arc<RwLock<Vec<…>>>` that
+    /// reflects the current set of trials.  Callers may acquire a read
+    /// lock for efficient, allocation-free access.
+    fn trials_arc(&self) -> &Arc<RwLock<Vec<CompletedTrial<V>>>>;
+
+    /// Atomically returns the next unique trial ID.
+    ///
+    /// Each call increments an internal counter so that consecutive
+    /// calls always produce distinct IDs.
+    fn next_trial_id(&self) -> u64;
+
+    /// Reload from an external source (e.g. a file written by another
+    /// process).  Returns `true` if the in-memory buffer was updated.
+    ///
+    /// The default implementation is a no-op that returns `false`.
+    fn refresh(&self) -> bool {
+        false
+    }
+}

src/study.rs

@@ -4,8 +4,8 @@ use core::any::Any;
 use core::fmt;
 #[cfg(feature = "async")]
 use core::future::Future;
+use core::marker::PhantomData;
 use core::ops::ControlFlow;
-use core::sync::atomic::{AtomicU64, Ordering};
 use core::time::Duration;
 use std::collections::{HashMap, VecDeque};
 use std::sync::Arc;
@@ -49,10 +49,8 @@ where
     sampler: Arc<dyn Sampler>,
     /// The pruner used to decide whether to stop trials early.
     pruner: Arc<dyn Pruner>,
-    /// Completed trials (wrapped in Arc for sharing with Trial).
-    completed_trials: Arc<RwLock<Vec<CompletedTrial<V>>>>,
-    /// Counter for generating unique trial IDs.
-    next_trial_id: AtomicU64,
+    /// Trial storage backend (default: [`MemoryStorage`](crate::storage::MemoryStorage)).
+    storage: Arc<dyn crate::storage::Storage<V>>,
     /// Optional factory for creating sampler-aware trials.
     /// Set automatically for `Study<f64>` so that `create_trial()` and all
     /// optimization methods use the sampler without requiring `_with_sampler` suffixes.
@@ -84,11 +82,35 @@ where
     #[must_use]
     pub fn new(direction: Direction) -> Self
     where
-        V: 'static,
+        V: Send + Sync + 'static,
     {
         Self::with_sampler(direction, RandomSampler::new())
     }
 
+    /// Returns a [`StudyBuilder`] for constructing a study with a fluent API.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use optimizer::prelude::*;
+    ///
+    /// let study: Study<f64> = Study::builder()
+    ///     .minimize()
+    ///     .sampler(TpeSampler::new())
+    ///     .pruner(NopPruner)
+    ///     .build();
+    /// ```
+    #[must_use]
+    pub fn builder() -> StudyBuilder<V> {
+        StudyBuilder {
+            direction: Direction::Minimize,
+            sampler: None,
+            pruner: None,
+            storage: None,
+            _marker: PhantomData,
+        }
+    }
+
     /// Creates a study that minimizes the objective value.
     ///
     /// This is a shorthand for `Study::with_sampler(Direction::Minimize, sampler)`.
@@ -109,7 +131,7 @@ where
     #[must_use]
     pub fn minimize(sampler: impl Sampler + 'static) -> Self
     where
-        V: 'static,
+        V: Send + Sync + 'static,
     {
         Self::with_sampler(Direction::Minimize, sampler)
     }
@@ -134,7 +156,7 @@ where
     #[must_use]
     pub fn maximize(sampler: impl Sampler + 'static) -> Self
     where
-        V: 'static,
+        V: Send + Sync + 'static,
     {
         Self::with_sampler(Direction::Maximize, sampler)
     }
@@ -158,40 +180,28 @@ where
     /// ```
     pub fn with_sampler(direction: Direction, sampler: impl Sampler + 'static) -> Self
     where
-        V: 'static,
+        V: Send + Sync + 'static,
     {
-        let sampler: Arc<dyn Sampler> = Arc::new(sampler);
-        let completed_trials = Arc::new(RwLock::new(Vec::new()));
-
-        let pruner: Arc<dyn Pruner> = Arc::new(NopPruner);
-
-        // For Study<f64>, set up a trial factory that provides sampler integration.
-        // This uses Any downcasting to check at runtime whether V = f64.
-        let trial_factory = Self::make_trial_factory(&sampler, &completed_trials, &pruner);
-
-        Self {
+        Self::with_sampler_and_storage(
             direction,
             sampler,
-            pruner,
-            completed_trials,
-            next_trial_id: AtomicU64::new(0),
-            trial_factory,
-            enqueued_params: Arc::new(Mutex::new(VecDeque::new())),
-        }
+            crate::storage::MemoryStorage::<V>::new(),
+        )
     }
 
     /// Builds a trial factory for sampler integration when `V = f64`.
     fn make_trial_factory(
         sampler: &Arc<dyn Sampler>,
-        completed_trials: &Arc<RwLock<Vec<CompletedTrial<V>>>>,
+        storage: &Arc<dyn crate::storage::Storage<V>>,
         pruner: &Arc<dyn Pruner>,
     ) -> Option<Arc<dyn Fn(u64) -> Trial + Send + Sync>>
     where
         V: 'static,
     {
-        // Try to downcast the completed_trials Arc to the f64 specialization.
+        // Try to downcast the storage's trial buffer to the f64 specialization.
         // This succeeds only when V = f64, enabling automatic sampler integration.
-        let any_ref: &dyn Any = completed_trials;
+        let trials_arc = storage.trials_arc();
+        let any_ref: &dyn Any = trials_arc;
         let f64_trials: Option<&Arc<RwLock<Vec<CompletedTrial<f64>>>>> = any_ref.downcast_ref();
 
         f64_trials.map(|trials| {
@@ -210,7 +220,35 @@ where
         })
     }
 
+    /// Creates a study with a custom sampler and storage backend.
+    ///
+    /// This is the most general constructor — all other constructors
+    /// delegate to this one.
+    pub fn with_sampler_and_storage(
+        direction: Direction,
+        sampler: impl Sampler + 'static,
+        storage: impl crate::storage::Storage<V> + 'static,
+    ) -> Self
+    where
+        V: 'static,
+    {
+        let sampler: Arc<dyn Sampler> = Arc::new(sampler);
+        let pruner: Arc<dyn Pruner> = Arc::new(NopPruner);
+        let storage: Arc<dyn crate::storage::Storage<V>> = Arc::new(storage);
+        let trial_factory = Self::make_trial_factory(&sampler, &storage, &pruner);
+
+        Self {
+            direction,
+            sampler,
+            pruner,
+            storage,
+            trial_factory,
+            enqueued_params: Arc::new(Mutex::new(VecDeque::new())),
+        }
+    }
+
     /// Returns the optimization direction.
+    #[must_use]
     pub fn direction(&self) -> Direction {
         self.direction
     }
@@ -254,19 +292,19 @@ where
         pruner: impl Pruner + 'static,
     ) -> Self
     where
-        V: 'static,
+        V: Send + Sync + 'static,
     {
         let sampler: Arc<dyn Sampler> = Arc::new(sampler);
         let pruner: Arc<dyn Pruner> = Arc::new(pruner);
-        let completed_trials = Arc::new(RwLock::new(Vec::new()));
-        let trial_factory = Self::make_trial_factory(&sampler, &completed_trials, &pruner);
+        let storage: Arc<dyn crate::storage::Storage<V>> =
+            Arc::new(crate::storage::MemoryStorage::<V>::new());
+        let trial_factory = Self::make_trial_factory(&sampler, &storage, &pruner);
 
         Self {
             direction,
             sampler,
             pruner,
-            completed_trials,
-            next_trial_id: AtomicU64::new(0),
+            storage,
             trial_factory,
             enqueued_params: Arc::new(Mutex::new(VecDeque::new())),
         }
@@ -277,8 +315,7 @@ where
         V: 'static,
     {
         self.sampler = Arc::new(sampler);
-        self.trial_factory =
-            Self::make_trial_factory(&self.sampler, &self.completed_trials, &self.pruner);
+        self.trial_factory = Self::make_trial_factory(&self.sampler, &self.storage, &self.pruner);
     }
 
     /// Sets a new pruner for the study.
@@ -291,11 +328,11 @@ where
         V: 'static,
     {
         self.pruner = Arc::new(pruner);
-        self.trial_factory =
-            Self::make_trial_factory(&self.sampler, &self.completed_trials, &self.pruner);
+        self.trial_factory = Self::make_trial_factory(&self.sampler, &self.storage, &self.pruner);
     }
 
     /// Returns a reference to the study's pruner.
+    #[must_use]
     pub fn pruner(&self) -> &dyn Pruner {
         &*self.pruner
     }
@@ -375,7 +412,7 @@ where
 
     /// Generates the next unique trial ID.
     pub(crate) fn next_trial_id(&self) -> u64 {
-        self.next_trial_id.fetch_add(1, Ordering::SeqCst)
+        self.storage.next_trial_id()
     }
 
     /// Creates a new trial with a unique ID.
@@ -400,7 +437,10 @@ where
     /// let trial2 = study.create_trial();
     /// assert_eq!(trial2.id(), 1);
     /// ```
+    #[must_use]
     pub fn create_trial(&self) -> Trial {
+        self.storage.refresh();
+
         let id = self.next_trial_id();
         let mut trial = if let Some(factory) = &self.trial_factory {
             factory(id)
@@ -455,7 +495,8 @@ where
         );
         completed.state = TrialState::Complete;
         completed.constraints = trial.constraint_values().to_vec();
-        self.completed_trials.write().push(completed);
+
+        self.storage.push(completed);
     }
 
     /// Records a failed trial with an error message.
@@ -509,6 +550,7 @@ where
     /// let value = x_val * x_val;
     /// study.tell(trial, Ok::<_, &str>(value));
     /// ```
+    #[must_use]
     pub fn ask(&self) -> Trial {
         self.create_trial()
     }
@@ -565,7 +607,8 @@ where
         );
         completed.state = TrialState::Pruned;
         completed.constraints = trial.constraint_values().to_vec();
-        self.completed_trials.write().push(completed);
+
+        self.storage.push(completed);
     }
 
     /// Returns an iterator over all completed trials.
@@ -592,11 +635,12 @@ where
     ///     println!("Trial {} has value {:?}", completed.id, completed.value);
     /// }
     /// ```
+    #[must_use]
     pub fn trials(&self) -> Vec<CompletedTrial<V>>
     where
         V: Clone,
     {
-        self.completed_trials.read().clone()
+        self.storage.trials_arc().read().clone()
     }
 
     /// Returns the number of completed trials.
@@ -618,13 +662,16 @@ where
     /// study.complete_trial(trial, 0.5);
     /// assert_eq!(study.n_trials(), 1);
     /// ```
+    #[must_use]
     pub fn n_trials(&self) -> usize {
-        self.completed_trials.read().len()
+        self.storage.trials_arc().read().len()
     }
 
     /// Returns the number of pruned trials.
+    #[must_use]
     pub fn n_pruned_trials(&self) -> usize {
-        self.completed_trials
+        self.storage
+            .trials_arc()
             .read()
             .iter()
             .filter(|t| t.state == TrialState::Pruned)
@@ -703,7 +750,7 @@ where
     where
         V: Clone,
     {
-        let trials = self.completed_trials.read();
+        let trials = self.storage.trials_arc().read();
         let direction = self.direction;
 
         let best = trials
@@ -763,11 +810,12 @@ where
     /// Only includes completed trials (not failed or pruned).
     ///
     /// If fewer than `n` completed trials exist, returns all of them.
+    #[must_use]
     pub fn top_trials(&self, n: usize) -> Vec<CompletedTrial<V>>
     where
         V: Clone,
     {
-        let trials = self.completed_trials.read();
+        let trials = self.storage.trials_arc().read();
         let direction = self.direction;
         let mut completed: Vec<_> = trials
             .iter()
@@ -848,7 +896,7 @@ where
                     #[cfg(feature = "tracing")]
                     {
                         tracing::info!(trial_id, "trial completed");
-                        let trials = self.completed_trials.read();
+                        let trials = self.storage.trials_arc().read();
                         if trials
                             .iter()
                             .filter(|t| t.state == TrialState::Complete)
@@ -876,7 +924,8 @@ where
 
         // Return error if no trials completed successfully
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -975,7 +1024,8 @@ where
 
         // Return error if no trials completed successfully
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1101,7 +1151,8 @@ where
 
         // Return error if no trials completed successfully
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1197,7 +1248,7 @@ where
                     #[cfg(feature = "tracing")]
                     {
                         tracing::info!(trial_id, "trial completed");
-                        let trials = self.completed_trials.read();
+                        let trials = self.storage.trials_arc().read();
                         if trials
                             .iter()
                             .filter(|t| t.state == TrialState::Complete)
@@ -1210,7 +1261,7 @@ where
                     }
 
                     // Get the just-completed trial for the callback
-                    let trials = self.completed_trials.read();
+                    let trials = self.storage.trials_arc().read();
                     let Some(completed) = trials.last() else {
                         return Err(crate::Error::Internal(
                             "completed trial not found after adding",
@@ -1243,7 +1294,8 @@ where
 
         // Return error if no trials completed successfully
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1328,7 +1380,8 @@ where
         }
 
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1420,7 +1473,7 @@ where
                     #[cfg(feature = "tracing")]
                     {
                         tracing::info!(trial_id, "trial completed");
-                        let trials = self.completed_trials.read();
+                        let trials = self.storage.trials_arc().read();
                         if trials
                             .iter()
                             .filter(|t| t.state == TrialState::Complete)
@@ -1432,7 +1485,7 @@ where
                         }
                     }
 
-                    let trials = self.completed_trials.read();
+                    let trials = self.storage.trials_arc().read();
                     let Some(completed) = trials.last() else {
                         return Err(crate::Error::Internal(
                             "completed trial not found after adding",
@@ -1461,7 +1514,8 @@ where
         }
 
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1519,7 +1573,8 @@ where
         }
 
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1608,7 +1663,8 @@ where
         }
 
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1716,7 +1772,8 @@ where
 
         // Return error if no trials completed successfully
         let has_complete = self
-            .completed_trials
+            .storage
+            .trials_arc()
             .read()
             .iter()
             .any(|t| t.state == TrialState::Complete);
@@ -1746,7 +1803,7 @@ where
     pub fn to_csv(&self, mut writer: impl std::io::Write) -> std::io::Result<()> {
         use std::collections::BTreeMap;
 
-        let trials = self.completed_trials.read();
+        let trials = self.storage.trials_arc().read();
 
         // Collect all unique parameter labels (sorted for deterministic column order).
         let mut param_columns: BTreeMap<ParamId, String> = BTreeMap::new();
@@ -1873,7 +1930,7 @@ where
     pub fn summary(&self) -> String {
         use fmt::Write;
 
-        let trials = self.completed_trials.read();
+        let trials = self.storage.trials_arc().read();
         let n_complete = trials
             .iter()
             .filter(|t| t.state == TrialState::Complete)
@@ -1920,6 +1977,7 @@ where
     ///
     /// This clones the internal trial list, so it is suitable for
     /// analysis and iteration but not for hot paths.
+    #[must_use]
     pub fn iter(&self) -> std::vec::IntoIter<CompletedTrial<V>> {
         self.trials().into_iter()
     }
@@ -1967,7 +2025,7 @@ where
         use crate::param::ParamValue;
         use crate::types::TrialState;
 
-        let trials = self.completed_trials.read();
+        let trials = self.storage.trials_arc().read();
         let complete: Vec<_> = trials
             .iter()
             .filter(|t| t.state == TrialState::Complete)
@@ -2039,7 +2097,6 @@ where
     /// # Errors
     ///
     /// Returns [`crate::Error::NoCompletedTrials`] if fewer than 2 trials have completed.
-    #[cfg(feature = "fanova")]
     pub fn fanova(&self) -> crate::Result<crate::fanova::FanovaResult> {
         self.fanova_with_config(&crate::fanova::FanovaConfig::default())
     }
@@ -2052,7 +2109,6 @@ where
     /// # Errors
     ///
     /// Returns [`crate::Error::NoCompletedTrials`] if fewer than 2 trials have completed.
-    #[cfg(feature = "fanova")]
     #[allow(clippy::cast_precision_loss)]
     pub fn fanova_with_config(
         &self,
@@ -2064,7 +2120,7 @@ where
         use crate::param::ParamValue;
         use crate::types::TrialState;
 
-        let trials = self.completed_trials.read();
+        let trials = self.storage.trials_arc().read();
         let complete: Vec<_> = trials
             .iter()
             .filter(|t| t.state == TrialState::Complete)
@@ -2169,6 +2225,7 @@ impl Study<f64> {
         since = "0.2.0",
         note = "use `create_trial()` instead — it now uses the sampler automatically for Study<f64>"
     )]
+    #[must_use]
     pub fn create_trial_with_sampler(&self) -> Trial {
         self.create_trial()
     }
@@ -2258,7 +2315,160 @@ impl Study<f64> {
     }
 }
 
-#[cfg(feature = "visualization")]
+impl<V: PartialOrd + Send + Sync + 'static> Study<V> {
+    /// Creates a study with a custom sampler, pruner, and storage backend.
+    pub fn with_sampler_pruner_and_storage(
+        direction: Direction,
+        sampler: impl Sampler + 'static,
+        pruner: impl Pruner + 'static,
+        storage: impl crate::storage::Storage<V> + 'static,
+    ) -> Self {
+        let sampler: Arc<dyn Sampler> = Arc::new(sampler);
+        let pruner: Arc<dyn Pruner> = Arc::new(pruner);
+        let storage: Arc<dyn crate::storage::Storage<V>> = Arc::new(storage);
+        let trial_factory = Self::make_trial_factory(&sampler, &storage, &pruner);
+
+        Self {
+            direction,
+            sampler,
+            pruner,
+            storage,
+            trial_factory,
+            enqueued_params: Arc::new(Mutex::new(VecDeque::new())),
+        }
+    }
+}
+
+/// A builder for constructing [`Study`] instances with a fluent API.
+///
+/// Created via [`Study::builder()`]. Collects sampler, pruner, direction,
+/// and storage options before constructing the study.
+///
+/// # Defaults
+///
+/// - Direction: [`Minimize`](Direction::Minimize)
+/// - Sampler: [`RandomSampler`]
+/// - Pruner: [`NopPruner`]
+/// - Storage: [`MemoryStorage`](crate::storage::MemoryStorage)
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::prelude::*;
+///
+/// let study: Study<f64> = Study::builder()
+///     .maximize()
+///     .sampler(TpeSampler::new())
+///     .pruner(MedianPruner::new(Direction::Maximize).n_warmup_steps(5))
+///     .build();
+///
+/// assert_eq!(study.direction(), Direction::Maximize);
+/// ```
+pub struct StudyBuilder<V: PartialOrd = f64> {
+    direction: Direction,
+    sampler: Option<Box<dyn Sampler>>,
+    pruner: Option<Box<dyn Pruner>>,
+    storage: Option<Box<dyn crate::storage::Storage<V>>>,
+    _marker: PhantomData<V>,
+}
+
+impl<V: PartialOrd> StudyBuilder<V> {
+    /// Sets the optimization direction to minimize.
+    #[must_use]
+    pub fn minimize(mut self) -> Self {
+        self.direction = Direction::Minimize;
+        self
+    }
+
+    /// Sets the optimization direction to maximize.
+    #[must_use]
+    pub fn maximize(mut self) -> Self {
+        self.direction = Direction::Maximize;
+        self
+    }
+
+    /// Sets the optimization direction.
+    #[must_use]
+    pub fn direction(mut self, direction: Direction) -> Self {
+        self.direction = direction;
+        self
+    }
+
+    /// Sets the sampler used for parameter suggestions.
+    #[must_use]
+    pub fn sampler(mut self, sampler: impl Sampler + 'static) -> Self {
+        self.sampler = Some(Box::new(sampler));
+        self
+    }
+
+    /// Sets the pruner used for early stopping of trials.
+    #[must_use]
+    pub fn pruner(mut self, pruner: impl Pruner + 'static) -> Self {
+        self.pruner = Some(Box::new(pruner));
+        self
+    }
+
+    /// Sets a custom storage backend.
+    #[must_use]
+    pub fn storage(mut self, storage: impl crate::storage::Storage<V> + 'static) -> Self {
+        self.storage = Some(Box::new(storage));
+        self
+    }
+
+    /// Builds the [`Study`] with the configured options.
+    #[must_use]
+    pub fn build(self) -> Study<V>
+    where
+        V: Send + Sync + 'static,
+    {
+        let sampler = self
+            .sampler
+            .unwrap_or_else(|| Box::new(RandomSampler::new()));
+        let pruner = self.pruner.unwrap_or_else(|| Box::new(NopPruner));
+        let storage = self
+            .storage
+            .unwrap_or_else(|| Box::new(crate::storage::MemoryStorage::<V>::new()));
+
+        let sampler: Arc<dyn Sampler> = Arc::from(sampler);
+        let pruner: Arc<dyn Pruner> = Arc::from(pruner);
+        let storage: Arc<dyn crate::storage::Storage<V>> = Arc::from(storage);
+        let trial_factory = Study::make_trial_factory(&sampler, &storage, &pruner);
+
+        Study {
+            direction: self.direction,
+            sampler,
+            pruner,
+            storage,
+            trial_factory,
+            enqueued_params: Arc::new(Mutex::new(VecDeque::new())),
+        }
+    }
+}
+
+#[cfg(feature = "journal")]
+impl<V> Study<V>
+where
+    V: PartialOrd + Send + Sync + serde::Serialize + serde::de::DeserializeOwned + 'static,
+{
+    /// Creates a study backed by a JSONL journal file.
+    ///
+    /// Any existing trials in the file are loaded into memory and the
+    /// trial ID counter is set to one past the highest stored ID. New
+    /// trials are written through to the file on completion.
+    ///
+    /// # Errors
+    ///
+    /// Returns a [`Storage`](crate::Error::Storage) error if loading fails.
+    pub fn with_journal(
+        direction: Direction,
+        sampler: impl Sampler + 'static,
+        path: impl AsRef<std::path::Path>,
+    ) -> crate::Result<Self> {
+        let storage = crate::storage::JournalStorage::<V>::open(path)?;
+        Ok(Self::with_sampler_and_storage(direction, sampler, storage))
+    }
+}
+
 impl Study<f64> {
     /// Generates an HTML report with interactive Plotly.js charts.
     ///
@@ -2326,11 +2536,13 @@ impl<V: PartialOrd + Clone + serde::Serialize> Study<V> {
     /// Returns an I/O error if the file cannot be created or written.
     pub fn save(&self, path: impl AsRef<std::path::Path>) -> std::io::Result<()> {
         let path = path.as_ref();
+        let trials = self.trials();
+        let next_trial_id = trials.iter().map(|t| t.id).max().map_or(0, |id| id + 1);
         let snapshot = StudySnapshot {
             version: 1,
             direction: self.direction,
-            trials: self.trials(),
-            next_trial_id: self.next_trial_id.load(Ordering::Relaxed),
+            trials,
+            next_trial_id,
             metadata: HashMap::new(),
         };
 
@@ -2382,7 +2594,7 @@ impl<V: PartialOrd + Clone + Default + serde::Serialize> Study<V> {
 }
 
 #[cfg(feature = "serde")]
-impl<V: PartialOrd + Clone + serde::de::DeserializeOwned + 'static> Study<V> {
+impl<V: PartialOrd + Send + Sync + Clone + serde::de::DeserializeOwned + 'static> Study<V> {
     /// Loads a study from a JSON file.
     ///
     /// The loaded study uses a `RandomSampler` by default. Call
@@ -2396,12 +2608,12 @@ impl<V: PartialOrd + Clone + serde::de::DeserializeOwned + 'static> Study<V> {
         let file = std::fs::File::open(path)?;
         let snapshot: StudySnapshot<V> = serde_json::from_reader(file)
             .map_err(|e| std::io::Error::new(std::io::ErrorKind::InvalidData, e))?;
-        let study = Study::new(snapshot.direction);
-        *study.completed_trials.write() = snapshot.trials;
-        study
-            .next_trial_id
-            .store(snapshot.next_trial_id, Ordering::Relaxed);
-        Ok(study)
+        let storage = crate::storage::MemoryStorage::with_trials(snapshot.trials);
+        Ok(Self::with_sampler_and_storage(
+            snapshot.direction,
+            RandomSampler::new(),
+            storage,
+        ))
     }
 }
 

tests/differential_evolution_tests.rs

@@ -0,0 +1,350 @@
+use optimizer::prelude::*;
+use optimizer::sampler::differential_evolution::{
+    DifferentialEvolutionSampler, DifferentialEvolutionStrategy,
+};
+
+#[test]
+fn sphere_function() {
+    let sampler = DifferentialEvolutionSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(200, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 1.0,
+        "sphere best value should be < 1.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn rosenbrock_function() {
+    let sampler = DifferentialEvolutionSampler::builder()
+        .population_size(20)
+        .seed(42)
+        .build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(400, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            let val = (1.0 - xv).powi(2) + 100.0 * (yv - xv * xv).powi(2);
+            Ok::<_, Error>(val)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 50.0,
+        "rosenbrock best value should be < 50.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn rastrigin_function() {
+    let sampler = DifferentialEvolutionSampler::builder()
+        .population_size(30)
+        .mutation_factor(0.7)
+        .crossover_rate(0.9)
+        .seed(42)
+        .build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.12, 5.12).name("x");
+    let y = FloatParam::new(-5.12, 5.12).name("y");
+
+    study
+        .optimize(500, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            let val = 20.0
+                + (xv * xv - 10.0 * (2.0 * std::f64::consts::PI * xv).cos())
+                + (yv * yv - 10.0 * (2.0 * std::f64::consts::PI * yv).cos());
+            Ok::<_, Error>(val)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    // Rastrigin is multimodal; just check reasonable convergence
+    assert!(
+        best.value < 10.0,
+        "rastrigin best value should be < 10.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn bounds_respected() {
+    let sampler = DifferentialEvolutionSampler::with_seed(123);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-2.0, 3.0).name("x");
+    let y = FloatParam::new(0.0, 10.0).name("y");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv + yv)
+        })
+        .unwrap();
+
+    for trial in study.trials() {
+        let xv: f64 = trial.get(&x).unwrap();
+        let yv: f64 = trial.get(&y).unwrap();
+        assert!((-2.0..=3.0).contains(&xv), "x = {xv} out of bounds [-2, 3]");
+        assert!((0.0..=10.0).contains(&yv), "y = {yv} out of bounds [0, 10]");
+    }
+}
+
+#[test]
+fn strategy_best1() {
+    let sampler = DifferentialEvolutionSampler::builder()
+        .strategy(DifferentialEvolutionStrategy::Best1)
+        .population_size(15)
+        .seed(42)
+        .build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(200, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 5.0,
+        "Best1 strategy should converge, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn strategy_current_to_best1() {
+    let sampler = DifferentialEvolutionSampler::builder()
+        .strategy(DifferentialEvolutionStrategy::CurrentToBest1)
+        .population_size(15)
+        .seed(42)
+        .build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(200, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 5.0,
+        "CurrentToBest1 strategy should converge, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn mixed_params_float_and_categorical() {
+    let sampler = DifferentialEvolutionSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let cat = CategoricalParam::new(vec!["a", "b", "c"]).name("cat");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let cv = cat.suggest(trial)?;
+            let penalty = match cv {
+                "a" => 0.0,
+                "b" => 1.0,
+                _ => 2.0,
+            };
+            Ok::<_, Error>(xv * xv + penalty)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 10.0,
+        "best value should be < 10.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn seeded_reproducibility() {
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    let run = |seed: u64| {
+        let sampler = DifferentialEvolutionSampler::with_seed(seed);
+        let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+        study
+            .optimize(50, |trial| {
+                let xv = x.suggest(trial)?;
+                let yv = y.suggest(trial)?;
+                Ok::<_, Error>(xv * xv + yv * yv)
+            })
+            .unwrap();
+        study.trials().iter().map(|t| t.value).collect::<Vec<_>>()
+    };
+
+    let results1 = run(42);
+    let results2 = run(42);
+    assert_eq!(results1, results2, "same seed should produce same results");
+}
+
+#[test]
+fn different_seeds_different_results() {
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    let run = |seed: u64| {
+        let sampler = DifferentialEvolutionSampler::with_seed(seed);
+        let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+        study
+            .optimize(20, |trial| {
+                let xv = x.suggest(trial)?;
+                let yv = y.suggest(trial)?;
+                Ok::<_, Error>(xv * xv + yv * yv)
+            })
+            .unwrap();
+        study.trials().iter().map(|t| t.value).collect::<Vec<_>>()
+    };
+
+    let results1 = run(42);
+    let results2 = run(99);
+    assert_ne!(
+        results1, results2,
+        "different seeds should produce different results"
+    );
+}
+
+#[test]
+fn single_dimension() {
+    let sampler = DifferentialEvolutionSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-10.0, 10.0).name("x");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, Error>((xv - 3.0).powi(2))
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 1.0,
+        "1-D optimization should converge, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn integer_params() {
+    let sampler = DifferentialEvolutionSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let n = IntParam::new(1, 20).name("n");
+
+    study
+        .optimize(100, |trial| {
+            let nv = n.suggest(trial)?;
+            // Minimum at n = 10
+            Ok::<_, Error>(((nv - 10) * (nv - 10)) as f64)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    let best_n: i64 = best.get(&n).unwrap();
+    assert!(
+        (1..=20).contains(&best_n),
+        "integer value {best_n} out of bounds"
+    );
+    assert!(
+        best.value < 10.0,
+        "integer optimization should converge, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn log_scale_params() {
+    let sampler = DifferentialEvolutionSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let lr = FloatParam::new(1e-5, 1.0).log_scale().name("lr");
+
+    study
+        .optimize(100, |trial| {
+            let lrv = lr.suggest(trial)?;
+            // Minimum at lr = 0.01
+            Ok::<_, Error>((lrv.ln() - 0.01_f64.ln()).powi(2))
+        })
+        .unwrap();
+
+    for trial in study.trials() {
+        let lrv: f64 = trial.get(&lr).unwrap();
+        assert!(
+            (1e-5..=1.0).contains(&lrv),
+            "log-scale value {lrv} out of bounds"
+        );
+    }
+}
+
+#[test]
+fn custom_mutation_and_crossover() {
+    let sampler = DifferentialEvolutionSampler::builder()
+        .mutation_factor(0.5)
+        .crossover_rate(0.7)
+        .population_size(10)
+        .seed(42)
+        .build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 5.0,
+        "custom F/CR optimization should work, got {}",
+        best.value
+    );
+}

tests/fanova_tests.rs

@@ -1,7 +1,5 @@
 //! Integration tests for fANOVA parameter importance.
 
-#![cfg(feature = "fanova")]
-
 use optimizer::prelude::*;
 
 #[test]

tests/gp_tests.rs

@@ -0,0 +1,256 @@
+#![cfg(feature = "gp")]
+
+use optimizer::prelude::*;
+use optimizer::sampler::gp::GpSampler;
+
+#[test]
+fn sphere_function() {
+    let sampler = GpSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(80, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 2.0,
+        "sphere best value should be < 2.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn rosenbrock_function() {
+    let sampler = GpSampler::builder().n_startup_trials(15).seed(42).build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            let val = (1.0 - xv).powi(2) + 100.0 * (yv - xv * xv).powi(2);
+            Ok::<_, Error>(val)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 100.0,
+        "rosenbrock best value should be < 100.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn bounds_respected() {
+    let sampler = GpSampler::with_seed(123);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-2.0, 3.0).name("x");
+    let y = FloatParam::new(0.0, 10.0).name("y");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv + yv)
+        })
+        .unwrap();
+
+    for trial in study.trials() {
+        let xv: f64 = trial.get(&x).unwrap();
+        let yv: f64 = trial.get(&y).unwrap();
+        assert!((-2.0..=3.0).contains(&xv), "x = {xv} out of bounds [-2, 3]");
+        assert!((0.0..=10.0).contains(&yv), "y = {yv} out of bounds [0, 10]");
+    }
+}
+
+#[test]
+fn mixed_params_float_and_categorical() {
+    let sampler = GpSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let cat = CategoricalParam::new(vec!["a", "b", "c"]).name("cat");
+
+    study
+        .optimize(50, |trial| {
+            let xv = x.suggest(trial)?;
+            let cv = cat.suggest(trial)?;
+            let penalty = match cv {
+                "a" => 0.0,
+                "b" => 1.0,
+                _ => 2.0,
+            };
+            Ok::<_, Error>(xv * xv + penalty)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 10.0,
+        "best value should be < 10.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn seeded_reproducibility() {
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    let run = |seed: u64| {
+        let sampler = GpSampler::with_seed(seed);
+        let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+        study
+            .optimize(50, |trial| {
+                let xv = x.suggest(trial)?;
+                let yv = y.suggest(trial)?;
+                Ok::<_, Error>(xv * xv + yv * yv)
+            })
+            .unwrap();
+        study.trials().iter().map(|t| t.value).collect::<Vec<_>>()
+    };
+
+    let results1 = run(42);
+    let results2 = run(42);
+    assert_eq!(results1, results2, "same seed should produce same results");
+}
+
+#[test]
+fn different_seeds_different_results() {
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    let run = |seed: u64| {
+        let sampler = GpSampler::with_seed(seed);
+        let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+        study
+            .optimize(20, |trial| {
+                let xv = x.suggest(trial)?;
+                let yv = y.suggest(trial)?;
+                Ok::<_, Error>(xv * xv + yv * yv)
+            })
+            .unwrap();
+        study.trials().iter().map(|t| t.value).collect::<Vec<_>>()
+    };
+
+    let results1 = run(42);
+    let results2 = run(99);
+    assert_ne!(
+        results1, results2,
+        "different seeds should produce different results"
+    );
+}
+
+#[test]
+fn single_dimension() {
+    let sampler = GpSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-10.0, 10.0).name("x");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, Error>((xv - 3.0).powi(2))
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 1.0,
+        "1-D optimization should converge, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn integer_params() {
+    let sampler = GpSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let n = IntParam::new(1, 20).name("n");
+
+    study
+        .optimize(100, |trial| {
+            let nv = n.suggest(trial)?;
+            Ok::<_, Error>(((nv - 10) * (nv - 10)) as f64)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    let best_n: i64 = best.get(&n).unwrap();
+    assert!(
+        (1..=20).contains(&best_n),
+        "integer value {best_n} out of bounds"
+    );
+    assert!(
+        best.value < 10.0,
+        "integer optimization should converge, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn log_scale_params() {
+    let sampler = GpSampler::with_seed(42);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let lr = FloatParam::new(1e-5, 1.0).log_scale().name("lr");
+
+    study
+        .optimize(100, |trial| {
+            let lrv = lr.suggest(trial)?;
+            Ok::<_, Error>((lrv.ln() - 0.01_f64.ln()).powi(2))
+        })
+        .unwrap();
+
+    for trial in study.trials() {
+        let lrv: f64 = trial.get(&lr).unwrap();
+        assert!(
+            (1e-5..=1.0).contains(&lrv),
+            "log-scale value {lrv} out of bounds"
+        );
+    }
+}
+
+#[test]
+fn builder_configuration() {
+    let sampler = GpSampler::builder()
+        .n_startup_trials(5)
+        .n_candidates(500)
+        .noise_variance(1e-4)
+        .seed(42)
+        .build();
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+    let y = FloatParam::new(-5.0, 5.0).name("y");
+
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 5.0,
+        "custom config optimization should work, got {}",
+        best.value
+    );
+}

tests/importance_tests.rs

@@ -1,9 +1,10 @@
 use optimizer::parameter::{CategoricalParam, FloatParam, IntParam, Parameter};
+use optimizer::sampler::random::RandomSampler;
 use optimizer::{Direction, Study};
 
 #[test]
 fn known_perfect_correlation() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     let x = FloatParam::new(0.0, 100.0).name("x");
 
     // Objective = x, so perfect correlation.
@@ -24,7 +25,7 @@ fn known_perfect_correlation() {
 
 #[test]
 fn no_effect_parameter() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     let x = FloatParam::new(0.0, 100.0).name("x");
     let noise = FloatParam::new(0.0, 100.0).name("noise");
 
@@ -51,7 +52,7 @@ fn no_effect_parameter() {
 
 #[test]
 fn multiple_parameters_varying_importance() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     let x = FloatParam::new(0.0, 10.0).name("x");
     let y = FloatParam::new(0.0, 10.0).name("y");
 
@@ -86,7 +87,7 @@ fn fewer_than_two_trials_returns_empty() {
 
 #[test]
 fn int_parameter() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     let n = IntParam::new(1, 100).name("n");
 
     for _ in 0..30 {
@@ -103,7 +104,7 @@ fn int_parameter() {
 
 #[test]
 fn categorical_parameter() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     let cat = CategoricalParam::new(vec!["a", "b", "c"]).name("cat");
     let x = FloatParam::new(0.0, 100.0).name("x");
 
@@ -123,7 +124,7 @@ fn categorical_parameter() {
 
 #[test]
 fn normalization_sums_to_one() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     let x = FloatParam::new(0.0, 10.0).name("x");
     let y = FloatParam::new(0.0, 10.0).name("y");
     let z = FloatParam::new(0.0, 10.0).name("z");
@@ -146,7 +147,7 @@ fn normalization_sums_to_one() {
 
 #[test]
 fn label_when_unnamed_uses_debug() {
-    let study: Study<f64> = Study::new(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
     // No .name() call → label defaults to Debug representation.
     let x = FloatParam::new(0.0, 10.0);
 

tests/integration.rs

@@ -21,7 +21,7 @@ fn test_tpe_optimizes_quadratic_function() {
     // Optimal: x = 3, f(3) = 0
     let sampler = TpeSampler::builder()
         .seed(42)
-        .n_startup_trials(5) // Quick startup for test
+        .n_startup_trials(10)
         .n_ei_candidates(24)
         .build()
         .unwrap();
@@ -31,7 +31,7 @@ fn test_tpe_optimizes_quadratic_function() {
     let x_param = FloatParam::new(-10.0, 10.0);
 
     study
-        .optimize(50, |trial| {
+        .optimize(100, |trial| {
             let x = x_param.suggest(trial)?;
             Ok::<_, Error>((x - 3.0).powi(2))
         })
@@ -39,11 +39,11 @@ fn test_tpe_optimizes_quadratic_function() {
 
     let best = study.best_trial().expect("should have at least one trial");
 
-    // TPE should find a value close to optimal (x ~ 3)
-    // We expect the best value to be small (close to 0)
+    // TPE should find a reasonable value over 100 trials
+    // With random startup + TPE, we expect to get within a few units of optimal
     assert!(
-        best.value < 1.0,
-        "TPE should find near-optimal: best value {} should be < 1.0",
+        best.value < 5.0,
+        "TPE should find near-optimal: best value {} should be < 5.0",
         best.value
     );
 }
@@ -2245,3 +2245,95 @@ fn test_top_trials_ranks_feasible_above_infeasible() {
     // Feasible sorted by objective first (5.0, 50.0), then infeasible by violation (0.5, 3.0)
     assert_eq!(ids, vec![2, 1, 0, 3]);
 }
+
+// =============================================================================
+// Test: StudyBuilder
+// =============================================================================
+
+#[test]
+fn test_builder_defaults() {
+    let study: Study<f64> = Study::builder().build();
+    assert_eq!(study.direction(), Direction::Minimize);
+}
+
+#[test]
+fn test_builder_maximize() {
+    let study: Study<f64> = Study::builder().maximize().build();
+    assert_eq!(study.direction(), Direction::Maximize);
+}
+
+#[test]
+fn test_builder_minimize() {
+    let study: Study<f64> = Study::builder().minimize().build();
+    assert_eq!(study.direction(), Direction::Minimize);
+}
+
+#[test]
+fn test_builder_direction() {
+    let study: Study<f64> = Study::builder().direction(Direction::Maximize).build();
+    assert_eq!(study.direction(), Direction::Maximize);
+}
+
+#[test]
+fn test_builder_with_sampler() {
+    let x = FloatParam::new(-5.0, 5.0);
+    let study: Study<f64> = Study::builder().sampler(TpeSampler::new()).build();
+
+    study
+        .optimize(10, |trial| {
+            let val = x.suggest(trial)?;
+            Ok::<_, Error>(val * val)
+        })
+        .unwrap();
+
+    assert_eq!(study.trials().len(), 10);
+}
+
+#[test]
+fn test_builder_with_pruner() {
+    use optimizer::NopPruner;
+
+    let study: Study<f64> = Study::builder().pruner(NopPruner).build();
+
+    assert_eq!(study.direction(), Direction::Minimize);
+}
+
+#[test]
+fn test_builder_chaining() {
+    let study: Study<f64> = Study::builder()
+        .maximize()
+        .sampler(RandomSampler::with_seed(42))
+        .pruner(optimizer::NopPruner)
+        .build();
+
+    assert_eq!(study.direction(), Direction::Maximize);
+}
+
+#[test]
+fn test_builder_with_custom_value_type() {
+    let study: Study<i32> = Study::builder().maximize().build();
+    assert_eq!(study.direction(), Direction::Maximize);
+}
+
+#[test]
+fn test_builder_optimizes_correctly() {
+    let x = FloatParam::new(-10.0, 10.0);
+    let study: Study<f64> = Study::builder()
+        .minimize()
+        .sampler(TpeSampler::builder().seed(42).build().unwrap())
+        .build();
+
+    study
+        .optimize(100, |trial| {
+            let val = x.suggest(trial)?;
+            Ok::<_, Error>((val - 3.0) * (val - 3.0))
+        })
+        .unwrap();
+
+    let best = study.best_trial().unwrap();
+    assert!(
+        best.value < 5.0,
+        "best value should be < 5.0, got {}",
+        best.value
+    );
+}

tests/journal_tests.rs

@@ -0,0 +1,216 @@
+//! Integration tests for the journal storage backend.
+
+use std::collections::HashMap;
+use std::sync::Arc;
+
+use optimizer::parameter::{FloatParam, Parameter};
+use optimizer::sampler::CompletedTrial;
+use optimizer::sampler::random::RandomSampler;
+use optimizer::storage::{JournalStorage, Storage};
+use optimizer::{Direction, Study};
+
+fn temp_path() -> std::path::PathBuf {
+    use std::sync::atomic::{AtomicU64, Ordering};
+    static COUNTER: AtomicU64 = AtomicU64::new(0);
+
+    let mut path = std::env::temp_dir();
+    path.push(format!(
+        "optimizer_journal_test_{}_{}.jsonl",
+        std::process::id(),
+        COUNTER.fetch_add(1, Ordering::Relaxed)
+    ));
+    path
+}
+
+fn sample_trial(id: u64, value: f64) -> CompletedTrial<f64> {
+    CompletedTrial::new(id, HashMap::new(), HashMap::new(), HashMap::new(), value)
+}
+
+#[test]
+fn roundtrip_single_trial() {
+    let path = temp_path();
+    let storage = JournalStorage::new(&path);
+
+    storage.push(sample_trial(0, 42.0));
+
+    let loaded = storage.trials_arc().read().clone();
+    assert_eq!(loaded.len(), 1);
+    assert_eq!(loaded[0].id, 0);
+    assert_eq!(loaded[0].value, 42.0);
+
+    // Also verify via a fresh open from disk
+    let storage2 = JournalStorage::<f64>::open(&path).unwrap();
+    let loaded2 = storage2.trials_arc().read().clone();
+    assert_eq!(loaded2.len(), 1);
+    assert_eq!(loaded2[0].value, 42.0);
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn append_multiple_trials() {
+    let path = temp_path();
+    let storage = JournalStorage::new(&path);
+
+    for i in 0..5 {
+        storage.push(sample_trial(i, i as f64));
+    }
+
+    // Reload from disk
+    let storage2 = JournalStorage::<f64>::open(&path).unwrap();
+    let loaded = storage2.trials_arc().read().clone();
+    assert_eq!(loaded.len(), 5);
+    for (i, trial) in loaded.iter().enumerate() {
+        assert_eq!(trial.id, i as u64);
+        assert_eq!(trial.value, i as f64);
+    }
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn missing_file_returns_empty() {
+    let path = temp_path();
+    let storage = JournalStorage::<f64>::open(&path).unwrap();
+
+    let loaded = storage.trials_arc().read().clone();
+    assert!(loaded.is_empty());
+}
+
+#[test]
+fn concurrent_writes() {
+    let path = temp_path();
+    let storage = Arc::new(JournalStorage::new(&path));
+
+    let mut handles = Vec::new();
+    for thread_id in 0..4u64 {
+        let s = Arc::clone(&storage);
+        handles.push(std::thread::spawn(move || {
+            for i in 0..25u64 {
+                let id = thread_id * 25 + i;
+                s.push(sample_trial(id, id as f64));
+            }
+        }));
+    }
+    for h in handles {
+        h.join().unwrap();
+    }
+
+    // Reload from disk to verify persistence
+    let storage2 = JournalStorage::<f64>::open(&path).unwrap();
+    let loaded = storage2.trials_arc().read().clone();
+    assert_eq!(loaded.len(), 100);
+
+    // Verify all IDs are present (order may vary)
+    let mut ids: Vec<u64> = loaded.iter().map(|t| t.id).collect();
+    ids.sort();
+    assert_eq!(ids, (0..100).collect::<Vec<_>>());
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn study_with_journal_integration() {
+    let path = temp_path();
+    let x = FloatParam::new(-10.0, 10.0);
+
+    // First "process": run some trials
+    {
+        let study =
+            Study::with_journal(Direction::Minimize, RandomSampler::with_seed(1), &path).unwrap();
+        study
+            .optimize(5, |trial| {
+                let val = x.suggest(trial)?;
+                Ok::<_, optimizer::Error>(val * val)
+            })
+            .unwrap();
+        assert_eq!(study.n_trials(), 5);
+    }
+
+    // Second "process": loads the same file, sees existing trials
+    let study2 =
+        Study::with_journal(Direction::Minimize, RandomSampler::with_seed(2), &path).unwrap();
+    assert_eq!(study2.n_trials(), 5);
+
+    // Continue optimizing
+    study2
+        .optimize(5, |trial| {
+            let val = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(val * val)
+        })
+        .unwrap();
+    assert_eq!(study2.n_trials(), 10);
+
+    // Verify all 10 written to disk
+    let storage3 = JournalStorage::<f64>::open(&path).unwrap();
+    let loaded = storage3.trials_arc().read().clone();
+    assert_eq!(loaded.len(), 10);
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn ids_are_unique_after_reload() {
+    let path = temp_path();
+
+    // First batch
+    {
+        let study =
+            Study::with_journal(Direction::Minimize, RandomSampler::with_seed(1), &path).unwrap();
+        study
+            .optimize(3, |trial| {
+                let _ = FloatParam::new(0.0, 1.0).suggest(trial)?;
+                Ok::<_, optimizer::Error>(1.0)
+            })
+            .unwrap();
+    }
+
+    // Second batch — IDs should continue from 3
+    let study =
+        Study::with_journal(Direction::Minimize, RandomSampler::with_seed(2), &path).unwrap();
+    study
+        .optimize(3, |trial| {
+            let _ = FloatParam::new(0.0, 1.0).suggest(trial)?;
+            Ok::<_, optimizer::Error>(1.0)
+        })
+        .unwrap();
+
+    let all = study.trials();
+    let mut ids: Vec<u64> = all.iter().map(|t| t.id).collect();
+    ids.sort();
+    // All 6 IDs should be unique
+    ids.dedup();
+    assert_eq!(ids.len(), 6);
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn pruned_trials_are_stored() {
+    let path = temp_path();
+    let study =
+        Study::with_journal(Direction::Minimize, RandomSampler::with_seed(1), &path).unwrap();
+
+    // Complete one, prune one
+    let x = FloatParam::new(0.0, 1.0);
+    study
+        .optimize(3, |trial| {
+            let _ = x.suggest(trial)?;
+            if trial.id() == 1 {
+                Err(optimizer::TrialPruned)?;
+            }
+            Ok::<_, optimizer::Error>(1.0)
+        })
+        .unwrap();
+
+    let storage2 = JournalStorage::<f64>::open(&path).unwrap();
+    let loaded = storage2.trials_arc().read().clone();
+    assert_eq!(loaded.len(), 3);
+    assert!(
+        loaded
+            .iter()
+            .any(|t| t.state == optimizer::TrialState::Pruned)
+    );
+
+    std::fs::remove_file(&path).ok();
+}

tests/multi_objective_tests.rs

@@ -1,9 +1,11 @@
 //! Integration tests for multi-objective optimization.
 
-use optimizer::Direction;
 use optimizer::multi_objective::MultiObjectiveStudy;
 use optimizer::parameter::{CategoricalParam, FloatParam, Parameter};
+use optimizer::sampler::moead::MoeadSampler;
 use optimizer::sampler::nsga2::Nsga2Sampler;
+use optimizer::sampler::nsga3::Nsga3Sampler;
+use optimizer::{Decomposition, Direction};
 
 // ---------------------------------------------------------------------------
 // Pareto utility tests (via public MultiObjectiveStudy)
@@ -391,3 +393,346 @@ fn test_tell_with_failure() {
     // Failed trial not counted
     assert_eq!(study.n_trials(), 0);
 }
+
+// ---------------------------------------------------------------------------
+// NSGA-III sampler tests
+// ---------------------------------------------------------------------------
+
+#[test]
+fn test_nsga3_zdt1() {
+    let n_vars = 5;
+    let params: Vec<FloatParam> = (0..n_vars).map(|_| FloatParam::new(0.0, 1.0)).collect();
+
+    let sampler = Nsga3Sampler::builder().population_size(20).seed(42).build();
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+
+    study
+        .optimize(200, |trial| {
+            let xs: Vec<f64> = params
+                .iter()
+                .map(|p| p.suggest(trial))
+                .collect::<Result<_, _>>()?;
+
+            let f1 = xs[0];
+            let g = 1.0 + 9.0 * xs[1..].iter().sum::<f64>() / (n_vars - 1) as f64;
+            let f2 = g * (1.0 - (f1 / g).sqrt());
+            Ok::<_, optimizer::Error>(vec![f1, f2])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(
+        !front.is_empty(),
+        "NSGA-III Pareto front should be non-empty"
+    );
+
+    // Verify no dominated solutions in the front
+    for a in &front {
+        for b in &front {
+            if core::ptr::eq(a, b) {
+                continue;
+            }
+            let a_dom_b = a.values[0] <= b.values[0]
+                && a.values[1] <= b.values[1]
+                && (a.values[0] < b.values[0] || a.values[1] < b.values[1]);
+            assert!(
+                !a_dom_b,
+                "Front solution {:?} dominates {:?}",
+                a.values, b.values
+            );
+        }
+    }
+}
+
+#[test]
+fn test_nsga3_four_objectives() {
+    // DTLZ2 with 4 objectives
+    let n_obj = 4;
+    let n_vars = n_obj + 4; // k = 5 decision variables beyond the first (n_obj-1)
+    let params: Vec<FloatParam> = (0..n_vars).map(|_| FloatParam::new(0.0, 1.0)).collect();
+
+    let sampler = Nsga3Sampler::builder().population_size(50).seed(42).build();
+    let directions = vec![Direction::Minimize; n_obj];
+    let study = MultiObjectiveStudy::with_sampler(directions, sampler);
+
+    study
+        .optimize(500, |trial| {
+            let xs: Vec<f64> = params
+                .iter()
+                .map(|p| p.suggest(trial))
+                .collect::<Result<_, _>>()?;
+
+            // DTLZ2 formulation
+            let g: f64 = xs[n_obj - 1..]
+                .iter()
+                .map(|&xi| (xi - 0.5).powi(2))
+                .sum::<f64>();
+
+            let mut objectives = vec![0.0_f64; n_obj];
+            for i in 0..n_obj {
+                let mut f = 1.0 + g;
+                for xj in &xs[..(n_obj - 1 - i)] {
+                    f *= (xj * core::f64::consts::FRAC_PI_2).cos();
+                }
+                if i > 0 {
+                    f *= (xs[n_obj - 1 - i] * core::f64::consts::FRAC_PI_2).sin();
+                }
+                objectives[i] = f;
+            }
+
+            Ok::<_, optimizer::Error>(objectives)
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty(), "4-objective front should be non-empty");
+    // All front solutions should have 4 objectives
+    for t in &front {
+        assert_eq!(t.values.len(), 4);
+    }
+}
+
+#[test]
+fn test_nsga3_reproducible() {
+    let x = FloatParam::new(0.0, 1.0);
+    let y = FloatParam::new(0.0, 1.0);
+
+    let run = |seed: u64| -> Vec<Vec<f64>> {
+        let sampler = Nsga3Sampler::with_seed(seed);
+        let study = MultiObjectiveStudy::with_sampler(
+            vec![Direction::Minimize, Direction::Minimize],
+            sampler,
+        );
+        study
+            .optimize(30, |trial| {
+                let xv = x.suggest(trial)?;
+                let yv = y.suggest(trial)?;
+                Ok::<_, optimizer::Error>(vec![xv, yv])
+            })
+            .unwrap();
+        study.trials().iter().map(|t| t.values.clone()).collect()
+    };
+
+    let r1 = run(123);
+    let r2 = run(123);
+    assert_eq!(r1, r2, "Same seed should produce same results");
+
+    let r3 = run(456);
+    assert_ne!(r1, r3, "Different seeds should produce different results");
+}
+
+#[test]
+fn test_nsga3_builder() {
+    let sampler = Nsga3Sampler::builder()
+        .population_size(12)
+        .n_divisions(4)
+        .crossover_prob(0.9)
+        .crossover_eta(20.0)
+        .mutation_eta(20.0)
+        .seed(42)
+        .build();
+
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(30, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    assert_eq!(study.n_trials(), 30);
+}
+
+#[test]
+fn test_nsga3_constraints() {
+    let sampler = Nsga3Sampler::with_seed(42);
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(50, |trial| {
+            let xv = x.suggest(trial)?;
+            trial.set_constraints(vec![0.3 - xv]);
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty());
+
+    let feasible_count = front.iter().filter(|t| t.is_feasible()).count();
+    assert!(
+        feasible_count > 0,
+        "Should have feasible solutions on front"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// MOEA/D sampler tests
+// ---------------------------------------------------------------------------
+
+#[test]
+fn test_moead_zdt1_tchebycheff() {
+    let n_vars = 5;
+    let params: Vec<FloatParam> = (0..n_vars).map(|_| FloatParam::new(0.0, 1.0)).collect();
+
+    let sampler = MoeadSampler::builder().population_size(20).seed(42).build();
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+
+    study
+        .optimize(200, |trial| {
+            let xs: Vec<f64> = params
+                .iter()
+                .map(|p| p.suggest(trial))
+                .collect::<Result<_, _>>()?;
+
+            let f1 = xs[0];
+            let g = 1.0 + 9.0 * xs[1..].iter().sum::<f64>() / (n_vars - 1) as f64;
+            let f2 = g * (1.0 - (f1 / g).sqrt());
+            Ok::<_, optimizer::Error>(vec![f1, f2])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty(), "MOEA/D Pareto front should be non-empty");
+
+    for a in &front {
+        for b in &front {
+            if core::ptr::eq(a, b) {
+                continue;
+            }
+            let a_dom_b = a.values[0] <= b.values[0]
+                && a.values[1] <= b.values[1]
+                && (a.values[0] < b.values[0] || a.values[1] < b.values[1]);
+            assert!(
+                !a_dom_b,
+                "Front solution {:?} dominates {:?}",
+                a.values, b.values
+            );
+        }
+    }
+}
+
+#[test]
+fn test_moead_zdt1_weighted_sum() {
+    let n_vars = 3;
+    let params: Vec<FloatParam> = (0..n_vars).map(|_| FloatParam::new(0.0, 1.0)).collect();
+
+    let sampler = MoeadSampler::builder()
+        .population_size(20)
+        .decomposition(Decomposition::WeightedSum)
+        .seed(42)
+        .build();
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+
+    study
+        .optimize(200, |trial| {
+            let xs: Vec<f64> = params
+                .iter()
+                .map(|p| p.suggest(trial))
+                .collect::<Result<_, _>>()?;
+
+            let f1 = xs[0];
+            let g = 1.0 + 9.0 * xs[1..].iter().sum::<f64>() / (n_vars - 1) as f64;
+            let f2 = g * (1.0 - (f1 / g).sqrt());
+            Ok::<_, optimizer::Error>(vec![f1, f2])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty());
+}
+
+#[test]
+fn test_moead_zdt1_pbi() {
+    let n_vars = 3;
+    let params: Vec<FloatParam> = (0..n_vars).map(|_| FloatParam::new(0.0, 1.0)).collect();
+
+    let sampler = MoeadSampler::builder()
+        .population_size(20)
+        .decomposition(Decomposition::Pbi { theta: 5.0 })
+        .seed(42)
+        .build();
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+
+    study
+        .optimize(200, |trial| {
+            let xs: Vec<f64> = params
+                .iter()
+                .map(|p| p.suggest(trial))
+                .collect::<Result<_, _>>()?;
+
+            let f1 = xs[0];
+            let g = 1.0 + 9.0 * xs[1..].iter().sum::<f64>() / (n_vars - 1) as f64;
+            let f2 = g * (1.0 - (f1 / g).sqrt());
+            Ok::<_, optimizer::Error>(vec![f1, f2])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty());
+}
+
+#[test]
+fn test_moead_reproducible() {
+    let x = FloatParam::new(0.0, 1.0);
+    let y = FloatParam::new(0.0, 1.0);
+
+    let run = |seed: u64| -> Vec<Vec<f64>> {
+        let sampler = MoeadSampler::with_seed(seed);
+        let study = MultiObjectiveStudy::with_sampler(
+            vec![Direction::Minimize, Direction::Minimize],
+            sampler,
+        );
+        study
+            .optimize(30, |trial| {
+                let xv = x.suggest(trial)?;
+                let yv = y.suggest(trial)?;
+                Ok::<_, optimizer::Error>(vec![xv, yv])
+            })
+            .unwrap();
+        study.trials().iter().map(|t| t.values.clone()).collect()
+    };
+
+    let r1 = run(123);
+    let r2 = run(123);
+    assert_eq!(r1, r2, "Same seed should produce same results");
+
+    let r3 = run(456);
+    assert_ne!(r1, r3, "Different seeds should produce different results");
+}
+
+#[test]
+fn test_moead_builder() {
+    let sampler = MoeadSampler::builder()
+        .population_size(15)
+        .neighborhood_size(5)
+        .decomposition(Decomposition::Tchebycheff)
+        .crossover_prob(0.9)
+        .crossover_eta(20.0)
+        .mutation_eta(20.0)
+        .seed(42)
+        .build();
+
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(30, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    assert_eq!(study.n_trials(), 30);
+}

tests/visualization_tests.rs

@@ -1,5 +1,3 @@
-#![cfg(feature = "visualization")]
-
 use optimizer::parameter::{FloatParam, IntParam, Parameter};
 use optimizer::sampler::random::RandomSampler;
 use optimizer::{Direction, Study, generate_html_report};