mixle.doe package

Design and analysis of computer experiments for mixle.

This package covers the full loop of reasoning about an expensive black-box model f(x) over a bounded input space:

  • Designs – space-filling (Latin hypercube, maximin, Sobol’/Halton, MaxPro maximum-projection), classical factorial / fractional-factorial / Plackett-Burman / response-surface / mixture, and optimal designs (D/A/I/G/E/c criteria).

  • Bayesian optimization – the single-point acquisitions (EI/PI/UCB/Thompson/knowledge-gradient), plus the research frontier: rigorous Monte-Carlo q-EI and local-penalization batch design, Max-value Entropy Search, trust-region BO (TuRBO) for high dimensions, constrained and multi-objective BO, and cost-aware multi-fidelity BO.

  • Active learning & optimal design – sequential design to learn a surrogate (ALM / ALC-IMSE) or model parameters (expected information gain, closed-form and nested-Monte-Carlo).

  • Analysis – global sensitivity (Sobol’/Morris/FAST/DGSM), forward uncertainty propagation, and Kennedy-O’Hagan calibration to field data. (These were the standalone mixle.uq package; folded in here they share this package’s quasi-Monte-Carlo sampling, GP surrogate, and kernels.)

The space-filling and classical design generators all return a plain (n, d) numpy matrix of input points scaled into the supplied per-dimension bounds:

>>> from mixle.doe import latin_hypercube
>>> x = latin_hypercube([(0.0, 1.0), (-2.0, 2.0)], n=8, seed=0)
>>> x.shape
(8, 2)
full_factorial(bounds, levels)[source]

Return a full-factorial grid design over bounds.

levels is the number of evenly-spaced levels per dimension (a scalar applied to all dimensions, or one entry per dimension, each >= 1). The result has prod(levels) rows in row-major (last axis fastest) order. A dimension with one level is placed at its midpoint.

Parameters:
Return type:

ndarray

halton_design(bounds, n, seed=None, *, scramble=True)[source]

Return an n-point Halton low-discrepancy design over bounds.

Halton is a quasi-random sequence (coprime van der Corput sequences per axis) that, like Sobol’, fills space more evenly than iid uniform sampling and works for any n (no power-of-two preference). scramble=True randomizes it reproducibly from seed.

Parameters:
Return type:

ndarray

latin_hypercube(bounds, n, seed=None, *, center=False)[source]

Return an n-point Latin hypercube design over bounds.

Each of the d axes is partitioned into n equal strata and exactly one sample falls in each stratum (the defining LHS property), with the per-axis stratum assignments independently permuted. With center=True each sample sits at its stratum midpoint; otherwise it is drawn uniformly within the stratum.

Parameters:
Return type:

ndarray

maximin_latin_hypercube(bounds, n, seed=None, *, trials=32)[source]

Return the best of trials Latin hypercube designs by the maximin criterion.

Generates trials independent LHS designs and keeps the one whose minimum pairwise (Euclidean, bound-normalized) distance is largest, a simple way to improve space-fillingness.

Parameters:
Return type:

ndarray

maxpro_design(bounds, n, seed=None, *, restarts=3, swaps=1500, maxiter=300)[source]

Return a maximum-projection (MaxPro) space-filling design (Joseph, Gul & Ba 2015).

MaxPro minimizes sum_{i<j} 1 / prod_k (x_ik - x_jk)^2, which forces the points to spread out in every subspace projection, not just the full space – the design of choice when only some inputs turn out to matter (factor screening / sloppy models), since the points stay space-filling even projected onto the active subset.

Built by the paper’s two-stage procedure, repeated from restarts Latin-hypercube starts: (1) a MaxProLHD coordinate-exchange phase (swaps within-column swaps) finds a good basin while keeping the LHS structure, then (2) continuous L-BFGS-B optimization (analytic gradient, maxiter steps) moves the points off the grid to the true MaxPro optimum – the refinement that cuts the criterion by far more than the swap phase alone. The best design over restarts is returned as an (n, d) array.

Parameters:
Return type:

ndarray

random_design(bounds, n, seed=None)[source]

Return n iid uniform points over bounds as an (n, d) array.

Parameters:
Return type:

ndarray

sobol_design(bounds, n, seed=None, *, scramble=True)[source]

Return an n-point Sobol’ low-discrepancy design over bounds.

Sobol’ is a quasi-random sequence whose discrepancy (deviation from perfectly even coverage) is far lower than iid uniform sampling, so it fills the space more evenly for moderate n. scramble=True applies Owen scrambling, which randomizes the sequence reproducibly from seed while preserving the low-discrepancy structure; scramble=False returns the raw deterministic sequence. Balance properties are best when n is a power of two.

Parameters:
Return type:

ndarray

fractional_factorial(bounds, generators, *, coded=False)[source]

Two-level fractional factorial 2**(k-p) from pyDOE-style generator strings.

generators names one column per factor as a product of base factors, e.g. "a b c ab ac" – a 2**(5-2) design whose factors d, e are deliberately aliased as d = ab, e = ac. Base factors are the distinct single letters; each token is the elementwise product of its letters, optionally negated with a leading -. The number of tokens must equal len(bounds).

Returns a (2**k, d) design (k = number of base factors) mapped into bounds – or the raw coded +/-1 matrix if coded=True.

Parameters:
Return type:

ndarray

plackett_burman(bounds, *, coded=False)[source]

Plackett-Burman two-level screening design for len(bounds) factors.

