diff --git a/README.md b/README.md
index 1d60eae..d0765c8 100644
--- a/README.md
+++ b/README.md
@@ -91,8 +91,8 @@ study.run(n_jobs=-1)
 
 # 4. Inspect results
 print(study.summary())
-front = study.front()  # non-dominated configs
-hv = study.front_hypervolume(ref=...)  # hypervolume indicator
+front = study.front("benchmark")  # non-dominated config indices
+hv = study.front_hypervolume("benchmark", ref_point)  # hypervolume indicator
 ```
 
 ### Protocols
@@ -156,7 +156,7 @@ factors = [
 ]
 
 grid = build_grid(factors, method="sobol", n_samples=1024)
-si = screen(run_fn, factors, method="morris", n_eval=3000)
+si = screen(run_fn, factors, method="morris", n_trajectories=3000)
 ```
 
 ### Execution
@@ -176,9 +176,9 @@ results = run_adaptive(world, scorer, factors, observables, n_trials=600)
 ```python
 from trade_study import extract_front, hypervolume, pareto_rank
 
-front_mask = extract_front(results.scores, directions)
+front_idx = extract_front(results.scores, directions)
 ranks = pareto_rank(results.scores, directions)
-hv = hypervolume(results.scores[front_mask], directions, ref=ref_point)
+hv = hypervolume(results.scores[front_idx], ref_point, directions)
 ```
 
 ### Stacking
diff --git a/docs/api/study.md b/docs/api/study.md
index 02975fc..bb4afbc 100644
--- a/docs/api/study.md
+++ b/docs/api/study.md
@@ -7,3 +7,7 @@ Multi-phase study orchestration.
 ::: trade_study.Study
 
 ::: trade_study.top_k_pareto_filter
+
+::: trade_study.weighted_sum_filter
+
+::: trade_study.feasibility_filter
diff --git a/docs/api/viz.md b/docs/api/viz.md
new file mode 100644
index 0000000..d0fdc99
--- /dev/null
+++ b/docs/api/viz.md
@@ -0,0 +1,11 @@
+# Visualization
+
+Plotting utilities for trade-study results.
+
+::: trade_study.plot_front
+
+::: trade_study.plot_parallel
+
+::: trade_study.plot_scores
+
+::: trade_study.plot_calibration
diff --git a/docs/assets/cstr_front.png b/docs/assets/cstr_front.png
index 1a04ace..482c2de 100644
Binary files a/docs/assets/cstr_front.png and b/docs/assets/cstr_front.png differ
diff --git a/docs/assets/cstr_parallel.png b/docs/assets/cstr_parallel.png
index 22d62a2..69ce15d 100644
Binary files a/docs/assets/cstr_parallel.png and b/docs/assets/cstr_parallel.png differ
diff --git a/docs/assets/cstr_selectivity.png b/docs/assets/cstr_selectivity.png
index be1d598..4542837 100644
Binary files a/docs/assets/cstr_selectivity.png and b/docs/assets/cstr_selectivity.png differ
diff --git a/docs/guide/bayesian.md b/docs/guide/bayesian.md
index 57fc472..313e0cd 100644
--- a/docs/guide/bayesian.md
+++ b/docs/guide/bayesian.md
@@ -79,9 +79,12 @@ metrics:
 
 - **CRPS** (Continuous Ranked Probability Score) — a proper scoring
   rule that rewards sharp, well-calibrated predictive distributions.
+- **Energy score** — a multivariate proper scoring rule
+  complementary to CRPS.
 - **Coverage (95 %)** — the fraction of test points whose true value
   falls inside the 95 % posterior predictive interval.
 - **RMSE** — root mean squared error of the posterior mean.
+- **MAE** — mean absolute error of the posterior mean.
 
 ```python
 --8<-- "examples/bayesian_study.py:world"
@@ -110,6 +113,16 @@ cost model for the number of observations:
 --8<-- "examples/bayesian_study.py:annotation"
 ```
 
