chore: release v0.8.1

Implement functional ANOVA decomposition behind the `fanova` feature
flag. A self-contained random forest is trained on trial data, then
marginal predictions are used to compute per-parameter main effects
and pairwise interaction effects, normalized to sum to 1.0.

Adds Study::fanova() / fanova_with_config(), FanovaResult, and
FanovaConfig. Includes unit and integration tests covering dominant
parameters, interaction detection, consistency with correlation-based
importance, and error handling.

.github/workflows/ci.yml

@@ -8,6 +8,7 @@ on:
 
 permissions:
   contents: read
+  pull-requests: write
 
 env:
   CARGO_TERM_COLOR: always
@@ -73,8 +74,14 @@ jobs:
       - uses: Swatinem/rust-cache@v2
       - name: Run sync example (ml_hyperparameter_tuning)
         run: cargo run --example ml_hyperparameter_tuning
+      - name: Run sync example (benchmark_convergence)
+        run: cargo run --example benchmark_convergence
+      - name: Run derive example (parameter_api)
+        run: cargo run --example parameter_api --features derive
       - 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
 
   docs:
     name: Docs
@@ -250,6 +257,99 @@ jobs:
       - name: Typos Check
         run: typos src/
 
+  bench-check:
+    name: Benchmark Smoke Test
+    if: github.event_name == 'push'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v6
+      - name: Install Rust
+        run: |
+          rustup override set stable
+          rustup update stable
+      - uses: Swatinem/rust-cache@v2
+      - name: Run benchmarks (smoke test)
+        run: cargo bench --all-features --bench samplers --bench optimization
+
+  bench-compare:
+    name: Benchmark Comparison
+    if: github.event_name == 'pull_request'
+    runs-on: ubuntu-latest
+    steps:
+      - name: Install Rust
+        run: |
+          rustup override set stable
+          rustup update stable
+
+      - name: Install critcmp
+        run: cargo install critcmp
+
+      # --- Run benchmarks on the BASE (target) branch ---
+      - uses: actions/checkout@v6
+        with:
+          ref: ${{ github.event.pull_request.base.sha }}
+          clean: false
+
+      - uses: Swatinem/rust-cache@v2
+        with:
+          key: bench-base
+
+      - name: Bench baseline
+        run: cargo bench --all-features --bench samplers --bench optimization -- --save-baseline base
+
+      # --- Run benchmarks on the HEAD (PR) branch ---
+      - uses: actions/checkout@v6
+        with:
+          ref: ${{ github.event.pull_request.head.sha }}
+          clean: false
+
+      - uses: Swatinem/rust-cache@v2
+        with:
+          key: bench-head
+
+      - name: Bench PR head
+        run: cargo bench --all-features --bench samplers --bench optimization -- --save-baseline head
+
+      # --- Compare and post comment ---
+      - name: Compare benchmarks
+        id: compare
+        run: |
+          EOF=$(dd if=/dev/urandom bs=15 count=1 status=none | base64)
+          echo "result<<$EOF" >> "$GITHUB_OUTPUT"
+          critcmp base head --color never >> "$GITHUB_OUTPUT"
+          echo "$EOF" >> "$GITHUB_OUTPUT"
+
+      - name: Find existing comment
+        uses: peter-evans/find-comment@v4
+        id: find-comment
+        with:
+          issue-number: ${{ github.event.pull_request.number }}
+          comment-author: github-actions[bot]
+          body-includes: "## Benchmark Comparison"
+
+      - name: Post or update PR comment
+        uses: peter-evans/create-or-update-comment@v5
+        with:
+          comment-id: ${{ steps.find-comment.outputs.comment-id }}
+          issue-number: ${{ github.event.pull_request.number }}
+          edit-mode: replace
+          body: |
+            ## Benchmark Comparison
+
+            **Base:** `${{ github.event.pull_request.base.sha }}` | **Head:** `${{ github.event.pull_request.head.sha }}`
+
+            <details>
+            <summary>Click to expand full results</summary>
+
+            ```
+            ${{ steps.compare.outputs.result }}
+            ```
+
+            </details>
+
+            > Benchmarks run on `ubuntu-latest` via [Criterion](https://github.com/bheisler/criterion.rs) + [critcmp](https://github.com/BurntSushi/critcmp).
+            > Results may vary due to shared CI runners. Look for consistent >5% changes.
+
   publish:
     name: Publish to crates.io
     runs-on: ubuntu-latest

Cargo.toml

@@ -3,7 +3,7 @@ members = ["optimizer-derive"]
 
 [package]
 name = "optimizer"
-version = "0.6.0"
+version = "0.8.1"
 edition = "2024"
 rust-version = "1.88"
 license = "MIT"
@@ -21,15 +21,36 @@ thiserror = "2"
 parking_lot = "0.12"
 tokio = { version = "1", features = ["sync", "rt-multi-thread"], optional = true }
 optimizer-derive = { version = "0.1.0", path = "optimizer-derive", optional = true }
+serde = { version = "1", features = ["derive"], optional = true }
+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 }
 
 [features]
 default = []
 async = ["dep:tokio"]
 derive = ["dep:optimizer-derive"]
+serde = ["dep:serde", "dep:serde_json"]
+tracing = ["dep:tracing"]
+sobol = ["dep:sobol_burley"]
+cma-es = ["dep:nalgebra"]
+visualization = []
+fanova = []
 
 [dev-dependencies]
 tokio = { version = "1", features = ["rt-multi-thread", "macros", "time"] }
 optimizer-derive = { version = "0.1.0", path = "optimizer-derive" }
+serde_json = "1"
+criterion = { version = "0.8", features = ["html_reports"] }
+
+[[bench]]
+name = "samplers"
+harness = false
+
+[[bench]]
+name = "optimization"
+harness = false
 
 [[example]]
 name = "async_api_optimization"
@@ -44,3 +65,8 @@ path = "examples/ml_hyperparameter_tuning.rs"
 name = "parameter_api"
 path = "examples/parameter_api.rs"
 required-features = ["derive"]
+
+[[example]]
+name = "visualization_demo"
+path = "examples/visualization_demo.rs"
+required-features = ["visualization"]

benches/optimization.rs

@@ -0,0 +1,112 @@
+#[allow(dead_code)]
+mod test_functions;
+
+use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main};
+use optimizer::parameter::Parameter;
+use optimizer::sampler::random::RandomSampler;
+use optimizer::sampler::tpe::TpeSampler;
+use optimizer::{FloatParam, Study};
+
+fn make_params(dims: usize) -> Vec<FloatParam> {
+    (0..dims)
+        .map(|i| FloatParam::new(-5.0, 5.0).name(format!("x{i}")))
+        .collect()
+}
+
+fn bench_tpe_sphere(c: &mut Criterion) {
+    let mut group = c.benchmark_group("tpe_sphere");
+    group.sample_size(10);
+
+    for dims in [2, 10, 50] {
+        let params = make_params(dims);
+        group.bench_with_input(BenchmarkId::new("dims", dims), &params, |b, params| {
+            b.iter(|| {
+                let study = Study::minimize(TpeSampler::builder().seed(42).build().unwrap());
+                study
+                    .optimize(100, |trial| {
+                        let x: Vec<f64> = params
+                            .iter()
+                            .map(|p| p.suggest(trial))
+                            .collect::<Result<_, _>>()
+                            .unwrap();
+                        Ok::<_, optimizer::Error>(test_functions::sphere(&x))
+                    })
+                    .unwrap();
+            });
+        });
+    }
+    group.finish();
+}
+
+fn bench_tpe_rosenbrock(c: &mut Criterion) {
+    let mut group = c.benchmark_group("tpe_rosenbrock");
+    group.sample_size(10);
+
+    for dims in [2, 10] {
+        let params = make_params(dims);
+        group.bench_with_input(BenchmarkId::new("dims", dims), &params, |b, params| {
+            b.iter(|| {
+                let study = Study::minimize(TpeSampler::builder().seed(42).build().unwrap());
+                study
+                    .optimize(100, |trial| {
+                        let x: Vec<f64> = params
+                            .iter()
+                            .map(|p| p.suggest(trial))
+                            .collect::<Result<_, _>>()
+                            .unwrap();
+                        Ok::<_, optimizer::Error>(test_functions::rosenbrock(&x))
+                    })
+                    .unwrap();
+            });
+        });
+    }
+    group.finish();
+}
+
+fn bench_random_vs_tpe(c: &mut Criterion) {
+    let mut group = c.benchmark_group("random_vs_tpe");
+    group.sample_size(10);
+    let params = make_params(5);
+
+    group.bench_function("random_5d", |b| {
+        b.iter(|| {
+            let study = Study::minimize(RandomSampler::with_seed(42));
+            study
+                .optimize(100, |trial| {
+                    let x: Vec<f64> = params
+                        .iter()
+                        .map(|p| p.suggest(trial))
+                        .collect::<Result<_, _>>()
+                        .unwrap();
+                    Ok::<_, optimizer::Error>(test_functions::sphere(&x))
+                })
+                .unwrap();
+        });
+    });
+
+    group.bench_function("tpe_5d", |b| {
+        b.iter(|| {
+            let study = Study::minimize(TpeSampler::builder().seed(42).build().unwrap());
+            study
+                .optimize(100, |trial| {
+                    let x: Vec<f64> = params
+                        .iter()
+                        .map(|p| p.suggest(trial))
+                        .collect::<Result<_, _>>()
+                        .unwrap();
+                    Ok::<_, optimizer::Error>(test_functions::sphere(&x))
+                })
+                .unwrap();
+        });
+    });
+
+    group.finish();
+}
+
+criterion_group!(
+    benches,
+    bench_tpe_sphere,
+    bench_tpe_rosenbrock,
+    bench_random_vs_tpe
+);
+criterion_main!(benches);

benches/samplers.rs

@@ -0,0 +1,118 @@
+use std::collections::HashMap;
+
+use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main};
+use optimizer::parameter::Parameter;
+use optimizer::sampler::Sampler;
+use optimizer::sampler::grid::GridSearchSampler;
+use optimizer::sampler::random::RandomSampler;
+use optimizer::sampler::tpe::TpeSampler;
+use optimizer::{CompletedTrial, FloatParam};
+
+/// Build a synthetic history of `n` completed trials over `dims` float parameters.
+fn build_history(n: usize, dims: usize) -> Vec<CompletedTrial<f64>> {
+    let params: Vec<FloatParam> = (0..dims)
+        .map(|i| FloatParam::new(-5.0, 5.0).name(format!("x{i}")))
+        .collect();
+
+    let mut history = Vec::with_capacity(n);
+    let sampler = RandomSampler::with_seed(42);
+
+    for trial_id in 0..n {
+        let id = trial_id as u64;
+        let mut param_values = HashMap::new();
+        let mut distributions = HashMap::new();
+        let mut param_labels = HashMap::new();
+        for p in &params {
+            let dist = p.distribution();
+            let val = sampler.sample(&dist, id, &history);
+            param_values.insert(p.id(), val);
+            distributions.insert(p.id(), dist);
+            param_labels.insert(p.id(), p.label());
+        }
+        // Use sphere function value as objective
+        let value: f64 = param_values
+            .values()
+            .map(|v| {
+                let optimizer::ParamValue::Float(f) = v else {
+                    unreachable!()
+                };
+                f * f
+            })
+            .sum();
+        history.push(CompletedTrial::new(
+            id,
+            param_values,
+            distributions,
+            param_labels,
+            value,
+        ));
+    }
+    history
+}
+
+fn bench_tpe_sample(c: &mut Criterion) {
+    let mut group = c.benchmark_group("tpe_sample");
+    let dist = FloatParam::new(-5.0, 5.0).distribution();
+    let tpe = TpeSampler::builder().seed(42).build().unwrap();
+
+    for history_size in [10, 100, 1000] {
+        let history = build_history(history_size, 2);
+        group.bench_with_input(
+            BenchmarkId::new("history", history_size),
+            &history,
+            |b, history| {
+                b.iter(|| tpe.sample(&dist, history.len() as u64, history));
+            },
+        );
+    }
+    group.finish();
+}
+
+fn bench_random_sample(c: &mut Criterion) {
+    let mut group = c.benchmark_group("random_sample");
+    let dist = FloatParam::new(-5.0, 5.0).distribution();
+    let sampler = RandomSampler::with_seed(42);
+
+    for history_size in [10, 100, 1000] {
+        let history = build_history(history_size, 2);
+        group.bench_with_input(
+            BenchmarkId::new("history", history_size),
+            &history,
+            |b, history| {
+                b.iter(|| sampler.sample(&dist, history.len() as u64, history));
+            },
+        );
+    }
+    group.finish();
+}
+
+fn bench_grid_sample(c: &mut Criterion) {
+    let mut group = c.benchmark_group("grid_sample");
+    let dist = FloatParam::new(-5.0, 5.0).distribution();
+    let history: Vec<CompletedTrial<f64>> = Vec::new();
+
+    for grid_points in [5, 10, 50] {
+        group.bench_with_input(
+            BenchmarkId::new("points", grid_points),
+            &grid_points,
+            |b, _| {
+                b.iter(|| {
+                    // Fresh sampler each iteration since grid tracks used points
+                    let sampler = GridSearchSampler::builder()
+                        .n_points_per_param(grid_points)
+                        .build();
+                    sampler.sample(&dist, 0, &history)
+                });
+            },
+        );
+    }
+    group.finish();
+}
+
+criterion_group!(
+    benches,
+    bench_tpe_sample,
+    bench_random_sample,
+    bench_grid_sample
+);
+criterion_main!(benches);

benches/test_functions.rs

@@ -0,0 +1,84 @@
+//! Standard optimization test functions for benchmarking.
+
+/// Sphere function: unimodal, convex. Global minimum f(0,...,0) = 0.
+pub fn sphere(x: &[f64]) -> f64 {
+    x.iter().map(|xi| xi * xi).sum()
+}
+
+/// Rosenbrock function: narrow valley. Global minimum f(1,...,1) = 0.
+pub fn rosenbrock(x: &[f64]) -> f64 {
+    x.windows(2)
+        .map(|w| 100.0 * (w[1] - w[0] * w[0]).powi(2) + (1.0 - w[0]).powi(2))
+        .sum()
+}
+
+/// Rastrigin function: highly multimodal. Global minimum f(0,...,0) = 0.
+pub fn rastrigin(x: &[f64]) -> f64 {
+    let n = x.len() as f64;
+    10.0 * n
+        + x.iter()
+            .map(|xi| xi * xi - 10.0 * (2.0 * std::f64::consts::PI * xi).cos())
+            .sum::<f64>()
+}
+
+/// Ackley function: nearly flat with a deep well. Global minimum f(0,...,0) = 0.
+pub fn ackley(x: &[f64]) -> f64 {
+    let n = x.len() as f64;
+    let sum_sq: f64 = x.iter().map(|xi| xi * xi).sum();
+    let sum_cos: f64 = x
+        .iter()
+        .map(|xi| (2.0 * std::f64::consts::PI * xi).cos())
+        .sum();
+    -20.0 * (-0.2 * (sum_sq / n).sqrt()).exp() - (sum_cos / n).exp() + 20.0 + std::f64::consts::E
+}
+
+/// Branin function (2D only). Three global minima with f* ≈ 0.397887.
+///
+/// # Panics
+///
+/// Panics if `x` does not have exactly 2 elements.
+pub fn branin(x: &[f64]) -> f64 {
+    assert!(x.len() == 2, "Branin requires exactly 2 dimensions");
+    let (x1, x2) = (x[0], x[1]);
+    let pi = std::f64::consts::PI;
+    let a = 1.0;
+    let b = 5.1 / (4.0 * pi * pi);
+    let c = 5.0 / pi;
+    let r = 6.0;
+    let s = 10.0;
+    let t = 1.0 / (8.0 * pi);
+    a * (x2 - b * x1 * x1 + c * x1 - r).powi(2) + s * (1.0 - t) * x1.cos() + s
+}
+
+/// Hartmann 6D function. Global minimum f* ≈ -3.3224.
+///
+/// # Panics
+///
+/// Panics if `x` does not have exactly 6 elements.
+pub fn hartmann6(x: &[f64]) -> f64 {
+    assert!(x.len() == 6, "Hartmann6 requires exactly 6 dimensions");
+
+    let alpha = [1.0, 1.2, 3.0, 3.2];
+    let a_matrix = [
+        [10.0, 3.0, 17.0, 3.5, 1.7, 8.0],
+        [0.05, 10.0, 17.0, 0.1, 8.0, 14.0],
+        [3.0, 3.5, 1.7, 10.0, 17.0, 8.0],
+        [17.0, 8.0, 0.05, 10.0, 0.1, 14.0],
+    ];
+    let p_matrix = [
+        [0.1312, 0.1696, 0.5569, 0.0124, 0.8283, 0.5886],
+        [0.2329, 0.4135, 0.8307, 0.3736, 0.1004, 0.9991],
+        [0.2348, 0.1451, 0.3522, 0.2883, 0.3047, 0.6650],
+        [0.4047, 0.8828, 0.8732, 0.5743, 0.1091, 0.0381],
+    ];
+
+    let mut result = 0.0;
+    for i in 0..4 {
+        let mut inner = 0.0;
+        for (j, xj) in x.iter().enumerate() {
+            inner += a_matrix[i][j] * (xj - p_matrix[i][j]).powi(2);
+        }
+        result -= alpha[i] * (-inner).exp();
+    }
+    result
+}

deny.toml

@@ -1,7 +1,11 @@
 [advisories]
 version = 2
 db-path = "~/.cargo/advisory-db"
-ignore = []
+ignore = [
+    # paste is unmaintained but it's a transitive dep via simba -> nalgebra
+    # with no fix available upstream yet
+    "RUSTSEC-2024-0436",
+]
 
 [licenses]
 version = 2
@@ -9,7 +13,6 @@ allow = [
     "MIT",
     "Apache-2.0",
     "Apache-2.0 WITH LLVM-exception",
-    "BSD-2-Clause",
     "Unicode-3.0",
 ]
 confidence-threshold = 0.8

examples/benchmark_convergence.rs

@@ -0,0 +1,124 @@
+use std::ops::ControlFlow;
+use std::time::Instant;
+
+use optimizer::parameter::Parameter;
+use optimizer::sampler::random::RandomSampler;
+use optimizer::sampler::tpe::TpeSampler;
+use optimizer::{FloatParam, Study};
+
+/// Standard optimization test functions.
+mod functions {
+    pub fn sphere(x: &[f64]) -> f64 {
+        x.iter().map(|xi| xi * xi).sum()
+    }
+
+    pub fn rosenbrock(x: &[f64]) -> f64 {
+        x.windows(2)
+            .map(|w| 100.0 * (w[1] - w[0] * w[0]).powi(2) + (1.0 - w[0]).powi(2))
+            .sum()
+    }
+
+    pub fn rastrigin(x: &[f64]) -> f64 {
+        let n = x.len() as f64;
+        10.0 * n
+            + x.iter()
+                .map(|xi| xi * xi - 10.0 * (2.0 * std::f64::consts::PI * xi).cos())
+                .sum::<f64>()
+    }
+}
+
+fn run_convergence(
+    name: &str,
+    sampler_name: &str,
+    study: Study<f64>,
+    params: &[FloatParam],
+    objective: fn(&[f64]) -> f64,
+    n_trials: usize,
+) {
+    let start = Instant::now();
+
+    study
+        .optimize_with_callback(
+            n_trials,
+            |trial| {
+                let x: Vec<f64> = params
+                    .iter()
+                    .map(|p| p.suggest(trial))
+                    .collect::<Result<_, _>>()
+                    .unwrap();
+                Ok::<_, optimizer::Error>(objective(&x))
+            },
+            |study, _trial| {
+                let elapsed = start.elapsed().as_millis();
+                let best = study.best_value().unwrap();
+                let n = study.n_trials();
+                println!("{n},{best},{elapsed},{sampler_name},{name}");
+                ControlFlow::Continue(())
+            },
+        )
+        .unwrap();
+}
+
+fn main() {
+    println!("trial,best_value,elapsed_ms,sampler,function");
+
+    let dims = 5;
+    let params: Vec<FloatParam> = (0..dims)
+        .map(|i| FloatParam::new(-5.0, 5.0).name(format!("x{i}")))
+        .collect();
+    let n_trials = 200;
+
+    // Sphere: Random vs TPE
+    run_convergence(
+        "sphere_5d",
+        "random",
+        Study::minimize(RandomSampler::with_seed(1)),
+        &params,
+        functions::sphere,
+        n_trials,
+    );
+    run_convergence(
+        "sphere_5d",
+        "tpe",
+        Study::minimize(TpeSampler::builder().seed(1).build().unwrap()),
+        &params,
+        functions::sphere,
+        n_trials,
+    );
+
+    // Rosenbrock: Random vs TPE
+    run_convergence(
+        "rosenbrock_5d",
+        "random",
+        Study::minimize(RandomSampler::with_seed(2)),
+        &params,
+        functions::rosenbrock,
+        n_trials,
+    );
+    run_convergence(
+        "rosenbrock_5d",
+        "tpe",
+        Study::minimize(TpeSampler::builder().seed(2).build().unwrap()),
+        &params,
+        functions::rosenbrock,
+        n_trials,
+    );
+
+    // Rastrigin: Random vs TPE
+    run_convergence(
+        "rastrigin_5d",
+        "random",
+        Study::minimize(RandomSampler::with_seed(3)),
+        &params,
+        functions::rastrigin,
+        n_trials,
+    );
+    run_convergence(
+        "rastrigin_5d",
+        "tpe",
+        Study::minimize(TpeSampler::builder().seed(3).build().unwrap()),
+        &params,
+        functions::rastrigin,
+        n_trials,
+    );
+}

examples/visualization_demo.rs

@@ -0,0 +1,45 @@
+use optimizer::prelude::*;
+
+fn main() {
+    // Multi-parameter optimization with TPE sampler.
+    let sampler = TpeSampler::builder().seed(42).build().unwrap();
+    let mut study: Study<f64> = Study::with_sampler(Direction::Minimize, sampler);
+    study.set_pruner(MedianPruner::new(Direction::Minimize));
+
+    let lr = FloatParam::new(1e-5, 1e-1)
+        .log_scale()
+        .name("learning_rate");
+    let n_layers = IntParam::new(1, 5).name("n_layers");
+    let dropout = FloatParam::new(0.0, 0.5).step(0.05).name("dropout");
+    let batch_size = CategoricalParam::new(vec![16, 32, 64, 128]).name("batch_size");
+
+    study
+        .optimize(80, |trial| {
+            let lr_val = lr.suggest(trial)?;
+            let layers = n_layers.suggest(trial)?;
+            let drop = dropout.suggest(trial)?;
+            let bs = batch_size.suggest(trial)?;
+
+            // Simulate training with intermediate reporting.
+            let mut loss = 1.0;
+            for epoch in 0..10 {
+                loss *= 0.7 + 0.3 * lr_val.ln().abs() / 12.0;
+                loss += drop * 0.05;
+                loss += (1.0 / bs as f64) * 0.1;
+                loss -= layers as f64 * 0.02;
+                trial.report(epoch, loss);
+                if trial.should_prune() {
+                    return Err(TrialPruned.into());
+                }
+            }
+
+            Ok::<_, Error>(loss)
+        })
+        .unwrap();
+
+    println!("{}", study.summary());
+
+    let path = "optimization_report.html";
+    generate_html_report(&study, path).unwrap();
+    println!("\nReport saved to {path}");
+}

optimizer-derive/Cargo.toml

@@ -13,4 +13,3 @@ proc-macro = true
 [dependencies]
 syn = { version = "2", features = ["full"] }
 quote = "1"
-proc-macro2 = "1"

src/distribution.rs

@@ -2,6 +2,7 @@
 
 /// Distribution for floating-point parameters.
 #[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub struct FloatDistribution {
     /// Lower bound (inclusive).
     pub low: f64,