Returns N runs where N is the smallest multiple of four that is at least d + 1 (so the design is saturated or near-saturated, ideal for screening many factors in few runs). For N a power of two the design is a Hadamard matrix (a resolution-III fractional factorial); for N in {12, 20, 24} a known cyclic generator is used; otherwise N is rounded up to the next power of two so a design always exists. Main effects are mutually orthogonal but aliased with two-factor interactions – use it to find the few large effects, then follow up with a fuller design.

Parameters:
Return type:

ndarray

central_composite(bounds, *, center=4, alpha='rotatable', coded=False)[source]

Central-composite design (CCD) for fitting a full second-order response surface.

A CCD stacks three parts: the 2**d two-level factorial corners (estimate linear and interaction terms), 2*d axial / star points at distance alpha on each axis (estimate the pure-quadratic curvature), and center replicates at the centre (estimate pure error and curvature). alpha sets the axial distance:

  • "rotatable" (default) – alpha = (2**d)**0.25, so prediction variance depends only on distance from the centre;

  • "orthogonal" – the value making the second-order terms orthogonal (depends on center);

  • "face" (a face-centred CCD / CCF) – alpha = 1, keeping every run inside the cube;

  • a positive float – used directly.

Returns (2**d + 2*d + center, d) rows mapped into bounds (or coded if coded=True).

Parameters:
Return type:

ndarray

box_behnken(bounds, *, center=None, coded=False)[source]

Box-Behnken response-surface design (3 levels per factor, no corner runs).

For every pair of factors it runs the four (+/-1, +/-1) combinations with all other factors at the centre, plus center centre replicates. Unlike a CCD it never sets all factors to an extreme at once (no cube corners), which is useful when those combinations are expensive or infeasible, and it needs only three levels per factor. Requires d >= 3.

Returns (4 * C(d, 2) + center, d) rows mapped into bounds (or coded if coded=True).

Parameters:
Return type:

ndarray

simplex_lattice(q, m)[source]

{q, m} simplex-lattice design for q mixture components.

Each component takes one of the m + 1 evenly spaced proportions 0, 1/m, ..., 1 and the proportions in a blend sum to one. The design is exactly the set of such blends – C(q+m-1, m) points – and supports a degree-m Scheffe mixture polynomial.

Returns a (C(q+m-1, m), q) array of blends (rows sum to 1).

Parameters:
Return type:

ndarray

simplex_centroid(q)[source]

Simplex-centroid design for q mixture components.

Runs the centroid of every non-empty subset of components: the q pure components, the C(q,2) binary 1/2:1/2 blends, … up to the overall centroid (all components at 1/q) – 2**q - 1 blends in total. Supports the special cubic mixture model.

Returns a (2**q - 1, q) array of blends (rows sum to 1).

Parameters:

q (int)

Return type:

ndarray

to_pseudocomponents(blends, lower)[source]

Map canonical simplex blends onto pseudo-components with per-component lower bounds.

With lower bounds l_i (summing to < 1), a constrained mixture’s feasible region is itself a smaller simplex; this maps a canonical blend x onto the real proportions a_i = l_i + (1 - sum l) * x_i (the standard L-pseudocomponent transform), so any simplex design above can be run inside the constrained region. Rows still sum to 1.

Parameters:
Return type:

ndarray

factorial_effects(design, y, *, interactions=True, coded=False)[source]

Estimate main effects and two-factor interactions from a two-level design.

Fits the linear model y ~ 1 + x_i (+ x_i x_j) in coded +/-1 units by least squares; the coefficients are half the classical effects. design is the (n, d) run matrix (the real factor levels, coded to +/-1 automatically – or pass coded=True if it is already +/-1), y the measured responses. Set interactions=False for a main-effects-only (e.g. screening) fit.

Parameters:
Return type:

FactorialEffects

class FactorialEffects(terms, coef, effects, intercept, residual_std)[source]

Bases: object

Estimated effects from a two-level factorial / fractional-factorial / Plackett-Burman design.

Parameters:
terms

term names ("intercept", "x0", "x0:x1", …).

Type:

list[str]

coef

least-squares regression coefficients in coded +/-1 units.

Type:

numpy.ndarray

effects

the classical effect per term – the change in mean response as a factor moves from its low to its high level, i.e. 2 * coef (the intercept entry is just the grand mean).

Type:

numpy.ndarray

intercept

the grand mean of the response.

Type:

float

residual_std

residual standard deviation when the design has spare runs (else None).

Type:

float | None

terms: list[str]
coef: ndarray
effects: ndarray
intercept: float
residual_std: float | None
as_dict()[source]

Map each non-intercept term to its effect.

Return type:

dict[str, float]

response_surface(x, y)[source]

Fit a full second-order (quadratic) response surface and analyse its stationary point.

Least-squares-fits y = b0 + sum b_i x_i + sum_{i<=j} b_{ij} x_i x_j to the design runs x (n, d) and responses y, then solves for the stationary point x* = -1/2 B^{-1} b and classifies it from the eigenvalues of the quadratic matrix B. Fit on the coded design for a well-conditioned model (the classic central-composite / Box-Behnken workflow).

Return type:

ResponseSurface

class ResponseSurface(coef, terms, b, B, stationary_point, eigenvalues, kind, residual_std)[source]

Bases: object

A fitted second-order response surface y = b0 + b'x + x'Bx and its canonical analysis.

Parameters:
coef

full coefficient vector (intercept, linears, then the upper-triangular second-order terms).

Type:

numpy.ndarray