+## Constraints
+
+A `Constraint` defines a hard feasibility bound on an observable.
+`feasibility_filter` builds a phase filter that drops any design
+violating the constraints:
+
+```python
+--8<-- "examples/bayesian_study.py:constraint"
+```
+
 ## Factor screening
 
 Before running the full grid, Morris screening identifies which
@@ -157,6 +170,56 @@ overconfident nor excessively vague.
 
 ![CRPS strip plot](../assets/bayesian_crps.png)
 
+## Quality indicators
+
+`hypervolume` measures the dominated volume of the Pareto front;
+`igd_plus` (inverted generational distance plus) measures how close
+the front is to a reference set.  Together they quantify front
+quality — useful for comparing design methods or grid sizes:
+
+```python
+--8<-- "examples/bayesian_study.py:igd_plus"
+```
+
+## Feasibility filtering
+
+Apply `Constraint` + `feasibility_filter` to keep only designs
+meeting hard bounds (here, `coverage_95 >= 0.90`):
+
+```python
+--8<-- "examples/bayesian_study.py:feasibility"
+```
+
+## Weighted-sum filtering
+
+`weighted_sum_filter` scalarises objectives into a single score
+via min-max normalisation and weighted combination, then keeps
+the top-K designs.  Useful when you want a simple ranking
+without full Pareto analysis:
+
+```python
+--8<-- "examples/bayesian_study.py:weighted_sum"
+```
+
+## Sobol grid
+
+`build_grid(method="sobol")` generates quasi-random points with
+better space-filling properties than LHS for low dimensions:
+
+```python
+--8<-- "examples/bayesian_study.py:sobol"
+```
+
+## Study workflow
+
+The `Study` class orchestrates multi-phase workflows.
+`Study.summary()` returns per-phase statistics and
+`Study.stack()` computes stacking weights directly:
+
+```python
+--8<-- "examples/bayesian_study.py:study_workflow"
+```
+
 ## Score-based stacking
 
 `stack_scores` finds simplex-constrained weights that minimise the
@@ -222,9 +285,9 @@ rather than a design factor:
 
 ## What to try next
 
-- Swap `method="morris"` for `method="sobol"` in `screen()` for
-  variance-based sensitivity indices.
-- Use `Constraint` + `feasibility_filter` to enforce
-  `coverage_95 >= 0.90` before stacking.
 - Try `stack_bayesian()` on models that expose log-likelihood
   (requires `arviz`).