@@ -15,6 +16,7 @@ pub struct FloatDistribution {
 
 /// Distribution for integer parameters.
 #[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub struct IntDistribution {
     /// Lower bound (inclusive).
     pub low: i64,
@@ -28,6 +30,7 @@ pub struct IntDistribution {
 
 /// Distribution for categorical parameters.
 #[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub struct CategoricalDistribution {
     /// Number of choices available.
     pub n_choices: usize,
@@ -35,6 +38,7 @@ pub struct CategoricalDistribution {
 
 /// Enum wrapping all parameter distribution types.
 #[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub enum Distribution {
     /// A floating-point distribution.
     Float(FloatDistribution),

src/error.rs

@@ -76,6 +76,15 @@ pub enum Error {
     #[error("trial was pruned")]
     TrialPruned,
 
+    /// Returned when the objective returns the wrong number of values.
+    #[error("objective dimension mismatch: expected {expected} values, got {got}")]
+    ObjectiveDimensionMismatch {
+        /// The expected number of objective values.
+        expected: usize,
+        /// The actual number of objective values returned.
+        got: usize,
+    },
+
     /// Returned when an internal invariant is violated.
     #[error("internal error: {0}")]
     Internal(&'static str),

src/fanova.rs

@@ -0,0 +1,536 @@
+//! fANOVA (functional ANOVA) parameter importance via random forest.
+//!
+//! Decomposes the variance of the objective function into contributions
+//! from individual parameters (main effects) and parameter interactions.
+//!
+//! The algorithm:
+//! 1. Fits a random forest to `(parameters) -> objective_value`
+//! 2. Applies functional ANOVA decomposition to the forest
+//! 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 {
+    /// Per-parameter importance (fraction of total variance explained).
+    /// Sorted by descending importance.
+    pub main_effects: Vec<(String, f64)>,
+    /// Pairwise interaction importance (fraction of total variance explained).
+    /// Sorted by descending importance.
+    pub interactions: Vec<((String, String), f64)>,
+}
+
+/// Configuration for fANOVA analysis.
+#[derive(Debug, Clone)]
+pub struct FanovaConfig {
+    /// Number of trees in the random forest (default: 64).
+    pub n_trees: usize,
+    /// Maximum depth of each tree. `None` for unlimited (default: `None`).
+    pub max_depth: Option<usize>,
+    /// Minimum samples required to split a node (default: 2).
+    pub min_samples_split: usize,
+    /// Minimum samples required in a leaf node (default: 1).
+    pub min_samples_leaf: usize,
+    /// Random seed for reproducibility (default: `Some(42)`).
+    pub seed: Option<u64>,
+}
+
+impl Default for FanovaConfig {
+    fn default() -> Self {
+        Self {
+            n_trees: 64,
+            max_depth: None,
+            min_samples_split: 2,
+            min_samples_leaf: 1,
+            seed: Some(42),
+        }
+    }
+}
+
+// --- Decision Tree ---
+
+/// A node in the regression tree (arena-allocated).
+#[derive(Debug, Clone)]
+enum TreeNode {
+    Leaf {
+        value: f64,
+        n_samples: usize,
+    },
+    Split {
+        feature: usize,
+        threshold: f64,
+        left: usize,
+        right: usize,
+        n_samples: usize,
+    },
+}
+
+/// A regression decision tree for fANOVA.
+#[derive(Debug, Clone)]
+struct DecisionTree {
+    nodes: Vec<TreeNode>,
+}
+
+impl DecisionTree {
+    /// Build a tree from the given data using the specified bootstrap indices.
+    fn build(
+        data: &[Vec<f64>],
+        targets: &[f64],
+        indices: &[usize],
+        config: &FanovaConfig,
+        rng: &mut StdRng,
+    ) -> Self {
+        let mut tree = Self { nodes: Vec::new() };
+        tree.build_node(data, targets, indices, 0, config, rng);
+        tree
+    }
+
+    #[allow(clippy::cast_precision_loss)]
+    fn build_node(
+        &mut self,
+        data: &[Vec<f64>],
+        targets: &[f64],
+        indices: &[usize],
+        depth: usize,
+        config: &FanovaConfig,
+        rng: &mut StdRng,
+    ) -> usize {
+        let n = indices.len();
+        let mean = indices.iter().map(|&i| targets[i]).sum::<f64>() / n as f64;
+
+        // Stopping conditions
+        if n < config.min_samples_split || config.max_depth.is_some_and(|d| depth >= d) {
+            let idx = self.nodes.len();
+            self.nodes.push(TreeNode::Leaf {
+                value: mean,
+                n_samples: n,
+            });
+            return idx;
+        }
+
+        // Pure node check (all targets identical)
+        #[allow(clippy::float_cmp)]
+        if indices.iter().all(|&i| targets[i] == targets[indices[0]]) {
+            let idx = self.nodes.len();
+            self.nodes.push(TreeNode::Leaf {
+                value: mean,
+                n_samples: n,
+            });
+            return idx;
+        }
+
+        let n_features = data[0].len();
+        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
+        let max_features = ((n_features as f64).sqrt().ceil() as usize)
+            .max(1)
+            .min(n_features);
+        let candidates = partial_shuffle(n_features, max_features, rng);
+
+        // Total variance at this node
+        let total_var: f64 = indices.iter().map(|&i| (targets[i] - mean).powi(2)).sum();
+        if total_var == 0.0 {
+            let idx = self.nodes.len();
+            self.nodes.push(TreeNode::Leaf {
+                value: mean,
+                n_samples: n,
+            });
+            return idx;
+        }
+
+        let mut best_score = f64::NEG_INFINITY;
+        let mut best_feature = 0;
+        let mut best_threshold = 0.0;
+
+        for &feat in &candidates {
+            let mut values: Vec<f64> = indices.iter().map(|&i| data[i][feat]).collect();
+            values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
+            values.dedup();
+
+            if values.len() < 2 {
+                continue;
+            }
+
+            for w in values.windows(2) {
+                let threshold = f64::midpoint(w[0], w[1]);
+                let (l_sum, l_sq, l_n, r_sum, r_sq, r_n) =
+                    split_stats(data, targets, indices, feat, threshold);
+
+                if l_n < config.min_samples_leaf || r_n < config.min_samples_leaf {
+                    continue;
+                }
+
+                let l_var = l_sq - l_sum * l_sum / l_n as f64;
+                let r_var = r_sq - r_sum * r_sum / r_n as f64;
+                let score = total_var - l_var - r_var;
+
+                if score > best_score {
+                    best_score = score;
+                    best_feature = feat;
+                    best_threshold = threshold;
+                }
+            }
+        }
+
+        if best_score <= 0.0 {
+            let idx = self.nodes.len();
+            self.nodes.push(TreeNode::Leaf {
+                value: mean,
+                n_samples: n,
+            });
+            return idx;
+        }
+
+        let (left_indices, right_indices): (Vec<usize>, Vec<usize>) = indices
+            .iter()
+            .partition(|&&i| data[i][best_feature] <= best_threshold);
+
+        if left_indices.is_empty() || right_indices.is_empty() {
+            let idx = self.nodes.len();
+            self.nodes.push(TreeNode::Leaf {
+                value: mean,
+                n_samples: n,
+            });
+            return idx;
+        }
+
+        // Reserve slot for this split node (placeholder replaced below)
+        let node_idx = self.nodes.len();
+        self.nodes.push(TreeNode::Leaf {
+            value: 0.0,
+            n_samples: 0,
+        });
+
+        let left = self.build_node(data, targets, &left_indices, depth + 1, config, rng);
+        let right = self.build_node(data, targets, &right_indices, depth + 1, config, rng);
+
+        self.nodes[node_idx] = TreeNode::Split {
+            feature: best_feature,
+            threshold: best_threshold,
+            left,
+            right,
+            n_samples: n,
+        };
+
+        node_idx
+    }
+
+    /// Compute marginal prediction for a given feature subset.
+    ///
+    /// Features in `subset` use values from `feature_values`.
+    /// Features not in `subset` are marginalized by weighting branches
+    /// proportionally to their training-data fractions.
+    fn marginal_predict(&self, subset: &[usize], feature_values: &[f64]) -> f64 {
+        self.marginal_predict_at(0, subset, feature_values)
+    }
+
+    #[allow(clippy::cast_precision_loss)]
+    fn marginal_predict_at(&self, idx: usize, subset: &[usize], vals: &[f64]) -> f64 {
+        match self.nodes[idx] {
+            TreeNode::Leaf { value, .. } => value,
+            TreeNode::Split {
+                feature,
+                threshold,
+                left,
+                right,
+                n_samples,
+            } => {
+                if subset.contains(&feature) {
+                    if vals[feature] <= threshold {
+                        self.marginal_predict_at(left, subset, vals)
+                    } else {
+                        self.marginal_predict_at(right, subset, vals)
+                    }
+                } else {
+                    let l_n = self.n_samples(left) as f64;
+                    let r_n = self.n_samples(right) as f64;
+                    let total = n_samples as f64;
+                    (l_n / total) * self.marginal_predict_at(left, subset, vals)
+                        + (r_n / total) * self.marginal_predict_at(right, subset, vals)
+                }
+            }
+        }
+    }
+
+    fn n_samples(&self, idx: usize) -> usize {
+        match self.nodes[idx] {
+            TreeNode::Leaf { n_samples, .. } | TreeNode::Split { n_samples, .. } => n_samples,
+        }
+    }
+}
+
+// --- 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> {
+    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);
+        indices.swap(i, j);
+    }
+    indices.truncate(k);
+    indices
+}
+
+/// Compute left/right split statistics for variance reduction.
+#[allow(clippy::cast_precision_loss)]
+fn split_stats(
+    data: &[Vec<f64>],
+    targets: &[f64],
+    indices: &[usize],
+    feature: usize,
+    threshold: f64,
+) -> (f64, f64, usize, f64, f64, usize) {
+    let (mut l_sum, mut l_sq, mut l_n) = (0.0, 0.0, 0usize);
+    let (mut r_sum, mut r_sq, mut r_n) = (0.0, 0.0, 0usize);
+
+    for &i in indices {
+        let y = targets[i];
+        if data[i][feature] <= threshold {
+            l_sum += y;
+            l_sq += y * y;
+            l_n += 1;
+        } else {
+            r_sum += y;
+            r_sq += y * y;
+            r_n += 1;
+        }
+    }
+
+    (l_sum, l_sq, l_n, r_sum, r_sq, r_n)
+}
+
+/// Population variance of a slice.
+#[allow(clippy::cast_precision_loss)]
+fn variance(values: &[f64]) -> f64 {
+    if values.is_empty() {
+        return 0.0;
+    }
+    let n = values.len() as f64;
+    let mean = values.iter().sum::<f64>() / n;
+    values.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / n
+}
+
+// --- Public API ---
+
+/// Run fANOVA analysis on pre-processed numerical data.
+///
+/// `data` is `n_samples` rows, each with `n_features` columns.
+/// `targets` has one entry per sample.
+/// `feature_names` maps feature index to human-readable name.
+#[allow(clippy::cast_precision_loss)]
+pub(crate) fn compute_fanova(
+    data: &[Vec<f64>],
+    targets: &[f64],
+    feature_names: &[String],
+    config: &FanovaConfig,
+) -> FanovaResult {
+    let n_samples = data.len();
+    let n_features = data[0].len();
+
+    let mut rng: StdRng = config
+        .seed
+        .map_or_else(rand::make_rng, StdRng::seed_from_u64);
+
+    // 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();
+            DecisionTree::build(data, targets, &bootstrap, config, &mut rng)
+        })
+        .collect();
+
+    // Compute main effects: V_j = Var[E[f | x_j]]
+    let main_var: Vec<f64> = (0..n_features)
+        .map(|j| {
+            let subset = [j];
+            let preds: Vec<f64> = (0..n_samples)
+                .map(|i| {
+                    trees
+                        .iter()
+                        .map(|t| t.marginal_predict(&subset, &data[i]))
+                        .sum::<f64>()
+                        / trees.len() as f64
+                })
+                .collect();
+            variance(&preds)
+        })
+        .collect();
+
+    // Compute pairwise interaction effects: V_{j,k} - V_j - V_k
+    let mut interactions: Vec<((String, String), f64)> = Vec::new();
+    for j in 0..n_features {
+        for k in (j + 1)..n_features {
+            let subset = [j, k];
+            let preds: Vec<f64> = (0..n_samples)
+                .map(|i| {
+                    trees
+                        .iter()
+                        .map(|t| t.marginal_predict(&subset, &data[i]))
+                        .sum::<f64>()
+                        / trees.len() as f64
+                })
+                .collect();
+            let joint = variance(&preds);
+            let interaction = (joint - main_var[j] - main_var[k]).max(0.0);
+            if interaction > 1e-10 {
+                interactions.push((
+                    (feature_names[j].clone(), feature_names[k].clone()),
+                    interaction,
+                ));
+            }
+        }
+    }
+
+    // Normalize so all importances sum to 1.0
+    let total: f64 =
+        main_var.iter().sum::<f64>() + interactions.iter().map(|(_, v)| *v).sum::<f64>();
+
+    let mut main_effects: Vec<(String, f64)> = feature_names
+        .iter()
+        .zip(&main_var)
+        .map(|(name, &v)| (name.clone(), if total > 0.0 { v / total } else { 0.0 }))
+        .collect();
+    main_effects.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(core::cmp::Ordering::Equal));
+
+    if total > 0.0 {
+        for entry in &mut interactions {
+            entry.1 /= total;
+        }
+    }
+    interactions.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(core::cmp::Ordering::Equal));
+
+    FanovaResult {
+        main_effects,
+        interactions,
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn single_dominant_parameter() {
+        // f(x, y) = x — only x matters
+        let mut rng = StdRng::seed_from_u64(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)])
+            .collect();
+        let targets: Vec<f64> = data.iter().map(|row| row[0]).collect();
+
+        let result = compute_fanova(
+            &data,
+            &targets,
+            &["x".into(), "y".into()],
+            &FanovaConfig::default(),
+        );
+
+        assert_eq!(result.main_effects[0].0, "x");
+        assert!(
+            result.main_effects[0].1 > 0.8,
+            "x importance = {}",
+            result.main_effects[0].1
+        );
+    }
+
+    #[test]
+    fn interaction_detection() {
+        // f(x, y) = x * y — both matter and interact
+        let mut rng = StdRng::seed_from_u64(0);
+        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)])
+            .collect();
+        let targets: Vec<f64> = data.iter().map(|row| row[0] * row[1]).collect();
+
+        let config = FanovaConfig {
+            n_trees: 128,
+            ..FanovaConfig::default()
+        };
+        let result = compute_fanova(&data, &targets, &["x".into(), "y".into()], &config);
+
+        assert!(
+            !result.interactions.is_empty(),
+            "should detect x*y interaction"
+        );
+        assert!(
+            result.interactions[0].1 > 0.05,
+            "interaction importance = {}",
+            result.interactions[0].1
+        );
+    }
+
+    #[test]
+    fn variance_computation() {
+        assert!((variance(&[1.0, 2.0, 3.0, 4.0, 5.0]) - 2.0).abs() < 1e-10);
+        assert!(variance(&[5.0, 5.0, 5.0]).abs() < 1e-10);
+        assert!(variance(&[]).abs() < 1e-10);
+    }
+
+    #[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 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),
+                ]
+            })
+            .collect();
+        let targets: Vec<f64> = data.iter().map(|r| 3.0 * r[0] + 0.1 * r[1]).collect();
+
+        let result = compute_fanova(
+            &data,
+            &targets,
+            &["x".into(), "y".into(), "z".into()],
+            &FanovaConfig::default(),
+        );
+
+        // x should be the most important
+        assert_eq!(result.main_effects[0].0, "x");
+        assert!(result.main_effects[0].1 > 0.5);
+
+        // z should have near-zero importance
+        let z_imp = result
+            .main_effects
+            .iter()
+            .find(|(name, _)| name == "z")
+            .map_or(0.0, |(_, v)| *v);
+        assert!(z_imp < 0.1, "z importance = {z_imp}");
+    }
+
+    #[test]
+    fn importances_sum_to_one() {
+        let mut rng = StdRng::seed_from_u64(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)])
+            .collect();
+        let targets: Vec<f64> = data.iter().map(|r| r[0] + r[1]).collect();
+
+        let result = compute_fanova(
+            &data,
+            &targets,
+            &["x".into(), "y".into()],
+            &FanovaConfig::default(),
+        );
+
+        let total: f64 = result.main_effects.iter().map(|(_, v)| *v).sum::<f64>()
+            + result.interactions.iter().map(|(_, v)| *v).sum::<f64>();
+        assert!(
+            (total - 1.0).abs() < 1e-10,
+            "importances should sum to 1.0, got {total}"
+        );
+    }
+}

src/importance.rs