terms

matching term names.

Type:

list[str]

b

linear coefficient vector (d,).

Type:

numpy.ndarray

B

symmetric (d, d) matrix of quadratic coefficients (cross terms split onto both halves).

Type:

numpy.ndarray

stationary_point

x* solving grad = b + 2 B x = 0 (least-squares if B is singular).

Type:

numpy.ndarray

eigenvalues

eigenvalues of B – all negative => the stationary point is a maximum, all positive => a minimum, mixed signs => a saddle (a ridge if some are ~0).

Type:

numpy.ndarray

kind

"maximum" / "minimum" / "saddle".

Type:

str

residual_std

residual standard deviation when the design has spare runs (else None).

Type:

float | None

coef: ndarray
terms: list[str]
b: ndarray
B: ndarray
stationary_point: ndarray
eigenvalues: ndarray
kind: str
residual_std: float | None
predict(x)[source]

Predict the response at points x (m, d) from the fitted surface.

Return type:

ndarray

gradient(x)[source]

The response gradient b + 2 B x at x – its direction is the path of steepest ascent.

Return type:

ndarray

design_diagnostics(design, model, *, ref=None)[source]

Quality diagnostics for a design under a model – “is this a good design to run?”.

Builds the model matrix F = model(design) and reports, relative to a hypothetical perfectly orthogonal design (where each is 1.0):

  • d_efficiency – det(M)**(1/p) / n: overall coefficient-estimation precision;

  • a_efficiency – p / (n * trace(M^-1)): average coefficient variance;

  • g_efficiency – p / (n * max prediction variance) over ref (or the design itself);

  • condition_number of M (large => near-collinear / fragile to fit);

  • max_correlation – the largest absolute pairwise correlation among the non-intercept model columns (the aliasing check; 0 for an orthogonal design).

model is a model-matrix function such as mixle.doe.optimal.polynomial_features(). Use it on the coded design for meaningful efficiencies.

Return type:

dict

class OptimizationResult(x, y)[source]

Bases: object

Common outcome of a model-based optimization run: the full evaluation history.

x is the (N, d) matrix of evaluated points and y the corresponding objective values (an (N,) vector for single-objective runs, an (N, M) matrix for multi-objective). Concrete result types extend this with their best-point / Pareto-front fields.

Parameters:
x: ndarray
y: ndarray
class BayesOptResult(x, y, best_x, best_y)[source]

Bases: OptimizationResult

Outcome of a Bayesian-optimization run.

Parameters:
best_x: ndarray
best_y: float
expected_improvement(mean, std, best, xi=0.0, *, maximize=False, **_)[source]

Return the expected-improvement acquisition at points with surrogate mean and std.

For minimization the improvement over the incumbent best is best - mean - xi; for maximization it is mean - best - xi. xi >= 0 trades exploration for exploitation. Points with zero predictive std get zero EI. Higher is better (maximized over candidates).

Parameters:
Return type:

ndarray

knowledge_gradient(mean, cov, noise=1.0e-6)[source]

Knowledge-gradient acquisition value of one observation at each candidate (Frazier et al. 2009).

Given the Gaussian-process posterior mean and joint cov over a candidate set, returns, for each candidate x, the expected increase in the best posterior mean after fantasizing one (noisy) observation there: KG(x) = E_y[max_x' mu_{n+1}(x')] - max_x' mu_n(x'). Maximizing KG is the one-step Bayes-optimal, look-ahead choice (it values information, not just immediate improvement), so it explores where an observation would most change the believed optimum. Assumes a maximization objective; negate mean for minimization. Computed exactly via the piecewise-linear epigraph of the fantasized posterior means, and >= 0 by construction.

Parameters:
Return type:

ndarray

propose_knowledge_gradient(x, y, bounds, *, maximize=False, n_candidates=512, seed=None, gp=None, fit_kwargs=None, noise=1.0e-6)[source]

Propose the next evaluation point by maximizing the knowledge gradient over a candidate set.

Fits the GP surrogate to (x, y), draws n_candidates Latin-hypercube points, evaluates the joint posterior, and returns the candidate with the largest knowledge_gradient() – the look-ahead Bayesian-optimization proposal. maximize selects the objective sense (the mean is negated for minimization).

Parameters:
Return type:

ndarray

log_expected_improvement(mean, std, best, xi=0.0, *, maximize=False, **_)[source]

Return the log expected-improvement acquisition – a numerically stable EI (Ament et al. 2023).

Mathematically log(EI), but computed so it stays finite and informative deep in the no-improvement tail where EI itself underflows to 0 (and log EI to -inf), keeping the optimizer’s ordering and gradients usable. Same argmax as expected_improvement(); points with zero predictive std get -inf. Higher is better. The z >= 0 branch is the direct well-conditioned form; the z < 0 tail uses the scaled complementary error function erfcx (Phi(z)/phi(z) = sqrt(pi/2) erfcx(-z/sqrt2)), which is bounded there and so never under/overflows.

Parameters:
Return type:

ndarray

probability_of_improvement(mean, std, best, *, maximize=False, xi=0.0, **_)[source]

Return the probability-of-improvement acquisition.

The probability that a candidate improves on the incumbent best by at least xi: P(f < best - xi) for minimization, P(f > best + xi) for maximization. Where the predictive std is zero the improvement is deterministic (1.0 if it improves, else 0.0). Higher is better.

Parameters:
Return type:

ndarray

upper_confidence_bound(mean, std, best=0.0, *, maximize=False, kappa=1.96, **_)[source]