+- Add more phases to the `Study` — e.g. a third phase that
+  re-evaluates with even higher fidelity.
+- Experiment with different `weighted_sum_filter` weight vectors
+  to explore different trade-off preferences.
diff --git a/docs/guide/cstr.md b/docs/guide/cstr.md
index c13e720..7f40530 100644
--- a/docs/guide/cstr.md
+++ b/docs/guide/cstr.md
@@ -99,6 +99,25 @@ operating parameters with physical bounds:
 --8<-- "examples/cstr_study.py:factors"
 ```
 
+## Constraints
+
+Hard feasibility bounds can be defined with `Constraint` objects.
+`feasibility_filter` builds a phase filter that drops designs
+violating any constraint:
+
+```python
+--8<-- "examples/cstr_study.py:constraint"
+```
+
+## Factor screening
+
+Before committing to a full grid sweep, Sobol screening identifies
+which factors influence the objectives most:
+
+```python
+--8<-- "examples/cstr_study.py:screening"
+```
+
 ## Build the study phases
 
 A `Study` chains multiple **Phases**.  Phase 1 explores the design space
@@ -154,9 +173,37 @@ conversion–selectivity trade-off.
 
 ## What to try next
 
-- Add a **screening** step with `screen()` to identify which factors
-  matter most before running the full grid.
-- Use `run_adaptive()` for Bayesian optimization instead of LHS.
 - Call `stack_bayesian()` on the Pareto front to compute model-averaging
   weights.
 - Save results with `save_results()` for later analysis.
+- Combine `feasibility_filter` with `top_k_pareto_filter` across
+  successive phases for constrained multi-objective optimisation.
+
+## Progress callbacks and parallel execution
+
+`run_grid` accepts a `callback` for progress reporting and `n_jobs`
+for parallel execution via joblib:
+
+```python
+--8<-- "examples/cstr_study.py:callback"
+```
+
+## Adaptive exploration
+
+When the design space is too large for a grid sweep, use
+`Phase(grid="adaptive")` for optuna-driven multi-objective
+Bayesian optimisation.  The `factors` argument on `Study` provides
+the parameter bounds:
+
+```python
+--8<-- "examples/cstr_study.py:adaptive"
+```
+
+## Feasibility-constrained phases
+
+Use `feasibility_filter` as a phase filter to enforce hard
+constraints before passing designs to subsequent phases:
+
+```python
+--8<-- "examples/cstr_study.py:feasibility"
+```
diff --git a/docs/index.md b/docs/index.md
index 42892d2..3d3fe63 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -29,8 +29,9 @@ design-of-experiments studies:
 | [Protocols](api/protocols.md) | Core types: `Observable`, `Direction`, `Scorer`, `Simulator`, etc. |
 | [Design](api/design.md) | `Factor`, `build_grid`, `screen`, `reduce_factors` |
 | [Runner](api/runner.md) | `run_grid`, `run_adaptive` |
-| [Study](api/study.md) | `Phase`, `Study`, `top_k_pareto_filter` |
+| [Study](api/study.md) | `Phase`, `Study`, `top_k_pareto_filter`, `weighted_sum_filter`, `feasibility_filter` |
 | [Scoring](api/scoring.md) | `score`, `coverage_curve` |
 | [Pareto](api/pareto.md) | `extract_front`, `pareto_rank`, `hypervolume`, `igd_plus` |
 | [Stacking](api/stacking.md) | `stack_bayesian`, `stack_scores`, `ensemble_predict` |
+| [Visualization](api/viz.md) | `plot_front`, `plot_parallel`, `plot_scores`, `plot_calibration` |
 | [I/O](api/io.md) | `load_results`, `save_results` |
diff --git a/examples/bayesian_study.py b/examples/bayesian_study.py
index 2dda582..52f8868 100644
--- a/examples/bayesian_study.py
+++ b/examples/bayesian_study.py
@@ -17,6 +17,7 @@
 
 from trade_study import (
     Annotation,
+    Constraint,
     Direction,
     Factor,
     FactorType,
@@ -27,6 +28,9 @@
     coverage_curve,
     ensemble_predict,
     extract_front,
+    feasibility_filter,
+    hypervolume,
+    igd_plus,
     load_results,
     plot_calibration,
     plot_front,
@@ -39,6 +43,7 @@
     screen,
     stack_scores,
     top_k_pareto_filter,
+    weighted_sum_filter,
 )
 
 ASSET_DIR = "docs/assets"
@@ -191,7 +196,7 @@ def score(
         observations: Any,
         config: dict[str, Any],
     ) -> dict[str, float]:
-        """Compute CRPS, 95 percent coverage, and RMSE on posterior mean.
+        """Compute CRPS, energy, 95 percent coverage, RMSE, and MAE.
 
         Args:
             truth: True test-set values (y_test).
@@ -199,12 +204,20 @@ def score(
             config: Hyperparameter dict (unused).
 
         Returns:
-            Scores for crps, coverage_95, and rmse.
+            Scores for crps, energy, coverage_95, rmse, and mae.
         """
-        crps = score("crps", observations["pred_samples"], truth)
+        crps_val = score("crps", observations["pred_samples"], truth)
+        energy_val = score("energy", observations["pred_samples"], truth)
         cov95 = score("coverage", observations["pred_samples"], truth, level=0.95)
-        rmse = score("rmse", observations["pred_mean"], truth)
-        return {"crps": crps, "coverage_95": cov95, "rmse": rmse}
+        rmse_val = score("rmse", observations["pred_mean"], truth)
+        mae_val = score("mae", observations["pred_mean"], truth)
+        return {
+            "crps": crps_val,
+            "energy": energy_val,
+            "coverage_95": cov95,
+            "rmse": rmse_val,
+            "mae": mae_val,
+        }
 
 
 # --8<-- [end:world]
@@ -215,8 +228,10 @@ def score(
 # --8<-- [start:observables]
 observables = [
     Observable("crps", Direction.MINIMIZE),
+    Observable("energy", Direction.MINIMIZE, weight=0.5),
     Observable("coverage_95", Direction.MAXIMIZE, weight=0.5),
     Observable("rmse", Direction.MINIMIZE),
+    Observable("mae", Direction.MINIMIZE, weight=0.3),
 ]
 # --8<-- [end:observables]
 
@@ -237,6 +252,16 @@ def score(
 )
 # --8<-- [end:annotation]
 
+# --8<-- [start:constraint]
+# Require at least 90% empirical coverage at the 95% nominal level
+min_coverage = Constraint(
+    name="min_coverage",
+    observable="coverage_95",
+    op=">=",
+    threshold=0.90,
+)
+# --8<-- [end:constraint]
+
 
 def _run_screening(
     world: BayesianRegressionSimulator,
@@ -270,19 +295,26 @@ def _print_front(results: Any, front_idx: Any) -> None:
     # --8<-- [start:results]
     print(f"\nPareto front: {len(front_idx)} / {results.scores.shape[0]} designs")
 
-    print(
-        f"\n{'prior_var':>10s}  {'noise':>6s}  {'n_obs':>5s}  "
-        f"{'CRPS':>6s}  {'Cov95':>6s}  {'RMSE':>6s}  {'Cost':>5s}"
+    header = (
+        f"{'prior_var':>10s}  {'noise':>6s}  {'n_obs':>5s}  "
+        f"{'CRPS':>6s}  {'Energy':>7s}  {'Cov95':>6s}  "
+        f"{'RMSE':>6s}  {'MAE':>6s}  {'Cost':>5s}"
     )
-    print("-" * 58)
+    print(f"\n{header}")
+    print("-" * len(header))
     for i in front_idx:
         cfg = results.configs[i]
-        crps_val, cov, rmse_val = results.scores[i]
+        crps_val = results.scores[i, 0]
+        energy_val = results.scores[i, 1]
+        cov = results.scores[i, 2]
+        rmse_val = results.scores[i, 3]
+        mae_val = results.scores[i, 4]
         cost = results.annotations[i, 0] if results.annotations is not None else 0.0
         print(
             f"{cfg['prior_var']:10.2f}  {cfg['noise_scale']:6.2f}  "
             f"{cfg['n_obs']:5d}  "
-            f"{crps_val:6.3f}  {cov:6.3f}  {rmse_val:6.3f}  {cost:5.1f}"
+            f"{crps_val:6.3f}  {energy_val:7.3f}  {cov:6.3f}  "
+            f"{rmse_val:6.3f}  {mae_val:6.3f}  {cost:5.1f}"
         )
     # --8<-- [end:results]
 
@@ -314,7 +346,7 @@ def _run_stacking(results: Any, front_idx: Any, world: Any) -> None:
 
     ens_mean = ensemble_predict(front_predictions, stacking_weights)
     ens_rmse = float(np.sqrt(np.mean((ens_mean - Y_TEST) ** 2)))
-    best_single = float(results.scores[front_idx, 2].min())
+    best_single = float(results.scores[front_idx, 3].min())  # rmse column
     print(f"\nBest single-model RMSE: {best_single:.4f}")
     print(f"Stacked ensemble RMSE:  {ens_rmse:.4f}")
     # --8<-- [end:stacking]
@@ -388,31 +420,165 @@ def _save_plots(results: Any, directions: Any, nominal: Any, empirical: Any) ->
     # --8<-- [end:plots]
 
 
+def _run_igd_plus(results: Any, front_idx: Any) -> None:
+    """Compute IGD+ relative to a synthetic ideal front."""
+    # --8<-- [start:igd_plus]
+    directions = [o.direction for o in observables]
+    weights = [o.weight for o in observables]
+
+    front_scores = results.scores[front_idx]
+
+    # Build a synthetic reference front from per-objective best values
+    n_obj = front_scores.shape[1]
+    ideal = np.tile(np.median(front_scores, axis=0), (n_obj, 1))
+    for j, d in enumerate(directions):
+        ideal[j, j] = (
+            front_scores[:, j].min()
+            if d == Direction.MINIMIZE
+            else front_scores[:, j].max()
+        )
+
+    ref_point = np.array([
+        front_scores[:, j].max()
+        if d == Direction.MINIMIZE
+        else front_scores[:, j].min()
+        for j, d in enumerate(directions)
+    ])
+    hv = hypervolume(front_scores, ref_point, directions, weights)
+    igd = igd_plus(front_scores, ideal, directions, weights)
+    print(f"\nHypervolume: {hv:.4f}")
+    print(f"IGD+:        {igd:.4f}")
+    # --8<-- [end:igd_plus]
+
+
+def _run_feasibility(results: Any) -> None:
+    """Demonstrate Constraint + feasibility_filter."""
+    # --8<-- [start:feasibility]
+    all_idx = np.arange(results.scores.shape[0])
+    feas_filter = feasibility_filter(constraints=[min_coverage])
+    kept = feas_filter(results, observables)
+    print("\nFeasibility filter (coverage_95 >= 0.90):")
+    print(f"  {len(all_idx)} total -> {len(kept)} feasible designs")
+    # --8<-- [end:feasibility]
+
+
+def _run_weighted_sum(results: Any) -> None:
+    """Demonstrate weighted_sum_filter for scalarisation-based selection."""
+    # --8<-- [start:weighted_sum]
+    ws_filter = weighted_sum_filter(
+        weights={"crps": 1.0, "coverage_95": 0.5, "rmse": 0.8},
+        k=8,
+    )
+    kept = ws_filter(results, observables)
+    print("\nWeighted-sum filter (top 8 by scalarised score):")
+    for i in kept[:5]:
+        cfg = results.configs[i]
+        print(
+            f"  prior_var={cfg['prior_var']:.2f}, "
+            f"noise={cfg['noise_scale']:.2f}, "
+            f"n_obs={cfg['n_obs']}"
+        )
+    if len(kept) > 5:
+        print(f"  ... and {len(kept) - 5} more")
+    # --8<-- [end:weighted_sum]
+
+
+def _run_sobol_grid(world: Any, scorer: Any) -> None:
+    """Demonstrate Sobol quasi-random grid generation."""
+    # --8<-- [start:sobol]
+    sobol_grid = build_grid(factors, method="sobol", n_samples=40, seed=0)
+    for cfg in sobol_grid:
+        cfg["n_obs"] = round(cfg["n_obs"])
+
+    print(f"\nSobol grid: {len(sobol_grid)} configurations")
+    sobol_results = run_grid(
+        world=world,
+        scorer=scorer,
+        grid=sobol_grid,
+        observables=observables,
+    )
+    directions = [o.direction for o in observables]
+    weights = [o.weight for o in observables]
+    sobol_front = extract_front(sobol_results.scores, directions, weights)
+    print(f"Sobol Pareto front: {len(sobol_front)} designs")
+    # --8<-- [end:sobol]
+
+
+def _run_study_workflow() -> None:
+    """Demonstrate Study with summary(), stack(), and feasibility_filter."""
+    # --8<-- [start:study_workflow]
+    world = BayesianRegressionSimulator()
+    scorer = BayesianRegressionScorer()
+
+    grid = build_grid(factors, method="lhs", n_samples=60, seed=42)
+    for cfg in grid:
+        cfg["n_obs"] = round(cfg["n_obs"])
+
+    study = Study(
+        world=world,
+        scorer=scorer,
+        observables=observables,
+        phases=[
+            Phase(
+                name="explore",
+                grid=grid,
+                filter_fn=feasibility_filter(
+                    constraints=[min_coverage],
+                ),
+            ),
+            Phase(name="rank", grid="carry"),
+        ],
+        annotations=[compute_cost],
+    )
+    study.run()
+
+    # Study.summary() — per-phase statistics
+    summary = study.summary()
+    for phase_name, stats in summary.items():
+        print(f"\n  Phase '{phase_name}':")
+        print(f"    trials={stats['n_trials']}, front={stats['n_front']}")
+        for obs_name, rng in stats["observable_ranges"].items():
+            print(f"    {obs_name}: [{rng['min']:.4f}, {rng['max']:.4f}]")
+
+    # Study.stack() — convenience method for score-based stacking
+    stacking_w = study.stack("rank", maximize=False)
+    print("Study.stack() weights (top entries):")
+    for i, w in enumerate(stacking_w):
+        if w > 0.01:
+            print(f"  design {i}: {w:.3f}")
+    # --8<-- [end:study_workflow]
+
+
 def _run_multi_fidelity() -> None:
-    """Demonstrate multi-fidelity via Phase-level world overrides."""
+    """Demonstrate multi-fidelity via Phase-level world and scorer overrides."""
     # --8<-- [start:multifidelity]
     # Cheap surrogate: only 50 posterior draws (fast, noisier scores)
     cheap_world = BayesianRegressionSimulator(n_samples=50)
     # Expensive model: 2000 posterior draws (slow, precise scores)
     expensive_world = BayesianRegressionSimulator(n_samples=2000)
 
+    # Phase-level scorer override: use a separate scorer for screening
+    screen_scorer = BayesianRegressionScorer()
+    validate_scorer = BayesianRegressionScorer()
+
     grid = build_grid(factors, method="lhs", n_samples=60, seed=42)
     for cfg in grid:
         cfg["n_obs"] = round(cfg["n_obs"])
 
     study = Study(
         world=expensive_world,
-        scorer=BayesianRegressionScorer(),
+        scorer=validate_scorer,
         observables=observables,
         phases=[
-            # Phase 1: screen 60 designs with the cheap surrogate
+            # Phase 1: screen with cheap world AND its own scorer
             Phase(
                 name="screen",
                 grid=grid,
                 world=cheap_world,
+                scorer=screen_scorer,
                 filter_fn=top_k_pareto_filter(k=10),
             ),
-            # Phase 2: validate top 10 with the expensive model
+            # Phase 2: validate top 10 with the expensive model + default scorer
             Phase(name="validate", grid="carry"),
         ],
         annotations=[compute_cost],
@@ -460,10 +626,15 @@ def main() -> None:
     front_idx = extract_front(results.scores, directions, weights)
 
     _print_front(results, front_idx)
+    _run_igd_plus(results, front_idx)
+    _run_feasibility(results)
+    _run_weighted_sum(results)
     _run_stacking(results, front_idx, world)
     nominal, empirical = _run_calibration(results, front_idx, world)
     _run_persistence(results)
     _save_plots(results, directions, nominal, empirical)
+    _run_sobol_grid(world, scorer)
+    _run_study_workflow()
     _run_multi_fidelity()
 
 
diff --git a/examples/cstr_study.py b/examples/cstr_study.py
index 846b521..dd31e06 100644
--- a/examples/cstr_study.py
+++ b/examples/cstr_study.py
@@ -19,6 +19,7 @@
 import numpy as np
 
 from trade_study import (
+    Constraint,
     Direction,
     Factor,
     FactorType,
@@ -28,9 +29,13 @@
     Study,
     build_grid,
     extract_front,
+    feasibility_filter,
+    igd_plus,
     plot_front,
     plot_parallel,
     plot_scores,
+    reduce_factors,
+    screen,
     top_k_pareto_filter,
 )
 
@@ -169,6 +174,22 @@ def score(
 ]
 # --8<-- [end:factors]
 
+# --8<-- [start:constraint]
+# Require conversion >= 0.50 and energy cost <= 100
+min_conversion = Constraint(
+    name="min_conversion",
+    observable="conversion",
+    op=">=",
+    threshold=0.50,
+)
+max_energy = Constraint(
+    name="max_energy",
+    observable="energy_cost",
+    op="<=",
+    threshold=100.0,
+)
+# --8<-- [end:constraint]
+
 
 # --8<-- [start:refine]
 def refine_grid(
@@ -235,6 +256,46 @@ def refine_grid(
 # --8<-- [end:phases]
 
 
+def _run_screening() -> None:
+    """Run Sobol screening to identify influential factors."""
+    # --8<-- [start:screening]
+    world = CSTRSimulator()
+    scorer = CSTRScorer()
+
+    # Screening requires continuous factors; define ranges matching the
+    # discrete levels so SALib can build its sample matrix.
+    screening_factors = [
+        Factor("temperature", FactorType.CONTINUOUS, bounds=(330.0, 390.0)),
+        Factor("residence_time", FactorType.CONTINUOUS, bounds=(20.0, 110.0)),
+        Factor("inlet_concentration", FactorType.CONTINUOUS, bounds=(0.5, 2.5)),
+        Factor("coolant_flow", FactorType.CONTINUOUS, bounds=(1.0, 3.0)),
+    ]
+
+    def run_fn(config: dict[str, Any]) -> dict[str, float]:
+        """Compose simulator + scorer for screening.
+
+        Returns:
+            Score dictionary from the scorer.
+        """
+        truth, obs = world.generate(config)
+        return scorer.score(truth, obs, config)
+
+    importance = screen(
+        run_fn,
+        screening_factors,
+        method="sobol",
+        n_trajectories=8,
+        seed=42,
+    )
+    print("Sobol screening:")
+    for name, vals in importance.items():
+        print(f"  {name}: {vals}")
+
+    reduced = reduce_factors(screening_factors, importance, threshold=0.01)
+    print(f"\nImportant factors: {[f.name for f in reduced]}")
+    # --8<-- [end:screening]
+
+
 def _plot_kinetics(plt: Any) -> None:
     """Plot conversion & selectivity vs temperature for several residence times."""
     temps = np.linspace(320, 400, 200)
@@ -327,8 +388,133 @@ def _print_results(study: Study) -> None:
     print(f"\nHypervolume: {hv:.2f}")
 
 
+def _print_quality(study: Study) -> None:
+    """Print IGD+ and per-phase summary statistics."""
+    directions = [o.direction for o in observables]
+    r2 = study.results("refinement")
+
+    # IGD+ relative to synthetic ideal
+    r2_front_idx = study.front("refinement")
+    front_scores = r2.scores[r2_front_idx]
+    n_obj = front_scores.shape[1]
+    ideal = np.tile(np.median(front_scores, axis=0), (n_obj, 1))
+    for j, d in enumerate(directions):
+        ideal[j, j] = (
+            front_scores[:, j].min()
+            if d == Direction.MINIMIZE
+            else front_scores[:, j].max()
+        )
+    igd = igd_plus(front_scores, ideal, directions)
+    print(f"IGD+:        {igd:.4f}")
+
+    # Study.summary()
+    summary = study.summary()
+    for phase_name, stats in summary.items():
+        print(f"\n  Phase '{phase_name}':")
+        print(f"    trials={stats['n_trials']}, front={stats['n_front']}")
+        for obs_name, rng in stats["observable_ranges"].items():
+            print(f"    {obs_name}: [{rng['min']:.4f}, {rng['max']:.4f}]")
+
+
+def _run_with_callback() -> None:
+    """Demonstrate run_grid with a progress callback and n_jobs."""
+    # --8<-- [start:callback]
+    from trade_study import TrialResult, run_grid
+
+    world = CSTRSimulator()
+    scorer = CSTRScorer()
+    grid = build_grid(factors, method="full")
+
+    completed = 0
+
+    def progress(_idx: int, total: int, _result: TrialResult) -> None:
+        """Print progress every 20 trials.
+
+        Args:
+            _idx: Current trial index (unused).
+            total: Total number of trials.
+            _result: The completed trial result (unused).
+        """
+        nonlocal completed
+        completed += 1
+        if completed % 20 == 0 or completed == total:
+            print(f"  Progress: {completed}/{total} trials done")
+
+    results = run_grid(
+        world=world,
+        scorer=scorer,
+        grid=grid,
+        observables=observables,
+        n_jobs=2,
+        callback=progress,
+    )
+    print(f"\nParallel run (n_jobs=2): {results.scores.shape[0]} trials")
+    # --8<-- [end:callback]
+
+
+def _run_adaptive() -> None:
+    """Demonstrate adaptive (optuna-driven) phase."""
+    # --8<-- [start:adaptive]
+    # Adaptive phase uses optuna NSGA-II to explore the design space
+    adaptive_phases = [
+        Phase(
+            name="adaptive_search",
+            grid="adaptive",
+            n_trials=50,
+            filter_fn=top_k_pareto_filter(15),
+        ),
+        Phase(name="validate", grid="carry"),
+    ]
+
+    study = Study(
+        world=CSTRSimulator(),
+        scorer=CSTRScorer(),
+        observables=observables,
+        phases=adaptive_phases,
+        factors=factors,
+    )
+    study.run()
+
+    r = study.results("validate")
+    directions = [o.direction for o in observables]
+    front = extract_front(r.scores, directions)
+    print(f"\nAdaptive study: {r.scores.shape[0]} validated, {len(front)} on front")
+    # --8<-- [end:adaptive]
+
+
+def _run_feasibility() -> None:
+    """Demonstrate feasibility_filter with constraints."""
+    # --8<-- [start:feasibility]
+    feas_phases = [
+        Phase(
+            name="sweep",
+            grid=build_grid(factors, method="full"),
+            filter_fn=feasibility_filter(
+                constraints=[min_conversion, max_energy],
+            ),
+        ),
+        Phase(name="feasible", grid="carry"),
+    ]
+
+    study = Study(
+        world=CSTRSimulator(),
+        scorer=CSTRScorer(),
+        observables=observables,
+        phases=feas_phases,
+    )
+    study.run()
+
+    sweep_r = study.results("sweep")
+    feas_r = study.results("feasible")
+    print("\nFeasibility filter:")
+    print(f"  {sweep_r.scores.shape[0]} total -> {feas_r.scores.shape[0]} feasible")
+    # --8<-- [end:feasibility]
+
+
 def main() -> None:
     """Run the CSTR trade study and print results."""
+    _run_screening()
+
     # --8<-- [start:run]
     study = Study(
         world=CSTRSimulator(),
@@ -341,6 +527,7 @@ def main() -> None:
 
     # --8<-- [start:results]
     _print_results(study)
+    _print_quality(study)
     # --8<-- [end:results]
 
     # --8<-- [start:plots]
@@ -376,6 +563,10 @@ def main() -> None:
     plt.close(fig_sel2)
     # --8<-- [end:plots]
 
+    _run_with_callback()
+    _run_adaptive()
+    _run_feasibility()
+
 
 if __name__ == "__main__":
     main()
diff --git a/mkdocs.yml b/mkdocs.yml
index 7d9c027..7a98342 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -74,4 +74,5 @@ nav:
       - Scoring: api/scoring.md
       - Pareto: api/pareto.md
       - Stacking: api/stacking.md
+      - Visualization: api/viz.md
       - I/O: api/io.md