@@ -0,0 +1,93 @@
+//! Parameter importance via Spearman rank correlation.
+
+/// Assign average ranks to a slice of `f64` values (handles ties).
+#[allow(clippy::cast_precision_loss, clippy::float_cmp)]
+pub(crate) fn rank(values: &[f64]) -> Vec<f64> {
+    let n = values.len();
+    let mut indexed: Vec<(usize, f64)> = values.iter().copied().enumerate().collect();
+    indexed.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(core::cmp::Ordering::Equal));
+
+    let mut ranks = vec![0.0; n];
+    let mut i = 0;
+    while i < n {
+        // Find the run of tied values.
+        let mut j = i + 1;
+        while j < n && indexed[j].1 == indexed[i].1 {
+            j += 1;
+        }
+        // Average rank for the tie group (1-based ranks).
+        let avg = (i + 1..=j).sum::<usize>() as f64 / (j - i) as f64;
+        for item in &indexed[i..j] {
+            ranks[item.0] = avg;
+        }
+        i = j;
+    }
+    ranks
+}
+
+/// Pearson correlation coefficient on two equal-length slices.
+#[allow(clippy::cast_precision_loss)]
+fn pearson(x: &[f64], y: &[f64]) -> f64 {
+    let n = x.len() as f64;
+    let mean_x = x.iter().sum::<f64>() / n;
+    let mean_y = y.iter().sum::<f64>() / n;
+
+    let mut cov = 0.0;
+    let mut var_x = 0.0;
+    let mut var_y = 0.0;
+    for (xi, yi) in x.iter().zip(y.iter()) {
+        let dx = xi - mean_x;
+        let dy = yi - mean_y;
+        cov += dx * dy;
+        var_x += dx * dx;
+        var_y += dy * dy;
+    }
+
+    let denom = (var_x * var_y).sqrt();
+    if denom == 0.0 { 0.0 } else { cov / denom }
+}
+
+/// Spearman rank correlation (Pearson on ranks).
+pub(crate) fn spearman(x: &[f64], y: &[f64]) -> f64 {
+    pearson(&rank(x), &rank(y))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn rank_no_ties() {
+        let ranks = rank(&[30.0, 10.0, 20.0]);
+        assert_eq!(ranks, vec![3.0, 1.0, 2.0]);
+    }
+
+    #[test]
+    fn rank_with_ties() {
+        let ranks = rank(&[10.0, 20.0, 20.0, 30.0]);
+        assert_eq!(ranks, vec![1.0, 2.5, 2.5, 4.0]);
+    }
+
+    #[test]
+    fn perfect_positive_correlation() {
+        let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let y = vec![2.0, 4.0, 6.0, 8.0, 10.0];
+        let r = spearman(&x, &y);
+        assert!((r - 1.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn perfect_negative_correlation() {
+        let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let y = vec![10.0, 8.0, 6.0, 4.0, 2.0];
+        let r = spearman(&x, &y);
+        assert!((r + 1.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn zero_variance_returns_zero() {
+        let x = vec![5.0, 5.0, 5.0];
+        let y = vec![1.0, 2.0, 3.0];
+        assert!(spearman(&x, &y).abs() < f64::EPSILON);
+    }
+}

src/lib.rs

@@ -17,6 +17,11 @@
 //! - **Random Search** - Simple random sampling for baseline comparisons
 //! - **TPE (Tree-Parzen Estimator)** - Bayesian optimization for efficient search
 //! - **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)
+//! - **BOHB** - Bayesian Optimization + `HyperBand` for budget-aware TPE sampling
+//! - **NSGA-II** - Non-dominated Sorting Genetic Algorithm II for multi-objective optimization
+//! - **MOTPE** - Multi-Objective Tree-Parzen Estimator for Bayesian multi-objective optimization
 //!
 //! Additional features include:
 //!
@@ -182,19 +187,58 @@
 //!
 //! - `async`: Enable async optimization methods (requires tokio)
 //! - `derive`: Enable `#[derive(Categorical)]` for enum parameters
+//! - `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
+//! - `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
+
+/// Emit a `tracing::info!` event when the `tracing` feature is enabled.
+/// No-op otherwise.
+#[cfg(feature = "tracing")]
+macro_rules! trace_info {
+    ($($arg:tt)*) => { tracing::info!($($arg)*) };
+}
+
+#[cfg(not(feature = "tracing"))]
+macro_rules! trace_info {
+    ($($arg:tt)*) => {};
+}
+
+/// Emit a `tracing::debug!` event when the `tracing` feature is enabled.
+/// No-op otherwise.
+#[cfg(feature = "tracing")]
+macro_rules! trace_debug {
+    ($($arg:tt)*) => { tracing::debug!($($arg)*) };
+}
+
+#[cfg(not(feature = "tracing"))]
+macro_rules! trace_debug {
+    ($($arg:tt)*) => {};
+}
 
 mod distribution;
 mod error;
+#[cfg(feature = "fanova")]
+mod fanova;
+mod importance;
 mod kde;
+pub mod multi_objective;
 mod param;
 pub mod parameter;
+pub mod pareto;
 pub mod pruner;
 pub mod sampler;
 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")]
 pub use optimizer_derive::Categorical;
 pub use param::ParamValue;
@@ -206,12 +250,23 @@ pub use pruner::{
     SuccessiveHalvingPruner, ThresholdPruner, WilcoxonPruner,
 };
 pub use sampler::CompletedTrial;
+pub use sampler::bohb::BohbSampler;
+#[cfg(feature = "cma-es")]
+pub use sampler::cma_es::CmaEsSampler;
 pub use sampler::grid::GridSearchSampler;
+pub use sampler::motpe::MotpeSampler;
+pub use sampler::nsga2::Nsga2Sampler;
 pub use sampler::random::RandomSampler;
+#[cfg(feature = "sobol")]
+pub use sampler::sobol::SobolSampler;
 pub use sampler::tpe::TpeSampler;
 pub use study::Study;
+#[cfg(feature = "serde")]
+pub use study::StudySnapshot;
 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.
 ///
@@ -223,6 +278,9 @@ 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;
     pub use crate::parameter::{
         BoolParam, Categorical, CategoricalParam, EnumParam, FloatParam, IntParam, Parameter,
@@ -232,10 +290,21 @@ pub mod prelude {
         SuccessiveHalvingPruner, ThresholdPruner,
     };
     pub use crate::sampler::CompletedTrial;
+    pub use crate::sampler::bohb::BohbSampler;
+    #[cfg(feature = "cma-es")]
+    pub use crate::sampler::cma_es::CmaEsSampler;
     pub use crate::sampler::grid::GridSearchSampler;
+    pub use crate::sampler::motpe::MotpeSampler;
+    pub use crate::sampler::nsga2::Nsga2Sampler;
     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 = "serde")]
+    pub use crate::study::StudySnapshot;
     pub use crate::trial::{AttrValue, Trial};
     pub use crate::types::Direction;
+    #[cfg(feature = "visualization")]
+    pub use crate::visualization::generate_html_report;
 }

src/multi_objective.rs

@@ -0,0 +1,412 @@
+//! Multi-objective optimization via a dedicated study type.
+//!
+//! [`MultiObjectiveStudy`] manages trials that return multiple objective
+//! values. It supports arbitrary numbers of objectives with per-objective
+//! directions (minimize or maximize). Use [`pareto_front()`](MultiObjectiveStudy::pareto_front)
+//! to retrieve the Pareto-optimal solutions.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::Direction;
+//! use optimizer::multi_objective::MultiObjectiveStudy;
+//! use optimizer::parameter::{FloatParam, Parameter};
+//!
+//! let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+//! let x = FloatParam::new(0.0, 1.0);
+//!
+//! study
+//!     .optimize(20, |trial| {
+//!         let xv = x.suggest(trial)?;
+//!         Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+//!     })
+//!     .unwrap();
+//!
+//! let front = study.pareto_front();
+//! assert!(!front.is_empty());
+//! ```
+
+use core::sync::atomic::{AtomicU64, Ordering};
+use std::collections::HashMap;
+use std::sync::Arc;
+
+use parking_lot::RwLock;
+
+use crate::distribution::Distribution;
+use crate::param::ParamValue;
+use crate::parameter::{ParamId, Parameter};
+use crate::pruner::NopPruner;
+use crate::sampler::random::RandomSampler;
+use crate::sampler::{CompletedTrial, Sampler};
+use crate::trial::{AttrValue, Trial};
+use crate::types::{Direction, TrialState};
+
+// ---------------------------------------------------------------------------
+// MultiObjectiveTrial
+// ---------------------------------------------------------------------------
+
+/// A completed trial with multiple objective values.
+#[derive(Clone, Debug)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
+pub struct MultiObjectiveTrial {
+    /// The unique identifier for this trial.
+    pub id: u64,
+    /// The sampled parameter values, keyed by parameter id.
+    pub params: HashMap<ParamId, ParamValue>,
+    /// The parameter distributions used, keyed by parameter id.
+    pub distributions: HashMap<ParamId, Distribution>,
+    /// Human-readable labels for parameters, keyed by parameter id.
+    pub param_labels: HashMap<ParamId, String>,
+    /// The objective values (one per objective).
+    pub values: Vec<f64>,
+    /// The state of the trial.
+    pub state: TrialState,
+    /// User-defined attributes stored during the trial.
+    pub user_attrs: HashMap<String, AttrValue>,
+    /// Constraint values for this trial (<=0.0 means feasible).
+    #[cfg_attr(feature = "serde", serde(default))]
+    pub constraints: Vec<f64>,
+}
+
+impl MultiObjectiveTrial {
+    /// Returns the typed value for the given parameter.
+    ///
+    /// Returns `None` if the parameter was not used in this trial.
+    ///
+    /// # Panics
+    ///
+    /// Panics if the stored value is incompatible with the parameter type.
+    pub fn get<P: Parameter>(&self, param: &P) -> Option<P::Value> {
+        self.params.get(&param.id()).map(|v| {
+            param
+                .cast_param_value(v)
+                .expect("parameter type mismatch: stored value incompatible with parameter")
+        })
+    }
+
+    /// Returns `true` if all constraints are satisfied (values <= 0.0).
+    ///
+    /// A trial with no constraints is considered feasible.
+    #[must_use]
+    pub fn is_feasible(&self) -> bool {
+        self.constraints.iter().all(|&c| c <= 0.0)
+    }
+
+    /// Gets a user attribute by key.
+    #[must_use]
+    pub fn user_attr(&self, key: &str) -> Option<&AttrValue> {
+        self.user_attrs.get(key)
+    }
+
+    /// Returns all user attributes.
+    #[must_use]
+    pub fn user_attrs(&self) -> &HashMap<String, AttrValue> {
+        &self.user_attrs
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MultiObjectiveSampler trait
+// ---------------------------------------------------------------------------
+
+/// Trait for samplers aware of multi-objective history.
+///
+/// Separate from [`Sampler`] because NSGA-II needs access to
+/// `&[MultiObjectiveTrial]` (with vector-valued objectives) and
+/// `&[Direction]` (one direction per objective).
+pub trait MultiObjectiveSampler: Send + Sync {
+    /// Samples a parameter value from the given distribution.
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[MultiObjectiveTrial],
+        directions: &[Direction],
+    ) -> ParamValue;
+}
+
+// ---------------------------------------------------------------------------
+// RandomMultiObjectiveSampler
+// ---------------------------------------------------------------------------
+
+/// Default MO sampler that delegates to [`RandomSampler`].
+pub(crate) struct RandomMultiObjectiveSampler(RandomSampler);
+
+impl RandomMultiObjectiveSampler {
+    pub(crate) fn new() -> Self {
+        Self(RandomSampler::new())
+    }
+}
+
+impl MultiObjectiveSampler for RandomMultiObjectiveSampler {
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        _history: &[MultiObjectiveTrial],
+        _directions: &[Direction],
+    ) -> ParamValue {
+        self.0.sample(distribution, trial_id, &[])
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MoSamplerBridge — bridges MultiObjectiveSampler to Sampler trait
+// ---------------------------------------------------------------------------
+
+/// Bridges a [`MultiObjectiveSampler`] to the [`Sampler`] trait so that
+/// `Trial::with_sampler()` can use it.
+struct MoSamplerBridge {
+    inner: Arc<dyn MultiObjectiveSampler>,
+    history: Arc<RwLock<Vec<MultiObjectiveTrial>>>,
+    directions: Vec<Direction>,
+}
+
+impl Sampler for MoSamplerBridge {
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        _history: &[CompletedTrial],
+    ) -> ParamValue {
+        let mo_history = self.history.read();
+        self.inner
+            .sample(distribution, trial_id, &mo_history, &self.directions)
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MultiObjectiveStudy
+// ---------------------------------------------------------------------------
+
+/// A study for multi-objective optimization.
+///
+/// Manages trials that return multiple objective values. Supports
+/// arbitrary numbers of objectives with independent minimize/maximize
+/// directions.
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::Direction;
+/// use optimizer::multi_objective::MultiObjectiveStudy;
+/// use optimizer::parameter::{FloatParam, Parameter};
+///
+/// // Bi-objective: minimize both
+/// let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+/// 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();
+///
+/// let front = study.pareto_front();
+/// assert!(!front.is_empty());
+/// ```
+pub struct MultiObjectiveStudy {
+    directions: Vec<Direction>,
+    sampler: Arc<dyn MultiObjectiveSampler>,
+    completed_trials: Arc<RwLock<Vec<MultiObjectiveTrial>>>,
+    next_trial_id: AtomicU64,
+}
+
+impl MultiObjectiveStudy {
+    /// Creates a new multi-objective study with the given directions.
+    ///
+    /// Uses a random sampler by default.
+    ///
+    /// # Arguments
+    ///
+    /// * `directions` - One direction per objective (minimize or maximize).
+    #[must_use]
+    pub fn new(directions: Vec<Direction>) -> Self {
+        Self {
+            directions,
+            sampler: Arc::new(RandomMultiObjectiveSampler::new()),
+            completed_trials: Arc::new(RwLock::new(Vec::new())),
+            next_trial_id: AtomicU64::new(0),
+        }
+    }
+
+    /// Creates a new study with a custom multi-objective sampler.
+    #[must_use]
+    pub fn with_sampler(
+        directions: Vec<Direction>,
+        sampler: impl MultiObjectiveSampler + 'static,
+    ) -> Self {
+        Self {
+            directions,
+            sampler: Arc::new(sampler),
+            completed_trials: Arc::new(RwLock::new(Vec::new())),
+            next_trial_id: AtomicU64::new(0),
+        }
+    }
+
+    /// Returns the optimization directions.
+    #[must_use]
+    pub fn directions(&self) -> &[Direction] {
+        &self.directions
+    }
+
+    /// Returns the number of objectives.
+    #[must_use]
+    pub fn n_objectives(&self) -> usize {
+        self.directions.len()
+    }
+
+    /// Returns the number of completed trials.
+    #[must_use]
+    pub fn n_trials(&self) -> usize {
+        self.completed_trials.read().len()
+    }
+
+    /// Returns all completed trials.
+    #[must_use]
+    pub fn trials(&self) -> Vec<MultiObjectiveTrial> {
+        self.completed_trials.read().clone()
+    }
+
+    /// Returns the Pareto-optimal trials (front 0).
+    #[must_use]
+    pub fn pareto_front(&self) -> Vec<MultiObjectiveTrial> {
+        let trials = self.completed_trials.read();
+        let complete: Vec<_> = trials
+            .iter()
+            .filter(|t| t.state == TrialState::Complete)
+            .collect();
+
+        if complete.is_empty() {
+            return Vec::new();
+        }
+
+        let values: Vec<Vec<f64>> = complete.iter().map(|t| t.values.clone()).collect();
+        let fronts = crate::pareto::fast_non_dominated_sort(&values, &self.directions);
+
+        if fronts.is_empty() {
+            return Vec::new();
+        }
+
+        fronts[0].iter().map(|&i| complete[i].clone()).collect()
+    }
+
+    /// Creates a new trial wired to the study's MO sampler.
+    fn create_trial(&self) -> Trial {
+        let id = self.next_trial_id.fetch_add(1, Ordering::SeqCst);
+
+        let bridge: Arc<dyn Sampler> = Arc::new(MoSamplerBridge {
+            inner: Arc::clone(&self.sampler),
+            history: Arc::clone(&self.completed_trials),
+            directions: self.directions.clone(),
+        });
+
+        // Dummy f64 history — the bridge ignores it.
+        let dummy_history: Arc<RwLock<Vec<CompletedTrial<f64>>>> =
+            Arc::new(RwLock::new(Vec::new()));
+
+        Trial::with_sampler(id, bridge, dummy_history, Arc::new(NopPruner))
+    }
+
+    /// Records a completed trial.
+    fn complete_trial(&self, mut trial: Trial, values: Vec<f64>) {
+        trial.set_complete();
+        let mo_trial = MultiObjectiveTrial {
+            id: trial.id(),
+            params: trial.params().clone(),
+            distributions: trial.distributions().clone(),
+            param_labels: trial.param_labels().clone(),
+            values,
+            state: TrialState::Complete,
+            user_attrs: trial.user_attrs().clone(),
+            constraints: trial.constraint_values().to_vec(),
+        };
+        self.completed_trials.write().push(mo_trial);
+    }
+
+    /// Records a failed trial (not stored in history).
+    fn fail_trial(trial: &mut Trial) {
+        trial.set_failed();
+    }
+
+    /// Request a new trial for the ask/tell interface.
+    ///
+    /// After creating the trial, suggest parameters on it, evaluate your
+    /// objective externally, then pass the trial back to [`tell()`](Self::tell).
+    pub fn ask(&self) -> Trial {
+        self.create_trial()
+    }
+
+    /// Report the result of a trial obtained from [`ask()`](Self::ask).
+    ///
+    /// Pass `Ok(values)` for a successful evaluation or `Err(reason)` for a failure.
+    ///
+    /// # Errors
+    ///
+    /// Returns `ObjectiveDimensionMismatch` if the number of values doesn't
+    /// match the number of directions.
+    pub fn tell(
+        &self,
+        mut trial: Trial,
+        result: core::result::Result<Vec<f64>, impl ToString>,
+    ) -> crate::Result<()> {
+        if let Ok(values) = result {
+            if values.len() != self.directions.len() {
+                return Err(crate::Error::ObjectiveDimensionMismatch {
+                    expected: self.directions.len(),
+                    got: values.len(),
+                });
+            }
+            self.complete_trial(trial, values);
+        } else {
+            Self::fail_trial(&mut trial);
+        }
+        Ok(())
+    }
+
+    /// Runs multi-objective optimization for `n_trials` trials.
+    ///
+    /// The objective function must return a `Vec<f64>` with one value per
+    /// objective.
+    ///
+    /// # Errors
+    ///
+    /// Returns `ObjectiveDimensionMismatch` if the objective returns the wrong
+    /// number of values. Returns `NoCompletedTrials` if all trials fail.
+    pub fn optimize<F, E>(&self, n_trials: usize, mut objective: F) -> crate::Result<()>
+    where
+        F: FnMut(&mut Trial) -> core::result::Result<Vec<f64>, E>,
+        E: ToString,
+    {
+        for _ in 0..n_trials {
+            let mut trial = self.create_trial();
+
+            match objective(&mut trial) {
+                Ok(values) => {
+                    if values.len() != self.directions.len() {
+                        return Err(crate::Error::ObjectiveDimensionMismatch {
+                            expected: self.directions.len(),
+                            got: values.len(),
+                        });
+                    }
+                    self.complete_trial(trial, values);
+                }
+                Err(_) => {
+                    Self::fail_trial(&mut trial);
+                }
+            }
+        }
+
+        let has_complete = self
+            .completed_trials
+            .read()
+            .iter()
+            .any(|t| t.state == TrialState::Complete);
+        if !has_complete {
+            return Err(crate::Error::NoCompletedTrials);
+        }
+
+        Ok(())
+    }
+}

src/param.rs

@@ -6,6 +6,7 @@
 /// For categorical parameters, the `Categorical` variant stores
 /// the index into the choices array.
 #[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub enum ParamValue {
     /// A floating-point parameter value.
     Float(f64),

src/parameter.rs

@@ -37,7 +37,8 @@ static NEXT_PARAM_ID: AtomicU64 = AtomicU64::new(0);
 ///
 /// Each parameter is assigned a unique `ParamId` at creation time. Cloning a parameter
 /// copies its `ParamId`, so clones refer to the same logical parameter.
-#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
+#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub struct ParamId(u64);
 
 impl ParamId {

src/pareto.rs

@@ -0,0 +1,593 @@
+//! Pareto front analysis utilities for multi-objective optimization.
+//!
+//! Provides functions for analyzing and working with Pareto fronts:
+//!
+//! - [`hypervolume`] — measure the quality of a Pareto front
+//! - [`non_dominated_sort`] — rank solutions into successive fronts
+//! - [`pareto_front_indices`] — filter to non-dominated solutions only
+//! - [`crowding_distance`] — measure diversity within a front
+//!
+//! Internally also provides fast non-dominated sorting (Deb et al., 2002)
+//! used by [`MultiObjectiveStudy::pareto_front()`](crate::MultiObjectiveStudy::pareto_front)
+//! and [`Nsga2Sampler`](crate::Nsga2Sampler).
+
+use crate::types::Direction;
+
+/// Returns `true` if solution `a` Pareto-dominates solution `b`.
+///
+/// A solution dominates another if it is at least as good in all objectives
+/// and strictly better in at least one, respecting the given directions.
+#[allow(clippy::module_name_repetitions)]
+pub(crate) fn dominates(a: &[f64], b: &[f64], directions: &[Direction]) -> bool {
+    debug_assert_eq!(a.len(), b.len());
+    debug_assert_eq!(a.len(), directions.len());
+
+    let mut strictly_better = false;
+    for ((&av, &bv), dir) in a.iter().zip(b.iter()).zip(directions.iter()) {
+        let better = match dir {
+            Direction::Minimize => av < bv,
+            Direction::Maximize => av > bv,
+        };
+        let worse = match dir {
+            Direction::Minimize => av > bv,
+            Direction::Maximize => av < bv,
+        };
+        if worse {
+            return false;
+        }
+        if better {
+            strictly_better = true;
+        }
+    }
+    strictly_better
+}
+
+/// Constrained dominance: feasible beats infeasible, among infeasible
+/// prefer lower total constraint violation, among feasible use Pareto dominance.
+pub(crate) fn constrained_dominates(
+    a_values: &[f64],
+    b_values: &[f64],
+    a_constraints: &[f64],
+    b_constraints: &[f64],
+    directions: &[Direction],
+) -> bool {
+    let a_feasible = a_constraints.iter().all(|&c| c <= 0.0);
+    let b_feasible = b_constraints.iter().all(|&c| c <= 0.0);
+
+    match (a_feasible, b_feasible) {
+        (true, false) => true,
+        (false, true) => false,
+        (false, false) => {
+            let a_violation: f64 = a_constraints.iter().map(|c| c.max(0.0)).sum();
+            let b_violation: f64 = b_constraints.iter().map(|c| c.max(0.0)).sum();
+            a_violation < b_violation
+        }
+        (true, true) => dominates(a_values, b_values, directions),
+    }
+}
+
+/// Fast non-dominated sorting (Deb et al., 2002).
+///
+/// Returns `Vec<Vec<usize>>` where `fronts[0]` is the Pareto front,
+/// each inner vec contains indices into `values`.
+///
+/// Complexity: O(M * N^2) where M = objectives, N = solutions.
+#[allow(clippy::cast_possible_truncation)]
+pub(crate) fn fast_non_dominated_sort(
+    values: &[Vec<f64>],
+    directions: &[Direction],
+) -> Vec<Vec<usize>> {
+    fast_non_dominated_sort_constrained(values, directions, &[])
+}
+
+/// Fast non-dominated sorting with constraint support.
+///
+/// `constraints` is either empty (no constraints) or has the same length
+/// as `values`, where each entry is the constraint vector for that solution.
+#[allow(clippy::cast_possible_truncation)]
+pub(crate) fn fast_non_dominated_sort_constrained(
+    values: &[Vec<f64>],
+    directions: &[Direction],
+    constraints: &[Vec<f64>],
+) -> Vec<Vec<usize>> {
+    let n = values.len();
+    if n == 0 {
+        return Vec::new();
+    }
+
+    let has_constraints = !constraints.is_empty();
+    let empty_constraints: Vec<f64> = Vec::new();
+
+    // S_p: set of solutions dominated by p
+    let mut dominated_by: Vec<Vec<usize>> = vec![Vec::new(); n];
+    // n_p: domination count for p
+    let mut domination_count: Vec<usize> = vec![0; n];
+
+    for i in 0..n {
+        for j in (i + 1)..n {
+            let (a_c, b_c) = if has_constraints {
+                (&constraints[i], &constraints[j])
+            } else {
+                (&empty_constraints, &empty_constraints)
+            };
+
+            let i_dom_j = if has_constraints {
+                constrained_dominates(&values[i], &values[j], a_c, b_c, directions)
+            } else {
+                dominates(&values[i], &values[j], directions)
+            };
+            let j_dom_i = if has_constraints {
+                constrained_dominates(&values[j], &values[i], b_c, a_c, directions)
+            } else {
+                dominates(&values[j], &values[i], directions)
+            };
+
+            if i_dom_j {
+                dominated_by[i].push(j);
+                domination_count[j] += 1;
+            } else if j_dom_i {
+                dominated_by[j].push(i);
+                domination_count[i] += 1;
+            }
+        }
+    }
+
+    let mut fronts: Vec<Vec<usize>> = Vec::new();
+    let mut current_front: Vec<usize> = (0..n).filter(|&i| domination_count[i] == 0).collect();
+
+    while !current_front.is_empty() {
+        let mut next_front: Vec<usize> = Vec::new();
+        for &p in &current_front {
+            for &q in &dominated_by[p] {
+                domination_count[q] -= 1;
+                if domination_count[q] == 0 {
+                    next_front.push(q);
+                }
+            }
+        }
+        fronts.push(current_front);
+        current_front = next_front;
+    }
+
+    fronts
+}
+
+/// Crowding distance for one front (index-based, internal API).
+///
+/// Boundary solutions get `f64::INFINITY`. Returns one distance value per
+/// solution in the front, in the same order as `front_indices`.
+#[allow(clippy::cast_precision_loss)]
+pub(crate) fn crowding_distance_indexed(front_indices: &[usize], values: &[Vec<f64>]) -> Vec<f64> {
+    let n = front_indices.len();
+    if n <= 2 {
+        return vec![f64::INFINITY; n];
+    }
+
+    let m = values[front_indices[0]].len(); // number of objectives
+    let mut distances = vec![0.0_f64; n];
+
+    // Helper to look up objective value for a front member.
+    let val = |front_pos: usize, obj: usize| -> f64 { values[front_indices[front_pos]][obj] };
+
+    for obj in 0..m {
+        // Sort front positions by this objective
+        let mut sorted: Vec<usize> = (0..n).collect();
+        sorted.sort_by(|&a, &b| {
+            val(a, obj)
+                .partial_cmp(&val(b, obj))
+                .unwrap_or(core::cmp::Ordering::Equal)
+        });
+
+        // Boundary solutions get infinity
+        distances[sorted[0]] = f64::INFINITY;
+        distances[sorted[n - 1]] = f64::INFINITY;
+
+        let range = val(sorted[n - 1], obj) - val(sorted[0], obj);
+        if range > 0.0 {
+            for i in 1..(n - 1) {
+                distances[sorted[i]] += (val(sorted[i + 1], obj) - val(sorted[i - 1], obj)) / range;
+            }
+        }
+    }
+
+    distances
+}
+
+// ---------------------------------------------------------------------------
+// Public API
+// ---------------------------------------------------------------------------
+
+/// Compute the hypervolume indicator of a Pareto front.
+///
+/// The hypervolume is the volume of the objective space dominated by
+/// the Pareto front and bounded by a reference point. Higher values
+/// indicate a better front.
+///
+/// Each entry in `front` is one solution's objective values.
+/// `reference_point` should be worse than all front members in every
+/// objective (e.g., the worst acceptable values).
+///
+/// # Panics
+///
+/// Panics (in debug) if dimensions of `front`, `reference_point`, and
+/// `directions` are inconsistent.
+#[must_use]
+#[allow(clippy::cast_precision_loss)]
+pub fn hypervolume(front: &[Vec<f64>], reference_point: &[f64], directions: &[Direction]) -> f64 {
+    if front.is_empty() {
+        return 0.0;
+    }
+    let d = reference_point.len();
+    debug_assert!(front.iter().all(|p| p.len() == d));
+    debug_assert_eq!(d, directions.len());
+
+    // Normalize to minimize-space (negate maximized objectives).
+    let normalized: Vec<Vec<f64>> = front
+        .iter()
+        .map(|p| {
+            p.iter()
+                .zip(directions)
+                .map(|(&v, dir)| match dir {
+                    Direction::Minimize => v,
+                    Direction::Maximize => -v,
+                })
+                .collect()
+        })
+        .collect();
+
+    let ref_norm: Vec<f64> = reference_point
+        .iter()
+        .zip(directions)
+        .map(|(&v, dir)| match dir {
+            Direction::Minimize => v,
+            Direction::Maximize => -v,
+        })
+        .collect();
+
+    // Keep only points strictly dominated by the reference point.
+    let filtered: Vec<Vec<f64>> = normalized
+        .into_iter()
+        .filter(|p| p.iter().zip(&ref_norm).all(|(&pv, &rv)| pv < rv))
+        .collect();
+
+    if filtered.is_empty() {
+        return 0.0;
+    }
+
+    hv_recursive(&filtered, &ref_norm)
+}
+
+/// Recursive hypervolume via slicing on the last objective.
+///
+/// All points are in minimize-space and dominated by `reference`.
+#[allow(clippy::cast_precision_loss)]
+fn hv_recursive(points: &[Vec<f64>], reference: &[f64]) -> f64 {
+    let d = reference.len();
+
+    // Base case: 1-D hypervolume is just the gap from the best point to ref.
+    if d == 1 {
+        let min_val = points.iter().map(|p| p[0]).fold(f64::INFINITY, f64::min);
+        return (reference[0] - min_val).max(0.0);
+    }
+
+    // Single point: hypervolume is the product of gaps.
+    if points.len() == 1 {
+        return points[0]
+            .iter()
+            .zip(reference)
+            .map(|(&p, &r)| (r - p).max(0.0))
+            .product();
+    }
+
+    // Sort by last objective ascending.
+    let mut sorted: Vec<&Vec<f64>> = points.iter().collect();
+    sorted.sort_by(|a, b| {
+        a[d - 1]
+            .partial_cmp(&b[d - 1])
+            .unwrap_or(core::cmp::Ordering::Equal)
+    });
+
+    let sub_ref: Vec<f64> = reference[..d - 1].to_vec();
+    let mut result = 0.0;
+
+    for i in 0..sorted.len() {
+        let height = if i + 1 < sorted.len() {
+            sorted[i + 1][d - 1] - sorted[i][d - 1]
+        } else {
+            reference[d - 1] - sorted[i][d - 1]
+        };
+
+        if height <= 0.0 {
+            continue;
+        }
+
+        // Project points[0..=i] onto the first d-1 dimensions and
+        // keep only the non-dominated subset.
+        let projected: Vec<Vec<f64>> = sorted[..=i].iter().map(|p| p[..d - 1].to_vec()).collect();
+        let non_dom = non_dominated_minimize(&projected);
+
+        if !non_dom.is_empty() {
+            result += height * hv_recursive(&non_dom, &sub_ref);
+        }
+    }
+
+    result
+}
+
+/// Return the non-dominated subset of `points` in minimize-space.
+fn non_dominated_minimize(points: &[Vec<f64>]) -> Vec<Vec<f64>> {
+    let mut result = Vec::new();
+    'outer: for (i, p) in points.iter().enumerate() {
+        for (j, q) in points.iter().enumerate() {
+            if i == j {
+                continue;
+            }
+            // Check if q dominates p (all <=, at least one <).
+            let mut all_leq = true;
+            let mut any_lt = false;
+            for (&qv, &pv) in q.iter().zip(p.iter()) {
+                if qv > pv {
+                    all_leq = false;
+                    break;
+                }
+                if qv < pv {
+                    any_lt = true;
+                }
+            }
+            if all_leq && any_lt {
+                continue 'outer;
+            }
+        }
+        result.push(p.clone());
+    }
+    result
+}
+
+/// Compute non-dominated sorting of a set of solutions.
+///
+/// Returns a vec of fronts, where `fronts[0]` is the Pareto front,
+/// `fronts[1]` is the next best, etc. Each inner vec contains indices
+/// into the original `solutions` slice.
+///
+/// Uses the fast non-dominated sorting algorithm from
+/// Deb et al. (2002) with O(M N²) complexity.
+#[must_use]
+pub fn non_dominated_sort(solutions: &[Vec<f64>], directions: &[Direction]) -> Vec<Vec<usize>> {
+    fast_non_dominated_sort(solutions, directions)
+}
+
+/// Filter solutions to return only non-dominated (Pareto-optimal) indices.
+///
+/// Equivalent to `non_dominated_sort(solutions, directions)[0]` but
+/// communicates the intent more clearly.
+#[must_use]
+pub fn pareto_front_indices(solutions: &[Vec<f64>], directions: &[Direction]) -> Vec<usize> {
+    let fronts = fast_non_dominated_sort(solutions, directions);
+    fronts.into_iter().next().unwrap_or_default()
+}
+
+/// Compute crowding distance for diversity measurement.
+///
+/// Returns one distance value per solution in `front` (same order).
+/// Boundary solutions (best/worst in any objective) receive
+/// [`f64::INFINITY`]. Interior solutions get a finite positive value
+/// proportional to the gap between their neighbors.
+///
+/// `directions` is accepted for API consistency but does not affect
+/// the result, since crowding distance measures spacing regardless of
+/// optimization direction.
+#[must_use]
+#[allow(clippy::cast_precision_loss, clippy::needless_range_loop)]
+pub fn crowding_distance(front: &[Vec<f64>], _directions: &[Direction]) -> Vec<f64> {
+    let n = front.len();
+    if n <= 2 {
+        return vec![f64::INFINITY; n];
+    }
+
+    let m = front[0].len();
+    let mut distances = vec![0.0_f64; n];
+
+    for obj in 0..m {
+        let mut sorted: Vec<usize> = (0..n).collect();
+        sorted.sort_by(|&a, &b| {
+            front[a][obj]
+                .partial_cmp(&front[b][obj])
+                .unwrap_or(core::cmp::Ordering::Equal)
+        });
+
+        distances[sorted[0]] = f64::INFINITY;
+        distances[sorted[n - 1]] = f64::INFINITY;
+
+        let range = front[sorted[n - 1]][obj] - front[sorted[0]][obj];
+        if range > 0.0 {
+            for i in 1..(n - 1) {
+                distances[sorted[i]] +=
+                    (front[sorted[i + 1]][obj] - front[sorted[i - 1]][obj]) / range;
+            }
+        }
+    }
+
+    distances
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_dominates_basic() {
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        assert!(dominates(&[1.0, 1.0], &[2.0, 2.0], &dirs));
+        assert!(!dominates(&[2.0, 2.0], &[1.0, 1.0], &dirs));
+        // Equal does not dominate
+        assert!(!dominates(&[1.0, 1.0], &[1.0, 1.0], &dirs));
+    }
+
+    #[test]
+    fn test_dominates_incomparable() {
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        assert!(!dominates(&[1.0, 3.0], &[3.0, 1.0], &dirs));
+        assert!(!dominates(&[3.0, 1.0], &[1.0, 3.0], &dirs));
+    }
+
+    #[test]
+    fn test_dominates_maximize() {
+        let dirs = [Direction::Maximize, Direction::Minimize];
+        // a = (5, 1) vs b = (3, 2): a is better in both
+        assert!(dominates(&[5.0, 1.0], &[3.0, 2.0], &dirs));
+        assert!(!dominates(&[3.0, 2.0], &[5.0, 1.0], &dirs));
+    }
+
+    #[test]
+    fn test_nds_known() {
+        let values = vec![
+            vec![1.0, 5.0], // front 0
+            vec![5.0, 1.0], // front 0
+            vec![3.0, 3.0], // front 0 (non-dominated)
+            vec![4.0, 4.0], // front 1 (dominated by #2)
+            vec![6.0, 6.0], // front 2
+        ];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let fronts = fast_non_dominated_sort(&values, &dirs);
+
+        assert_eq!(fronts.len(), 3);
+        let mut f0 = fronts[0].clone();
+        f0.sort_unstable();
+        assert_eq!(f0, vec![0, 1, 2]);
+        assert_eq!(fronts[1], vec![3]);
+        assert_eq!(fronts[2], vec![4]);
+    }
+
+    #[test]
+    fn test_crowding_indexed_boundaries() {
+        let values = vec![vec![1.0, 5.0], vec![3.0, 3.0], vec![5.0, 1.0]];
+        let front = vec![0, 1, 2];
+        let cd = crowding_distance_indexed(&front, &values);
+        assert!(cd[0].is_infinite());
+        assert!(cd[2].is_infinite());
+        assert!(cd[1].is_finite());
+        assert!(cd[1] > 0.0);
+    }
+
+    // ---- Public API tests ----
+
+    #[test]
+    fn test_hypervolume_2d_minimize() {
+        // Front: (1,3), (2,2), (3,1) with ref (4,4) — all minimize
+        let front = vec![vec![1.0, 3.0], vec![2.0, 2.0], vec![3.0, 1.0]];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let hv = hypervolume(&front, &[4.0, 4.0], &dirs);
+        // Strip 1: x=[1,2), h=4-3=1 → area=1
+        // Strip 2: x=[2,3), h=4-2=2 → area=2
+        // Strip 3: x=[3,4], h=4-1=3 → area=3
+        // Total = 6
+        assert!((hv - 6.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_hypervolume_2d_maximize() {
+        // Front: (3,1), (2,2), (1,3) with ref (0,0) — all maximize
+        let front = vec![vec![3.0, 1.0], vec![2.0, 2.0], vec![1.0, 3.0]];
+        let dirs = [Direction::Maximize, Direction::Maximize];
+        let hv = hypervolume(&front, &[0.0, 0.0], &dirs);
+        // In negate-space: points become (-3,-1),(-2,-2),(-1,-3), ref=(0,0)
+        // Same geometry as minimize test above → area = 6
+        assert!((hv - 6.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_hypervolume_single_point() {
+        let front = vec![vec![1.0, 1.0]];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let hv = hypervolume(&front, &[3.0, 3.0], &dirs);
+        // Rectangle: (3-1) * (3-1) = 4
+        assert!((hv - 4.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_hypervolume_empty_front() {
+        let front: Vec<Vec<f64>> = vec![];
+        let dirs = [Direction::Minimize];
+        assert!(hypervolume(&front, &[1.0], &dirs).abs() < f64::EPSILON);
+    }
+
+    #[test]
+    fn test_hypervolume_point_at_ref() {
+        // Point not strictly better than ref → contributes nothing
+        let front = vec![vec![5.0, 5.0]];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let hv = hypervolume(&front, &[5.0, 5.0], &dirs);
+        assert!(hv.abs() < f64::EPSILON);
+    }
+
+    #[test]
+    fn test_hypervolume_3d() {
+        // Single point in 3D: (1,1,1) with ref (2,2,2)
+        let front = vec![vec![1.0, 1.0, 1.0]];
+        let dirs = [
+            Direction::Minimize,
+            Direction::Minimize,
+            Direction::Minimize,
+        ];
+        let hv = hypervolume(&front, &[2.0, 2.0, 2.0], &dirs);
+        assert!((hv - 1.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_non_dominated_sort_public() {
+        let values = vec![
+            vec![1.0, 5.0],
+            vec![5.0, 1.0],
+            vec![3.0, 3.0],
+            vec![4.0, 4.0],
+        ];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let fronts = non_dominated_sort(&values, &dirs);
+        assert_eq!(fronts.len(), 2);
+        let mut f0 = fronts[0].clone();
+        f0.sort_unstable();
+        assert_eq!(f0, vec![0, 1, 2]);
+        assert_eq!(fronts[1], vec![3]);
+    }
+
+    #[test]
+    fn test_pareto_front_indices_basic() {
+        let values = vec![
+            vec![1.0, 5.0],
+            vec![5.0, 1.0],
+            vec![3.0, 3.0],
+            vec![4.0, 4.0],
+        ];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let mut idx = pareto_front_indices(&values, &dirs);
+        idx.sort_unstable();
+        assert_eq!(idx, vec![0, 1, 2]);
+    }
+
+    #[test]
+    fn test_pareto_front_indices_empty() {
+        let values: Vec<Vec<f64>> = vec![];
+        let dirs = [Direction::Minimize];
+        assert!(pareto_front_indices(&values, &dirs).is_empty());
+    }
+
+    #[test]
+    fn test_crowding_distance_public() {
+        let front = vec![vec![1.0, 5.0], vec![3.0, 3.0], vec![5.0, 1.0]];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let cd = crowding_distance(&front, &dirs);
+        assert!(cd[0].is_infinite());
+        assert!(cd[2].is_infinite());
+        assert!(cd[1].is_finite());
+        assert!(cd[1] > 0.0);
+    }
+
+    #[test]
+    fn test_crowding_distance_single_point() {
+        let front = vec![vec![2.0, 3.0]];
+        let dirs = [Direction::Minimize, Direction::Minimize];
+        let cd = crowding_distance(&front, &dirs);
+        assert_eq!(cd.len(), 1);
+        assert!(cd[0].is_infinite());
+    }
+}

src/sampler/bohb.rs

@@ -0,0 +1,691 @@
+//! BOHB (Bayesian Optimization + `HyperBand`) sampler.
+//!
+//! BOHB combines TPE's model-guided sampling with Hyperband's budget-aware
+//! evaluation. Instead of building one global TPE model, BOHB conditions
+//! its TPE model on trials evaluated at a specific budget level, giving
+//! better-calibrated proposals for each rung of the Hyperband schedule.
+//!
+//! # How it works
+//!
+//! 1. Compute all Hyperband rung steps (budget levels) from the config.
+//! 2. On each `sample()` call, scan the history's `intermediate_values`
+//!    to find the **largest budget level** with enough observations
+//!    (`>= min_points_in_model`).
+//! 3. Build a filtered history where each trial's `value` is replaced
+//!    with its intermediate value at that budget level.
+//! 4. Delegate to an internal [`TpeSampler`] for the actual sampling.
+//! 5. Fall back to random sampling if no budget level has enough data.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::sampler::bohb::BohbSampler;
+//! use optimizer::{Direction, Study};
+//!
+//! let bohb = BohbSampler::new();
+//! let pruner = bohb.matching_pruner(Direction::Minimize);
+//! let study: Study<f64> = Study::with_sampler_and_pruner(Direction::Minimize, bohb, pruner);
+//! ```
+//!
+//! Using the builder for custom configuration:
+//!
+//! ```
+//! use optimizer::sampler::bohb::BohbSampler;
+//!
+//! let bohb = BohbSampler::builder()
+//!     .min_resource(1)
+//!     .max_resource(81)
+//!     .reduction_factor(3)
+//!     .min_points_in_model(10)
+//!     .seed(42)
+//!     .build()
+//!     .unwrap();
+//! ```
+
+use crate::distribution::Distribution;
+use crate::error::Result;
+use crate::param::ParamValue;
+use crate::pruner::HyperbandPruner;
+use crate::sampler::tpe::TpeSampler;
+use crate::sampler::{CompletedTrial, Sampler};
+use crate::types::Direction;
+
+/// A BOHB sampler that combines TPE with Hyperband budget awareness.
+///
+/// BOHB filters trial history by budget level before delegating to TPE,
+/// so the surrogate model is conditioned on trials evaluated at the same
+/// resource level. This produces better-calibrated parameter proposals
+/// than using a single global model across all budgets.
+///
+/// Use [`BohbSampler::matching_pruner`] to create a [`HyperbandPruner`]
+/// with matching parameters.
+pub struct BohbSampler {
+    min_resource: u64,
+    max_resource: u64,
+    reduction_factor: u64,
+    min_points_in_model: usize,
+    tpe: TpeSampler,
+}
+
+impl BohbSampler {
+    /// Creates a new BOHB sampler with default settings.
+    ///
+    /// Defaults:
+    /// - `min_resource`: 1
+    /// - `max_resource`: 81
+    /// - `reduction_factor`: 3
+    /// - `min_points_in_model`: 10
+    /// - TPE: default settings
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            min_resource: 1,
+            max_resource: 81,
+            reduction_factor: 3,
+            min_points_in_model: 10,
+            tpe: TpeSampler::new(),
+        }
+    }
+
+    /// Creates a builder for configuring a BOHB sampler.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use optimizer::sampler::bohb::BohbSampler;
+    ///
+    /// let sampler = BohbSampler::builder()
+    ///     .min_resource(1)
+    ///     .max_resource(27)
+    ///     .reduction_factor(3)
+    ///     .min_points_in_model(5)
+    ///     .seed(42)
+    ///     .build()
+    ///     .unwrap();
+    /// ```
+    #[must_use]
+    pub fn builder() -> BohbSamplerBuilder {
+        BohbSamplerBuilder::new()
+    }
+
+    /// Creates a [`HyperbandPruner`] with matching Hyperband parameters.
+    ///
+    /// This ensures the pruner's budget schedule is consistent with the
+    /// budget levels used by BOHB for model conditioning.
+    #[must_use]
+    pub fn matching_pruner(&self, direction: Direction) -> HyperbandPruner {
+        HyperbandPruner::new()
+            .min_resource(self.min_resource)
+            .max_resource(self.max_resource)
+            .reduction_factor(self.reduction_factor)
+            .direction(direction)
+    }
+
+    /// Compute all unique budget levels (rung steps) across all Hyperband brackets.
+    ///
+    /// Returns sorted ascending.
+    #[allow(
+        clippy::cast_precision_loss,
+        clippy::cast_possible_truncation,
+        clippy::cast_sign_loss
+    )]
+    fn all_budget_levels(&self) -> Vec<u64> {
+        let eta = self.reduction_factor as f64;
+        let ratio = self.max_resource as f64 / self.min_resource as f64;
+        let s_max = (ratio.ln() / eta.ln()).floor() as u64;
+
+        let mut levels = Vec::new();
+        for bracket in 0..=s_max {
+            let exponent = s_max.saturating_sub(bracket);
+            let min_resource_bracket =
+                (self.max_resource as f64 / eta.powi(exponent as i32)).ceil() as u64;
+
+            let mut rung: u32 = 0;
+            while let Some(power) = self.reduction_factor.checked_pow(rung) {
+                let step = min_resource_bracket.saturating_mul(power);
+                if step > self.max_resource {
+                    break;
+                }
+                levels.push(step);
+                rung += 1;
+            }
+        }
+
+        levels.sort_unstable();
+        levels.dedup();
+        levels
+    }
+
+    /// Build a filtered history for a specific budget level.
+    ///
+    /// For each trial that has an intermediate value at the given budget step,
+    /// creates a new `CompletedTrial` with `value` replaced by the intermediate
+    /// value at that step.
+    fn filter_history_for_budget(history: &[CompletedTrial], budget: u64) -> Vec<CompletedTrial> {
+        history
+            .iter()
+            .filter_map(|trial| {
+                trial
+                    .intermediate_values
+                    .iter()
+                    .find(|(step, _)| *step == budget)
+                    .map(|(_, iv)| CompletedTrial {
+                        id: trial.id,
+                        params: trial.params.clone(),
+                        distributions: trial.distributions.clone(),
+                        param_labels: trial.param_labels.clone(),
+                        value: *iv,
+                        intermediate_values: trial.intermediate_values.clone(),
+                        state: trial.state,
+                        user_attrs: trial.user_attrs.clone(),
+                        constraints: trial.constraints.clone(),
+                    })
+            })
+            .collect()
+    }
+}
+
+impl Default for BohbSampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Sampler for BohbSampler {
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[CompletedTrial],
+    ) -> ParamValue {
+        // Find the largest budget level with enough observations
+        let levels = self.all_budget_levels();
+
+        for &budget in levels.iter().rev() {
+            let count = history
+                .iter()
+                .filter(|t| {
+                    t.intermediate_values
+                        .iter()
+                        .any(|(step, _)| *step == budget)
+                })
+                .count();
+
+            if count >= self.min_points_in_model {
+                let filtered = Self::filter_history_for_budget(history, budget);
+                return self.tpe.sample(distribution, trial_id, &filtered);
+            }
+        }
+
+        // Not enough data at any budget level: delegate to TPE with empty history
+        // which triggers its uniform-random startup behavior.
+        self.tpe.sample(distribution, trial_id, &[])
+    }
+}
+
+/// Builder for configuring a [`BohbSampler`].
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::bohb::BohbSamplerBuilder;
+///
+/// let sampler = BohbSamplerBuilder::new()
+///     .min_resource(1)
+///     .max_resource(81)
+///     .reduction_factor(3)
+///     .gamma(0.15)
+///     .seed(42)
+///     .build()
+///     .unwrap();
+/// ```
+pub struct BohbSamplerBuilder {
+    min_resource: u64,
+    max_resource: u64,
+    reduction_factor: u64,
+    min_points_in_model: usize,
+    tpe_builder: crate::sampler::tpe::TpeSamplerBuilder,
+}
+
+impl BohbSamplerBuilder {
+    /// Creates a new builder with default settings.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            min_resource: 1,
+            max_resource: 81,
+            reduction_factor: 3,
+            min_points_in_model: 10,
+            tpe_builder: crate::sampler::tpe::TpeSamplerBuilder::new(),
+        }
+    }
+
+    /// Sets the minimum resource (budget) per trial.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `r` is 0.
+    #[must_use]
+    pub fn min_resource(mut self, r: u64) -> Self {
+        assert!(r > 0, "min_resource must be > 0, got {r}");
+        self.min_resource = r;
+        self
+    }
+
+    /// Sets the maximum resource (budget) per trial.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `r` is 0.
+    #[must_use]
+    pub fn max_resource(mut self, r: u64) -> Self {
+        assert!(r > 0, "max_resource must be > 0, got {r}");
+        self.max_resource = r;
+        self
+    }
+
+    /// Sets the reduction factor (eta).
+    ///
+    /// # Panics
+    ///
+    /// Panics if `eta` is less than 2.
+    #[must_use]
+    pub fn reduction_factor(mut self, eta: u64) -> Self {
+        assert!(eta >= 2, "reduction_factor must be >= 2, got {eta}");
+        self.reduction_factor = eta;
+        self
+    }
+
+    /// Sets the minimum number of observations at a budget level before
+    /// BOHB uses TPE instead of random sampling.
+    #[must_use]
+    pub fn min_points_in_model(mut self, n: usize) -> Self {
+        self.min_points_in_model = n;
+        self
+    }
+
+    /// Sets a fixed gamma value for the internal TPE sampler.
+    #[must_use]
+    pub fn gamma(mut self, gamma: f64) -> Self {
+        self.tpe_builder = self.tpe_builder.gamma(gamma);
+        self
+    }
+
+    /// Sets a custom gamma strategy for the internal TPE sampler.
+    #[must_use]
+    pub fn gamma_strategy<G: crate::sampler::tpe::GammaStrategy + 'static>(
+        mut self,
+        strategy: G,
+    ) -> Self {
+        self.tpe_builder = self.tpe_builder.gamma_strategy(strategy);
+        self
+    }
+
+    /// Sets the number of EI candidates for the internal TPE sampler.
+    #[must_use]
+    pub fn n_ei_candidates(mut self, n: usize) -> Self {
+        self.tpe_builder = self.tpe_builder.n_ei_candidates(n);
+        self
+    }
+
+    /// Sets a fixed KDE bandwidth for the internal TPE sampler.
+    #[must_use]
+    pub fn kde_bandwidth(mut self, bandwidth: f64) -> Self {
+        self.tpe_builder = self.tpe_builder.kde_bandwidth(bandwidth);
+        self
+    }
+
+    /// Sets a seed for reproducible sampling.
+    #[must_use]
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.tpe_builder = self.tpe_builder.seed(seed);
+        self
+    }
+
+    /// Builds the configured [`BohbSampler`].
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if the TPE configuration is invalid (e.g. gamma
+    /// not in (0, 1) or bandwidth not positive).
+    pub fn build(self) -> Result<BohbSampler> {
+        let tpe = self.tpe_builder.build()?;
+        Ok(BohbSampler {
+            min_resource: self.min_resource,
+            max_resource: self.max_resource,
+            reduction_factor: self.reduction_factor,
+            min_points_in_model: self.min_points_in_model,
+            tpe,
+        })
+    }
+}
+
+impl Default for BohbSamplerBuilder {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+#[cfg(test)]
+#[allow(clippy::cast_precision_loss)]
+mod tests {
+    use std::collections::HashMap;
+
+    use super::*;
+    use crate::distribution::{FloatDistribution, IntDistribution};
+    use crate::parameter::ParamId;
+    use crate::types::TrialState;
+
+    fn make_trial_with_intermediates(
+        id: u64,
+        value: f64,
+        params: Vec<(ParamId, ParamValue, Distribution)>,
+        intermediate_values: Vec<(u64, f64)>,
+    ) -> CompletedTrial {
+        let mut param_map = HashMap::new();
+        let mut dist_map = HashMap::new();
+        for (param_id, pv, dist) in params {
+            param_map.insert(param_id, pv);
+            dist_map.insert(param_id, dist);
+        }
+        CompletedTrial {
+            id,
+            params: param_map,
+            distributions: dist_map,
+            param_labels: HashMap::new(),
+            value,
+            intermediate_values,
+            state: TrialState::Complete,
+            user_attrs: HashMap::new(),
+            constraints: Vec::new(),
+        }
+    }
+
+    #[test]
+    fn budget_levels_default() {
+        let bohb = BohbSampler::new();
+        let levels = bohb.all_budget_levels();
+        // With min=1, max=81, eta=3:
+        // bracket 0: 1, 3, 9, 27, 81
+        // bracket 1: 3, 9, 27, 81
+        // bracket 2: 9, 27, 81
+        // bracket 3: 27, 81
+        // bracket 4: 81
+        // Unique sorted: [1, 3, 9, 27, 81]
+        assert_eq!(levels, vec![1, 3, 9, 27, 81]);
+    }
+
+    #[test]
+    fn budget_levels_eta2() {
+        let bohb = BohbSampler::builder()
+            .min_resource(1)
+            .max_resource(16)
+            .reduction_factor(2)
+            .build()
+            .unwrap();
+        let levels = bohb.all_budget_levels();
+        // s_max = floor(ln(16)/ln(2)) = 4
+        // bracket 0: 1, 2, 4, 8, 16
+        // bracket 1: 2, 4, 8, 16
+        // bracket 2: 4, 8, 16
+        // bracket 3: 8, 16
+        // bracket 4: 16
+        // Unique sorted: [1, 2, 4, 8, 16]
+        assert_eq!(levels, vec![1, 2, 4, 8, 16]);
+    }
+
+    #[test]
+    fn filter_history_selects_correct_budget() {
+        let x_id = ParamId::new();
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        let history = vec![
+            make_trial_with_intermediates(
+                0,
+                0.5,
+                vec![(x_id, ParamValue::Float(0.3), dist.clone())],
+                vec![(1, 0.9), (3, 0.7), (9, 0.5)],
+            ),
+            make_trial_with_intermediates(
+                1,
+                0.4,
+                vec![(x_id, ParamValue::Float(0.6), dist.clone())],
+                vec![(1, 0.8), (3, 0.4)],
+            ),
+            make_trial_with_intermediates(
+                2,
+                0.3,
+                vec![(x_id, ParamValue::Float(0.1), dist.clone())],
+                vec![(1, 0.7)],
+            ),
+        ];
+
+        // Budget 3: trials 0 and 1 have intermediate values at step 3
+        let filtered = BohbSampler::filter_history_for_budget(&history, 3);
+        assert_eq!(filtered.len(), 2);
+        assert!((filtered[0].value - 0.7).abs() < f64::EPSILON);
+        assert!((filtered[1].value - 0.4).abs() < f64::EPSILON);
+
+        // Budget 9: only trial 0
+        let filtered = BohbSampler::filter_history_for_budget(&history, 9);
+        assert_eq!(filtered.len(), 1);
+        assert!((filtered[0].value - 0.5).abs() < f64::EPSILON);
+
+        // Budget 27: nobody
+        let filtered = BohbSampler::filter_history_for_budget(&history, 27);
+        assert!(filtered.is_empty());
+    }
+
+    #[test]
+    fn matching_pruner_has_same_params() {
+        let bohb = BohbSampler::builder()
+            .min_resource(2)
+            .max_resource(64)
+            .reduction_factor(4)
+            .build()
+            .unwrap();
+        let pruner = bohb.matching_pruner(Direction::Minimize);
+
+        // We can't directly inspect HyperbandPruner fields, but we can
+        // verify it was created without panicking with the same params.
+        // The pruner's rung steps should match BOHB's budget levels.
+        // Just verify it doesn't panic.
+        drop(pruner);
+    }
+
+    #[test]
+    fn fallback_to_random_when_insufficient_data() {
+        let bohb = BohbSampler::builder()
+            .min_points_in_model(10)
+            .seed(42)
+            .build()
+            .unwrap();
+
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        // Only 3 trials with intermediate values (< min_points_in_model=10)
+        let x_id = ParamId::new();
+        let history: Vec<CompletedTrial> = (0..3)
+            .map(|i| {
+                make_trial_with_intermediates(
+                    i,
+                    i as f64,
+                    vec![(x_id, ParamValue::Float(i as f64 / 3.0), dist.clone())],
+                    vec![(1, i as f64)],
+                )
+            })
+            .collect();
+
+        // Should not panic, should sample within bounds
+        for trial_id in 0..20 {
+            let val = bohb.sample(&dist, trial_id, &history);
+            if let ParamValue::Float(v) = val {
+                assert!((0.0..=1.0).contains(&v));
+            } else {
+                panic!("Expected Float");
+            }
+        }
+    }
+
+    #[test]
+    fn uses_budget_level_when_enough_data() {
+        let bohb = BohbSampler::builder()
+            .min_points_in_model(5)
+            .seed(42)
+            .build()
+            .unwrap();
+
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 10.0,
+            log_scale: false,
+            step: None,
+        });
+
+        // Create 20 trials with intermediate values at budget 1.
+        // Good trials have x near 2.0, bad trials have x far from 2.0.
+        let x_id = ParamId::new();
+        let history: Vec<CompletedTrial> = (0..20)
+            .map(|i| {
+                let x = i as f64 / 2.0;
+                let iv_at_1 = (x - 2.0).powi(2);
+                make_trial_with_intermediates(
+                    i,
+                    iv_at_1, // final value same as intermediate for simplicity
+                    vec![(x_id, ParamValue::Float(x), dist.clone())],
+                    vec![(1, iv_at_1)],
+                )
+            })
+            .collect();
+
+        // Should use TPE on filtered history at budget 1
+        let val = bohb.sample(&dist, 100, &history);
+        if let ParamValue::Float(v) = val {
+            assert!((0.0..=10.0).contains(&v), "Value {v} out of bounds");
+        } else {
+            panic!("Expected Float");
+        }
+    }
+
+    #[test]
+    fn prefers_largest_budget_level() {
+        let bohb = BohbSampler::builder()
+            .min_resource(1)
+            .max_resource(9)
+            .reduction_factor(3)
+            .min_points_in_model(3)
+            .seed(42)
+            .build()
+            .unwrap();
+
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 10.0,
+            log_scale: false,
+            step: None,
+        });
+
+        // Budget levels: [1, 3, 9]
+        assert_eq!(bohb.all_budget_levels(), vec![1, 3, 9]);
+
+        // Create 5 trials with intermediates at budget 1 and 3
+        let x_id = ParamId::new();
+        let history: Vec<CompletedTrial> = (0..5)
+            .map(|i| {
+                let x = i as f64;
+                make_trial_with_intermediates(
+                    i,
+                    x,
+                    vec![(x_id, ParamValue::Float(x), dist.clone())],
+                    vec![(1, x * 2.0), (3, x)],
+                )
+            })
+            .collect();
+
+        // Budget 3 has 5 observations (>= 3), budget 9 has 0.
+        // BOHB should pick budget 3 (largest with enough data).
+        // The filtered history at budget 3 has values [0, 1, 2, 3, 4].
+        let filtered_3 = BohbSampler::filter_history_for_budget(&history, 3);
+        assert_eq!(filtered_3.len(), 5);
+        let filtered_9 = BohbSampler::filter_history_for_budget(&history, 9);
+        assert_eq!(filtered_9.len(), 0);
+
+        // Should sample successfully
+        let val = bohb.sample(&dist, 100, &history);
+        assert!(matches!(val, ParamValue::Float(_)));
+    }
+
+    #[test]
+    fn builder_validates_tpe_params() {
+        // Invalid gamma
+        let result = BohbSampler::builder().gamma(1.5).build();
+        assert!(result.is_err());
+
+        // Invalid bandwidth
+        let result = BohbSampler::builder().kde_bandwidth(-1.0).build();
+        assert!(result.is_err());
+    }
+
+    #[test]
+    #[should_panic(expected = "min_resource must be > 0")]
+    fn builder_rejects_zero_min_resource() {
+        let _ = BohbSampler::builder().min_resource(0);
+    }
+
+    #[test]
+    #[should_panic(expected = "max_resource must be > 0")]
+    fn builder_rejects_zero_max_resource() {
+        let _ = BohbSampler::builder().max_resource(0);
+    }
+
+    #[test]
+    #[should_panic(expected = "reduction_factor must be >= 2")]
+    fn builder_rejects_small_reduction_factor() {
+        let _ = BohbSampler::builder().reduction_factor(1);
+    }
+
+    #[test]
+    fn int_distribution_works() {
+        let bohb = BohbSampler::builder()
+            .min_points_in_model(3)
+            .seed(42)
+            .build()
+            .unwrap();
+
+        let dist = Distribution::Int(IntDistribution {
+            low: 0,
+            high: 100,
+            log_scale: false,
+            step: None,
+        });
+
+        let x_id = ParamId::new();
+        let history: Vec<CompletedTrial> = (0..10)
+            .map(|i| {
+                make_trial_with_intermediates(
+                    i,
+                    i as f64,
+                    vec![(x_id, ParamValue::Int(i.cast_signed() * 10), dist.clone())],
+                    vec![(1, i as f64)],
+                )
+            })
+            .collect();
+
+        let val = bohb.sample(&dist, 100, &history);
+        if let ParamValue::Int(v) = val {
+            assert!((0..=100).contains(&v));
+        } else {
+            panic!("Expected Int");
+        }
+    }
+}

src/sampler/mod.rs

@@ -1,7 +1,14 @@
 //! Sampler trait and implementations for parameter sampling.
 
+pub mod bohb;
+#[cfg(feature = "cma-es")]
+pub mod cma_es;
 pub mod grid;
+pub mod motpe;
+pub mod nsga2;
 pub mod random;
+#[cfg(feature = "sobol")]
+pub mod sobol;
 pub mod tpe;
 
 use std::collections::HashMap;
@@ -18,6 +25,7 @@ use crate::types::TrialState;
 /// parameter values, their distributions, and the objective value returned
 /// by the objective function.
 #[derive(Clone, Debug)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub struct CompletedTrial<V = f64> {
     /// The unique identifier for this trial.
     pub id: u64,
@@ -35,6 +43,9 @@ pub struct CompletedTrial<V = f64> {
     pub state: TrialState,
     /// User-defined attributes stored during the trial.
     pub user_attrs: HashMap<String, AttrValue>,
+    /// Constraint values for this trial (<=0.0 means feasible).
+    #[cfg_attr(feature = "serde", serde(default))]
+    pub constraints: Vec<f64>,
 }
 
 impl<V> CompletedTrial<V> {
@@ -55,6 +66,7 @@ impl<V> CompletedTrial<V> {
             intermediate_values: Vec::new(),
             state: TrialState::Complete,
             user_attrs: HashMap::new(),
+            constraints: Vec::new(),
         }
     }
 
@@ -77,6 +89,7 @@ impl<V> CompletedTrial<V> {
             intermediate_values,
             state: TrialState::Complete,
             user_attrs,
+            constraints: Vec::new(),
         }
     }
 
@@ -121,6 +134,14 @@ impl<V> CompletedTrial<V> {
         })
     }
 
+    /// Returns `true` if all constraints are satisfied (values <= 0.0).
+    ///
+    /// A trial with no constraints is considered feasible.
+    #[must_use]
+    pub fn is_feasible(&self) -> bool {
+        self.constraints.iter().all(|&c| c <= 0.0)
+    }
+
     /// Gets a user attribute by key.
     #[must_use]
     pub fn user_attr(&self, key: &str) -> Option<&AttrValue> {

src/sampler/motpe.rs

@@ -0,0 +1,932 @@
+//! Multi-Objective Tree-Parzen Estimator (MOTPE) sampler.
+//!
+//! Extends TPE to handle multi-objective optimization by using Pareto
+//! non-dominated sorting to define "good" vs "bad" trial regions for
+//! the KDE models, replacing the single-objective gamma-based split.
+//!
+//! # Algorithm
+//!
+//! In single-objective TPE, trials are sorted by value and split at a
+//! gamma percentile into good/bad groups. MOTPE replaces this with:
+//!
+//! 1. Compute non-dominated sorting on all completed trials
+//! 2. Use the Pareto front (rank 0) as "good" trials
+//! 3. Use dominated trials as "bad" trials
+//! 4. Build KDE l(x) from good, g(x) from bad
+//! 5. Sample candidates and score by l(x)/g(x)
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::Direction;
+//! use optimizer::multi_objective::MultiObjectiveStudy;
+//! use optimizer::parameter::{FloatParam, Parameter};
+//! use optimizer::sampler::motpe::MotpeSampler;
+//!
+//! let sampler = MotpeSampler::builder().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();
+//!
+//! let front = study.pareto_front();
+//! assert!(!front.is_empty());
+//! ```
+
+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};
+
+/// Multi-Objective TPE (MOTPE) sampler for multi-objective Bayesian optimization.
+///
+/// Uses Pareto non-dominated sorting to split completed trials into
+/// "good" (non-dominated, rank 0) and "bad" (dominated) groups, then
+/// fits kernel density estimators to each group and samples new points
+/// that maximize l(x)/g(x).
+///
+/// During the startup phase (fewer than `n_startup_trials` completed),
+/// MOTPE falls back to random sampling.
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::Direction;
+/// use optimizer::multi_objective::MultiObjectiveStudy;
+/// use optimizer::parameter::{FloatParam, Parameter};
+/// use optimizer::sampler::motpe::MotpeSampler;
+///
+/// let sampler = MotpeSampler::builder()
+///     .n_startup_trials(10)
+///     .n_ei_candidates(24)
+///     .seed(42)
+///     .build();
+///
+/// let study =
+///     MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+/// ```
+pub struct MotpeSampler {
+    /// Number of trials before MOTPE kicks in (uses random sampling before this).
+    n_startup_trials: usize,
+    /// Number of candidate samples to evaluate when selecting the next point.
+    n_ei_candidates: usize,
+    /// Optional fixed bandwidth for KDE. If None, uses Scott's rule.
+    kde_bandwidth: Option<f64>,
+    /// Thread-safe RNG for sampling.
+    rng: Mutex<StdRng>,
+}
+
+impl MotpeSampler {
+    /// Creates a new MOTPE sampler with default settings.
+    ///
+    /// Defaults:
+    /// - `n_startup_trials`: 11
+    /// - `n_ei_candidates`: 24
+    /// - `kde_bandwidth`: None (Scott's rule)
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            n_startup_trials: 11,
+            n_ei_candidates: 24,
+            kde_bandwidth: None,
+            rng: Mutex::new(rand::make_rng()),
+        }
+    }
+
+    /// Creates a new MOTPE sampler with a fixed seed.
+    #[must_use]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            n_startup_trials: 11,
+            n_ei_candidates: 24,
+            kde_bandwidth: None,
+            rng: Mutex::new(StdRng::seed_from_u64(seed)),
+        }
+    }
+
+    /// Creates a builder for configuring a MOTPE sampler.
+    #[must_use]
+    pub fn builder() -> MotpeSamplerBuilder {
+        MotpeSamplerBuilder::new()
+    }
+
+    /// Splits trials into good (non-dominated) and bad (dominated) groups
+    /// using Pareto non-dominated sorting.
+    fn split_trials<'a>(
+        history: &'a [MultiObjectiveTrial],
+        directions: &[Direction],
+    ) -> (Vec<&'a MultiObjectiveTrial>, Vec<&'a MultiObjectiveTrial>) {
+        let complete: Vec<(usize, &MultiObjectiveTrial)> = history
+            .iter()
+            .enumerate()
+            .filter(|(_, t)| t.state == TrialState::Complete)
+            .collect();
+
+        if complete.is_empty() {
+            return (vec![], vec![]);
+        }
+
+        let values: Vec<Vec<f64>> = complete.iter().map(|(_, t)| t.values.clone()).collect();
+        let constraints: Vec<Vec<f64>> = complete
+            .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(&values, directions, &constraints)
+        } else {
+            pareto::fast_non_dominated_sort(&values, directions)
+        };
+
+        if fronts.is_empty() {
+            return (vec![], vec![]);
+        }
+
+        // Front 0 = good (non-dominated), everything else = bad
+        let good: Vec<&MultiObjectiveTrial> = fronts[0].iter().map(|&i| complete[i].1).collect();
+
+        let bad: Vec<&MultiObjectiveTrial> = fronts[1..]
+            .iter()
+            .flatten()
+            .map(|&i| complete[i].1)
+            .collect();
+
+        (good, bad)
+    }
+
+    /// Samples uniformly from a distribution (used during startup phase).
+    #[allow(
+        clippy::cast_possible_truncation,
+        clippy::cast_precision_loss,
+        clippy::unused_self
+    )]
+    fn sample_uniform(distribution: &Distribution, rng: &mut StdRng) -> 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))
+            }
+        }
+    }
+
+    /// Samples using TPE for float distributions.
+    #[allow(clippy::too_many_arguments)]
+    fn sample_tpe_float(
+        &self,
+        low: f64,
+        high: f64,
+        log_scale: bool,
+        step: Option<f64>,
+        good_values: Vec<f64>,
+        bad_values: Vec<f64>,
+        rng: &mut StdRng,
+    ) -> f64 {
+        // Transform to internal space (log space if needed)
+        let (internal_low, internal_high, good_internal, bad_internal) = if log_scale {
+            let i_low = low.ln();
+            let i_high = high.ln();
+            let g: Vec<f64> = good_values.iter().map(|&v| v.ln()).collect();
+            let b: Vec<f64> = bad_values.iter().map(|&v| v.ln()).collect();
+            (i_low, i_high, g, b)
+        } else {
+            (low, high, good_values, bad_values)
+        };
+
+        // Fit KDEs to good and bad groups
+        let l_kde = match self.kde_bandwidth {
+            Some(bw) => KernelDensityEstimator::with_bandwidth(good_internal, bw),
+            None => KernelDensityEstimator::new(good_internal),
+        };
+        let g_kde = match self.kde_bandwidth {
+            Some(bw) => KernelDensityEstimator::with_bandwidth(bad_internal, bw),
+            None => KernelDensityEstimator::new(bad_internal),
+        };
+
+        // 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);
+        };
+
+        // Generate candidates from l(x) and select the one with best l(x)/g(x)
+        let mut best_candidate = internal_low;
+        let mut best_ratio = f64::NEG_INFINITY;
+
+        for _ in 0..self.n_ei_candidates {
+            let candidate = l_kde.sample(rng).clamp(internal_low, internal_high);
+
+            let l_density = l_kde.pdf(candidate);
+            let g_density = g_kde.pdf(candidate);
+
+            let ratio = if g_density < f64::EPSILON {
+                if l_density > f64::EPSILON {
+                    f64::INFINITY
+                } else {
+                    0.0
+                }
+            } else {
+                l_density / g_density
+            };
+
+            if ratio > best_ratio {
+                best_ratio = ratio;
+                best_candidate = candidate;
+            }
+        }
+
+        // Transform back from internal space
+        let mut value = if log_scale {
+            best_candidate.exp()
+        } else {
+            best_candidate
+        };
+
+        // Apply step constraint if present
+        if let Some(step) = step {
+            let k = ((value - low) / step).round();
+            value = low + k * step;
+        }
+
+        value.clamp(low, high)
+    }
+
+    /// Samples using TPE for integer distributions.
+    #[allow(
+        clippy::too_many_arguments,
+        clippy::cast_precision_loss,
+        clippy::cast_possible_truncation
+    )]
+    fn sample_tpe_int(
+        &self,
+        low: i64,
+        high: i64,
+        log_scale: bool,
+        step: Option<i64>,
+        good_values: &[i64],
+        bad_values: &[i64],
+        rng: &mut StdRng,
+    ) -> 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();
+
+        let float_value = self.sample_tpe_float(
+            low as f64,
+            high as f64,
+            log_scale,
+            step.map(|s| s as f64),
+            good_floats,
+            bad_floats,
+            rng,
+        );
+
+        let int_value = float_value.round() as i64;
+        let int_value = if let Some(step) = step {
+            let k = ((int_value - low) as f64 / step as f64).round() as i64;
+            low + k * step
+        } else {
+            int_value
+        };
+
+        int_value.clamp(low, high)
+    }
+
+    /// Samples using TPE for categorical distributions.
+    #[allow(clippy::cast_precision_loss)]
+    fn sample_tpe_categorical(
+        n_choices: usize,
+        good_indices: &[usize],
+        bad_indices: &[usize],
+        rng: &mut StdRng,
+    ) -> usize {
+        let mut good_counts = vec![0usize; n_choices];
+        let mut bad_counts = vec![0usize; n_choices];
+
+        for &idx in good_indices {
+            if idx < n_choices {
+                good_counts[idx] += 1;
+            }
+        }
+        for &idx in bad_indices {
+            if idx < n_choices {
+                bad_counts[idx] += 1;
+            }
+        }
+
+        // Laplace smoothing
+        let good_total = good_indices.len() as f64 + n_choices as f64;
+        let bad_total = bad_indices.len() as f64 + n_choices as f64;
+
+        let mut weights = vec![0.0f64; n_choices];
+        for i in 0..n_choices {
+            let l_prob = (good_counts[i] as f64 + 1.0) / good_total;
+            let g_prob = (bad_counts[i] as f64 + 1.0) / bad_total;
+            weights[i] = l_prob / g_prob;
+        }
+
+        // Sample proportionally to weights
+        let total_weight: f64 = weights.iter().sum();
+        let threshold = rng.random::<f64>() * total_weight;
+
+        let mut cumulative = 0.0;
+        for (i, &w) in weights.iter().enumerate() {
+            cumulative += w;
+            if cumulative >= threshold {
+                return i;
+            }
+        }
+
+        n_choices - 1
+    }
+}
+
+impl Default for MotpeSampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl MultiObjectiveSampler for MotpeSampler {
+    #[allow(clippy::too_many_lines)]
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        _trial_id: u64,
+        history: &[MultiObjectiveTrial],
+        directions: &[Direction],
+    ) -> ParamValue {
+        let mut rng = self.rng.lock();
+
+        // Fall back to random sampling during startup phase
+        let n_complete = history
+            .iter()
+            .filter(|t| t.state == TrialState::Complete)
+            .count();
+        if n_complete < self.n_startup_trials {
+            return Self::sample_uniform(distribution, &mut rng);
+        }
+
+        // Split trials into good (Pareto front) and bad (dominated)
+        let (good_trials, bad_trials) = Self::split_trials(history, directions);
+
+        if good_trials.is_empty() || bad_trials.is_empty() {
+            return Self::sample_uniform(distribution, &mut rng);
+        }
+
+        match distribution {
+            Distribution::Float(d) => {
+                let good_values: Vec<f64> = good_trials
+                    .iter()
+                    .flat_map(|t| t.params.values())
+                    .filter_map(|v| match v {
+                        ParamValue::Float(f) => Some(*f),
+                        _ => None,
+                    })
+                    .filter(|&v| v >= d.low && v <= d.high)
+                    .collect();
+
+                let bad_values: Vec<f64> = bad_trials
+                    .iter()
+                    .flat_map(|t| t.params.values())
+                    .filter_map(|v| match v {
+                        ParamValue::Float(f) => Some(*f),
+                        _ => None,
+                    })
+                    .filter(|&v| v >= d.low && v <= d.high)
+                    .collect();
+
+                if good_values.is_empty() || bad_values.is_empty() {
+                    return Self::sample_uniform(distribution, &mut rng);
+                }
+
+                let value = self.sample_tpe_float(
+                    d.low,
+                    d.high,
+                    d.log_scale,
+                    d.step,
+                    good_values,
+                    bad_values,
+                    &mut rng,
+                );
+                ParamValue::Float(value)
+            }
+            Distribution::Int(d) => {
+                let good_values: Vec<i64> = good_trials
+                    .iter()
+                    .flat_map(|t| t.params.values())
+                    .filter_map(|v| match v {
+                        ParamValue::Int(i) => Some(*i),
+                        _ => None,
+                    })
+                    .filter(|&v| v >= d.low && v <= d.high)
+                    .collect();
+
+                let bad_values: Vec<i64> = bad_trials
+                    .iter()
+                    .flat_map(|t| t.params.values())
+                    .filter_map(|v| match v {
+                        ParamValue::Int(i) => Some(*i),
+                        _ => None,
+                    })
+                    .filter(|&v| v >= d.low && v <= d.high)
+                    .collect();
+
+                if good_values.is_empty() || bad_values.is_empty() {
+                    return Self::sample_uniform(distribution, &mut rng);
+                }
+
+                let value = self.sample_tpe_int(
+                    d.low,
+                    d.high,
+                    d.log_scale,
+                    d.step,
+                    &good_values,
+                    &bad_values,
+                    &mut rng,
+                );
+                ParamValue::Int(value)
+            }
+            Distribution::Categorical(d) => {
+                let good_indices: Vec<usize> = good_trials
+                    .iter()
+                    .flat_map(|t| t.params.values())
+                    .filter_map(|v| match v {
+                        ParamValue::Categorical(i) => Some(*i),
+                        _ => None,
+                    })
+                    .filter(|&i| i < d.n_choices)
+                    .collect();
+
+                let bad_indices: Vec<usize> = bad_trials
+                    .iter()
+                    .flat_map(|t| t.params.values())
+                    .filter_map(|v| match v {
+                        ParamValue::Categorical(i) => Some(*i),
+                        _ => None,
+                    })
+                    .filter(|&i| i < d.n_choices)
+                    .collect();
+
+                if good_indices.is_empty() || bad_indices.is_empty() {
+                    return Self::sample_uniform(distribution, &mut rng);
+                }
+
+                let index = Self::sample_tpe_categorical(
+                    d.n_choices,
+                    &good_indices,
+                    &bad_indices,
+                    &mut rng,
+                );
+                ParamValue::Categorical(index)
+            }
+        }
+    }
+}
+
+/// Builder for configuring a [`MotpeSampler`].
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::motpe::MotpeSamplerBuilder;
+///
+/// let sampler = MotpeSamplerBuilder::new()
+///     .n_startup_trials(15)
+///     .n_ei_candidates(32)
+///     .seed(42)
+///     .build();
+/// ```
+#[derive(Debug, Clone)]
+pub struct MotpeSamplerBuilder {
+    n_startup_trials: usize,
+    n_ei_candidates: usize,
+    kde_bandwidth: Option<f64>,
+    seed: Option<u64>,
+}
+
+impl MotpeSamplerBuilder {
+    /// Creates a new builder with default settings.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            n_startup_trials: 11,
+            n_ei_candidates: 24,
+            kde_bandwidth: None,
+            seed: None,
+        }
+    }
+
+    /// Sets the number of startup trials before MOTPE sampling begins.
+    #[must_use]
+    pub fn n_startup_trials(mut self, n: usize) -> Self {
+        self.n_startup_trials = n;
+        self
+    }
+
+    /// Sets the number of EI candidates to evaluate per sample.
+    #[must_use]
+    pub fn n_ei_candidates(mut self, n: usize) -> Self {
+        self.n_ei_candidates = n;
+        self
+    }
+
+    /// Sets a fixed bandwidth for the kernel density estimator.
+    ///
+    /// By default, Scott's rule is used for automatic bandwidth selection.
+    #[must_use]
+    pub fn kde_bandwidth(mut self, bandwidth: f64) -> Self {
+        self.kde_bandwidth = Some(bandwidth);
+        self
+    }
+
+    /// Sets a seed for reproducible sampling.
+    #[must_use]
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.seed = Some(seed);
+        self
+    }
+
+    /// Builds the configured [`MotpeSampler`].
+    #[must_use]
+    pub fn build(self) -> MotpeSampler {
+        let rng = match self.seed {
+            Some(s) => StdRng::seed_from_u64(s),
+            None => rand::make_rng(),
+        };
+
+        MotpeSampler {
+            n_startup_trials: self.n_startup_trials,
+            n_ei_candidates: self.n_ei_candidates,
+            kde_bandwidth: self.kde_bandwidth,
+            rng: Mutex::new(rng),
+        }
+    }
+}
+
+impl Default for MotpeSamplerBuilder {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+#[cfg(test)]
+#[allow(
+    clippy::similar_names,
+    clippy::cast_sign_loss,
+    clippy::cast_precision_loss
+)]
+mod tests {
+    use std::collections::HashMap;
+
+    use super::*;
+    use crate::distribution::{CategoricalDistribution, FloatDistribution, IntDistribution};
+    use crate::parameter::ParamId;
+
+    fn create_mo_trial(
+        id: u64,
+        values: Vec<f64>,
+        params: Vec<(ParamId, ParamValue, Distribution)>,
+    ) -> MultiObjectiveTrial {
+        let mut param_map = HashMap::new();
+        let mut dist_map = HashMap::new();
+        let label_map = HashMap::new();
+        for (param_id, pv, dist) in params {
+            param_map.insert(param_id, pv);
+            dist_map.insert(param_id, dist);
+        }
+        MultiObjectiveTrial {
+            id,
+            params: param_map,
+            distributions: dist_map,
+            param_labels: label_map,
+            values,
+            state: TrialState::Complete,
+            user_attrs: HashMap::new(),
+            constraints: Vec::new(),
+        }
+    }
+
+    #[test]
+    fn test_motpe_startup_random_sampling() {
+        let sampler = MotpeSampler::with_seed(42);
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+        let directions = [Direction::Minimize, Direction::Minimize];
+
+        // With no history, should use random sampling
+        let history: Vec<MultiObjectiveTrial> = vec![];
+        for _ in 0..50 {
+            let value = sampler.sample(&dist, 0, &history, &directions);
+            if let ParamValue::Float(v) = value {
+                assert!((0.0..=1.0).contains(&v));
+            } else {
+                panic!("Expected Float value");
+            }
+        }
+    }
+
+    #[test]
+    fn test_motpe_split_pareto() {
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+        let directions = [Direction::Minimize, Direction::Minimize];
+        let x_id = ParamId::new();
+
+        // Create trials: Pareto front = {(0.1, 0.9), (0.5, 0.5), (0.9, 0.1)}
+        // Dominated = {(0.6, 0.8), (0.8, 0.7)}
+        let history = vec![
+            create_mo_trial(
+                0,
+                vec![0.1, 0.9],
+                vec![(x_id, ParamValue::Float(0.1), dist.clone())],
+            ),
+            create_mo_trial(
+                1,
+                vec![0.5, 0.5],
+                vec![(x_id, ParamValue::Float(0.5), dist.clone())],
+            ),
+            create_mo_trial(
+                2,
+                vec![0.9, 0.1],
+                vec![(x_id, ParamValue::Float(0.9), dist.clone())],
+            ),
+            create_mo_trial(
+                3,
+                vec![0.6, 0.8],
+                vec![(x_id, ParamValue::Float(0.6), dist.clone())],
+            ),
+            create_mo_trial(
+                4,
+                vec![0.8, 0.7],
+                vec![(x_id, ParamValue::Float(0.8), dist.clone())],
+            ),
+        ];
+
+        let (good, bad) = MotpeSampler::split_trials(&history, &directions);
+        assert_eq!(good.len(), 3, "Pareto front should have 3 members");
+        assert_eq!(bad.len(), 2, "2 dominated trials");
+    }
+
+    #[test]
+    fn test_motpe_samples_float() {
+        let sampler = MotpeSampler::builder()
+            .n_startup_trials(5)
+            .n_ei_candidates(24)
+            .seed(42)
+            .build();
+
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+        let directions = [Direction::Minimize, Direction::Minimize];
+        let x_id = ParamId::new();
+
+        // Build history where values near 0.3 are on the Pareto front
+        let mut history = Vec::new();
+        for i in 0..20 {
+            let x = f64::from(i) / 20.0;
+            // Pareto front: f1 = (x - 0.3)^2, f2 = (x - 0.3)^2 + 0.1
+            // Best solutions cluster around x = 0.3
+            let f1 = (x - 0.3).powi(2);
+            let f2 = (x - 0.7).powi(2);
+            history.push(create_mo_trial(
+                i as u64,
+                vec![f1, f2],
+                vec![(x_id, ParamValue::Float(x), dist.clone())],
+            ));
+        }
+
+        // MOTPE should produce values within [0, 1]
+        for i in 0..50 {
+            let value = sampler.sample(&dist, 100 + i, &history, &directions);
+            if let ParamValue::Float(v) = value {
+                assert!((0.0..=1.0).contains(&v), "Value {v} out of range");
+            } else {
+                panic!("Expected Float value");
+            }
+        }
+    }
+
+    #[test]
+    fn test_motpe_int_sampling() {
+        let sampler = MotpeSampler::builder().n_startup_trials(5).seed(42).build();
+
+        let dist = Distribution::Int(IntDistribution {
+            low: 0,
+            high: 100,
+            log_scale: false,
+            step: None,
+        });
+        let directions = [Direction::Minimize, Direction::Minimize];
+        let x_id = ParamId::new();
+
+        let mut history = Vec::new();
+        for i in 0..20 {
+            let x = i * 5;
+            let f1 = ((x as f64) - 30.0).powi(2);
+            let f2 = ((x as f64) - 70.0).powi(2);
+            history.push(create_mo_trial(
+                i as u64,
+                vec![f1, f2],
+                vec![(x_id, ParamValue::Int(x), dist.clone())],
+            ));
+        }
+
+        for i in 0..50 {
+            let value = sampler.sample(&dist, 100 + i, &history, &directions);
+            if let ParamValue::Int(v) = value {
+                assert!((0..=100).contains(&v), "Value {v} out of range");
+            } else {
+                panic!("Expected Int value");
+            }
+        }
+    }
+
+    #[test]
+    fn test_motpe_categorical_sampling() {
+        let sampler = MotpeSampler::builder().n_startup_trials(5).seed(42).build();
+
+        let dist = Distribution::Categorical(CategoricalDistribution { n_choices: 3 });
+        let directions = [Direction::Minimize, Direction::Minimize];
+        let cat_id = ParamId::new();
+
+        // Category 1 is on the Pareto front, others are dominated
+        let mut history = Vec::new();
+        for i in 0..15 {
+            let category = i % 3;
+            let (f1, f2) = match category {
+                0 => (0.8, 0.8), // dominated
+                1 => (0.1, 0.9), // Pareto front
+                2 => (0.9, 0.1), // Pareto front
+                _ => unreachable!(),
+            };
+            history.push(create_mo_trial(
+                i as u64,
+                vec![f1, f2],
+                vec![(
+                    cat_id,
+                    ParamValue::Categorical(category as usize),
+                    dist.clone(),
+                )],
+            ));
+        }
+
+        let mut counts = vec![0usize; 3];
+        for i in 0..200 {
+            let value = sampler.sample(&dist, 100 + i, &history, &directions);
+            if let ParamValue::Categorical(idx) = value {
+                assert!(idx < 3, "Category {idx} out of range");
+                counts[idx] += 1;
+            } else {
+                panic!("Expected Categorical value");
+            }
+        }
+
+        // Categories 1 and 2 (on Pareto front) should dominate category 0
+        assert!(
+            counts[1] + counts[2] > counts[0],
+            "Pareto-front categories should be sampled more: {counts:?}"
+        );
+    }
+
+    #[test]
+    fn test_motpe_reproducibility() {
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+        let directions = [Direction::Minimize, Direction::Minimize];
+        let x_id = ParamId::new();
+
+        let history: Vec<MultiObjectiveTrial> = (0..20)
+            .map(|i| {
+                let x = f64::from(i) / 20.0;
+                create_mo_trial(
+                    i as u64,
+                    vec![x, 1.0 - x],
+                    vec![(x_id, ParamValue::Float(x), dist.clone())],
+                )
+            })
+            .collect();
+
+        let sampler1 = MotpeSampler::builder()
+            .seed(12345)
+            .n_startup_trials(5)
+            .build();
+        let sampler2 = MotpeSampler::builder()
+            .seed(12345)
+            .n_startup_trials(5)
+            .build();
+
+        for i in 0..10 {
+            let v1 = sampler1.sample(&dist, i, &history, &directions);
+            let v2 = sampler2.sample(&dist, i, &history, &directions);
+            assert_eq!(v1, v2, "Samples should be identical with same seed");
+        }
+    }
+
+    #[test]
+    fn test_motpe_with_study() {
+        use crate::multi_objective::MultiObjectiveStudy;
+        use crate::parameter::{FloatParam, Parameter};
+
+        let sampler = MotpeSampler::builder().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::<_, crate::Error>(vec![xv, 1.0 - xv])
+            })
+            .unwrap();
+
+        let front = study.pareto_front();
+        assert!(!front.is_empty(), "Should have Pareto-optimal solutions");
+
+        // All front solutions should have values summing to ~1.0
+        for trial in &front {
+            let sum: f64 = trial.values.iter().sum();
+            assert!(
+                (sum - 1.0).abs() < 0.01,
+                "Pareto front values should sum to ~1.0, got {sum}"
+            );
+        }
+    }
+
+    #[test]
+    fn test_motpe_builder_defaults() {
+        let sampler = MotpeSamplerBuilder::new().build();
+        assert_eq!(sampler.n_startup_trials, 11);
+        assert_eq!(sampler.n_ei_candidates, 24);
+        assert!(sampler.kde_bandwidth.is_none());
+    }
+
+    #[test]
+    fn test_motpe_builder_custom() {
+        let sampler = MotpeSamplerBuilder::new()
+            .n_startup_trials(20)
+            .n_ei_candidates(48)
+            .kde_bandwidth(0.5)
+            .seed(99)
+            .build();
+        assert_eq!(sampler.n_startup_trials, 20);
+        assert_eq!(sampler.n_ei_candidates, 48);
+        assert_eq!(sampler.kde_bandwidth, Some(0.5));
+    }
+}