Return the confidence-bound acquisition (UCB for maximization, LCB for minimization).

Returns a merit that is maximized over candidates: the optimistic bound mean + kappa * std when maximizing, and kappa * std - mean when minimizing (so picking the largest merit selects the most promising low-objective point). kappa >= 0 trades exploration for exploitation; best is ignored. Higher is better.

Parameters:
Return type:

ndarray

thompson_sampling(mean, std, best=0.0, *, maximize=False, rng=None, **_)[source]

Thompson-sampling acquisition: one marginal posterior draw N(mean, std) per candidate.

Returns a merit (maximized over candidates) equal to the drawn value when maximizing and its negation when minimizing, so the selected point is the optimum of the sampled objective. A randomized, exploration-aware acquisition – repeated proposals explore competing optima in proportion to posterior probability, with no exploration knob to tune. This is the cheap marginal variant (independent per-candidate draws, ignoring the GP’s cross-candidate correlation); pass an rng (via acq_kwargs) for reproducible proposals. best is ignored.

Parameters:
Return type:

ndarray

register_acquisition(name, fn, aliases=())[source]

Register an acquisition fn under name (and any aliases).

fn is called as fn(mean, std, best, *, maximize, **params) and must return a merit array that the proposal loop maximizes over candidates. This is the extension point for new acquisitions – registering is all that is needed, no edits to propose_next/minimize.

Parameters:

  • name (str)

  • fn (Acquisition)

  • aliases (tuple[str, ...])

Return type:

None

available_acquisitions()[source]

Return the sorted names (and aliases) of all registered acquisitions.

Return type:

list[str]

minimize(objective, bounds, n_init=5, n_iter=15, seed=None, *, maximize=False, xi=0.0, acq='ei', acq_kwargs=None, n_candidates=512, fit_kwargs=None)[source]

Run sequential GP Bayesian optimization of a scalar objective over bounds.

Seeds with an n_init-point Latin-hypercube design, then runs n_iter acquisition-driven steps using acq ("ei" by default; also "pi" / "ucb" or any registered acquisition). Minimizes by default; set maximize=True to maximize. objective takes a (d,) point and returns a float.

Parameters:
Return type:

BayesOptResult

propose_next(x, y, bounds, n_candidates=512, seed=None, *, maximize=False, xi=0.0, acq='ei', acq_kwargs=None, gp=None, fit_kwargs=None, return_acquisition=False)[source]

Propose the next point to evaluate by maximizing an acquisition function.

Fits a GP to (x, y), scores n_candidates Latin-hypercube points by the acq acquisition ("ei" / "pi" / "ucb" or any registered name / callable), and returns the best candidate (a (d,) array), optionally with its merit. xi is forwarded to acquisitions that use it (EI, PI); per-acquisition parameters such as kappa go in acq_kwargs.

Parameters:
Return type:

ndarray | tuple[ndarray, float]

propose_batch(x, y, bounds, q, n_candidates=512, seed=None, *, maximize=False, xi=0.0, acq='ei', acq_kwargs=None, fit_kwargs=None)[source]

Propose a batch of q points to evaluate together, via the kriging-believer heuristic.

Each pick maximizes the acquisition; the GP posterior mean at the chosen point is then appended as a fantasized observation and the surrogate refit, so the next pick is steered away from it. Returns a (q, d) array. This needs no true objective evaluations between picks, so it suits parallel/asynchronous experiment campaigns.

Parameters:
Return type:

ndarray

optimal_design(bounds, n, *, candidates=None, model=None, criterion='D', n_candidates=256, n_restarts=5, max_iter=100, ref=None, seed=None)[source]

Return an n-point optimal design selected from a candidate pool by Fedorov exchange.

The pool is either generated as a Sobol design of n_candidates points over per-dimension bounds, or supplied directly as an (P, d) candidates array (in which case bounds may be None). model is a model-matrix function (default polynomial_features() degree 1, i.e. linear with intercept); criterion selects the optimality merit ("D" / "A" / "I" or any registered name / callable). The best design over n_restarts random starts is returned as an (n, d) array. For "I" optimality, prediction variance is averaged over ref (a model matrix) when given, else over the candidate pool.

Raises if n is below the number of model parameters (the information matrix would be singular).

Parameters:
Return type:

ndarray

polynomial_features(degree=1, *, bias=True)[source]

Return a model-matrix function building polynomial features up to degree (with interactions).

The returned f(X) maps an (n, d) point array to an (n, p) model matrix: an optional intercept column, then every monomial prod(x_j) over index multisets of size 1..degree (so degree=1 is linear, degree=2 is a full quadratic response surface including the cross terms x_i x_j).

Parameters:
Return type:

Callable[[ndarray], ndarray]

d_criterion(info, *, ref=None)[source]

D-optimality merit: log det M (-inf if M is singular). Higher is better.

Parameters:
Return type:

float

a_criterion(info, *, ref=None)[source]

A-optimality merit: -trace(M^{-1}) (-inf if singular). Higher is better.

Parameters:
Return type:

float

i_criterion(info, *, ref=None)[source]

I-optimality merit: -mean prediction variance over ref (-inf if singular).

The prediction variance at a reference row g is g M^{-1} g; this returns the negative mean over the reference model matrix ref so larger is better. Falls back to A-optimality when no reference set is supplied.

Parameters:
Return type:

float

g_criterion(info, *, ref=None)[source]

G-optimality merit: -max prediction variance over ref (-inf if singular).

