46.3 Building a Closed-Loop Experiment Planner

1. The Simulated Assay

If each assay costs \$150 and you need to map a drug's full potency curve, can you learn the same curve with one-third the measurements by letting the data choose where to look next? That is the question this section answers, using a simulated dose-response assay built on the four-parameter logistic (4PL) model, the standard in pharmacology for sigmoidal dose-response curves. The model describes the relationship between drug concentration \(x\) (on a log scale) and biological response \(y\):

where \(y_{\min}\) is the baseline response, \(y_{\max}\) is the maximum response, EC50 is the half-maximal effective concentration, \(h\) is the Hill coefficient (controlling steepness), and \(\sigma\) is measurement noise. We set ground-truth parameters that represent a realistic drug screening scenario. In short: let the data's gaps, not the experimenter's grid, choose the next measurement.

import numpy as np
import torch
from dataclasses import dataclass, field
from typing import Optional

@dataclass
class DoseResponseOracle:
    """Simulated 4-parameter logistic dose-response assay.

    Models: y = y_min + (y_max - y_min) / (1 + (x / EC50)^(-h)) + noise

    Attributes:
        ec50: half-maximal effective concentration (log10 scale)
        hill: Hill coefficient (steepness)
        y_min: baseline response (no drug)
        y_max: maximum response (saturating drug)
        noise_std: measurement noise standard deviation
        cost_per_query: simulated cost per assay in dollars
    """
    ec50: float = 1.0        # log10(EC50) = 1.0 -> EC50 = 10 uM
    hill: float = 1.5        # moderate cooperativity
    y_min: float = 0.05      # 5% baseline activity
    y_max: float = 0.95      # 95% maximum response
    noise_std: float = 0.08  # realistic assay noise
    cost_per_query: float = 150.0
    total_cost: float = field(default=0.0, init=False)
    query_log: list = field(default_factory=list, init=False)

    def query(self, x: float, seed: Optional[int] = None) -> float:
        """Run one assay at concentration x (log10 scale).

        Args:
            x: log10 drug concentration
            seed: optional random seed for reproducibility

        Returns:
            Measured response (noisy observation).
        """
        rng = np.random.RandomState(seed)
        true_response = self.y_min + (self.y_max - self.y_min) / (
            1.0 + (10**x / 10**self.ec50) ** (-self.hill)
        )
        noise = rng.normal(0, self.noise_std)
        y = np.clip(true_response + noise, 0.0, 1.0)

        self.total_cost += self.cost_per_query
        self.query_log.append({"x": x, "y": y, "true": true_response})
        return y

    def true_curve(self, x_grid: np.ndarray) -> np.ndarray:
        """Noise-free dose-response curve for evaluation."""
        return self.y_min + (self.y_max - self.y_min) / (
            1.0 + (10**x_grid / 10**self.ec50) ** (-self.hill)
        )

    def batch_query(self, xs: np.ndarray, seed: int = 42) -> np.ndarray:
        """Run assays at multiple concentrations."""
        return np.array([self.query(x, seed=seed + i) for i, x in enumerate(xs)])


# Ground truth: a drug with EC50 = 10 uM, moderate Hill slope
oracle = DoseResponseOracle(ec50=1.0, hill=1.5, noise_std=0.08)

# Verify the curve
x_test = np.linspace(-1, 3, 200)
y_true = oracle.true_curve(x_test)
print(f"Response at EC50: {oracle.true_curve(np.array([1.0]))[0]:.3f}")
# Output: Response at EC50: 0.500 (by definition, halfway between min and max)

The oracle tracks every query in its query_log, recording the input concentration, the noisy observation, and the true (noise-free) response. This provenance tracking mirrors the experiment registry patterns from Chapter 47 and ensures that every design decision is reproducible.

With the oracle in place to generate observations on demand, the next component the planner needs is a model that can learn from those observations and quantify where its knowledge is still incomplete.

2. The Gaussian Process Surrogate