src/sampler/nsga2.rs

@@ -0,0 +1,763 @@
+//! NSGA-II (Non-dominated Sorting Genetic Algorithm II) sampler.
+//!
+//! Implements multi-objective optimization using non-dominated sorting,
+//! crowding distance, SBX crossover, and polynomial mutation.
+//!
+//! # Examples
+//!
+//! ```
+//! use optimizer::Direction;
+//! use optimizer::multi_objective::MultiObjectiveStudy;
+//! use optimizer::parameter::{FloatParam, Parameter};
+//! use optimizer::sampler::nsga2::Nsga2Sampler;
+//!
+//! let sampler = Nsga2Sampler::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)?;
+//!         Ok::<_, optimizer::Error>(vec![xv * xv, (xv - 1.0).powi(2)])
+//!     })
+//!     .unwrap();
+//! ```
+
+use std::collections::HashMap;
+
+use parking_lot::Mutex;
+use rand::rngs::StdRng;
+use rand::{RngExt, SeedableRng};
+
+use crate::distribution::Distribution;
+use crate::multi_objective::MultiObjectiveTrial;
+use crate::param::ParamValue;
+use crate::pareto;
+use crate::types::Direction;
+
+/// NSGA-II sampler for multi-objective optimization.
+///
+/// Provides non-dominated sorting, crowding distance selection,
+/// SBX crossover, and polynomial mutation.
+pub struct Nsga2Sampler {
+    state: Mutex<Nsga2State>,
+}
+
+impl Nsga2Sampler {
+    /// Creates a new NSGA-II sampler with a random seed.
+    #[must_use]
+    pub fn new() -> Self {
+        Self {
+            state: Mutex::new(Nsga2State::new(Nsga2Config::default(), None)),
+        }
+    }
+
+    /// Creates a new NSGA-II sampler with a fixed seed.
+    #[must_use]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            state: Mutex::new(Nsga2State::new(Nsga2Config::default(), Some(seed))),
+        }
+    }
+
+    /// Creates a builder for configuring an `Nsga2Sampler`.
+    #[must_use]
+    pub fn builder() -> Nsga2SamplerBuilder {
+        Nsga2SamplerBuilder::default()
+    }
+}
+
+impl Default for Nsga2Sampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+/// Builder for [`Nsga2Sampler`].
+#[derive(Debug, Clone, Default)]
+pub struct Nsga2SamplerBuilder {
+    population_size: Option<usize>,
+    crossover_prob: Option<f64>,
+    crossover_eta: Option<f64>,
+    mutation_eta: Option<f64>,
+    seed: Option<u64>,
+}
+
+impl Nsga2SamplerBuilder {
+    /// Sets the population size. Default: `4 + floor(3 * ln(n_params))`, minimum 4.
+    #[must_use]
+    pub fn population_size(mut self, size: usize) -> Self {
+        self.population_size = Some(size);
+        self
+    }
+
+    /// Sets the crossover probability. Default: 0.9.
+    #[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 [`Nsga2Sampler`].
+    #[must_use]
+    pub fn build(self) -> Nsga2Sampler {
+        let config = Nsga2Config {
+            user_population_size: self.population_size,
+            crossover_prob: self.crossover_prob.unwrap_or(0.9),
+            crossover_eta: self.crossover_eta.unwrap_or(20.0),
+            mutation_eta: self.mutation_eta.unwrap_or(20.0),
+        };
+        Nsga2Sampler {
+            state: Mutex::new(Nsga2State::new(config, self.seed)),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Internal types
+// ---------------------------------------------------------------------------
+
+#[derive(Clone, Debug)]
+struct Nsga2Config {
+    user_population_size: Option<usize>,
+    crossover_prob: f64,
+    crossover_eta: f64,
+    mutation_eta: f64,
+}
+
+impl Default for Nsga2Config {
+    fn default() -> Self {
+        Self {
+            user_population_size: None,
+            crossover_prob: 0.9,
+            crossover_eta: 20.0,
+            mutation_eta: 20.0,
+        }
+    }
+}
+
+/// 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,
+    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,
+            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,
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// MultiObjectiveSampler implementation
+// ---------------------------------------------------------------------------
+
+impl crate::multi_objective::MultiObjectiveSampler for Nsga2Sampler {
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        history: &[MultiObjectiveTrial],
+        directions: &[Direction],
+    ) -> 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),
+        }
+    }
+}
+
+/// 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(
+    state: &mut Nsga2State,
+    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);
+        }
+        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;
+    }
+
+    // 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;
+}
+
+// ---------------------------------------------------------------------------
+// NSGA-II generation algorithm
+// ---------------------------------------------------------------------------
+
+/// Performs NSGA-II selection: non-dominated sort + crowding distance,
+/// then selects `pop_size` parents from the population.
+fn nsga2_select(
+    state: &mut Nsga2State,
+    population: &[&MultiObjectiveTrial],
+    directions: &[Direction],
+) -> (Vec<Vec<ParamValue>>, Vec<usize>, Vec<f64>) {
+    let pop_size = state.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();
+    let has_constraints = constraints.iter().any(|c| !c.is_empty());
+
+    let fronts = if has_constraints {
+        pareto::fast_non_dominated_sort_constrained(&values, directions, &constraints)
+    } else {
+        pareto::fast_non_dominated_sort(&values, directions)
+    };
+
+    let n = population.len();
+    let mut rank = vec![0_usize; n];
+    let mut crowding = vec![0.0_f64; n];
+
+    for (front_rank, front) in fronts.iter().enumerate() {
+        let cd = pareto::crowding_distance_indexed(front, &values);
+        for (i, &idx) in front.iter().enumerate() {
+            rank[idx] = front_rank;
+            crowding[idx] = cd[i];
+        }
+    }
+
+    let mut selected: Vec<usize> = Vec::with_capacity(pop_size);
+    for front in &fronts {
+        if selected.len() + front.len() <= pop_size {
+            selected.extend_from_slice(front);
+        } else {
+            let remaining = pop_size - selected.len();
+            let mut front_sorted: Vec<usize> = front.clone();
+            front_sorted.sort_by(|&a, &b| {
+                crowding[b]
+                    .partial_cmp(&crowding[a])
+                    .unwrap_or(core::cmp::Ordering::Equal)
+            });
+            selected.extend_from_slice(&front_sorted[..remaining]);
+            break;
+        }
+    }
+
+    while selected.len() < pop_size {
+        selected.push(state.rng.random_range(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))
+        .collect();
+
+    let sel_rank: Vec<usize> = selected.iter().map(|&i| rank[i]).collect();
+    let sel_crowding: Vec<f64> = selected.iter().map(|&i| crowding[i]).collect();
+
+    (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;
+
+    if population.len() < 2 {
+        return (0..pop_size)
+            .map(|_| {
+                let params = state
+                    .dimensions
+                    .iter()
+                    .map(|d| sample_random(&mut state.rng, &d.distribution))
+                    .collect();
+                Candidate { params }
+            })
+            .collect();
+    }
+
+    let (parents, sel_rank, sel_crowding) = nsga2_select(state, population, directions);
+
+    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 (mut child1, mut child2) = crossover(
+            &mut state.rng,
+            &parents[p1],
+            &parents[p2],
+            &state.dimensions,
+            state.config.crossover_prob,
+            state.config.crossover_eta,
+        );
+
+        mutate(
+            &mut state.rng,
+            &mut child1,
+            &state.dimensions,
+            state.config.mutation_eta,
+        );
+        mutate(
+            &mut state.rng,
+            &mut child2,
+            &state.dimensions,
+            state.config.mutation_eta,
+        );
+
+        offspring.push(Candidate { params: child1 });
+        if offspring.len() < pop_size {
+            offspring.push(Candidate { params: child2 });
+        }
+    }
+
+    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);
+
+    if ranks[a] < ranks[b] {
+        a
+    } else if ranks[b] < ranks[a] {
+        b
+    } else if crowding[a] >= crowding[b] {
+        a
+    } else {
+        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/sobol.rs

@@ -0,0 +1,420 @@
+//! Quasi-random sampler using Sobol low-discrepancy sequences.
+
+use parking_lot::Mutex;
+use sobol_burley::sample;
+
+use crate::distribution::Distribution;
+use crate::param::ParamValue;
+use crate::sampler::{CompletedTrial, Sampler};
+
+/// Internal state for tracking the dimension counter within a trial.
+struct SobolState {
+    /// The `trial_id` of the current trial (used to reset dimension counter).
+    current_trial: u64,
+    /// Next Sobol dimension to use for the current trial.
+    next_dimension: u32,
+}
+
+/// Quasi-random sampler using Sobol low-discrepancy sequences.
+///
+/// Provides better uniform coverage of the parameter space than
+/// [`RandomSampler`](super::random::RandomSampler). Useful as a baseline or
+/// for the startup phase of model-based samplers.
+///
+/// Unlike random sampling, Sobol sequences are deterministic and fill the
+/// space more evenly, reducing the number of trials needed to adequately
+/// cover the search space.
+///
+/// Each trial uses a different Sobol sequence index, and each parameter
+/// within a trial maps to a different Sobol dimension. Parameters must be
+/// suggested in the same order across trials for consistent dimension
+/// assignment.
+///
+/// Sobol sequences are most effective in moderate dimensions (up to ~20).
+/// For very high dimensions, the uniformity advantage diminishes.
+///
+/// Requires the `sobol` feature flag.
+///
+/// # Examples
+///
+/// ```
+/// use optimizer::sampler::sobol::SobolSampler;
+/// use optimizer::{Direction, Study};
+///
+/// let study: Study<f64> = Study::with_sampler(Direction::Minimize, SobolSampler::new());
+/// ```
+pub struct SobolSampler {
+    seed: u32,
+    state: Mutex<SobolState>,
+}
+
+impl SobolSampler {
+    /// Creates a new Sobol sampler with a default seed of 0.
+    #[must_use]
+    pub fn new() -> Self {
+        Self::with_seed(0)
+    }
+
+    /// Creates a new Sobol sampler with the given seed.
+    ///
+    /// Different seeds produce statistically independent Sobol sequences.
+    /// Using the same seed will produce the same sequence of sampled values.
+    #[must_use]
+    #[allow(clippy::cast_possible_truncation)]
+    pub fn with_seed(seed: u64) -> Self {
+        Self {
+            seed: seed as u32,
+            state: Mutex::new(SobolState {
+                current_trial: u64::MAX,
+                next_dimension: 0,
+            }),
+        }
+    }
+}
+
+impl Default for SobolSampler {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Sampler for SobolSampler {
+    #[allow(clippy::cast_possible_truncation, clippy::cast_precision_loss)]
+    fn sample(
+        &self,
+        distribution: &Distribution,
+        trial_id: u64,
+        _history: &[CompletedTrial],
+    ) -> ParamValue {
+        let mut state = self.state.lock();
+
+        // Reset dimension counter when a new trial starts.
+        if state.current_trial != trial_id {
+            state.current_trial = trial_id;
+            state.next_dimension = 0;
+        }
+
+        let dimension = state.next_dimension;
+        state.next_dimension = dimension + 1;
+
+        // Use trial_id as the Sobol sequence index.
+        let index = trial_id as u32;
+
+        // Generate a quasi-random point in [0, 1).
+        let point = f64::from(sample(index, dimension, self.seed));
+
+        map_point_to_distribution(point, distribution)
+    }
+}
+
+/// Maps a uniform [0, 1) point to a value within the given distribution.
+#[allow(
+    clippy::cast_possible_truncation,
+    clippy::cast_precision_loss,
+    clippy::cast_sign_loss
+)]
+fn map_point_to_distribution(point: f64, 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();
+                (log_low + point * (log_high - log_low)).exp()
+            } else if let Some(step) = d.step {
+                let n_steps = ((d.high - d.low) / step).floor() as i64;
+                let k = (point * (n_steps + 1) as f64).floor() as i64;
+                let k = k.min(n_steps);
+                d.low + (k as f64) * step
+            } else {
+                d.low + point * (d.high - d.low)
+            };
+            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 = (log_low + point * (log_high - log_low)).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 = (point * (n_steps + 1) as f64).floor() as i64;
+                let k = k.min(n_steps);
+                d.low + k * step
+            } else {
+                let range = d.high - d.low + 1;
+                let k = (point * range as f64).floor() as i64;
+                (d.low + k).min(d.high)
+            };
+            ParamValue::Int(value)
+        }
+        Distribution::Categorical(d) => {
+            let index = (point * d.n_choices as f64).floor() as usize;
+            let index = index.min(d.n_choices - 1);
+            ParamValue::Categorical(index)
+        }
+    }
+}
+
+#[cfg(test)]
+#[allow(
+    clippy::cast_possible_truncation,
+    clippy::cast_precision_loss,
+    clippy::cast_sign_loss
+)]
+mod tests {
+    use super::*;
+    use crate::distribution::{CategoricalDistribution, FloatDistribution, IntDistribution};
+
+    #[test]
+    fn float_within_bounds() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Float(FloatDistribution {
+            low: -5.0,
+            high: 5.0,
+            log_scale: false,
+            step: None,
+        });
+
+        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 float_log_scale_within_bounds() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Float(FloatDistribution {
+            low: 1e-5,
+            high: 1.0,
+            log_scale: true,
+            step: None,
+        });
+
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Float(v) = value {
+                assert!(
+                    (1e-5..=1.0).contains(&v),
+                    "value {v} out of bounds at trial {i}"
+                );
+            } else {
+                panic!("Expected Float value");
+            }
+        }
+    }
+
+    #[test]
+    fn float_step_respects_grid() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: Some(0.25),
+        });
+
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Float(v) = value {
+                assert!((0.0..=1.0).contains(&v), "value {v} out of bounds");
+                let k = (v / 0.25).round() as i64;
+                let expected = k as f64 * 0.25;
+                assert!((v - expected).abs() < 1e-10, "value {v} not on step grid");
+            } else {
+                panic!("Expected Float value");
+            }
+        }
+    }
+
+    #[test]
+    fn int_within_bounds() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Int(IntDistribution {
+            low: 0,
+            high: 10,
+            log_scale: false,
+            step: None,
+        });
+
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Int(v) = value {
+                assert!(
+                    (0..=10).contains(&v),
+                    "value {v} out of bounds at trial {i}"
+                );
+            } else {
+                panic!("Expected Int value");
+            }
+        }
+    }
+
+    #[test]
+    fn int_log_scale_within_bounds() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Int(IntDistribution {
+            low: 1,
+            high: 1000,
+            log_scale: true,
+            step: None,
+        });
+
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Int(v) = value {
+                assert!(
+                    (1..=1000).contains(&v),
+                    "value {v} out of bounds at trial {i}"
+                );
+            } else {
+                panic!("Expected Int value");
+            }
+        }
+    }
+
+    #[test]
+    fn int_step_respects_grid() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Int(IntDistribution {
+            low: 0,
+            high: 10,
+            log_scale: false,
+            step: Some(2),
+        });
+
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Int(v) = value {
+                assert!((0..=10).contains(&v), "value {v} out of bounds");
+                assert!(v % 2 == 0, "value {v} not on step grid");
+            } else {
+                panic!("Expected Int value");
+            }
+        }
+    }
+
+    #[test]
+    fn categorical_within_bounds() {
+        let sampler = SobolSampler::with_seed(42);
+        let dist = Distribution::Categorical(CategoricalDistribution { n_choices: 5 });
+
+        for i in 0..100 {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Categorical(idx) = value {
+                assert!(idx < 5, "index {idx} out of bounds at trial {i}");
+            } else {
+                panic!("Expected Categorical value");
+            }
+        }
+    }
+
+    #[test]
+    fn deterministic_with_same_seed() {
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        let sampler1 = SobolSampler::with_seed(42);
+        let sampler2 = SobolSampler::with_seed(42);
+
+        for i in 0..20 {
+            let v1 = sampler1.sample(&dist, i, &[]);
+            let v2 = sampler2.sample(&dist, i, &[]);
+            assert_eq!(v1, v2, "mismatch at trial {i}");
+        }
+    }
+
+    #[test]
+    fn different_seeds_produce_different_sequences() {
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        let sampler1 = SobolSampler::with_seed(0);
+        let sampler2 = SobolSampler::with_seed(12345);
+
+        let mut any_different = false;
+        for i in 0..20 {
+            let v1 = sampler1.sample(&dist, i, &[]);
+            let v2 = sampler2.sample(&dist, i, &[]);
+            if v1 != v2 {
+                any_different = true;
+                break;
+            }
+        }
+        assert!(
+            any_different,
+            "different seeds should produce different sequences"
+        );
+    }
+
+    #[test]
+    fn better_coverage_than_random() {
+        // Sobol sequence should cover [0,1] more uniformly than random.
+        // We measure this by checking that the Sobol samples fill all
+        // 10 equal-width bins with only 20 samples.
+        let sampler = SobolSampler::with_seed(0);
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        let n_bins = 10;
+        let n_samples = 20;
+        let mut bins = vec![0u32; n_bins];
+
+        for i in 0..n_samples {
+            let value = sampler.sample(&dist, i, &[]);
+            if let ParamValue::Float(v) = value {
+                let bin = ((v * n_bins as f64).floor() as usize).min(n_bins - 1);
+                bins[bin] += 1;
+            }
+        }
+
+        let filled_bins = bins.iter().filter(|&&c| c > 0).count();
+        assert!(
+            filled_bins >= 8,
+            "Expected at least 8/10 bins filled, got {filled_bins}: {bins:?}"
+        );
+    }
+
+    #[test]
+    fn multi_parameter_uses_different_dimensions() {
+        // When sampling multiple parameters per trial, each parameter
+        // should get a different Sobol dimension, producing different values.
+        let sampler = SobolSampler::with_seed(0);
+        let dist = Distribution::Float(FloatDistribution {
+            low: 0.0,
+            high: 1.0,
+            log_scale: false,
+            step: None,
+        });
+
+        // Sample two parameters for trial 0.
+        let v1 = sampler.sample(&dist, 0, &[]);
+        let v2 = sampler.sample(&dist, 0, &[]);
+
+        // They should differ (different Sobol dimensions for the same index).
+        assert_ne!(
+            v1, v2,
+            "multi-parameter samples should use different dimensions"
+        );
+    }
+}