Where I-optimality minimizes the average prediction variance, G-optimality minimizes its worst case over the reference region, so the fitted surface has a bounded error everywhere. Returns the negative maximum of g M^{-1} g over the reference rows (falls back to the largest coefficient variance max diag(M^{-1}) when no reference set is given).

Parameters:
Return type:

float

e_criterion(info, *, ref=None)[source]

E-optimality merit: the smallest eigenvalue of M (higher is better).

Maximizing the minimum eigenvalue of the information matrix shrinks the variance along the worst-determined parameter contrast, so no direction in coefficient space is left poorly estimated. 0 for a singular (rank-deficient) design.

Parameters:
Return type:

float

c_criterion(c)[source]

Return a c-optimality criterion targeting the linear combination c'.beta of coefficients.

c-optimality minimizes the variance of a specific quantity of interest c'.beta (e.g. a contrast or a prediction at one point). The returned criterion has merit -c' M^{-1} c (-inf if singular), so it plugs straight into optimal_design() or register_criterion().

Parameters:

c (ndarray)

Return type:

Criterion

register_criterion(name, fn, aliases=())[source]

Register an optimality criterion fn under name (and any aliases).

fn is called as fn(info, *, ref) with the information matrix M = F.T @ F and must return a merit that optimal_design() maximizes. This is the extension point for new criteria – registering is all that is needed, no edits to the exchange loop.

Parameters:
Return type:

None

available_criteria()[source]

Return the sorted names (and aliases) of all registered optimality criteria.

Return type:

list[str]

class ConstrainedBayesOptResult(x, y, best_x, best_y, c, feasible)[source]

Bases: BayesOptResult

Outcome of a constrained Bayesian-optimization run.

c holds the (N, K) observed constraint values (feasible rows have all entries <= 0) and feasible is the corresponding boolean mask. best_x / best_y are the best feasible point; if no feasible point was found they fall back to the least-infeasible observation.

Parameters:
c: ndarray
feasible: ndarray
probability_of_feasibility(mean, std)[source]

Return the per-point probability that all constraints are satisfied (c_k <= 0).

mean and std are (n, K) posterior predictive moments of the K constraint surrogates. Returns an (n,) array, the product over constraints of Phi(-mean_k / std_k). Where a constraint’s std is zero the feasibility is deterministic (1.0 if mean <= 0).

Parameters:
Return type:

ndarray

propose_next_constrained(x, y, c, bounds, n_candidates=512, seed=None, *, maximize=False, xi=0.0, acq='ei', acq_kwargs=None, fit_kwargs=None, return_acquisition=False)[source]

Propose the next point under inequality constraints c_k(x) <= 0.

Fits a GP to the objective (x, y) and one GP per constraint column of c (an (N, K) array), then maximizes the feasibility-weighted acquisition over n_candidates Latin-hypercube points. Until a feasible observation exists the acquisition factor is held at 1 so the search targets feasibility; afterwards the incumbent is the best feasible objective. Returns the chosen (d,) point, optionally with its merit.

Parameters:
Return type:

ndarray | tuple[ndarray, float]

constrained_minimize(objective, constraints, bounds, n_init=5, n_iter=15, seed=None, *, maximize=False, xi=0.0, acq='ei', acq_kwargs=None, n_candidates=512, fit_kwargs=None)[source]

Constrained GP Bayesian optimization of objective subject to constraints over bounds.

Each callable in constraints maps a (d,) point to a scalar that is feasible when <= 0. Seeds with an n_init-point Latin-hypercube design, then runs n_iter feasibility-weighted acquisition steps. Minimizes the objective by default; returns the best feasible point (or the least-infeasible one if none feasible) along with the full evaluation history.

Parameters:
Return type:

ConstrainedBayesOptResult

class MultiObjectiveResult(x, y, pareto_mask, pareto_x, pareto_y)[source]

Bases: OptimizationResult

Outcome of a multi-objective Bayesian-optimization run.

y is the (N, M) matrix of observed objective vectors (all minimized); pareto_mask flags the non-dominated rows, and pareto_x / pareto_y are those points and their objective vectors.

Parameters:
pareto_mask: ndarray
pareto_x: ndarray
pareto_y: ndarray
pareto_mask(y)[source]

Return a boolean mask of the non-dominated rows of y (an (N, M) minimization objective).

Row i is dominated when some other row is <= it on every objective and strictly < on at least one; the mask is True for the rows that survive (the Pareto-optimal set).

Parameters:

y (Any)

Return type:

ndarray

multi_minimize(objectives, bounds, n_init=10, n_iter=20, seed=None, *, rho=0.05, n_candidates=512, fit_kwargs=None)[source]

Multi-objective GP Bayesian optimization of objectives over bounds (ParEGO).

Each callable in objectives maps a (d,) point to a scalar; all are minimized. Seeds with an n_init-point Latin-hypercube design, then runs n_iter steps, each drawing a random Tchebycheff weighting, scalarizing the observed objectives, and taking one GP-EI step on that scalar. Returns the full evaluation history and the Pareto-optimal subset.

Parameters:
Return type:

MultiObjectiveResult

class BayesianOptimizer(bounds, *, acq='ei', acq_kwargs=None, maximize=False, n_init=None, xi=0.0, n_candidates=512, fit_kwargs=None, seed=None)[source]

Bases: object

Stateful ask-tell wrapper around the GP Bayesian-optimization loop.