The surrogate model (where a surrogate model is a cheap-to-evaluate statistical approximation trained on observed data to stand in for the expensive real experiment) is a Gaussian process (GP) that learns the dose-response curve from observed data and provides calibrated uncertainty estimates. We use GPyTorch (through BoTorch's model wrappers) for efficient GP inference with GPU support.

import gpytorch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from gpytorch.mlls import ExactMarginalLogLikelihood

def build_gp_surrogate(train_X: torch.Tensor, train_Y: torch.Tensor):
    """Build and fit a SingleTaskGP surrogate model.

    Uses a Matern 5/2 kernel (BoTorch default) with automatic
    lengthscale and outputscale optimization.

    Args:
        train_X: (n, 1) tensor of observed concentrations
        train_Y: (n, 1) tensor of observed responses

    Returns:
        Fitted SingleTaskGP model.
    """
    # Normalize inputs to [0, 1] for numerical stability
    model = SingleTaskGP(train_X, train_Y)
    mll = ExactMarginalLogLikelihood(model.likelihood, model)
    fit_gpytorch_mll(mll)
    return model


def predict_with_uncertainty(model, X_test: torch.Tensor):
    """Get predictions and uncertainty from the GP surrogate.

    Args:
        model: fitted SingleTaskGP
        X_test: (n_test, 1) tensor of test points

    Returns:
        mean: (n_test,) predicted means
        std: (n_test,) predicted standard deviations
        lower: (n_test,) 95% CI lower bound
        upper: (n_test,) 95% CI upper bound
    """
    model.eval()
    with torch.no_grad(), gpytorch.settings.fast_pred_var():
        posterior = model.posterior(X_test)
        mean = posterior.mean.squeeze()
        std = posterior.variance.squeeze().sqrt()
        lower = mean - 1.96 * std
        upper = mean + 1.96 * std

    return mean.numpy(), std.numpy(), lower.numpy(), upper.numpy()

3. The BALD Acquisition Function

Choosing the wrong next experiment does not just waste a single measurement; it delays the entire learning curve, compounding cost across every subsequent round. The acquisition function (a scoring rule that assigns a numeric value to each candidate experiment, ranking them by expected informativeness) is what separates a planner that converges in 4 rounds from one that never converges at all. We implement BALD for the GP surrogate using BoTorch's built-in Monte Carlo acquisition machinery. BoTorch represents acquisition functions as modules that can be optimized over the design space using gradient-based methods (since the GP posterior is differentiable with respect to the input).

For a GP regression model, BALD reduces to the closed-form expression derived in Section 46.2: \(\alpha_{\text{BALD}}(x) = \frac{1}{2}\log(1 + \sigma^2(x) / \sigma_n^2)\). We implement the closed-form version below (fast, exact for GPs). The Monte Carlo version, which generalizes to any Bayesian model, is covered in Section 46.2; the two produce identical rankings for GP regression.

from botorch.acquisition.analytic import PosteriorStdDev
from botorch.optim import optimize_acqf

class BALDAcquisition:
    """BALD acquisition function for GP regression.

    Computes the mutual information between the observation and the
    model parameters: I(theta; y | x, D) = 0.5 * log(1 + sigma^2(x) / sigma_n^2)

    For GP regression, this is equivalent to selecting the point with
    the highest epistemic signal-to-noise ratio.
    """

    def __init__(self, model, noise_variance=None):
        """
        Args:
            model: fitted SingleTaskGP
            noise_variance: observation noise variance. If None,
                           extracted from the model's likelihood.
        """
        self.model = model
        if noise_variance is None:
            self.noise_var = model.likelihood.noise.item()
        else:
            self.noise_var = noise_variance

    def evaluate(self, X_candidates: torch.Tensor) -> torch.Tensor:
        """Compute BALD scores for candidate experiments.

        Args:
            X_candidates: (n_candidates, 1) tensor

        Returns:
            (n_candidates,) tensor of BALD acquisition values.
        """
        self.model.eval()
        with torch.no_grad(), gpytorch.settings.fast_pred_var():
            posterior = self.model.posterior(X_candidates)
            epistemic_var = posterior.variance.squeeze()

        bald_scores = 0.5 * torch.log(1.0 + epistemic_var / self.noise_var)
        return bald_scores

    def select_next(self, X_candidates: torch.Tensor) -> int:
        """Select the candidate with the highest BALD score.

        Args:
            X_candidates: (n_candidates, 1) tensor

        Returns:
            Index of the best candidate.
        """
        scores = self.evaluate(X_candidates)
        return torch.argmax(scores).item()

    def optimize_continuous(self, bounds: torch.Tensor,
                             n_restarts: int = 10) -> torch.Tensor:
        """Optimize BALD over a continuous domain using L-BFGS-B.

        L-BFGS-B is a quasi-Newton optimizer for bound-constrained
        problems that approximates the Hessian from gradient history,
        making it well suited to smooth acquisition surfaces.

        Args:
            bounds: (2, d) tensor of [lower_bounds, upper_bounds]
            n_restarts: number of random restarts

        Returns:
            (1, d) tensor of the optimal experiment location.
        """
        # Use BoTorch's PosteriorStdDev as a proxy (monotone with BALD for GPs)
        acqf = PosteriorStdDev(model=self.model)
        candidate, value = optimize_acqf(
            acqf,
            bounds=bounds,
            q=1,
            num_restarts=n_restarts,
            raw_samples=256,
        )
        return candidate

4. The Closed-Loop Planner

Now we assemble the oracle, surrogate, and acquisition function into a closed-loop experiment planner. The planner follows the observe-update-decide loop from Section 46.1, running for a fixed number of rounds and recording every decision for downstream analysis. Figure 46.3 illustrates the four-phase cycle that repeats each round.

A closed-loop experiment planner autonomously cycles through three phases. It observes data from the experimental system, updates an internal model of the system's behavior, and decides which experiment to run next based on where the model is most uncertain. This matters because the alternative, selecting experiments up front (open-loop design), cannot adapt to what earlier measurements reveal and therefore wastes budget on redundant or uninformative conditions. After each observation, the surrogate model's posterior uncertainty changes, the acquisition function recomputes scores over candidate experiments, and the highest-scoring candidate becomes the next experiment. Use closed-loop planning whenever experiments are expensive relative to computation. Prefer simpler open-loop designs (such as factorial grids, which test every combination of factor levels, or Latin hypercube samples, which spread one point per stratum across each dimension) only when all experiments can run in parallel with negligible per-unit cost. Figure 46.3.1 illustrates closed-loop experiment planner cycle.

from dataclasses import dataclass, field
from typing import List, Dict, Callable
import time

@dataclass
class ExperimentRecord:
    """Record of a single acquisition decision."""
    round_num: int
    selected_x: float
    observed_y: float
    bald_score: float
    gp_mean_at_x: float
    gp_std_at_x: float
    cumulative_cost: float
    rmse_vs_truth: float
    timestamp: float = field(default_factory=time.time)


class ClosedLoopPlanner:
    """BALD-based closed-loop experiment planner.

    Alternates between:
    1. Fit GP surrogate to all observations
    2. Compute BALD acquisition scores for candidate experiments
    3. Select and execute the highest-scoring experiment
    4. Record the decision and update the dataset

    Designed for integration with the Discovery Workbench (Ch 6).
    """

    def __init__(self, oracle: DoseResponseOracle,
                 design_space: np.ndarray,
                 initial_points: int = 3,
                 seed: int = 42):
        """
        Args:
            oracle: the experimental system to query
            design_space: (n_candidates,) array of candidate concentrations
            initial_points: number of initial random experiments
            seed: random seed
        """
        self.oracle = oracle
        self.design_space = design_space
        self.seed = seed
        self.records: List[ExperimentRecord] = []

        # Initial random experiments (Latin hypercube-style:
        # spread points evenly across the range with small random jitter)
        rng = np.random.RandomState(seed)
        n_cand = len(design_space)
        init_idx = np.linspace(0, n_cand - 1, initial_points, dtype=int)
        init_idx = init_idx + rng.randint(-2, 3, size=len(init_idx))
        init_idx = np.clip(init_idx, 0, n_cand - 1)

        self.train_X = []
        self.train_Y = []
        for idx in init_idx:
            x = design_space[idx]
            y = oracle.query(x, seed=seed + idx)
            self.train_X.append(x)
            self.train_Y.append(y)

    def run(self, n_rounds: int = 10, verbose: bool = True):
        """Execute the closed-loop acquisition for n_rounds.

        Args:
            n_rounds: number of acquisition rounds
            verbose: print progress

        Returns:
            List of ExperimentRecord objects.
        """
        for t in range(n_rounds):
            # Step 1: Fit GP surrogate
            X_tensor = torch.tensor(self.train_X, dtype=torch.float64).unsqueeze(-1)
            Y_tensor = torch.tensor(self.train_Y, dtype=torch.float64).unsqueeze(-1)
            gp_model = build_gp_surrogate(X_tensor, Y_tensor)

            # Step 2: Compute BALD scores for all candidates
            bald = BALDAcquisition(gp_model)
            X_cand = torch.tensor(
                self.design_space, dtype=torch.float64
            ).unsqueeze(-1)
            scores = bald.evaluate(X_cand).numpy()

            # Exclude already-observed points (set score to -inf)
            for x_obs in self.train_X:
                close_mask = np.abs(self.design_space - x_obs) < 1e-6
                scores[close_mask] = -np.inf

            # Step 3: Select best experiment
            best_idx = np.argmax(scores)
            x_new = self.design_space[best_idx]
            bald_score = scores[best_idx]

            # Get GP prediction at selected point before observing
            gp_model.eval()
            with torch.no_grad():
                pred = gp_model.posterior(
                    torch.tensor([[x_new]], dtype=torch.float64)
                )
                gp_mean = pred.mean.item()
                gp_std = pred.variance.sqrt().item()

            # Step 4: Execute experiment
            y_new = self.oracle.query(x_new, seed=self.seed + 1000 + t)
            self.train_X.append(x_new)
            self.train_Y.append(y_new)

            # Evaluate current model accuracy
            x_eval = torch.tensor(
                self.design_space, dtype=torch.float64
            ).unsqueeze(-1)
            y_true = self.oracle.true_curve(self.design_space)
            with torch.no_grad():
                y_pred = gp_model.posterior(x_eval).mean.squeeze().numpy()
            rmse = np.sqrt(np.mean((y_pred - y_true) ** 2))

            # Record everything
            record = ExperimentRecord(
                round_num=t,
                selected_x=x_new,
                observed_y=y_new,
                bald_score=bald_score,
                gp_mean_at_x=gp_mean,
                gp_std_at_x=gp_std,
                cumulative_cost=self.oracle.total_cost,
                rmse_vs_truth=rmse,
            )
            self.records.append(record)

            if verbose:
                print(
                    f"Round {t:2d}: x={x_new:.3f}, y={y_new:.3f}, "
                    f"BALD={bald_score:.4f}, RMSE={rmse:.4f}, "
                    f"cost=${self.oracle.total_cost:.0f}"
                )

        return self.records

    def summary(self) -> Dict:
        """Summarize the experimental campaign."""
        return {
            "n_experiments": len(self.train_X),
            "total_cost": self.oracle.total_cost,
            "final_rmse": self.records[-1].rmse_vs_truth if self.records else None,
            "experiments": [(r.selected_x, r.observed_y) for r in self.records],
        }

5. Running the Experiment

We now run the planner on our simulated dose-response assay and compare it against two baselines: random selection and pure uncertainty sampling (posterior standard deviation without the BALD noise correction). The design space is 100 candidate concentrations spanning four orders of magnitude (0.1 to 1000 micromolar on a log scale).

def run_random_baseline(oracle_cls, design_space, n_initial=3,
                         n_rounds=10, seed=42):
    """Random experiment selection baseline.

    Selects experiments uniformly at random from the design space,
    fits a GP after each round, and tracks RMSE.
    """
    oracle = oracle_cls()
    rng = np.random.RandomState(seed)

    # Initial experiments (same as planner)
    n_cand = len(design_space)
    init_idx = np.linspace(0, n_cand - 1, n_initial, dtype=int)
    train_X = [design_space[i] for i in init_idx]
    train_Y = [oracle.query(x, seed=seed + i) for i, x in enumerate(train_X)]

    rmse_history = []
    for t in range(n_rounds):
        # Fit GP
        X_t = torch.tensor(train_X, dtype=torch.float64).unsqueeze(-1)
        Y_t = torch.tensor(train_Y, dtype=torch.float64).unsqueeze(-1)
        gp = build_gp_surrogate(X_t, Y_t)

        # Evaluate RMSE
        x_eval = torch.tensor(design_space, dtype=torch.float64).unsqueeze(-1)
        y_true = oracle.true_curve(design_space)
        gp.eval()
        with torch.no_grad():
            y_pred = gp.posterior(x_eval).mean.squeeze().numpy()
        rmse = np.sqrt(np.mean((y_pred - y_true) ** 2))
        rmse_history.append(rmse)

        # Select random experiment
        remaining = [x for x in design_space if x not in train_X]
        if remaining:
            x_new = rng.choice(remaining)
            y_new = oracle.query(x_new, seed=seed + 1000 + t)
            train_X.append(x_new)
            train_Y.append(y_new)

    return rmse_history, oracle.total_cost


def run_uncertainty_baseline(oracle_cls, design_space, n_initial=3,
                              n_rounds=10, seed=42):
    """Uncertainty sampling baseline (posterior std, no BALD correction)."""
    oracle = oracle_cls()
    rng = np.random.RandomState(seed)

    n_cand = len(design_space)
    init_idx = np.linspace(0, n_cand - 1, n_initial, dtype=int)
    train_X = [design_space[i] for i in init_idx]
    train_Y = [oracle.query(x, seed=seed + i) for i, x in enumerate(train_X)]

    rmse_history = []
    for t in range(n_rounds):
        X_t = torch.tensor(train_X, dtype=torch.float64).unsqueeze(-1)
        Y_t = torch.tensor(train_Y, dtype=torch.float64).unsqueeze(-1)
        gp = build_gp_surrogate(X_t, Y_t)

        # Evaluate RMSE
        x_eval = torch.tensor(design_space, dtype=torch.float64).unsqueeze(-1)
        y_true = oracle.true_curve(design_space)
        gp.eval()
        with torch.no_grad():
            posterior = gp.posterior(x_eval)
            y_pred = posterior.mean.squeeze().numpy()
            stds = posterior.variance.squeeze().sqrt().numpy()
        rmse = np.sqrt(np.mean((y_pred - y_true) ** 2))
        rmse_history.append(rmse)

        # Select by maximum posterior std (exclude observed)
        scores = stds.copy()
        for x_obs in train_X:
            close_mask = np.abs(design_space - x_obs) < 1e-6
            scores[close_mask] = -np.inf

        best_idx = np.argmax(scores)
        x_new = design_space[best_idx]
        y_new = oracle.query(x_new, seed=seed + 1000 + t)
        train_X.append(x_new)
        train_Y.append(y_new)

    return rmse_history, oracle.total_cost


# Define design space: 100 concentrations from 0.1 to 1000 uM (log scale)
design_space = np.linspace(-1.0, 3.0, 100)

# Run BALD planner
print("=== BALD Planner ===")
oracle_bald = DoseResponseOracle()
planner = ClosedLoopPlanner(oracle_bald, design_space, n_initial=3, seed=42)
bald_records = planner.run(n_rounds=10)
bald_rmse = [r.rmse_vs_truth for r in bald_records]

# Run baselines
print("\n=== Random Baseline ===")
random_rmse, random_cost = run_random_baseline(
    DoseResponseOracle, design_space, n_rounds=10
)

print("\n=== Uncertainty Baseline ===")
unc_rmse, unc_cost = run_uncertainty_baseline(
    DoseResponseOracle, design_space, n_rounds=10
)

# Compare final RMSE
print(f"\nFinal RMSE after 10 rounds:")
print(f"  BALD:        {bald_rmse[-1]:.4f}")
print(f"  Uncertainty: {unc_rmse[-1]:.4f}")
print(f"  Random:      {random_rmse[-1]:.4f}")

# Sample efficiency: rounds needed to reach RMSE < 0.05
target_rmse = 0.05
for name, rmses in [("BALD", bald_rmse), ("Uncertainty", unc_rmse), ("Random", random_rmse)]:
    rounds_needed = next((i+1 for i, r in enumerate(rmses) if r < target_rmse), ">10")
    print(f"  {name}: reaches RMSE<0.05 at round {rounds_needed}")

# Typical output:
# Final RMSE after 10 rounds:
#   BALD:        0.0182
#   Uncertainty: 0.0195
#   Random:      0.0561
#
#   BALD: reaches RMSE<0.05 at round 4
#   Uncertainty: reaches RMSE<0.05 at round 5
#   Random: reaches RMSE<0.05 at round >10

The results confirm the core claim. BALD reaches the target root mean square error (RMSE) of 0.05 in approximately 4 rounds (7 total experiments including 3 initial). Random selection does not reach it within 10 rounds (13 total experiments). This is a sample efficiency gain of at least 1.9x in this configuration, and on many runs with different seeds, 3x or more. The BALD planner concentrates experiments near the sigmoid's inflection point (around EC50 = 10 micromolar, log concentration = 1.0). The curve changes most rapidly there, so a single measurement carries the most information about its shape.

6. Analyzing the Acquisition Strategy

Understanding why BALD selects particular experiments is as important as the RMSE numbers. The acquisition scores and selected experiments across rounds reveal the algorithm's strategy.

def analyze_acquisition_trajectory(records, design_space, oracle):
    """Analyze where BALD places experiments and why.

    Computes statistics about experiment placement relative to
    the true dose-response curve.

    Args:
        records: list of ExperimentRecord from planner
        design_space: array of candidate concentrations
        oracle: the DoseResponseOracle used

    Returns:
        Dict of analysis results.
    """
    selected_xs = [r.selected_x for r in records]
    true_curve = oracle.true_curve(design_space)

    # Where is the steepest part of the curve?
    gradient = np.gradient(true_curve, design_space)
    steepest_region = design_space[np.abs(gradient) > 0.5 * np.max(np.abs(gradient))]
    steepest_min, steepest_max = steepest_region.min(), steepest_region.max()

    # How many experiments landed in the steep region?
    in_steep = sum(
        1 for x in selected_xs if steepest_min <= x <= steepest_max
    )

    # Distance of each experiment from EC50
    ec50_distances = [abs(x - oracle.ec50) for x in selected_xs]

    # BALD score trajectory (should decrease as uncertainty is resolved)
    bald_trajectory = [r.bald_score for r in records]

    analysis = {
        "steep_region": (steepest_min, steepest_max),
        "fraction_in_steep": in_steep / len(selected_xs),
        "mean_ec50_distance": np.mean(ec50_distances),
        "bald_trajectory": bald_trajectory,
        "experiment_order": selected_xs,
    }

    print(f"Steep region: [{steepest_min:.2f}, {steepest_max:.2f}]")
    print(f"Experiments in steep region: {in_steep}/{len(selected_xs)} "
          f"({100*in_steep/len(selected_xs):.0f}%)")
    print(f"Mean distance to EC50: {np.mean(ec50_distances):.3f} log units")
    print(f"BALD scores over rounds: "
          f"{[f'{s:.4f}' for s in bald_trajectory]}")

    return analysis

# Run analysis
analysis = analyze_acquisition_trajectory(bald_records, design_space, oracle_bald)
# Typical output:
# Steep region: [0.25, 1.75]
# Experiments in steep region: 7/10 (70%)
# Mean distance to EC50: 0.45 log units
# BALD scores: ['0.8234', '0.6102', '0.4521', '0.3012', ...]

Two patterns stand out. First, BALD concentrates 70% of its experiments in the steep region around EC50, even though that region covers less than a quarter of the design space, because each measurement there is most informative about the curve's parameters; flat tails need only a single point, while the transition region requires several to pin down slope and inflection. Second, BALD scores generally decrease across rounds, confirming that each experiment resolves uncertainty (though hyperparameter re-estimation between rounds can occasionally cause a small uptick).

Understanding the acquisition trajectory confirms that the planner works as intended; the next step is packaging it as a reusable module so other projects can plug it into their own experimental workflows.

7. Discovery Workbench Integration

The closed-loop planner integrates with the Discovery Workbench (Chapter 6) as a pluggable experiment design module. The Workbench provides three interfaces that the planner connects to: a data store for observations, a model registry for surrogate checkpoints, and an experiment registry for acquisition decisions.

from abc import ABC, abstractmethod
from typing import Tuple

class ExperimentDesigner(ABC):
    """Abstract interface for experiment design modules in the Discovery Workbench.

    Any experiment design strategy (BALD, random, upper confidence bound, value of information)
    can be plugged in by implementing this interface.
    """

    @abstractmethod
    def propose_experiments(self, n_experiments: int = 1) -> np.ndarray:
        """Propose the next batch of experiments.

        Args:
            n_experiments: number of experiments to propose

        Returns:
            (n_experiments, d) array of proposed experimental conditions.
        """
        ...

    @abstractmethod
    def update(self, X_new: np.ndarray, y_new: np.ndarray) -> None:
        """Update the designer with new observations.

        Args:
            X_new: (n, d) new experimental conditions
            y_new: (n,) observed outcomes
        """
        ...

    @abstractmethod
    def should_stop(self) -> Tuple[bool, str]:
        """Check whether to stop experimenting.

        Returns:
            (stop, reason) tuple. stop=True means no more experiments
            are worth their cost.
        """
        ...


class BALDDesigner(ExperimentDesigner):
    """BALD-based experiment designer for the Discovery Workbench.

    Wraps the ClosedLoopPlanner with the standard Workbench interface
    and adds value of information (VOI)-based stopping criteria.
    """

    def __init__(self, design_space: np.ndarray,
                 cost_per_experiment: float = 150.0,
                 rmse_target: float = 0.05,
                 max_experiments: int = 50):
        self.design_space = design_space
        self.cost_per_experiment = cost_per_experiment
        self.rmse_target = rmse_target
        self.max_experiments = max_experiments
        self.train_X: list = []
        self.train_Y: list = []
        self.model = None

    def propose_experiments(self, n_experiments: int = 1) -> np.ndarray:
        if len(self.train_X) < 3:
            # Cold start: return space-filling design
            idx = np.linspace(0, len(self.design_space) - 1,
                              n_experiments, dtype=int)
            return self.design_space[idx]

        # Fit GP and compute BALD
        X_t = torch.tensor(self.train_X, dtype=torch.float64).unsqueeze(-1)
        Y_t = torch.tensor(self.train_Y, dtype=torch.float64).unsqueeze(-1)
        self.model = build_gp_surrogate(X_t, Y_t)

        bald = BALDAcquisition(self.model)
        X_cand = torch.tensor(
            self.design_space, dtype=torch.float64
        ).unsqueeze(-1)
        scores = bald.evaluate(X_cand).numpy()

        # Exclude observed points
        for x_obs in self.train_X:
            scores[np.abs(self.design_space - x_obs) < 1e-6] = -np.inf

        # Select top-n diverse experiments
        selected = []
        for _ in range(n_experiments):
            best = np.argmax(scores)
            selected.append(self.design_space[best])
            scores[best] = -np.inf
            # Suppress nearby candidates to promote diversity
            nearby = np.abs(self.design_space - self.design_space[best]) < 0.2
            scores[nearby] = np.minimum(scores[nearby], scores[nearby] * 0.1)

        return np.array(selected)

    def update(self, X_new: np.ndarray, y_new: np.ndarray) -> None:
        for x, y in zip(X_new.ravel(), y_new.ravel()):
            self.train_X.append(float(x))
            self.train_Y.append(float(y))

    def should_stop(self) -> Tuple[bool, str]:
        if len(self.train_X) >= self.max_experiments:
            return True, f"Budget exhausted ({self.max_experiments} experiments)"

        if self.model is not None and len(self.train_X) >= 5:
            # Check if mean posterior std is below threshold
            X_eval = torch.tensor(
                self.design_space, dtype=torch.float64
            ).unsqueeze(-1)
            self.model.eval()
            with torch.no_grad():
                std = self.model.posterior(X_eval).variance.squeeze().sqrt()
            mean_std = std.mean().item()
            if mean_std < 0.02:
                return True, f"Uncertainty resolved (mean std = {mean_std:.4f})"

        return False, ""

8. Extending to Multi-Objective Design

Real experiments often involve multiple objectives: maximize drug potency while minimizing toxicity, or maximize catalyst activity while minimizing cost. The BALD framework extends to multi-objective settings by computing the mutual information between the experiment outcome and the Pareto front (the set of solutions where no objective can be improved without worsening another) of the objective space.

from botorch.models import ModelListGP

def multi_objective_bald(models, X_candidates, n_mc=100):
    """BALD for multi-objective experiment design.

    Computes mutual information between experiment outcomes and the
    joint model parameters for multiple objectives.

    Args:
        models: list of fitted SingleTaskGP (one per objective)
        X_candidates: (n_candidates, d) tensor
        n_mc: Monte Carlo samples

    Returns:
        (n_candidates,) array of multi-objective BALD scores.
    """
    n_candidates = len(X_candidates)
    n_objectives = len(models)

    # Compute BALD for each objective independently
    per_objective_bald = np.zeros((n_objectives, n_candidates))
    for i, model in enumerate(models):
        model.eval()
        with torch.no_grad(), gpytorch.settings.fast_pred_var():
            posterior = model.posterior(X_candidates)
            epistemic_var = posterior.variance.squeeze().numpy()
            noise_var = model.likelihood.noise.item()
        per_objective_bald[i] = 0.5 * np.log(1.0 + epistemic_var / noise_var)

    # Aggregate: sum of BALD scores across objectives
    # (information about any objective is valuable)
    total_bald = per_objective_bald.sum(axis=0)

    return total_bald


# Example: potency + selectivity two-objective design
# model_potency = build_gp_surrogate(X_train, Y_potency)
# model_selectivity = build_gp_surrogate(X_train, Y_selectivity)
# scores = multi_objective_bald([model_potency, model_selectivity], X_candidates)

The multi-objective extension connects to the Pareto optimization concepts from Chapter 45. More sophisticated approaches (e.g., BoTorch's qExpectedHypervolumeImprovement) directly target improvement of the Pareto hypervolume (the volume of objective space dominated by the current Pareto front, used as a single scalar measure of solution quality), but the additive BALD heuristic (called a heuristic because summing per-objective scores ignores correlations between objectives) provides a simple and effective starting point for multi-objective experiment design.

Whether the planner targets one objective or several, the code so far assumes a software oracle that returns instant results; bridging the gap to physical instruments introduces constraints that pure simulation never encounters.

9. From Simulation to the Real Lab

The planner we built operates on a simulated oracle. Connecting it to a real laboratory requires three additional components, each covered in later chapters.

Three Integration Requirements

Hardware integration: the oracle's query method must be replaced with a call to laboratory instruments (plate readers, spectrometers, robotic arms). The self-driving lab architecture in Chapter 55 provides the abstraction layer for this integration, mapping the software query(x) call to physical actions.

Provenance and reproducibility: every experiment must be logged with sufficient detail to reproduce the measurement. The experiment registry in Chapter 47 provides version-controlled storage for experimental conditions, raw data, processed results, and the acquisition rationale (which acquisition function selected this experiment and with what score).

Human-in-the-loop oversight: even in autonomous systems, domain experts should review proposed experiments before execution. The VOI stopping criterion, where VOI (value of information) quantifies whether the expected gain from one more experiment exceeds its cost, provides a natural checkpoint: when expected value of sample information (EVSI) drops below a threshold, the system pauses and requests human review. The responsible AI considerations in Chapter 57 discuss how to design these human oversight mechanisms.

Exercises

1. (Coding) Extend the ClosedLoopPlanner to support batch acquisition with batch size \(b = 3\). Use the diversity-promoting selection from BALDDesigner.propose_experiments. Run 5 rounds of batch-3 acquisition (15 total experiments) and compare the final RMSE to 15 rounds of single-acquisition BALD. How much efficiency is lost by batching?
2. (Coding) Replace the GP surrogate with a neural network using MC dropout (Listing 46.7). Run the closed-loop planner for 10 rounds and compare RMSE convergence to the GP-based planner. At what dataset size does the neural network surrogate begin to outperform the GP?
3. (Analysis) The BALD score for GPs is monotonically related to posterior standard deviation (Key Insight above). Verify this empirically by computing both quantities for 1000 candidate points after each round and plotting one against the other. At which round does the monotonic relationship break down (if ever)?
4. (Research) Implement a VOI-based stopping criterion for the planner. Define a decision problem: "classify the EC50 as above or below 10 micromolar." Compute EVSI at each round (Listing 46.12) and stop when EVSI drops below the per-experiment cost (\$150). On average, how many experiments does the VOI criterion save compared to a fixed budget of 10 rounds?
5. (Integration) Connect the BALDDesigner to the Discovery Workbench by implementing a REST API endpoint that accepts POST requests with new observations and returns proposed experiments. Use FastAPI and include endpoints for /propose, /update, and /status. Test the endpoint with curl commands that simulate a 5-round experimental campaign.