src/trial.rs

@@ -15,6 +15,7 @@ use crate::types::TrialState;
 
 /// A user attribute value that can be stored on a trial.
 #[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub enum AttrValue {
     /// A floating-point attribute.
     Float(f64),
@@ -88,6 +89,8 @@ pub struct Trial {
     user_attrs: HashMap<String, AttrValue>,
     /// Pre-filled parameter values from enqueue (used instead of sampling).
     fixed_params: HashMap<ParamId, ParamValue>,
+    /// Constraint values for this trial (<=0.0 means feasible).
+    constraint_values: Vec<f64>,
 }
 
 impl core::fmt::Debug for Trial {
@@ -104,6 +107,7 @@ impl core::fmt::Debug for Trial {
             .field("has_pruner", &self.pruner.is_some())
             .field("user_attrs", &self.user_attrs)
             .field("fixed_params", &self.fixed_params)
+            .field("constraint_values", &self.constraint_values)
             .finish()
     }
 }
@@ -143,6 +147,7 @@ impl Trial {
             pruner: None,
             user_attrs: HashMap::new(),
             fixed_params: HashMap::new(),
+            constraint_values: Vec::new(),
         }
     }
 
@@ -174,6 +179,7 @@ impl Trial {
             pruner: Some(pruner),
             user_attrs: HashMap::new(),
             fixed_params: HashMap::new(),
+            constraint_values: Vec::new(),
         }
     }
 