ask proposes the next point (or a batch), tell records evaluated observations, and best returns the incumbent. Minimizes by default; set maximize=True to maximize. The acquisition is selected by acq ("ei" / "pi" / "ucb" or any registered name / callable) with per-acquisition parameters in acq_kwargs.

Parameters:

  • bounds (Bounds)

  • acq (str | Acquisition)

  • acq_kwargs (dict[str, Any] | None)

  • maximize (bool)

  • n_init (int | None)

  • xi (float)

  • n_candidates (int)

  • fit_kwargs (dict[str, Any] | None)

  • seed (int | RandomState | None)

property x: ndarray

Return the observed points as an (N, d) array.

property y: ndarray

Return the observed objective values as an (N,) array.

property n_observations: int

Return the number of recorded observations.

property best: BayesOptResult

Return the incumbent (best observed point) as a BayesOptResult.

ask(q=1)[source]

Return the next point to evaluate as a (d,) array, or a (q, d) batch when q > 1.

The first n_init points come from a space-filling design; subsequent points are GP acquisition proposals (a kriging-believer batch when q > 1).

Parameters:

q (int)

Return type:

ndarray

tell(x, y)[source]

Record one or more evaluated observations; returns self for chaining.

x is a (d,) point or (m, d) batch and y the matching scalar or (m,) values.

Parameters:
Return type:

BayesianOptimizer

monte_carlo_qei(mean, cov, best, *, maximize=False, samples=512, seed=0)[source]

Monte-Carlo multi-point Expected Improvement of a batch with joint posterior N(mean, cov).

Draws samples joint posterior realizations of the q batch points and averages the batch improvement over the incumbent best – max(best - min_i f_i, 0) for minimization, or max(max_i f_i - best, 0) for maximization. For q = 1 this reduces to ordinary EI.

Parameters:
Return type:

float

propose_qei_batch(x, y, bounds, q, *, n_candidates=256, mc_samples=256, maximize=False, seed=None, gp=None, fit_kwargs=None)[source]

Propose a q-point batch by greedy Monte-Carlo q-EI under the joint GP posterior.

Fits the GP to (x, y), then builds the batch one point at a time: each new point is the Latin-hypercube candidate maximizing the q-EI of {batch so far} + candidate (evaluated with common random numbers so the greedy comparison is fair). Because the joint posterior is used, an already-chosen point lowers the marginal value of nearby candidates, so the batch self-diversifies without any fantasized observations. Returns a (q, d) array.

Parameters:
Return type:

ndarray

propose_local_penalization(x, y, bounds, q, *, n_candidates=512, maximize=False, seed=None, gp=None, fit_kwargs=None)[source]

Propose a q-point batch by local penalization (Gonzalez et al. 2016).

Picks points sequentially from a single GP fit (no refitting): each pick maximizes the expected improvement multiplied by a soft exclusion factor around every pending pick. The exclusion radius is set from a Lipschitz estimate L of the objective (the largest posterior-mean gradient over the candidates) and the gap to the incumbent, so the penalty is principled rather than a fixed distance. Cheaper than q-EI for large q (one GP fit, closed-form penalties). Returns a (q, d) array.

Parameters:
Return type:

ndarray

max_value_entropy_search(mean, std, max_samples, *, maximize=True)[source]

Max-value Entropy Search acquisition at candidates with posterior mean / std.

Given samples max_samples of the global optimum value y*, returns the per-candidate mutual information I(y; y*) = (1/M) sum_m [ gamma_m phi(gamma_m)/(2 Phi(gamma_m)) - log Phi(gamma_m) ] with gamma_m = (y*_m - mu)/sd (maximization; for minimization the sense is flipped by the caller). Higher is better – it favors uncertain candidates near the believed optimum.

Parameters:
Return type:

ndarray

sample_max_values(mean, std, n_samples=64, *, seed=0)[source]

Sample plausible global-max values y* via the Gumbel approximation (Wang & Jegelka 2017).

The CDF of the maximum over the candidate cloud is P(max <= y) = prod_i Phi((y - mu_i)/sd_i); a Gumbel is fit to its 25/50/75 percentiles (found by bisection) and sampled. Returns an (n_samples,) array of y* draws (never below the best posterior mean).

Parameters:
Return type:

ndarray

propose_mes(x, y, bounds, *, n_candidates=512, max_samples=64, maximize=False, seed=None, gp=None, fit_kwargs=None)[source]

Propose the next point by Max-value Entropy Search.

Fits the GP to (x, y), samples max_samples optimum values via sample_max_values(), scores Latin-hypercube candidates by max_value_entropy_search(), and returns the maximizer as a (d,) array. maximize selects the optimization sense (default minimize, matching the rest of the BO layer).

Parameters:
Return type:

ndarray

turbo_minimize(objective, bounds, *, n_init=None, max_evals=100, batch_size=1, maximize=False, n_candidates=None, seed=None, fit_kwargs=None)[source]

Optimize a black-box objective over bounds with TuRBO (trust-region BO).

Starts from a Latin-hypercube design of n_init points (default 2*d), then repeatedly fits a GP on the normalized data, draws Thompson candidates inside the current trust region, evaluates the best batch_size of them, and resizes the region by success/failure. On collapse it restarts from a new design. Runs until max_evals objective calls. Returns {'x', 'y', 'X', 'Y', 'n_restarts'} with the best point/value and the full history.

Parameters:
Return type:

dict[str, Any]

class TrustRegion(dim, length=0.8, length_min=0.0078125, length_max=1.6, success_tol=3, failure_tol=0, _success=0, _failure=0)[source]

Bases: object

Expand/shrink state of a TuRBO trust region (side length in normalized [0, 1]^d coordinates).

The length doubles after success_tol consecutive improving batches and halves after failure_tol consecutive non-improving ones; collapsed is True once it drops below length_min and the caller should restart. failure_tol defaults to the dimension (Eriksson 2019).

Parameters:
dim: int
length: float = 0.8
length_min: float = 0.0078125
length_max: float = 1.6
success_tol: int = 3
failure_tol: int = 0
property collapsed: bool
update(improved)[source]

Record whether the latest batch improved the incumbent, and resize the region.

Parameters:

improved (bool)

Return type:

None

alm_scores(gp, x, y, candidates)[source]

Active Learning MacKay scores: the GP posterior predictive variance at each candidate.

Parameters:
Return type:

ndarray

alc_scores(gp, x, y, candidates, reference)[source]

Active Learning Cohn / IMSE scores: integrated posterior-variance reduction per candidate.

Adding candidate c reduces the posterior variance at a reference point r by cov_post(r, c)^2 / var_post(c); this returns the sum over the reference set (the negative change in integrated MSE), so the maximizer is the most globally informative next point.

Parameters:
Return type:

ndarray

propose_active_learning(x, y, bounds, *, method='alc', n_candidates=512, n_reference=256, seed=None, gp=None, fit_kwargs=None)[source]

Propose the next active-learning point (method='alc' IMSE, or 'alm' max variance).

Parameters:
Return type:

ndarray

active_learning_design(objective, bounds, *, n_init=None, max_evals=40, method='alc', seed=None, fit_kwargs=None)[source]

Sequentially place points to maximize surrogate accuracy, returning the design and its responses.

Starts from a Latin-hypercube design, then repeatedly fits the GP and adds the most informative point (by method) until max_evals evaluations. Returns {'X', 'Y'} – a design tailored to learn objective well everywhere, not just near an optimum.

Parameters:
Return type:

dict[str, Any]

expected_information_gain_linear(model_matrix, *, noise=1.0, prior_cov=None)[source]

Exact expected information gain of a linear-Gaussian design y = F.theta + eps (= Bayesian D-opt).

For a Gaussian prior theta ~ N(0, Sigma0) and observation noise eps ~ N(0, noise^2 I), the mutual information between the data and theta is 0.5 * log det(I + noise^-2 Sigma0 F^T F). model_matrix is the (n, p) design matrix F (e.g. from mixle.doe.optimal.polynomial_features()). Higher EIG = a more informative design.

Parameters:
Return type:

float

expected_information_gain_nmc(prior_sampler, log_likelihood, simulate, *, n_outer=256, n_inner=256, seed=None)[source]

Nested-Monte-Carlo expected information gain for a general (nonlinear) design.