@@ -261,7 +267,11 @@ impl Trial {
             return false;
         };
         let history_guard = history.read();
-        pruner.should_prune(self.id, step, &self.intermediate_values, &history_guard)
+        let prune = pruner.should_prune(self.id, step, &self.intermediate_values, &history_guard);
+        if prune {
+            trace_info!(trial_id = self.id, step, "pruner recommends stopping");
+        }
+        prune
     }
 
     /// Returns all intermediate values reported so far.
@@ -287,6 +297,20 @@ impl Trial {
         &self.user_attrs
     }
 
+    /// Sets constraint values for this trial.
+    ///
+    /// Each value represents a constraint; a value <= 0.0 means the constraint
+    /// is satisfied (feasible). A value > 0.0 means the constraint is violated.
+    pub fn set_constraints(&mut self, values: Vec<f64>) {
+        self.constraint_values = values;
+    }
+
+    /// Returns the constraint values for this trial.
+    #[must_use]
+    pub fn constraint_values(&self) -> &[f64] {
+        &self.constraint_values
+    }
+
     /// Sets the trial state to Complete.
     pub(crate) fn set_complete(&mut self) {
         self.state = TrialState::Complete;
@@ -366,6 +390,13 @@ impl Trial {
 
         let result = param.cast_param_value(&value)?;
 
+        trace_debug!(
+            trial_id = self.id,
+            param = %param.label(),
+            value = %value,
+            "parameter sampled"
+        );
+
         // Store distribution, value, and label
         self.distributions.insert(param_id, distribution);
         self.params.insert(param_id, value);

src/types.rs

@@ -2,6 +2,7 @@
 
 /// The direction of optimization.
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub enum Direction {
     /// Minimize the objective value.
     Minimize,
@@ -11,6 +12,7 @@ pub enum Direction {
 
 /// The state of a trial in its lifecycle.
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
+#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
 pub enum TrialState {
     /// The trial is currently running.
     Running,

src/visualization.rs

@@ -0,0 +1,458 @@
+//! HTML report generation for optimization visualization.
+//!
+//! Generates self-contained HTML files with embedded Plotly.js charts
+//! for offline visualization of optimization results.
+
+use core::fmt::Write as _;
+use std::collections::BTreeMap;
+use std::path::Path;
+
+use crate::distribution::Distribution;
+use crate::param::ParamValue;
+use crate::parameter::ParamId;
+use crate::sampler::CompletedTrial;
+use crate::types::{Direction, TrialState};
+
+/// Generate an HTML report with interactive Plotly.js charts.
+///
+/// Creates a self-contained HTML file at `path` containing:
+/// - **Optimization history**: Objective value vs trial number with best-so-far line
+/// - **Slice plots**: Objective value vs each parameter (1D scatter)
+/// - **Parallel coordinates**: Multi-parameter relationship view
+/// - **Trial timeline**: Duration index of each trial (horizontal bar)
+/// - **Intermediate values**: Learning curves per trial (if pruning data available)
+/// - **Parameter importance**: Bar chart (if enough completed trials)
+///
+/// # Errors
+///
+/// Returns an I/O error if the file cannot be created or written.
+pub fn generate_html_report(
+    study: &crate::Study<f64>,
+    path: impl AsRef<Path>,
+) -> std::io::Result<()> {
+    let trials = study.trials();
+    let direction = study.direction();
+    let importance = study.param_importance();
+
+    let html = build_html(&trials, direction, &importance);
+    std::fs::write(path, html)
+}
+
+fn build_html(
+    trials: &[CompletedTrial<f64>],
+    direction: Direction,
+    importance: &[(String, f64)],
+) -> String {
+    let mut html = String::with_capacity(8192);
+
+    let dir_label = match direction {
+        Direction::Minimize => "Minimize",
+        Direction::Maximize => "Maximize",
+    };
+
+    // Collect parameter metadata.
+    let param_info = collect_param_info(trials);
+    let has_intermediate = trials.iter().any(|t| !t.intermediate_values.is_empty());
+
+    let _ = write!(
+        html,
+        r#"<!DOCTYPE html>
+<html lang="en">
+<head>
+<meta charset="utf-8">
+<meta name="viewport" content="width=device-width, initial-scale=1">
+<title>Optimization Report</title>
+<script src="https://cdn.plot.ly/plotly-2.35.2.min.js"></script>
+<style>
+  * {{ margin: 0; padding: 0; box-sizing: border-box; }}
+  body {{ font-family: -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, sans-serif;
+         background: #f5f6fa; color: #2c3e50; padding: 24px; }}
+  h1 {{ text-align: center; margin-bottom: 8px; font-size: 1.8em; }}
+  .subtitle {{ text-align: center; color: #7f8c8d; margin-bottom: 24px; }}
+  .chart {{ background: #fff; border-radius: 8px; box-shadow: 0 2px 8px rgba(0,0,0,0.08);
+            margin-bottom: 24px; padding: 16px; }}
+  .chart-title {{ font-size: 1.1em; font-weight: 600; margin-bottom: 8px; }}
+</style>
+</head>
+<body>
+<h1>Optimization Report</h1>
+<p class="subtitle">{dir_label} &middot; {n} trials</p>
+"#,
+        n = trials.len(),
+    );
+
+    // Optimization history chart.
+    html.push_str("<div class=\"chart\"><div class=\"chart-title\">Optimization History</div><div id=\"history\"></div></div>\n");
+    write_history_chart(&mut html, trials, direction);
+
+    // Slice plots.
+    if !param_info.is_empty() {
+        html.push_str("<div class=\"chart\"><div class=\"chart-title\">Slice Plots</div><div id=\"slices\"></div></div>\n");
+        write_slice_charts(&mut html, trials, &param_info);
+    }
+
+    // Parallel coordinates.
+    if param_info.len() >= 2 {
+        html.push_str("<div class=\"chart\"><div class=\"chart-title\">Parallel Coordinates</div><div id=\"parcoords\"></div></div>\n");
+        write_parallel_coordinates(&mut html, trials, &param_info, direction);
+    }
+
+    // Parameter importance.
+    if !importance.is_empty() {
+        html.push_str("<div class=\"chart\"><div class=\"chart-title\">Parameter Importance</div><div id=\"importance\"></div></div>\n");
+        write_importance_chart(&mut html, importance);
+    }
+
+    // Trial timeline.
+    html.push_str("<div class=\"chart\"><div class=\"chart-title\">Trial Timeline</div><div id=\"timeline\"></div></div>\n");
+    write_timeline_chart(&mut html, trials);
+
+    // Intermediate values.
+    if has_intermediate {
+        html.push_str("<div class=\"chart\"><div class=\"chart-title\">Intermediate Values</div><div id=\"intermediate\"></div></div>\n");
+        write_intermediate_chart(&mut html, trials);
+    }
+
+    html.push_str("</body>\n</html>\n");
+    html
+}
+
+/// Metadata about each parameter seen across trials.
+struct ParamMeta {
+    label: String,
+    #[allow(dead_code)]
+    dist: Option<Distribution>,
+}
+
+/// Collect parameter labels and distributions across all trials.
+fn collect_param_info(trials: &[CompletedTrial<f64>]) -> BTreeMap<ParamId, ParamMeta> {
+    let mut info: BTreeMap<ParamId, ParamMeta> = BTreeMap::new();
+    for trial in trials {
+        for &id in trial.params.keys() {
+            info.entry(id).or_insert_with(|| {
+                let label = trial
+                    .param_labels
+                    .get(&id)
+                    .cloned()
+                    .unwrap_or_else(|| id.to_string());
+                let dist = trial.distributions.get(&id).cloned();
+                ParamMeta { label, dist }
+            });
+        }
+    }
+    info
+}
+
+// ---------------------------------------------------------------------------
+// Chart generators
+// ---------------------------------------------------------------------------
+
+fn write_history_chart(html: &mut String, trials: &[CompletedTrial<f64>], direction: Direction) {
+    let complete: Vec<_> = trials
+        .iter()
+        .filter(|t| t.state == TrialState::Complete)
+        .collect();
+    if complete.is_empty() {
+        return;
+    }
+
+    let mut ids = Vec::with_capacity(complete.len());
+    let mut vals = Vec::with_capacity(complete.len());
+    let mut best_vals = Vec::with_capacity(complete.len());
+    let mut best = complete[0].value;
+    for t in &complete {
+        ids.push(t.id);
+        vals.push(t.value);
+        best = match direction {
+            Direction::Minimize => best.min(t.value),
+            Direction::Maximize => best.max(t.value),
+        };
+        best_vals.push(best);
+    }
+
+    let _ = write!(
+        html,
+        r##"<script>
+Plotly.newPlot("history", [
+  {{ x: {ids:?}, y: {vals:?}, mode: "markers", name: "Objective", type: "scatter",
+     marker: {{ color: "#3498db", size: 6 }} }},
+  {{ x: {ids:?}, y: {best_vals:?}, mode: "lines", name: "Best so far", type: "scatter",
+     line: {{ color: "#e74c3c", width: 2 }} }}
+], {{ xaxis: {{ title: "Trial" }}, yaxis: {{ title: "Objective Value" }},
+     margin: {{ t: 10 }}, legend: {{ x: 1, xanchor: "right", y: 1 }} }},
+   {{ responsive: true }});
+</script>
+"##,
+    );
+}
+
+fn write_slice_charts(
+    html: &mut String,
+    trials: &[CompletedTrial<f64>],
+    param_info: &BTreeMap<ParamId, ParamMeta>,
+) {
+    let complete: Vec<_> = trials
+        .iter()
+        .filter(|t| t.state == TrialState::Complete)
+        .collect();
+    if complete.is_empty() {
+        return;
+    }
+
+    let n_params = param_info.len();
+    let cols = if n_params <= 2 { n_params } else { 2 };
+    let rows = n_params.div_ceil(cols);
+
+    // Build subplot titles and data.
+    let mut subplot_titles = Vec::new();
+    let mut traces = String::new();
+    for (i, (id, meta)) in param_info.iter().enumerate() {
+        subplot_titles.push(format!("\"{}\"", escape_js(&meta.label)));
+
+        let mut x_vals = Vec::new();
+        let mut y_vals = Vec::new();
+        for t in &complete {
+            if let Some(pv) = t.params.get(id) {
+                x_vals.push(param_value_to_f64(pv));
+                y_vals.push(t.value);
+            }
+        }
+
+        let subplot_idx = i + 1;
+        let xa = if subplot_idx == 1 {
+            "x".to_string()
+        } else {
+            format!("x{subplot_idx}")
+        };
+        let ya = if subplot_idx == 1 {
+            "y".to_string()
+        } else {
+            format!("y{subplot_idx}")
+        };
+        let _ = write!(
+            traces,
+            r##"{{ x: {x_vals:?}, y: {y_vals:?}, mode: "markers", type: "scatter",
+               xaxis: "{xa}", yaxis: "{ya}",
+               marker: {{ color: "#3498db", size: 5 }}, showlegend: false }},"##,
+        );
+    }
+
+    let _ = write!(
+        html,
+        r#"<script>
+Plotly.newPlot("slices", [{traces}],
+  {{ grid: {{ rows: {rows}, columns: {cols}, pattern: "independent" }},
+     annotations: [{annotations}],
+     margin: {{ t: 30 }}, showlegend: false }},
+  {{ responsive: true }});
+</script>
+"#,
+        annotations = build_subplot_annotations(&subplot_titles, rows, cols),
+    );
+}
+
+fn write_parallel_coordinates(
+    html: &mut String,
+    trials: &[CompletedTrial<f64>],
+    param_info: &BTreeMap<ParamId, ParamMeta>,
+    direction: Direction,
+) {
+    let complete: Vec<_> = trials
+        .iter()
+        .filter(|t| t.state == TrialState::Complete)
+        .collect();
+    if complete.is_empty() {
+        return;
+    }
+
+    let mut dimensions = String::new();
+
+    // Add objective value as the first dimension.
+    let obj_vals: Vec<f64> = complete.iter().map(|t| t.value).collect();
+    let _ = write!(
+        dimensions,
+        r#"{{ label: "Objective", values: {obj_vals:?} }},"#,
+    );
+
+    // Add each parameter as a dimension.
+    for (id, meta) in param_info {
+        let vals: Vec<f64> = complete
+            .iter()
+            .map(|t| t.params.get(id).map_or(f64::NAN, param_value_to_f64))
+            .collect();
+        let _ = write!(
+            dimensions,
+            r#"{{ label: "{label}", values: {vals:?} }},"#,
+            label = escape_js(&meta.label),
+        );
+    }
+
+    // Color by objective value: green = good, red = bad.
+    let (cmin, cmax) = min_max(&obj_vals);
+    let colorscale = match direction {
+        Direction::Minimize => r##"[[0,"#2ecc71"],[1,"#e74c3c"]]"##,
+        Direction::Maximize => r##"[[0,"#e74c3c"],[1,"#2ecc71"]]"##,
+    };
+
+    let _ = write!(
+        html,
+        r#"<script>
+Plotly.newPlot("parcoords", [{{
+  type: "parcoords",
+  line: {{ color: {obj_vals:?}, colorscale: {colorscale},
+           cmin: {cmin}, cmax: {cmax}, showscale: true }},
+  dimensions: [{dimensions}]
+}}], {{ margin: {{ t: 10 }} }}, {{ responsive: true }});
+</script>
+"#,
+    );
+}
+
+fn write_importance_chart(html: &mut String, importance: &[(String, f64)]) {
+    let names: Vec<_> = importance.iter().map(|(n, _)| format!("\"{n}\"")).collect();
+    let values: Vec<f64> = importance.iter().map(|(_, v)| *v).collect();
+
+    let _ = write!(
+        html,
+        r##"<script>
+Plotly.newPlot("importance", [{{
+  x: {values:?}, y: [{names}], type: "bar", orientation: "h",
+  marker: {{ color: "#9b59b6" }}
+}}], {{ xaxis: {{ title: "Importance (|Spearman correlation|)" }},
+       yaxis: {{ automargin: true }}, margin: {{ t: 10, l: 120 }} }},
+   {{ responsive: true }});
+</script>
+"##,
+        names = names.join(","),
+    );
+}
+
+fn write_timeline_chart(html: &mut String, trials: &[CompletedTrial<f64>]) {
+    let mut ids = Vec::with_capacity(trials.len());
+    let mut colors = Vec::with_capacity(trials.len());
+    let mut labels = Vec::with_capacity(trials.len());
+
+    for t in trials {
+        ids.push(format!("\"Trial {}\"", t.id));
+        let (color, label) = match t.state {
+            TrialState::Complete => ("#2ecc71", "Complete"),
+            TrialState::Pruned => ("#f39c12", "Pruned"),
+            TrialState::Failed => ("#e74c3c", "Failed"),
+            TrialState::Running => ("#3498db", "Running"),
+        };
+        colors.push(format!("\"{color}\""));
+        labels.push(format!("\"{label}\""));
+    }
+
+    // Use trial index as a proxy for duration (no wallclock data available).
+    let indices: Vec<usize> = (0..trials.len()).collect();
+
+    let _ = write!(
+        html,
+        r#"<script>
+Plotly.newPlot("timeline", [{{
+  y: [{ids}], x: {indices:?}, type: "bar", orientation: "h",
+  text: [{labels}], textposition: "auto",
+  marker: {{ color: [{colors}] }}
+}}], {{ xaxis: {{ title: "Trial Index" }}, yaxis: {{ automargin: true, autorange: "reversed" }},
+       margin: {{ t: 10, l: 80 }}, showlegend: false }},
+   {{ responsive: true }});
+</script>
+"#,
+        ids = ids.join(","),
+        colors = colors.join(","),
+        labels = labels.join(","),
+    );
+}
+
+fn write_intermediate_chart(html: &mut String, trials: &[CompletedTrial<f64>]) {
+    let trials_with_iv: Vec<_> = trials
+        .iter()
+        .filter(|t| !t.intermediate_values.is_empty())
+        .collect();
+    if trials_with_iv.is_empty() {
+        return;
+    }
+
+    let mut traces = String::new();
+    for t in &trials_with_iv {
+        let steps: Vec<u64> = t.intermediate_values.iter().map(|(s, _)| *s).collect();
+        let values: Vec<f64> = t.intermediate_values.iter().map(|(_, v)| *v).collect();
+        let color = match t.state {
+            TrialState::Pruned => "#f39c12",
+            _ => "#3498db",
+        };
+        let _ = write!(
+            traces,
+            r#"{{ x: {steps:?}, y: {values:?}, mode: "lines+markers", name: "Trial {id}",
+               line: {{ color: "{color}", width: 1 }}, marker: {{ size: 3 }} }},"#,
+            id = t.id,
+        );
+    }
+
+    let _ = write!(
+        html,
+        r#"<script>
+Plotly.newPlot("intermediate", [{traces}],
+  {{ xaxis: {{ title: "Step" }}, yaxis: {{ title: "Intermediate Value" }},
+     margin: {{ t: 10 }}, showlegend: true }},
+  {{ responsive: true }});
+</script>
+"#,
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Helpers
+// ---------------------------------------------------------------------------
+
+#[allow(clippy::cast_precision_loss)]
+fn param_value_to_f64(pv: &ParamValue) -> f64 {
+    match *pv {
+        ParamValue::Float(v) => v,
+        ParamValue::Int(v) => v as f64,
+        ParamValue::Categorical(v) => v as f64,
+    }
+}
+
+fn escape_js(s: &str) -> String {
+    s.replace('\\', "\\\\")
+        .replace('"', "\\\"")
+        .replace('\n', "\\n")
+}
+
+fn min_max(vals: &[f64]) -> (f64, f64) {
+    let mut mn = f64::INFINITY;
+    let mut mx = f64::NEG_INFINITY;
+    for &v in vals {
+        if v < mn {
+            mn = v;
+        }
+        if v > mx {
+            mx = v;
+        }
+    }
+    (mn, mx)
+}
+
+/// Build Plotly annotation objects to act as subplot titles.
+#[allow(clippy::cast_precision_loss)]
+fn build_subplot_annotations(titles: &[String], rows: usize, cols: usize) -> String {
+    let mut anns = Vec::new();
+    for (i, title) in titles.iter().enumerate() {
+        let row = i / cols;
+        let col = i % cols;
+        // Compute x/y anchor in paper coordinates.
+        let x = if cols == 1 {
+            0.5
+        } else {
+            (f64::from(u32::try_from(col).unwrap_or(0))) / (cols as f64 - 1.0)
+        };
+        let y = 1.0 - (f64::from(u32::try_from(row).unwrap_or(0))) / (rows as f64).max(1.0) + 0.02;
+        anns.push(format!(
+            r#"{{ text: {title}, x: {x:.3}, y: {y:.3}, xref: "paper", yref: "paper",
+               showarrow: false, font: {{ size: 12 }} }}"#,
+        ));
+    }
+    anns.join(",")
+}

tests/bohb_integration.rs

@@ -0,0 +1,133 @@
+//! Integration tests for the BOHB sampler.
+
+#![allow(
+    clippy::cast_sign_loss,
+    clippy::cast_precision_loss,
+    clippy::cast_possible_truncation
+)]
+
+use optimizer::parameter::{FloatParam, Parameter};
+use optimizer::sampler::bohb::BohbSampler;
+use optimizer::{Direction, Error, Study, TrialPruned};
+
+#[test]
+fn bohb_converges_on_quadratic() {
+    let bohb = BohbSampler::builder()
+        .min_resource(1)
+        .max_resource(9)
+        .reduction_factor(3)
+        .min_points_in_model(5)
+        .seed(42)
+        .build()
+        .unwrap();
+
+    let pruner = bohb.matching_pruner(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler_and_pruner(Direction::Minimize, bohb, pruner);
+
+    let x_param = FloatParam::new(-10.0, 10.0);
+
+    study
+        .optimize(60, |trial| {
+            let x = x_param.suggest(trial)?;
+
+            // Report intermediate values at budget steps 1, 3, 9
+            let obj = (x - 3.0).powi(2);
+            // Simulate budget-based evaluation with noise decreasing at higher budgets
+            trial.report(1, obj + 5.0);
+            trial.report(3, obj + 1.0);
+            trial.report(9, obj);
+
+            Ok::<_, Error>(obj)
+        })
+        .expect("optimization should succeed");
+
+    let best = study.best_trial().expect("should have trials");
+    assert!(
+        best.value < 10.0,
+        "BOHB should find a reasonable solution, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn bohb_with_pruning() {
+    let bohb = BohbSampler::builder()
+        .min_resource(1)
+        .max_resource(27)
+        .reduction_factor(3)
+        .min_points_in_model(3)
+        .seed(123)
+        .build()
+        .unwrap();
+
+    let pruner = bohb.matching_pruner(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler_and_pruner(Direction::Minimize, bohb, pruner);
+
+    let x_param = FloatParam::new(-5.0, 5.0);
+
+    study
+        .optimize(40, |trial| {
+            let x = x_param.suggest(trial)?;
+            let obj = x * x;
+
+            // Report at each rung step and check for pruning
+            for &step in &[1u64, 3, 9, 27] {
+                let noisy_obj = obj + 10.0 / step as f64;
+                trial.report(step, noisy_obj);
+
+                if trial.should_prune() {
+                    return Err(TrialPruned.into());
+                }
+            }
+
+            Ok::<_, Error>(obj)
+        })
+        .expect("optimization should succeed");
+
+    // Verify we have completed trials
+    let best = study.best_trial().expect("should have at least one trial");
+    assert!(
+        best.value < 25.0,
+        "best value {} should be reasonable",
+        best.value
+    );
+}
+
+#[test]
+fn bohb_uses_budget_conditioned_history() {
+    // Verify that BOHB conditions on budget level by testing that samples
+    // are influenced by intermediate values, not just final values.
+    let bohb = BohbSampler::builder()
+        .min_resource(1)
+        .max_resource(9)
+        .reduction_factor(3)
+        .min_points_in_model(3)
+        .seed(42)
+        .build()
+        .unwrap();
+
+    let pruner = bohb.matching_pruner(Direction::Minimize);
+    let study: Study<f64> = Study::with_sampler_and_pruner(Direction::Minimize, bohb, pruner);
+
+    let x_param = FloatParam::new(0.0, 10.0);
+
+    study
+        .optimize(30, |trial| {
+            let x = x_param.suggest(trial)?;
+            // Intermediate values that guide optimization toward x=2
+            trial.report(1, (x - 2.0).powi(2) + 1.0);
+            trial.report(3, (x - 2.0).powi(2) + 0.5);
+            trial.report(9, (x - 2.0).powi(2));
+
+            Ok::<_, Error>((x - 2.0).powi(2))
+        })
+        .expect("optimization should succeed");
+
+    let best = study.best_trial().unwrap();
+    let best_x: f64 = best.get(&x_param).unwrap();
+    // Should find x reasonably close to 2.0
+    assert!(
+        (best_x - 2.0).abs() < 5.0,
+        "BOHB should explore near x=2, got x={best_x}"
+    );
+}

tests/cma_es_tests.rs

@@ -0,0 +1,259 @@
+#![cfg(feature = "cma-es")]
+
+use optimizer::prelude::*;
+use optimizer::sampler::cma_es::CmaEsSampler;
+
+#[test]
+fn sphere_function() {
+    let sampler = CmaEsSampler::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 = CmaEsSampler::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(300, |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();
+    // Rosenbrock minimum is 0 at (1, 1); we just check reasonable convergence
+    assert!(
+        best.value < 50.0,
+        "rosenbrock best value should be < 50.0, got {}",
+        best.value
+    );
+}
+
+#[test]
+fn bounds_respected() {
+    let sampler = CmaEsSampler::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 = CmaEsSampler::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();
+    // Should find a reasonable value
+    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 = CmaEsSampler::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 = CmaEsSampler::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 = CmaEsSampler::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 = CmaEsSampler::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 = CmaEsSampler::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_population_size_and_sigma() {
+    let sampler = CmaEsSampler::builder()
+        .sigma0(1.0)
+        .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 config optimization should work, got {}",
+        best.value
+    );
+}

tests/export_tests.rs

@@ -0,0 +1,230 @@
+use optimizer::parameter::{FloatParam, IntParam, Parameter};
+use optimizer::sampler::random::RandomSampler;
+use optimizer::{Direction, Study};
+
+#[test]
+fn csv_empty_study_produces_header_only() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let mut buf = Vec::new();
+    study.to_csv(&mut buf).unwrap();
+    let csv = String::from_utf8(buf).unwrap();
+    assert_eq!(csv, "trial_id,value,state\n");
+}
+
+#[test]
+fn csv_includes_all_trial_data() {
+    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 = IntParam::new(1, 5).name("y");
+
+    study
+        .optimize(3, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv + yv as f64)
+        })
+        .unwrap();
+
+    let mut buf = Vec::new();
+    study.to_csv(&mut buf).unwrap();
+    let csv = String::from_utf8(buf).unwrap();
+
+    let lines: Vec<&str> = csv.lines().collect();
+    // Header + 3 data rows.
+    assert_eq!(lines.len(), 4);
+
+    // Header should contain our parameter names.
+    let header = lines[0];
+    assert!(header.starts_with("trial_id,value,state"));
+    assert!(header.contains("x"));
+    assert!(header.contains("y"));
+
+    // Each data row should have the right number of columns.
+    let n_cols = header.split(',').count();
+    for line in &lines[1..] {
+        assert_eq!(line.split(',').count(), n_cols);
+    }
+
+    // All rows should have "Complete" state.
+    for line in &lines[1..] {
+        assert!(line.contains("Complete"));
+    }
+}
+
+#[test]
+fn csv_handles_pruned_trials() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    // First trial: complete
+    let mut trial = study.create_trial();
+    let _ = x.suggest(&mut trial).unwrap();
+    study.complete_trial(trial, 1.0);
+
+    // Second trial: pruned
+    let mut trial = study.create_trial();
+    let _ = x.suggest(&mut trial).unwrap();
+    study.prune_trial(trial);
+
+    let mut buf = Vec::new();
+    study.to_csv(&mut buf).unwrap();
+    let csv = String::from_utf8(buf).unwrap();
+
+    let lines: Vec<&str> = csv.lines().collect();
+    assert_eq!(lines.len(), 3); // header + 2 data rows
+
+    // Pruned trial should have empty value.
+    let pruned_line = lines[2];
+    assert!(pruned_line.contains("Pruned"));
+    // The value field (second column) should be empty.
+    let cols: Vec<&str> = pruned_line.split(',').collect();
+    assert_eq!(cols[2], "Pruned");
+    assert_eq!(cols[1], ""); // empty value for pruned
+}
+
+#[test]
+fn csv_handles_different_parameter_sets() {
+    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");
+
+    // First trial: only x
+    let mut trial = study.create_trial();
+    let xv = x.suggest(&mut trial).unwrap();
+    study.complete_trial(trial, xv);
+
+    // Second trial: only y
+    let mut trial = study.create_trial();
+    let yv = y.suggest(&mut trial).unwrap();
+    study.complete_trial(trial, yv);
+
+    let mut buf = Vec::new();
+    study.to_csv(&mut buf).unwrap();
+    let csv = String::from_utf8(buf).unwrap();
+
+    let lines: Vec<&str> = csv.lines().collect();
+    assert_eq!(lines.len(), 3);
+
+    // Both x and y columns should exist.
+    let header = lines[0];
+    assert!(header.contains("x"));
+    assert!(header.contains("y"));
+
+    // Each row has the right column count (missing params are empty).
+    let n_cols = header.split(',').count();
+    for line in &lines[1..] {
+        assert_eq!(line.split(',').count(), n_cols);
+    }
+}
+
+#[test]
+fn csv_output_is_parseable() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let lr = FloatParam::new(0.001, 0.1).name("learning_rate");
+    let layers = IntParam::new(1, 5).name("n_layers");
+
+    study
+        .optimize(5, |trial| {
+            let l = lr.suggest(trial)?;
+            let n = layers.suggest(trial)?;
+            Ok::<_, optimizer::Error>(l * n as f64)
+        })
+        .unwrap();
+
+    let mut buf = Vec::new();
+    study.to_csv(&mut buf).unwrap();
+    let csv = String::from_utf8(buf).unwrap();
+
+    // Parse each row: every value field should be a valid f64 for complete trials.
+    let lines: Vec<&str> = csv.lines().collect();
+    for line in &lines[1..] {
+        let cols: Vec<&str> = line.split(',').collect();
+        // trial_id should be a number
+        cols[0].parse::<u64>().unwrap();
+        // value should be parseable as f64
+        cols[1].parse::<f64>().unwrap();
+        // state should be a known value
+        assert!(["Complete", "Pruned", "Failed", "Running"].contains(&cols[2]));
+    }
+}
+
+#[test]
+fn export_csv_writes_file() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+    study
+        .optimize(3, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv * xv)
+        })
+        .unwrap();
+
+    let dir = std::env::temp_dir().join("optimizer_export_test");
+    std::fs::create_dir_all(&dir).unwrap();
+    let path = dir.join("test_export.csv");
+
+    study.export_csv(&path).unwrap();
+
+    let contents = std::fs::read_to_string(&path).unwrap();
+    assert!(contents.starts_with("trial_id,value,state"));
+    assert!(contents.lines().count() == 4); // header + 3 rows
+
+    // Clean up.
+    let _ = std::fs::remove_dir_all(&dir);
+}
+
+#[cfg(feature = "serde")]
+#[test]
+fn export_json_writes_file() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+    study
+        .optimize(3, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv * xv)
+        })
+        .unwrap();
+
+    let dir = std::env::temp_dir().join("optimizer_json_export_test");
+    std::fs::create_dir_all(&dir).unwrap();
+    let path = dir.join("test_export.json");
+
+    study.export_json(&path).unwrap();
+
+    let contents = std::fs::read_to_string(&path).unwrap();
+    let parsed: serde_json::Value = serde_json::from_str(&contents).unwrap();
+    let arr = parsed.as_array().unwrap();
+    assert_eq!(arr.len(), 3);
+
+    // Each entry should have the expected fields.
+    for entry in arr {
+        assert!(entry.get("id").is_some());
+        assert!(entry.get("value").is_some());
+        assert!(entry.get("state").is_some());
+        assert!(entry.get("params").is_some());
+    }
+
+    // Clean up.
+    let _ = std::fs::remove_dir_all(&dir);
+}
+
+#[test]
+fn csv_includes_user_attributes() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(2, |trial| {
+            let xv = x.suggest(trial)?;
+            trial.set_user_attr("training_time_secs", 45.2);
+            Ok::<_, optimizer::Error>(xv * xv)
+        })
+        .unwrap();
+
+    let mut buf = Vec::new();
+    study.to_csv(&mut buf).unwrap();
+    let csv = String::from_utf8(buf).unwrap();
+
+    let header = csv.lines().next().unwrap();
+    assert!(header.contains("training_time_secs"));
+}