Estimates EIG = E_{theta, y}[ log p(y|theta) - log E_{theta'}[p(y|theta')] ] (Ryan 2003): draw outer theta_i ~ prior and y_i ~ p(y|theta_i) via simulate; the inner expectation is a mean over n_inner prior draws of exp(log_likelihood(theta', y_i)). prior_sampler(rng, n) returns (n, k) parameter draws; log_likelihood(thetas, y) returns a log-density per row of thetas at the single observation y; simulate(theta, rng) draws one y given theta.

Parameters:
Return type:

float

multi_fidelity_minimize(objective, bounds, *, fidelities=(0.5, 1.0), costs=None, target=None, n_init=None, max_cost=40.0, n_candidates=256, maximize=False, seed=None, fit_kwargs=None)[source]

Cost-aware multi-fidelity Bayesian optimization of objective(x, s).

objective(x, s) returns the response at input x and fidelity s (one of fidelities); the largest fidelity (or target) is the true objective. costs is the per-fidelity evaluation cost (default: the fidelity value itself). The loop fits a GP over [x, s], proposes x by Expected Improvement at the target fidelity, then evaluates at the fidelity maximizing target-variance reduction per unit cost, until the cumulative cost reaches max_cost. Returns {'x', 'y', 'X', 'Y', 'cost'} – the best target-fidelity point and the full augmented history.

Parameters:
Return type:

dict[str, Any]

sobol_indices(func, bounds, n=4096, *, seed=0, names=None)[source]

First- and total-order Sobol sensitivity indices (Saltelli sampling, Jansen estimators).

Parameters:

  • func (Callable[[ndarray], ndarray]) – a vectorized model f(X) -> y mapping an (m, d) array of inputs to an (m,) array of scalar outputs.

  • bounds (Sequence[tuple[float, float]]) – [(lo, hi), ...] per input – the input is taken uniform on the box.

  • n (int) – base sample size; the total number of model evaluations is n * (d + 2).

  • seed (int) – RNG seed for the (Sobol) base samples.

  • names (Sequence[str] | None) – optional input names for the returned dict.

Returns:

{'S1': (d,), 'ST': (d,), 'names': [...], 'var': float} – first-order S1[i] is the fraction of output variance from input i alone; total-order ST[i] includes all interactions involving i (so ST[i] - S1[i] measures i’s interaction strength, and ST[i] ~ 0 means input i can be fixed).

Return type:

dict[str, Any]

morris_screening(func, bounds, *, trajectories=20, levels=4, seed=0, names=None)[source]

Morris elementary-effects screening – a cheap first factor ranking.

Walks trajectories one-factor-at-a-time paths on a levels-grid; the mean absolute elementary effect mu_star[i] ranks influence and the spread sigma[i] flags nonlinearity/interactions. Cost: trajectories * (d + 1) evaluations – far fewer than Sobol, for an initial screen.

Parameters:
Return type:

dict[str, Any]

fast_indices(func, bounds, n=600, *, harmonics=6, seed=0, names=None)[source]

First-order sensitivity indices via Random Balance Designs FAST (RBD-FAST).

A Fourier alternative to sobol_indices() for the first-order indices: every input is driven along the same triangle-wave search curve but under an independent random permutation, the model is evaluated once over the n points, and for each input the output – reordered along that input’s curve – has its variance concentrated at the base frequency’s first harmonics harmonics. The ratio of that power to the total is the first-order index (with the Tarantola bias correction). Cost is a single batch of n evaluations, independent of dimension.

Returns {'S1': (d,), 'names': [...], 'var': float}.

Parameters:
Return type:

dict[str, Any]

dgsm(func, bounds, n=1024, *, seed=0, rel_step=1.0e-4, names=None)[source]

Derivative-based global sensitivity measures (DGSM): mean squared partial derivatives.

nu[i] = E[(df/dx_i)^2] over the input box, estimated by central finite differences at n low-discrepancy points. Unlike the first-order Sobol index, a DGSM is nonzero whenever an input matters anywhere – including purely through interactions – so it is a cheap, robust screen that upper-bounds the total Sobol index (Sobol & Kucherenko, via the Poincare inequality: ST[i] <= (L_i / pi)^2 * nu[i] / Var(y) for a uniform input of width L_i). The reported importance is L_i^2 * nu[i] normalized to sum to one – a dimensionless influence ranking.

Returns {'nu': (d,), 'importance': (d,), 'names': [...]}.

Parameters:
Return type:

dict[str, Any]

propagate(func, mean, cov=None, *, n=10000, method='montecarlo', quantiles=(0.05, 0.5, 0.95), seed=0)[source]

Propagate a Gaussian input N(mean, cov) through func to output statistics.

Parameters:

  • func (Callable[[ndarray], ndarray]) – a vectorized model f(X) -> y mapping (m, d) inputs to (m,) (or (m, k)) outputs.

  • mean (ndarray) – length-d input mean. cov: (d, d) input covariance (defaults to identity).

  • n (int) – Monte Carlo sample size (ignored by the unscented method).

  • method (str) – 'montecarlo' (sample + summarize, gives quantiles) or 'unscented' (sigma-point moment propagation, mean + std only). Looked up through the propagator registry.

  • quantiles (tuple[float, ...]) – output quantiles to report (Monte Carlo only).

  • cov (ndarray | None)

  • seed (int)

Returns:

{'mean', 'std', 'quantiles' (mc only), 'samples' (mc only)}.

Return type:

dict[str, Any]

register_propagator(name)[source]

Decorator registering a propagation method under name for propagate().

The decorated callable receives (func, mean, cov, *, n, quantiles, seed) (mean/cov already coerced to float arrays) and returns the output-statistics dict.

Parameters:

name (str)

Return type:

Callable[[Callable[[…], dict[str, Any]]], Callable[[…], dict[str, Any]]]

unscented_transform(func, mean, cov, *, alpha=1e-3, beta=2.0, kappa=0.0)[source]

Propagate (mean, cov) through func with the unscented (sigma-point) transform.

Returns the output (mean, covariance). Exact for an affine func; for a nonlinear one it captures the mean and covariance to second order with only 2d+1 evaluations.

Parameters:
Return type:

tuple[ndarray, ndarray]

calibrate(simulator, x, y, theta0, *, discrepancy=True, discrepancy_lengthscale=None, seed=0, max_iter=300)[source]

Calibrate simulator(x, theta) to field data (x, y) with a GP discrepancy term.

Maximizes the marginal likelihood of the residual r(theta) = y - eta(x, theta) under a GP + noise model, over theta and the discrepancy amplitude + noise. discrepancy=False drops the GP (plain nonlinear least squares) – useful to show the bias the discrepancy removes.

The discrepancy correlation length is fixed (discrepancy_lengthscale, default 10% of the input domain) rather than fitted: this is the standard resolution of the Kennedy-O’Hagan theta/delta identifiability problem – a short discrepancy length forces the GP to model only local model-form error, leaving the smooth global trend to the parametric simulator so theta stays identifiable. Set it to the scale of model error you expect.

Parameters:

  • simulator (Callable[[ndarray, ndarray], ndarray]) – eta(x, theta) -> predictions (vectorized over the rows of x).

  • x (ndarray) – field inputs and observations.

  • y (ndarray) – field inputs and observations.

  • theta0 (Sequence[float]) – initial calibration parameters (its length sets the parameter count).

  • discrepancy (bool) – include the GP discrepancy term (the Kennedy-O’Hagan model).

  • discrepancy_lengthscale (float | None) – fixed GP correlation length (default: 10% of the input domain).

  • seed (int)

  • max_iter (int)

Return type:

KOCalibration

class KOCalibration(theta, ls, amp, noise, simulator, x, y)[source]

Bases: object

Result of calibrate(): the fitted parameters, discrepancy GP, and a calibrated predictor.

predict(x_new, *, with_discrepancy=True)[source]

Calibrated prediction at x_new: simulator at the fitted theta, plus the GP discrepancy (the bias-corrected estimate of reality) unless with_discrepancy=False (the pure simulator).

Parameters:
Return type:

ndarray

Submodules