tests/fanova_tests.rs

@@ -0,0 +1,88 @@
+//! Integration tests for fANOVA parameter importance.
+
+#![cfg(feature = "fanova")]
+
+use optimizer::prelude::*;
+
+#[test]
+fn fanova_dominant_parameter() {
+    // f(x, y) = x^2 — x should dominate
+    let x = FloatParam::new(0.0, 10.0).name("x");
+    let y = FloatParam::new(0.0, 10.0).name("y");
+
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    study
+        .optimize(50, |trial| {
+            let xv = x.suggest(trial)?;
+            let _yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * xv)
+        })
+        .unwrap();
+
+    let result = study.fanova().unwrap();
+    assert_eq!(result.main_effects[0].0, "x");
+    assert!(
+        result.main_effects[0].1 > 0.7,
+        "x importance = {}",
+        result.main_effects[0].1
+    );
+}
+
+#[test]
+fn fanova_interaction() {
+    // f(x, y) = x * y — both matter and interact
+    let x = FloatParam::new(0.0, 10.0).name("x");
+    let y = FloatParam::new(0.0, 10.0).name("y");
+
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(7));
+    study
+        .optimize(100, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(xv * yv)
+        })
+        .unwrap();
+
+    let config = FanovaConfig {
+        n_trees: 128,
+        ..FanovaConfig::default()
+    };
+    let result = study.fanova_with_config(&config).unwrap();
+
+    // Should detect interaction
+    assert!(
+        !result.interactions.is_empty(),
+        "should detect x*y interaction"
+    );
+}
+
+#[test]
+fn fanova_consistent_with_correlation() {
+    // f(x, y) = 3*x + 0.5*y — x should rank higher in both methods
+    let x = FloatParam::new(0.0, 10.0).name("x");
+    let y = FloatParam::new(0.0, 10.0).name("y");
+
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(99));
+    study
+        .optimize(80, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, Error>(3.0 * xv + 0.5 * yv)
+        })
+        .unwrap();
+
+    let corr = study.param_importance();
+    let fanova = study.fanova().unwrap();
+
+    // Both methods should rank x above y
+    assert_eq!(corr[0].0, "x", "correlation should rank x first");
+    assert_eq!(fanova.main_effects[0].0, "x", "fanova should rank x first");
+}
+
+#[test]
+fn fanova_too_few_trials() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    let result = study.fanova();
+    assert!(result.is_err(), "should error with no trials");
+}

tests/importance_tests.rs

@@ -0,0 +1,166 @@
+use optimizer::parameter::{CategoricalParam, FloatParam, IntParam, Parameter};
+use optimizer::{Direction, Study};
+
+#[test]
+fn known_perfect_correlation() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(0.0, 100.0).name("x");
+
+    // Objective = x, so perfect correlation.
+    for _ in 0..30 {
+        let mut trial = study.ask();
+        let xv = x.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(xv));
+    }
+
+    let importance = study.param_importance();
+    assert_eq!(importance.len(), 1);
+    assert_eq!(importance[0].0, "x");
+    assert!(
+        (importance[0].1 - 1.0).abs() < 1e-10,
+        "single param should be 1.0"
+    );
+}
+
+#[test]
+fn no_effect_parameter() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(0.0, 100.0).name("x");
+    let noise = FloatParam::new(0.0, 100.0).name("noise");
+
+    // Objective depends only on x; noise is unused in objective.
+    for _ in 0..50 {
+        let mut trial = study.ask();
+        let xv = x.suggest(&mut trial).unwrap();
+        let _nv = noise.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(xv));
+    }
+
+    let importance = study.param_importance();
+    assert_eq!(importance.len(), 2);
+    // x should have much higher importance than noise.
+    let x_score = importance.iter().find(|(l, _)| l == "x").unwrap().1;
+    let noise_score = importance.iter().find(|(l, _)| l == "noise").unwrap().1;
+    assert!(
+        x_score > noise_score,
+        "x ({x_score}) should outrank noise ({noise_score})"
+    );
+    // x should dominate
+    assert!(x_score > 0.7, "x importance {x_score} should be dominant");
+}
+
+#[test]
+fn multiple_parameters_varying_importance() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(0.0, 10.0).name("x");
+    let y = FloatParam::new(0.0, 10.0).name("y");
+
+    // Objective = 10*x + 0.01*y → x should be far more important.
+    for _ in 0..50 {
+        let mut trial = study.ask();
+        let xv = x.suggest(&mut trial).unwrap();
+        let yv = y.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(10.0 * xv + 0.01 * yv));
+    }
+
+    let importance = study.param_importance();
+    assert_eq!(importance.len(), 2);
+    assert_eq!(importance[0].0, "x", "x should rank first");
+    assert!(importance[0].1 > importance[1].1);
+}
+
+#[test]
+fn fewer_than_two_trials_returns_empty() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    // 0 trials
+    assert!(study.param_importance().is_empty());
+
+    // 1 trial
+    let x = FloatParam::new(0.0, 1.0).name("x");
+    let mut trial = study.ask();
+    let xv = x.suggest(&mut trial).unwrap();
+    study.tell(trial, Ok::<_, &str>(xv));
+    assert!(study.param_importance().is_empty());
+}
+
+#[test]
+fn int_parameter() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let n = IntParam::new(1, 100).name("n");
+
+    for _ in 0..30 {
+        let mut trial = study.ask();
+        let nv = n.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(nv as f64));
+    }
+
+    let importance = study.param_importance();
+    assert_eq!(importance.len(), 1);
+    assert_eq!(importance[0].0, "n");
+    assert!(importance[0].1 > 0.9);
+}
+
+#[test]
+fn categorical_parameter() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let cat = CategoricalParam::new(vec!["a", "b", "c"]).name("cat");
+    let x = FloatParam::new(0.0, 100.0).name("x");
+
+    // Objective depends only on x; categorical is random noise.
+    for _ in 0..50 {
+        let mut trial = study.ask();
+        let _c = cat.suggest(&mut trial).unwrap();
+        let xv = x.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(xv));
+    }
+
+    let importance = study.param_importance();
+    assert_eq!(importance.len(), 2);
+    let x_score = importance.iter().find(|(l, _)| l == "x").unwrap().1;
+    assert!(x_score > 0.5, "x should dominate over categorical noise");
+}
+
+#[test]
+fn normalization_sums_to_one() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    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");
+
+    for _ in 0..50 {
+        let mut trial = study.ask();
+        let xv = x.suggest(&mut trial).unwrap();
+        let yv = y.suggest(&mut trial).unwrap();
+        let zv = z.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(xv + 0.5 * yv + 0.1 * zv));
+    }
+
+    let importance = study.param_importance();
+    let sum: f64 = importance.iter().map(|(_, s)| *s).sum();
+    assert!(
+        (sum - 1.0).abs() < 1e-10,
+        "scores should sum to 1.0, got {sum}"
+    );
+}
+
+#[test]
+fn label_when_unnamed_uses_debug() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    // No .name() call → label defaults to Debug representation.
+    let x = FloatParam::new(0.0, 10.0);
+
+    for _ in 0..10 {
+        let mut trial = study.ask();
+        let xv = x.suggest(&mut trial).unwrap();
+        study.tell(trial, Ok::<_, &str>(xv));
+    }
+
+    let importance = study.param_importance();
+    assert_eq!(importance.len(), 1);
+    assert!(
+        importance[0].0.starts_with("FloatParam"),
+        "expected Debug label, got {:?}",
+        importance[0].0
+    );
+}

tests/integration.rs

@@ -1892,3 +1892,356 @@ fn test_enqueue_counted_in_n_trials() {
     // All 5 trials count, including the 2 enqueued ones
     assert_eq!(study.n_trials(), 5);
 }
+
+// =============================================================================
+// Test: Study summary and Display
+// =============================================================================
+
+#[test]
+fn test_summary_with_completed_trials() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(1));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(5, |trial| {
+            let val = x.suggest(trial)?;
+            Ok::<_, Error>(val * val)
+        })
+        .unwrap();
+
+    let summary = study.summary();
+    assert!(summary.contains("Minimize"));
+    assert!(summary.contains("5 trials"));
+    assert!(summary.contains("Best value:"));
+    assert!(summary.contains("x = "));
+}
+
+#[test]
+fn test_summary_no_completed_trials() {
+    let study: Study<f64> = Study::new(Direction::Maximize);
+    let summary = study.summary();
+    assert!(summary.contains("Maximize"));
+    assert!(summary.contains("0 trials"));
+    assert!(!summary.contains("Best value:"));
+}
+
+#[test]
+fn test_summary_with_pruned_trials() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(1));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    // Manually create some complete and pruned trials
+    for _ in 0..3 {
+        let mut trial = study.create_trial();
+        let val = x.suggest(&mut trial).unwrap();
+        study.complete_trial(trial, val);
+    }
+    for _ in 0..2 {
+        let mut trial = study.create_trial();
+        let _ = x.suggest(&mut trial).unwrap();
+        study.prune_trial(trial);
+    }
+
+    let summary = study.summary();
+    // Should show breakdown when there are pruned trials
+    if study.n_pruned_trials() > 0 {
+        assert!(summary.contains("complete"));
+        assert!(summary.contains("pruned"));
+    }
+}
+
+#[test]
+fn test_display_matches_summary() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(1));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(3, |trial| {
+            let val = x.suggest(trial)?;
+            Ok::<_, Error>(val)
+        })
+        .unwrap();
+
+    assert_eq!(format!("{study}"), study.summary());
+}
+
+// =============================================================================
+// Tests: optimize_with_retries
+// =============================================================================
+
+#[test]
+fn test_retries_successful_trials_not_retried() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x_param = FloatParam::new(0.0, 10.0);
+    let call_count = std::cell::Cell::new(0u32);
+
+    study
+        .optimize_with_retries(5, 3, |trial| {
+            let x = x_param.suggest(trial)?;
+            call_count.set(call_count.get() + 1);
+            Ok::<_, Error>(x * x)
+        })
+        .unwrap();
+
+    // All trials succeed on first try — exactly 5 calls
+    assert_eq!(call_count.get(), 5);
+    assert_eq!(study.n_trials(), 5);
+}
+
+#[test]
+fn test_retries_failed_trials_retried_up_to_max() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x_param = FloatParam::new(0.0, 10.0);
+    let call_count = std::cell::Cell::new(0u32);
+
+    let result = study.optimize_with_retries(1, 3, |trial| {
+        let _ = x_param.suggest(trial).unwrap();
+        call_count.set(call_count.get() + 1);
+        Err::<f64, _>("always fails")
+    });
+
+    // 1 initial attempt + 3 retries = 4 total calls
+    assert_eq!(call_count.get(), 4);
+    // No trials completed
+    assert!(matches!(result, Err(Error::NoCompletedTrials)));
+}
+
+#[test]
+fn test_retries_permanently_failed_after_exhaustion() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x_param = FloatParam::new(0.0, 10.0);
+
+    let result = study.optimize_with_retries(3, 2, |trial| {
+        let _ = x_param.suggest(trial).unwrap();
+        Err::<f64, _>("transient error")
+    });
+
+    assert!(
+        matches!(result, Err(Error::NoCompletedTrials)),
+        "all trials should permanently fail"
+    );
+    assert_eq!(
+        study.n_trials(),
+        0,
+        "no completed trials should be recorded"
+    );
+}
+
+#[test]
+fn test_retries_uses_same_parameters() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x_param = FloatParam::new(0.0, 10.0);
+    let seen_values = std::cell::RefCell::new(Vec::new());
+    let call_count = std::cell::Cell::new(0u32);
+
+    study
+        .optimize_with_retries(1, 2, |trial| {
+            let x = x_param.suggest(trial).map_err(|e| e.to_string())?;
+            seen_values.borrow_mut().push(x);
+            call_count.set(call_count.get() + 1);
+            // Fail first two attempts, succeed on third
+            if call_count.get() < 3 {
+                Err::<f64, _>("transient".to_string())
+            } else {
+                Ok(x * x)
+            }
+        })
+        .unwrap();
+
+    let values = seen_values.borrow();
+    assert_eq!(values.len(), 3, "should be called 3 times (1 + 2 retries)");
+    // All three calls should have gotten the same parameter value
+    assert_eq!(values[0], values[1]);
+    assert_eq!(values[1], values[2]);
+}
+
+#[test]
+fn test_retries_n_trials_counts_unique_configs() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x_param = FloatParam::new(0.0, 10.0);
+    let call_count = std::cell::Cell::new(0u32);
+
+    study
+        .optimize_with_retries(3, 2, |trial| {
+            let x = x_param.suggest(trial).map_err(|e| e.to_string())?;
+            call_count.set(call_count.get() + 1);
+            // Fail first attempt of each config, succeed on retry
+            if call_count.get() % 2 == 1 {
+                Err::<f64, _>("transient".to_string())
+            } else {
+                Ok(x * x)
+            }
+        })
+        .unwrap();
+
+    // 3 unique configs, each needing 2 calls = 6 total calls
+    assert_eq!(call_count.get(), 6);
+    // But only 3 completed trials
+    assert_eq!(study.n_trials(), 3);
+}
+
+#[test]
+fn test_retries_with_zero_max_retries_same_as_optimize() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x_param = FloatParam::new(0.0, 10.0);
+    let call_count = std::cell::Cell::new(0u32);
+
+    study
+        .optimize_with_retries(5, 0, |trial| {
+            let x = x_param.suggest(trial)?;
+            call_count.set(call_count.get() + 1);
+            Ok::<_, Error>(x * x)
+        })
+        .unwrap();
+
+    assert_eq!(call_count.get(), 5);
+    assert_eq!(study.n_trials(), 5);
+}
+
+// =============================================================================
+// Tests: IntoIterator for &Study
+// =============================================================================
+
+#[test]
+fn test_into_iterator_iterates_all_trials() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x_param = FloatParam::new(0.0, 10.0);
+
+    for _ in 0..5 {
+        let mut trial = study.create_trial();
+        let x = x_param.suggest(&mut trial).unwrap();
+        study.complete_trial(trial, x * x);
+    }
+
+    let mut count = 0;
+    for trial in &study {
+        assert_eq!(trial.state, optimizer::TrialState::Complete);
+        count += 1;
+    }
+    assert_eq!(count, 5);
+}
+
+#[test]
+fn test_into_iterator_empty_study() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    let count = (&study).into_iter().count();
+    assert_eq!(count, 0);
+}
+
+#[test]
+fn test_into_iterator_preserves_insertion_order() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    for i in 0..3 {
+        let trial = study.create_trial();
+        study.complete_trial(trial, f64::from(i));
+    }
+
+    let ids: Vec<u64> = (&study).into_iter().map(|t| t.id).collect();
+    assert_eq!(ids, vec![0, 1, 2]);
+}
+
+// =============================================================================
+// Tests: Constraint handling
+// =============================================================================
+
+#[test]
+fn test_is_feasible_all_satisfied() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let mut trial = study.create_trial();
+    trial.set_constraints(vec![-1.0, 0.0, -0.5]);
+    study.complete_trial(trial, 1.0);
+
+    let completed = study.best_trial().unwrap();
+    assert!(completed.is_feasible());
+}
+
+#[test]
+fn test_is_feasible_one_violated() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let mut trial = study.create_trial();
+    trial.set_constraints(vec![-1.0, 0.5, -0.5]);
+    study.complete_trial(trial, 1.0);
+
+    let completed = study.best_trial().unwrap();
+    assert!(!completed.is_feasible());
+}
+
+#[test]
+fn test_is_feasible_empty_constraints() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let trial = study.create_trial();
+    study.complete_trial(trial, 1.0);
+
+    let completed = study.best_trial().unwrap();
+    assert!(completed.is_feasible());
+}
+
+#[test]
+fn test_best_trial_prefers_feasible() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    // Infeasible trial with better objective
+    let mut trial1 = study.create_trial();
+    trial1.set_constraints(vec![1.0]);
+    study.complete_trial(trial1, 0.1);
+
+    // Feasible trial with worse objective
+    let mut trial2 = study.create_trial();
+    trial2.set_constraints(vec![-1.0]);
+    study.complete_trial(trial2, 100.0);
+
+    let best = study.best_trial().unwrap();
+    assert_eq!(best.id, 1); // feasible trial wins
+    assert_eq!(best.value, 100.0);
+}
+
+#[test]
+fn test_best_trial_feasible_by_objective() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    // Feasible, worse objective
+    let mut trial1 = study.create_trial();
+    trial1.set_constraints(vec![-1.0]);
+    study.complete_trial(trial1, 10.0);
+
+    // Feasible, better objective
+    let mut trial2 = study.create_trial();
+    trial2.set_constraints(vec![-0.5]);
+    study.complete_trial(trial2, 2.0);
+
+    let best = study.best_trial().unwrap();
+    assert_eq!(best.id, 1); // lower objective wins among feasible
+    assert_eq!(best.value, 2.0);
+}
+
+#[test]
+fn test_top_trials_ranks_feasible_above_infeasible() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    // Infeasible, low violation
+    let mut t0 = study.create_trial();
+    t0.set_constraints(vec![0.5]);
+    study.complete_trial(t0, 1.0);
+
+    // Feasible, worst objective among feasible
+    let mut t1 = study.create_trial();
+    t1.set_constraints(vec![-1.0]);
+    study.complete_trial(t1, 50.0);
+
+    // Feasible, best objective among feasible
+    let mut t2 = study.create_trial();
+    t2.set_constraints(vec![-0.1]);
+    study.complete_trial(t2, 5.0);
+
+    // Infeasible, high violation
+    let mut t3 = study.create_trial();
+    t3.set_constraints(vec![3.0]);
+    study.complete_trial(t3, 0.5);
+
+    let top = study.top_trials(4);
+    let ids: Vec<u64> = top.iter().map(|t| t.id).collect();
+    // 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]);
+}

tests/multi_objective_tests.rs

@@ -0,0 +1,393 @@
+//! Integration tests for multi-objective optimization.
+
+use optimizer::Direction;
+use optimizer::multi_objective::MultiObjectiveStudy;
+use optimizer::parameter::{CategoricalParam, FloatParam, Parameter};
+use optimizer::sampler::nsga2::Nsga2Sampler;
+
+// ---------------------------------------------------------------------------
+// Pareto utility tests (via public MultiObjectiveStudy)
+// ---------------------------------------------------------------------------
+
+#[test]
+fn test_basic_two_objective_random() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    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();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty(), "Pareto front should be non-empty");
+
+    // Verify no solution in the front dominates another
+    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_dimension_mismatch_error() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    let x = FloatParam::new(0.0, 1.0);
+
+    let result = study.optimize(1, |trial| {
+        let xv = x.suggest(trial)?;
+        // Return wrong number of values
+        Ok::<_, optimizer::Error>(vec![xv])
+    });
+
+    assert!(result.is_err());
+    let err = result.unwrap_err();
+    assert!(
+        matches!(
+            err,
+            optimizer::Error::ObjectiveDimensionMismatch {
+                expected: 2,
+                got: 1
+            }
+        ),
+        "Expected ObjectiveDimensionMismatch, got: {err}"
+    );
+}
+
+#[test]
+fn test_ask_tell() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Maximize]);
+    let x = FloatParam::new(0.0, 10.0);
+
+    for _ in 0..10 {
+        let mut trial = study.ask();
+        let xv = x.suggest(&mut trial).unwrap();
+        study
+            .tell(trial, Ok::<_, &str>(vec![xv, 10.0 - xv]))
+            .unwrap();
+    }
+
+    assert_eq!(study.n_trials(), 10);
+    let front = study.pareto_front();
+    assert!(!front.is_empty());
+}
+
+#[test]
+fn test_ask_tell_dimension_mismatch() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    let trial = study.ask();
+    let result = study.tell(trial, Ok::<_, &str>(vec![1.0, 2.0, 3.0]));
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_n_trials_counting() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    assert_eq!(study.n_trials(), 0);
+
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(5, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    assert_eq!(study.n_trials(), 5);
+}
+
+#[test]
+fn test_three_objectives() {
+    let study = MultiObjectiveStudy::new(vec![
+        Direction::Minimize,
+        Direction::Minimize,
+        Direction::Maximize,
+    ]);
+    let x = FloatParam::new(0.0, 1.0);
+    let y = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(30, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, optimizer::Error>(vec![xv, yv, 1.0 - xv - yv])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    assert!(!front.is_empty());
+    assert_eq!(study.n_objectives(), 3);
+}
+
+#[test]
+fn test_directions_accessor() {
+    let dirs = vec![Direction::Minimize, Direction::Maximize];
+    let study = MultiObjectiveStudy::new(dirs.clone());
+    assert_eq!(study.directions(), &dirs);
+    assert_eq!(study.n_objectives(), 2);
+}
+
+#[test]
+fn test_trials_accessor() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(3, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    let trials = study.trials();
+    assert_eq!(trials.len(), 3);
+    for t in &trials {
+        assert_eq!(t.values.len(), 2);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// NSGA-II sampler tests
+// ---------------------------------------------------------------------------
+
+#[test]
+fn test_nsga2_zdt1() {
+    // ZDT1 benchmark: minimize both objectives
+    let n_vars = 5;
+    let params: Vec<FloatParam> = (0..n_vars).map(|_| FloatParam::new(0.0, 1.0)).collect();
+
+    let sampler = Nsga2Sampler::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(), "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_nsga2_with_seed_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 = Nsga2Sampler::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_nsga2_builder() {
+    let sampler = Nsga2Sampler::builder()
+        .population_size(10)
+        .crossover_prob(0.8)
+        .crossover_eta(15.0)
+        .mutation_eta(25.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_nsga2_categorical_params() {
+    let sampler = Nsga2Sampler::with_seed(42);
+    let study =
+        MultiObjectiveStudy::with_sampler(vec![Direction::Minimize, Direction::Minimize], sampler);
+
+    let x = FloatParam::new(0.0, 1.0);
+    let cat = CategoricalParam::new(vec!["a", "b", "c"]);
+
+    study
+        .optimize(30, |trial| {
+            let xv = x.suggest(trial)?;
+            let cv = cat.suggest(trial)?;
+            let bonus = match cv {
+                "a" => 0.0,
+                "b" => 0.5,
+                _ => 1.0,
+            };
+            Ok::<_, optimizer::Error>(vec![xv + bonus, 1.0 - xv])
+        })
+        .unwrap();
+
+    assert_eq!(study.n_trials(), 30);
+    let front = study.pareto_front();
+    assert!(!front.is_empty());
+}
+
+#[test]
+fn test_nsga2_constraints() {
+    let sampler = Nsga2Sampler::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)?;
+            // Constraint: x >= 0.3 (i.e. 0.3 - x <= 0)
+            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());
+
+    // Check that feasible solutions exist on the front
+    let feasible_count = front.iter().filter(|t| t.is_feasible()).count();
+    assert!(
+        feasible_count > 0,
+        "Should have feasible solutions on front"
+    );
+}
+
+#[test]
+fn test_multi_objective_trial_get() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(5, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(vec![xv, 10.0 - xv])
+        })
+        .unwrap();
+
+    let front = study.pareto_front();
+    for t in &front {
+        let xv: f64 = t.get(&x).unwrap();
+        assert!((0.0..=10.0).contains(&xv));
+    }
+}
+
+#[test]
+fn test_multi_objective_trial_is_feasible() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(10, |trial| {
+            let xv = x.suggest(trial)?;
+            trial.set_constraints(vec![0.5 - xv]); // feasible if x >= 0.5
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    let trials = study.trials();
+    for t in &trials {
+        let xv = t.values[0];
+        if xv >= 0.5 {
+            assert!(t.is_feasible());
+        } else {
+            assert!(!t.is_feasible());
+        }
+    }
+}
+
+#[test]
+fn test_multi_objective_trial_user_attrs() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(3, |trial| {
+            let xv = x.suggest(trial)?;
+            trial.set_user_attr("iteration", 42_i64);
+            Ok::<_, optimizer::Error>(vec![xv, 1.0 - xv])
+        })
+        .unwrap();
+
+    let trials = study.trials();
+    for t in &trials {
+        assert!(t.user_attr("iteration").is_some());
+    }
+}
+
+#[test]
+fn test_tell_with_failure() {
+    let study = MultiObjectiveStudy::new(vec![Direction::Minimize, Direction::Minimize]);
+
+    let trial = study.ask();
+    study
+        .tell(trial, Err::<Vec<f64>, _>("evaluation failed"))
+        .unwrap();
+
+    // Failed trial not counted
+    assert_eq!(study.n_trials(), 0);
+}

tests/serde_tests.rs

@@ -0,0 +1,295 @@
+#![cfg(feature = "serde")]
+
+use std::collections::HashMap;
+
+use optimizer::parameter::{FloatParam, IntParam, Parameter};
+use optimizer::sampler::CompletedTrial;
+use optimizer::{Direction, ParamValue, Study, StudySnapshot, TrialState};
+
+#[test]
+fn round_trip_save_load() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(-10.0, 10.0).name("x");
+    let n = IntParam::new(1, 100).name("n");
+
+    study
+        .optimize(5, |trial| {
+            let x_val = x.suggest(trial)?;
+            let n_val = n.suggest(trial)?;
+            Ok::<_, optimizer::Error>(x_val * x_val + n_val as f64)
+        })
+        .unwrap();
+
+    let dir = tempdir();
+    let path = dir.join("study.json");
+
+    study.save(&path).unwrap();
+    let loaded: Study<f64> = Study::load(&path).unwrap();
+
+    assert_eq!(loaded.direction(), study.direction());
+    assert_eq!(loaded.n_trials(), study.n_trials());
+
+    let orig_trials = study.trials();
+    let loaded_trials = loaded.trials();
+
+    for (orig, loaded) in orig_trials.iter().zip(loaded_trials.iter()) {
+        assert_eq!(orig.id, loaded.id);
+        assert!((orig.value - loaded.value).abs() < 1e-10);
+        assert_eq!(orig.state, loaded.state);
+        assert_eq!(orig.params.len(), loaded.params.len());
+        assert_eq!(orig.distributions, loaded.distributions);
+        assert_eq!(orig.param_labels, loaded.param_labels);
+    }
+
+    // Clean up
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn json_output_is_human_readable() {
+    let study: Study<f64> = Study::new(Direction::Maximize);
+    let x = FloatParam::new(0.0, 1.0).name("x");
+
+    study
+        .optimize(2, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v)
+        })
+        .unwrap();
+
+    let dir = tempdir();
+    let path = dir.join("study.json");
+    study.save(&path).unwrap();
+
+    let contents = std::fs::read_to_string(&path).unwrap();
+
+    // Verify it's pretty-printed JSON with recognizable fields
+    assert!(contents.contains("\"version\""));
+    assert!(contents.contains("\"direction\""));
+    assert!(contents.contains("\"trials\""));
+    assert!(contents.contains("\"next_trial_id\""));
+    assert!(contents.contains("\"Maximize\""));
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn round_trip_empty_study() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    let dir = tempdir();
+    let path = dir.join("empty.json");
+
+    study.save(&path).unwrap();
+    let loaded: Study<f64> = Study::load(&path).unwrap();
+
+    assert_eq!(loaded.direction(), Direction::Minimize);
+    assert_eq!(loaded.n_trials(), 0);
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn snapshot_version_field_is_present() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    let dir = tempdir();
+    let path = dir.join("version.json");
+    study.save(&path).unwrap();
+
+    let contents = std::fs::read_to_string(&path).unwrap();
+    let snapshot: StudySnapshot<f64> = serde_json::from_str(&contents).unwrap();
+
+    assert_eq!(snapshot.version, 1);
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn completed_trial_serde_round_trip() {
+    let trial = CompletedTrial::new(42, HashMap::new(), HashMap::new(), HashMap::new(), 2.78);
+
+    let json = serde_json::to_string(&trial).unwrap();
+    let deserialized: CompletedTrial<f64> = serde_json::from_str(&json).unwrap();
+
+    assert_eq!(deserialized.id, 42);
+    assert_eq!(deserialized.value, 2.78);
+    assert_eq!(deserialized.state, TrialState::Complete);
+}
+
+#[test]
+fn param_value_serde_round_trip() {
+    let values = vec![
+        ParamValue::Float(1.23),
+        ParamValue::Int(42),
+        ParamValue::Categorical(2),
+    ];
+
+    for val in &values {
+        let json = serde_json::to_string(val).unwrap();
+        let deserialized: ParamValue = serde_json::from_str(&json).unwrap();
+        assert_eq!(&deserialized, val);
+    }
+}
+
+#[test]
+fn direction_serde_round_trip() {
+    let min_json = serde_json::to_string(&Direction::Minimize).unwrap();
+    let max_json = serde_json::to_string(&Direction::Maximize).unwrap();
+
+    assert_eq!(
+        serde_json::from_str::<Direction>(&min_json).unwrap(),
+        Direction::Minimize
+    );
+    assert_eq!(
+        serde_json::from_str::<Direction>(&max_json).unwrap(),
+        Direction::Maximize
+    );
+}
+
+#[test]
+fn round_trip_preserves_trial_id_counter() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(0.0, 1.0);
+
+    study
+        .optimize(10, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v)
+        })
+        .unwrap();
+
+    let dir = tempdir();
+    let path = dir.join("counter.json");
+    study.save(&path).unwrap();
+
+    let loaded: Study<f64> = Study::load(&path).unwrap();
+
+    // Creating a new trial should use an ID >= 10
+    let trial = loaded.create_trial();
+    assert!(trial.id() >= 10);
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn checkpoint_file_created_at_interval() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(-10.0, 10.0).name("x");
+
+    let dir = tempdir();
+    let checkpoint = dir.join("checkpoint.json");
+
+    study
+        .optimize_with_checkpoint(10, 3, &checkpoint, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v * v)
+        })
+        .unwrap();
+
+    // Checkpoint should exist (written at trials 3, 6, 9)
+    assert!(checkpoint.exists(), "checkpoint file was not created");
+
+    // Load it and verify it's valid
+    let loaded: Study<f64> = Study::load(&checkpoint).unwrap();
+    // Last checkpoint was at trial 9, so it should have 9 trials
+    assert_eq!(loaded.n_trials(), 9);
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn checkpoint_overwrites_previous() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(0.0, 1.0);
+
+    let dir = tempdir();
+    let checkpoint = dir.join("checkpoint.json");
+
+    study
+        .optimize_with_checkpoint(6, 3, &checkpoint, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v)
+        })
+        .unwrap();
+
+    // The checkpoint at trial 6 should overwrite the one from trial 3
+    let loaded: Study<f64> = Study::load(&checkpoint).unwrap();
+    assert_eq!(loaded.n_trials(), 6);
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn resume_from_checkpoint_continues_trial_ids() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(-5.0, 5.0).name("x");
+
+    let dir = tempdir();
+    let checkpoint = dir.join("resume.json");
+
+    // Run 10 trials with checkpointing
+    study
+        .optimize_with_checkpoint(10, 5, &checkpoint, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v * v)
+        })
+        .unwrap();
+
+    // Load and continue
+    let loaded: Study<f64> = Study::load(&checkpoint).unwrap();
+    assert_eq!(loaded.n_trials(), 10);
+
+    let remaining = 15 - loaded.n_trials();
+    loaded
+        .optimize(remaining, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v * v)
+        })
+        .unwrap();
+
+    assert_eq!(loaded.n_trials(), 15);
+
+    // Verify no duplicate trial IDs
+    let trials = loaded.trials();
+    let mut ids: Vec<u64> = trials.iter().map(|t| t.id).collect();
+    ids.sort_unstable();
+    ids.dedup();
+    assert_eq!(ids.len(), 15, "duplicate trial IDs found");
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+#[test]
+fn atomic_write_no_temp_file_left_behind() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+    let x = FloatParam::new(0.0, 1.0);
+
+    let dir = tempdir();
+    let checkpoint = dir.join("atomic.json");
+
+    study
+        .optimize_with_checkpoint(3, 3, &checkpoint, |trial| {
+            let v = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(v)
+        })
+        .unwrap();
+
+    // The temp file should have been renamed, not left behind
+    let tmp_path = dir.join(".atomic.json.tmp");
+    assert!(!tmp_path.exists(), "temp file was not cleaned up");
+    assert!(checkpoint.exists(), "checkpoint file was not created");
+
+    std::fs::remove_dir_all(&dir).ok();
+}
+
+/// Helper to create a unique temporary directory.
+fn tempdir() -> std::path::PathBuf {
+    use std::sync::atomic::{AtomicU64, Ordering};
+    static COUNTER: AtomicU64 = AtomicU64::new(0);
+    let id = COUNTER.fetch_add(1, Ordering::Relaxed);
+    let dir =
+        std::env::temp_dir().join(format!("optimizer_serde_test_{}_{id}", std::process::id()));
+    std::fs::create_dir_all(&dir).unwrap();
+    dir
+}

tests/test_functions.rs

@@ -0,0 +1,45 @@
+#[path = "../benches/test_functions.rs"]
+mod test_functions;
+
+use test_functions::*;
+
+const TOL: f64 = 1e-10;
+
+#[test]
+fn sphere_at_optimum() {
+    assert!(sphere(&[0.0, 0.0]).abs() < TOL);
+    assert!(sphere(&[0.0; 10]).abs() < TOL);
+}
+
+#[test]
+fn rosenbrock_at_optimum() {
+    assert!(rosenbrock(&[1.0, 1.0]).abs() < TOL);
+    assert!(rosenbrock(&[1.0; 5]).abs() < TOL);
+}
+
+#[test]
+fn rastrigin_at_optimum() {
+    assert!(rastrigin(&[0.0, 0.0]).abs() < TOL);
+    assert!(rastrigin(&[0.0; 10]).abs() < TOL);
+}
+
+#[test]
+fn ackley_at_optimum() {
+    assert!(ackley(&[0.0, 0.0]).abs() < 1e-8);
+    assert!(ackley(&[0.0; 10]).abs() < 1e-8);
+}
+
+#[test]
+fn branin_at_optimum() {
+    let target = 0.397_887_357_729_738_1;
+    let val = branin(&[std::f64::consts::PI, 2.275]);
+    assert!((val - target).abs() < 1e-3);
+}
+
+#[test]
+fn hartmann6_at_optimum() {
+    let target = -3.3224;
+    let x_opt = [0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573];
+    let val = hartmann6(&x_opt);
+    assert!((val - target).abs() < 0.01);
+}

tests/visualization_tests.rs

@@ -0,0 +1,182 @@
+#![cfg(feature = "visualization")]
+
+use optimizer::parameter::{FloatParam, IntParam, Parameter};
+use optimizer::sampler::random::RandomSampler;
+use optimizer::{Direction, Study, generate_html_report};
+
+#[test]
+fn html_report_creates_file() {
+    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 = IntParam::new(1, 5).name("y");
+
+    study
+        .optimize(10, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv + yv as f64)
+        })
+        .unwrap();
+
+    let path = std::env::temp_dir().join("test_report_creates_file.html");
+    generate_html_report(&study, &path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+    assert!(content.contains("<!DOCTYPE html>"));
+    assert!(content.contains("plotly"));
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn html_report_contains_all_chart_sections() {
+    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(-5.0, 5.0).name("y");
+
+    study
+        .optimize(20, |trial| {
+            let xv = x.suggest(trial)?;
+            let yv = y.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv * xv + yv * yv)
+        })
+        .unwrap();
+
+    let path = std::env::temp_dir().join("test_report_all_charts.html");
+    generate_html_report(&study, &path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+
+    // Should contain all chart divs.
+    assert!(content.contains("id=\"history\""));
+    assert!(content.contains("id=\"slices\""));
+    assert!(content.contains("id=\"parcoords\""));
+    assert!(content.contains("id=\"importance\""));
+    assert!(content.contains("id=\"timeline\""));
+
+    // Should contain chart titles.
+    assert!(content.contains("Optimization History"));
+    assert!(content.contains("Slice Plots"));
+    assert!(content.contains("Parallel Coordinates"));
+    assert!(content.contains("Parameter Importance"));
+    assert!(content.contains("Trial Timeline"));
+
+    // Should show direction and trial count.
+    assert!(content.contains("Minimize"));
+    assert!(content.contains("20 trials"));
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn html_report_empty_study() {
+    let study: Study<f64> = Study::new(Direction::Minimize);
+
+    let path = std::env::temp_dir().join("test_report_empty.html");
+    generate_html_report(&study, &path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+    assert!(content.contains("<!DOCTYPE html>"));
+    assert!(content.contains("0 trials"));
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn html_report_single_param_no_parcoords() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(5, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv * xv)
+        })
+        .unwrap();
+
+    let path = std::env::temp_dir().join("test_report_single_param.html");
+    generate_html_report(&study, &path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+
+    // Should have slice plot but not parallel coordinates (needs >= 2 params).
+    assert!(content.contains("id=\"slices\""));
+    assert!(!content.contains("id=\"parcoords\""));
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn html_report_maximize_direction() {
+    let study: Study<f64> = Study::with_sampler(Direction::Maximize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(5, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv)
+        })
+        .unwrap();
+
+    let path = std::env::temp_dir().join("test_report_maximize.html");
+    generate_html_report(&study, &path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+    assert!(content.contains("Maximize"));
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn export_html_convenience_method() {
+    let study: Study<f64> = Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(5, |trial| {
+            let xv = x.suggest(trial)?;
+            Ok::<_, optimizer::Error>(xv * xv)
+        })
+        .unwrap();
+
+    let path = std::env::temp_dir().join("test_export_html.html");
+    study.export_html(&path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+    assert!(content.contains("<!DOCTYPE html>"));
+    assert!(content.contains("id=\"history\""));
+
+    std::fs::remove_file(&path).ok();
+}
+
+#[test]
+fn html_report_with_intermediate_values() {
+    use optimizer::pruner::MedianPruner;
+
+    let mut study: Study<f64> =
+        Study::with_sampler(Direction::Minimize, RandomSampler::with_seed(42));
+    study.set_pruner(MedianPruner::new(Direction::Minimize));
+    let x = FloatParam::new(0.0, 10.0).name("x");
+
+    study
+        .optimize(10, |trial| {
+            let xv = x.suggest(trial)?;
+            for step in 0..5 {
+                let val = xv * xv + step as f64;
+                trial.report(step, val);
+                if trial.should_prune() {
+                    return Err(optimizer::TrialPruned.into());
+                }
+            }
+            Ok::<_, optimizer::Error>(xv * xv)
+        })
+        .unwrap();
+
+    let path = std::env::temp_dir().join("test_report_intermediate.html");
+    generate_html_report(&study, &path).unwrap();
+
+    let content = std::fs::read_to_string(&path).unwrap();
+    assert!(content.contains("id=\"intermediate\""));
+    assert!(content.contains("Intermediate Values"));
+
+    std::fs::remove_file(&path).ok();
+}