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), 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.

Each restart’s continuous refinement is only kept if the optimizer reports success and returns a finite objective/gradient/point; otherwise that restart falls back to its swap-refined discrete design instead, which is always a legitimate (if less-refined) MaxPro candidate – see _maxpro_swap(). A restart can therefore never contribute an unverified, possibly-divergent optimizer iterate to the result (MXR-080-0175).

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

class DistillationDesign(indices, scores, task_counts=<factory>, modality_counts=<factory>, candidate_scores=<factory>, metadata=<factory>)[source]

Bases: object

Selected pool indices and diagnostics for a distillation design.

indices point into the exact candidate pool supplied to the selector. scores are the sequential merits assigned to the chosen candidates. candidate_scores preserves the base uncertainty/preference score before diversity, coverage, and cost terms are applied. metadata records the target coverage and weights needed to audit or reproduce the design.

Parameters:
distillation_design(features, n, *, task_labels=None, modalities=None, uncertainty=None, preference=None, cost=None, task_weights=None, modality_weights=None, reference_features=None, eligible=None, uncertainty_weight=1.0, diversity_weight=1.0, task_coverage_weight=1.0, modality_coverage_weight=1.0, preference_weight=1.0, cost_weight=1.0, seed=None)[source]

Select a pool subset for teacher labeling or student distillation.

features is an (N, d) embedding/feature matrix for candidate examples. task_labels and modalities may be one tag per row or a sequence of tags per row. uncertainty may be (N,) or (N, T); the latter is averaged across the task tags present on each candidate. Higher uncertainty/preference and better coverage increase merit; higher cost lowers it.

Parameters:
Return type:

DistillationDesign

multitask_distillation_design(*args, **kwargs)[source]

Alias for distillation_design() that reads naturally at multi-task call sites.

Parameters:
Return type:

DistillationDesign

cross_modal_distillation_design(modality_features, n, *, task_labels=None, uncertainty=None, preference=None, cost=None, required_modalities=None, min_modalities=2, alignment_weight=1.0, seed=None, **kwargs)[source]

Select paired cross-modal records for distillation/alignment training.

modality_features maps modality name to an (N, d_m) feature matrix. Rows with non-finite values are treated as missing for that modality. The selector fuses standardized modality features, tags each row by available modalities, and adds an alignment-disagreement preference for modality pairs that share the same embedding width. By default only rows with at least two modalities are eligible.

Parameters:
Return type:

DistillationDesign

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).

Rejects a generators spec that collapses a factor to a constant column (repeated letters cancelling entirely, e.g. "aa") or that aliases two of the named factors with each other (two tokens whose reduced letter-words are identical or exact opposites, e.g. tokens "a" and "-a") – both silently produced a degenerate design with no diagnostic before (MXR-080-0179). This does not affect the intentional aliasing a fractional design is built on: a compound token like "ab" is fine and expected, it just cannot repeat another factor’s own word. Use generator_alias_structure() to see which factors are aliased with which interaction.

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

generator_alias_structure(generators, *, max_order=2)[source]

Return signed defining words and complete alias chains through max_order.

The result is a JSON-friendly record with:

  • factor_labels mapping every public positional factor name to its source token;

  • base_factors mapping each latent generator letter to its named factor and sign;

  • generators expressing every compound generator using those public factor names;

  • defining_relations and the complete generated defining_group as signed words; and

  • alias_sets for the identity and every effect through the requested interaction order.

A word record has an integer sign (+1 or -1) and a factors list. For example, "a b ab" defines I = +x0:x1:x2; through order two, x0 is aliased with +x1:x2 and x2 with +x0:x1. Negated base or compound tokens propagate their signs through every relation and alias. This supplies the actual confounding structure rather than only repeating the compound generator fragments (MXR-080-1477).

Every latent base letter must have one named singleton token so aliases can be mapped unambiguously to the positional xN labels used by mixle.doe.analysis. The generator syntax is otherwise validated exactly like fractional_factorial().

Parameters:
Return type:

dict[str, object]

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, in this row order (see central_composite_point_kinds() to recover which output row is which): 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, finite float – used directly.

bounds defines the factorial cube (-1/+1 map to each factor’s low/high): the factorial-corner and centre rows are therefore always exactly within bounds. Axial points sit at coded distance alpha on a single axis, and for "rotatable"/"orthogonal" that distance is normally greater than 1 – a mathematical requirement of the rotatability/orthogonality property, not a clamp to bounds – so those rows routinely fall outside bounds on purpose. Only alpha="face" (or a numeric alpha <= 1) guarantees every row stays within bounds (MXR-080-0179); use central_composite_point_kinds() to tell the row kinds apart programmatically when that distinction matters to the caller.

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

Parameters:
Return type:

ndarray

central_composite_point_kinds(bounds, *, center=4)[source]

Row-kind labels for central_composite()’s output, in the same row order.

Returns a length-(2**d + 2*d + center) string array of "factorial" / "axial" / "center" entries. "factorial" and "center" rows are always exactly within bounds; "axial" rows are guaranteed within bounds only when central_composite() was called with alpha="face" (or a numeric alpha <= 1) – see its docstring for why rotatable/ orthogonal axial points routinely extend past bounds on purpose (MXR-080-0179). Pass the same center used in the matching central_composite() call so the lengths line up, e.g.:

design = central_composite(bounds, center=4, alpha="rotatable")
kinds = central_composite_point_kinds(bounds, center=4)
in_bounds_rows = design[kinds != "axial"]
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. Every row is always within bounds (no analogue of a CCD’s rotatable axial points here). 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 (finite, non-negative, summing to < 1), a constrained mixture’s feasible region is itself a smaller simplex; this maps a canonical blend x – itself a finite (n, q) array of genuine simplex points (non-negative, rows summing to 1) – 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.

Raises:

ValueError – lower is not a finite, non-negative, 1-D vector summing to less than 1 (i.e. one that itself admits at least one feasible simplex point), or blends is not a finite 2-D array of shape (n, len(lower)) whose rows are genuine canonical simplex points (non-negative, summing to 1). Invalid input is rejected outright instead of silently mapped to an out-of-simplex result (MXR-080-0180): e.g. blend [2, -1] with lower bounds [0.1, 0.1] used to silently return [1.7, -0.7], a negative “proportion”.

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 coded, in which case every entry must be exactly -1.0 or 1.0), y the measured responses. Set interactions=False for a main-effects-only (e.g. screening) fit.

Rejects designs that cannot support the requested model outright: non-finite values, fewer runs than model parameters, and (for coded=True) levels other than +/-1. A design that has enough runs but still cannot separate every term – aliasing, including a factor that never varies – is not rejected: the aliased terms come back as nan in coef / effects, with the one combined quantity the data does support reported in estimable_contrasts (see FactorialEffects). That automatic resolution only covers the classical case where the aliasing is exact term-for-term confounding, as in any standard fractional-factorial / Plackett-Burman alias structure; a design that is rank-deficient for some other reason (e.g. too few runs and not a clean fractional design) raises instead of guessing which contrasts might be estimable.

Parameters:
Return type:

FactorialEffects

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

Bases: object

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

Only estimable quantities are reported as effects. When the design cannot separate two or more terms – classical aliasing/confounding, including a factor that never varies – those individual terms have no unique estimate: their entries in coef / effects are nan. aliases says which terms are aliased with which; estimable_contrasts gives the one well-defined combined quantity the data actually supports for each such group.

Parameters:
terms

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

Type:

list[str]

coef

least-squares regression coefficients in coded +/-1 units (nan for a term that is aliased with another and so has no individual estimate).

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 (nan if the intercept itself is aliased, e.g. with a factor that never varies).

Type:

float

residual_std

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

Type:

float | None

se

standard error of each entry of effects. None throughout when the design has no spare degrees of freedom to estimate a residual variance (e.g. a saturated design); individual entries are nan wherever effects is (aliased terms).

Type:

numpy.ndarray | None

aliases

for each term, the other term names it is aliased with (empty list if the term is uniquely estimable).

Type:

dict[str, list[str]]

estimable_contrasts

for each alias group of two or more terms, keyed by its member term names joined by their signed defining relation (for example "x0-x1" for opposite columns), the (effect, se) of the exact combined quantity that is estimable – empty when the design has no aliasing.

Type:

dict[str, tuple[float, float | None]]

as_dict()[source]

Map each non-intercept term to its effect (nan where aliased – see estimable_contrasts).

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 the stationarity equation B x* = -b/2 for the stationary point and classifies it from the eigenvalues of the quadratic matrix B. When B is singular a solution only exists if -b/2 lies in its range, which is checked explicitly via the residual of the stationarity equation – a least-squares “solution” that does not actually zero the gradient is reported as kind="no_stationary_point" rather than a fake stationary point (see ResponseSurface). 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, model_rank, n_parameters, degrees_of_freedom, estimable)[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 the stationarity equation B x* = -b/2, or None when no such point exists – i.e. the gradient b + 2 B x is never zero anywhere, which can only happen when B is singular (kind == "no_stationary_point"; see below).

Type:

numpy.ndarray | None

eigenvalues

eigenvalues of B – all negative => a maximum, all positive => a minimum, mixed signs => a saddle; any eigenvalue ~0 (with a genuine stationary point) => a ridge, a whole line/subspace of equally-stationary points, of which stationary_point is one representative (the minimum-norm) point.

Type:

numpy.ndarray

kind

"maximum" / "minimum" / "saddle" / "ridge" / "no_stationary_point". The last case is distinct from a ridge: a ridge is a genuine (degenerate) stationary subspace, while "no_stationary_point" means the stationarity equation has no solution at all (e.g. the exact surface y = x0**2 + x1, whose gradient’s second component is always exactly 1).

Type:

str

residual_std

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

Type:

float | None

model_rank

rank of the fitted quadratic model matrix.

Type:

int

n_parameters

number of quadratic coefficients requested.

Type:

int

degrees_of_freedom

residual degrees of freedom, n_runs - n_parameters.

Type:

int

estimable

whether every requested coefficient is identified. Construction refuses a non-estimable fit, so returned surfaces always report True.

Type:

bool

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; None when the requested model is rank-deficient), plus a finite effective_condition_number within the estimable subspace and explicit rank / full_rank fields;

  • 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.

ref – when given – is a reference design: raw candidate points in the same (m, d) raw-factor-column form as design (e.g. a finer grid over the same region, or a wider region to probe extrapolation risk), not a pre-built model matrix. It is passed through the same model used to build the design’s own model matrix before its prediction variance is computed, so g_efficiency compares like with like; a ref whose raw column count does not match design’s raises. g_efficiency is guaranteed to fall on the documented (0, 1] scale only for the self-referential case (ref=None, or ref equal to design itself) – a reference set concentrated in a lower-variance region than the design’s own points (e.g. just its centre) can legitimately score above 1.0.

The raw design and every transformed model matrix must be nonempty, two-dimensional, finite, and preserve the run axis. A rank-zero model is rejected. A partially estimable model reports zero full-model efficiencies, condition_number=None, and its finite effective condition number rather than publishing NaN/Inf quality evidence.

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:
class BayesOptResult(x, y, best_x, best_y, n_evaluations, failed_evaluations, stopped_reason)[source]

Bases: OptimizationResult

Outcome of a Bayesian-optimization run.

Parameters:
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 (near-)zero predictive std get the DETERMINISTIC LIMIT max(improve, 0) – a point-mass posterior has no uncertainty to average over, so the “expected” improvement is exactly the guaranteed one, not 0. Higher is better (maximized over candidates).

Raises ValueError if std is negative/non-finite, mean/std shapes disagree, mean/best are non-finite, or xi is negative/non-finite.

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).

Raises ValueError if the surrogate’s predicted mean/covariance don’t match the candidate contract (wrong shape, non-finite, asymmetric covariance) or if every candidate’s knowledge-gradient value is non-finite (see _validate_prediction() / _select_index()).

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 (near-)zero predictive std get the deterministic limit log(max(improve, 0)) – -inf when the point-mass outcome does not improve on best, and the log of the guaranteed improvement when it does (see expected_improvement()). 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.

Raises ValueError if std is negative/non-finite, mean/std shapes disagree, mean/best are non-finite, or xi is negative/non-finite.

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). A NEGATIVE std is rejected outright (see _validate_acquisition_inputs()) rather than treated as this same deterministic case – it is an invalid input, not a zero-variance one. Higher is better.

Raises ValueError if std is negative/non-finite, mean/std shapes disagree, mean/best are non-finite, or xi is negative/non-finite.

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.

Raises ValueError if std is negative/non-finite, mean/std shapes disagree, mean/best are non-finite, or kappa is negative/non-finite.

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 low-cost marginal variant (independent per-candidate draws, ignoring the GP’s cross-candidate correlation). The proposal loop (_propose_one(), and so propose_next() / propose_batch() / BayesianOptimizer) automatically threads its own candidate-generation rng through to this acquisition, so a seeded caller gets reproducible draws with no extra wiring; pass rng explicitly via acq_kwargs only to override that with a different generator. Called directly with no rng at all, this draws from a fresh, unseeded generator (matching _as_rng’s own None convention) – reproducible only through the proposal loop’s threading above. best is ignored.

Raises ValueError if std is negative/non-finite, mean/std shapes disagree, mean/best are non-finite, or rng is neither None nor a numpy.random RandomState/Generator.

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.

A non-finite objective response is recorded in failed_evaluations and stops the run with stopped_reason='objective_failed'. It is never admitted to the fitted x/y history.

The GP surrogate needs the optional torch extra (pip install "mixle[torch]"). On a base install that is checked up front, before the first objective evaluation – never after the n_init design has already spent expensive black-box calls.

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.

Raises ValueError if the surrogate’s predicted mean/covariance don’t match the candidate contract (wrong shape, non-finite, asymmetric covariance) or if every candidate’s acquisition merit is non-finite (validated at the shared surrogate boundary, see _propose_one()).

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.

Raises ValueError if the surrogate’s predicted mean/covariance or fantasized single-point prediction don’t match the candidate contract, or if every candidate’s acquisition merit is non-finite at some step (see _validate_prediction() / _select_index()).

Parameters:
Return type:

ndarray

class VerifiableOracle(name, tier, score_fn, fidelity=None, timeout=None, max_outstanding_timeouts=4)[source]

Bases: object

A callable candidate -> OracleResult that declares its verifiability tier and fidelity.

score_fn does the actual verification (wrap a simulator, an executable check, a held-out ground-truth lookup, or a real measurement pipeline; mixle.task.toolcall’s ToolCaller is the same “external check as a callable” shape for tool calls). Construction raises if tier is not one of VERIFIABILITY_TIERS – “self-graded by a model” is not a valid tier and is rejected here, not silently accepted and discovered later.

score_fn may optionally accept a cancel_event: threading.Event keyword argument (detected once, at construction, via introspection – see _accepts_keyword). A side-effecting score_fn that wraps something cancellable (a subprocess, an HTTP request) should periodically check cancel_event.is_set() and tear its underlying call down promptly once set. This is a best-effort, cooperative cancellation boundary, not a forced kill – Python cannot forcibly terminate a thread. Calls that ignore cancellation remain quarantined, count against max_outstanding_timeouts, and eventually open a circuit that prevents additional worker creation. Their outcomes are captured in late_results when they finish. See _call_with_timeout.

Parameters:

  • name (str)

  • tier (str)

  • score_fn (Callable[[Any], OracleResult])

  • fidelity (str | None)

  • timeout (float | None)

  • max_outstanding_timeouts (int)

property late_results: tuple[LateOracleResult, ...]

Immutable snapshot of late outcomes accumulated so far.

property outstanding_timeouts: int

Number of timed-out workers that have not cooperatively stopped or completed.

class OracleResult(score, receipt=<factory>, cost=1.0, abstained=False, passed=True, valid=True, call_id=None, call_started=True)[source]

Bases: object

One candidate’s verification outcome: its score, a receipt of how it was scored, and its cost.

Construction validates its own shape (MXR-080-0189): a real, finite score unless abstained=True, a finite nonnegative cost, and a dict receipt. This is the boundary a score_fn’s raw, untrusted return value crosses (see VerifiableOracle._validate_returned), so a wrong type, a NaN/inf that leaked in from some computation, or a negative “cost” is caught immediately here rather than silently reaching a GP fit or a cost report.

abstained marks a result that is NOT a genuine observation of the oracle’s true objective at the scored candidate – for example, VerifiableOracle’s timeout abstention. An abstained result’s score is a placeholder (-inf by convention, never a claim about the true objective) kept only so the receipted history stays complete; consumers such as optimize_under_oracle() gate on abstained and must never feed one to a fit as if it were ground truth.

Parameters:
class DesignCandidate(x, result, phase=None)[source]

Bases: object

One proposed-and-verified candidate: the point tried and what the oracle said about it.

Parameters:

  • x (ndarray)

  • result (OracleResult)

  • phase (str | None)

class DesignRun(oracle_name, oracle_tier, oracle_fidelity, _oracle=None)[source]

Bases: object

The full receipted history of a design loop: every candidate tried, and the oracle’s identity.

Parameters:

  • oracle_name (str)

  • oracle_tier (str)

  • oracle_fidelity (str | None)

  • _oracle (VerifiableOracle | None)

property history: tuple[DesignCandidate, ...]

Immutable chronological snapshot of the append-only candidate ledger.

append(candidate)[source]

Append one immutable candidate record; existing ledger entries can never be replaced.

Parameters:

candidate (DesignCandidate)

Return type:

None

property oracle_calls: int

Return the number of candidates attempted against the oracle (abstentions included).

property candidate_attempts: int

All proposed candidates, including circuit-breaker abstentions that started no external call.

property genuine_history: tuple[DesignCandidate, ...]

Return history filtered to genuine (non-abstained) observations.

An abstention (see OracleResult.abstained – e.g. a timeout) is not a real observation of the oracle’s objective, so it is excluded here; use history directly for the full, unfiltered receipted record, abstentions included.

property best: DesignCandidate

Return the highest-scoring GENUINE candidate in the run history.

Abstained candidates (see OracleResult.abstained) never win: an abstention is not a real observation of the oracle’s objective, and surfacing one as “the best candidate found” would report a fabricated result as if it were verified. Raises if every candidate abstained (no genuine observation exists to report), distinct from the empty-history case.

scores()[source]

Return the run’s oracle scores in chronological order, abstentions included (not hidden).

Return type:

ndarray

report()[source]

Named receipt of the run: which oracle, at what tier/fidelity, the best candidate found.

provisional_total_cost is the cost known when calls returned. Late outcomes linked by immutable call IDs are reconciled into settled_late_cost and settled_total_cost. total_cost is only populated once no timed-out call remains outstanding or unresolved.

Return type:

dict[str, Any]

optimize_under_oracle(oracle, bounds, *, n_init=5, n_iter=15, seed=None, **bo_kwargs)[source]

Run a propose-verify-refit design loop under a fixed oracle budget.

The loop proposes candidates, verifies each one with oracle, keeps the receipted history, refits the proposal model on every GENUINE observation, and repeats under an n_init + n_iter budget of oracle calls.

oracle=None raises immediately with the explicit refusal (“no verifiable objective; cannot optimize”) – the hard precondition checked before any candidate is proposed. Continuous/low-dimensional bounds only (see module docstring); the proposal model is BayesianOptimizer, maximizing the oracle’s score.

An abstained result (OracleResult.abstained – e.g. a timeout; see VerifiableOracle) is never told to the proposal model as if it were a real observation (MXR-080-0189): a timeout means the oracle did not answer, not that the true score at that candidate is -inf, and training a fit on a fabricated -inf “observation” can corrupt it. The abstained candidate is still appended to run.history – the receipted record keeps every attempt, genuine or not – it is simply excluded from what the optimizer learns from. BayesianOptimizer.ask is explicitly designed to tolerate this (its initial-design dispensing is gated on points asked, not points told), so skipping tell for an abstention does not desynchronize or corrupt the proposal model’s state.

If the entire initial phase abstains, however, the model has no observation at all to fit an acquisition step on, and ask raises rather than guessing. The adaptive phase is therefore not entered in that case: the run stops and is returned with its receipts intact, so an all-abstention run still yields its audit report (status="no_verified_result") instead of dying inside the proposal model (MXR-080-1488).

Parameters:
Return type:

DesignRun

amplify_and_capture(oracle, bounds, *, n_init=5, n_iter=10, candidate_pool_size=200, degree=2, significance=0.05, seed=None)[source]

Round 1: search the oracle for a budget of n_init + n_iter calls (the amplified teacher).

The search must beat a budget-matched random baseline – also n_init + n_iter oracle-verified draws, not one – on a one-sided permutation test of “round 1’s best score exceeds the baseline’s best score by more than a fair reshuffling of the pooled 2 * budget scores would typically produce”, at the prespecified significance level. A single best-of-N-vs-one-draw comparison is a multiple-comparisons trap (see the module docstring): matching the baseline’s budget to the search’s, and testing the two best-of-budget scores against the null distribution obtained by repeatedly reshuffling which of the 2 * budget oracle-verified scores “belong” to round 1 versus the baseline (rather than a bare, un-tested “is my number bigger” check), is what makes beats_baseline a calibrated claim instead of a counter that trends toward “always passes” as the budget grows. If the search does not clear the gate, this returns the explicit stopped_early=True result with nothing distilled.

Otherwise: fit StudentTeacher from round 1’s history; round 2 uses the student to rank a large, FRESH candidate pool (never round 1’s own points – see the round 2 comment below) and spends the same oracle-call budget verifying only the top-ranked fresh candidates – student-guided, not blind, but every accepted score is still oracle-verified. Runs mixle.task.collapse.collapse_monitor() over the two rounds.

Parameters:
Return type:

AmplifyReport

class AmplifyReport(round1, round2, baseline, baseline_p_value, significance, beats_baseline, round2_beats_round1, incumbent_best_score, collapse, student, stopped_early, reason, seeds, round1_effective_n, baseline_effective_n, round2_effective_n)[source]

Bases: object

Receipt for a two-round amplify-and-capture run.

round1 is the initial oracle-driven search. baseline is a random draw matched to round 1’s exact oracle-call budget (not a single draw – see amplify_and_capture), so every oracle call either side spent is accounted for in its own receipted AmplificationRound. baseline_p_value is the one-sided permutation-test p-value for “round 1’s best score exceeds the baseline’s best score by more than chance”, tested against the prespecified significance level to produce beats_baseline. round2 is present only when round 1 clears that gate, and a captured student is used to propose the matched-budget follow-up batch from a FRESH candidate pool (round 1’s own points are never re-offered to round 2, so round2_beats_round1 measures generalization to unseen candidates, not re-verification of already-known ones). incumbent_best_score is bookkeeping only – the best score retained across whichever rounds ran – and, unlike round 1’s era of this report, plays no part in computing round2_beats_round1 or collapse. collapse records the trajectory-level collapse verdict when both rounds run.

Parameters:

  • round1 (AmplificationRound)

  • round2 (AmplificationRound | None)

  • baseline (AmplificationRound)

  • baseline_p_value (float | None)

  • significance (float)

  • beats_baseline (bool)

  • round2_beats_round1 (bool)

  • incumbent_best_score (float | None)

  • collapse (CollapseVerdict | None)

  • student (StudentTeacher | None)

  • stopped_early (bool)

  • reason (str | None)

  • seeds (AmplificationSeeds)

  • round1_effective_n (int)

  • baseline_effective_n (int)

  • round2_effective_n (int | None)

class AmplificationRound(run, best_score, xs)[source]

Bases: object

One oracle-verified amplification round and its best observed score.

Parameters:
class StudentTeacher(coef, degree)[source]

Bases: object

Captured regression surrogate used to rank follow-up design candidates.

The object is callable as student(x) -> predicted_score so it can be used like other teacher/task-model callables, but its role here is only proposal ranking. Accepted scores still come from the verifiable oracle.

Parameters:
fit_student(run, *, degree=2)[source]

Fit StudentTeacher from run’s GENUINE (non-abstained) oracle-verified history.

Fits against run.genuine_history, not run.history: an abstained candidate’s placeholder -inf score (OracleResult.abstained – e.g. a timeout) is not a real observation of the oracle’s objective, and least-squares regression has no tolerance for a non-finite target at all – unlike the permutation test in amplify_and_capture, where a lone -inf only nudges a p-value, here a SINGLE abstained candidate anywhere in run.history drives every coefficient in coef to nan (confirmed: a clean fit against 8 well-posed genuine points, then adding one abstained candidate to the same history turns the entire coefficient vector to nan), silently turning the captured student that guides round 2’s candidate ranking into a useless, uniformly-nan predictor with no visible error at the call site.

Parameters:

  • run (DesignRun)

  • degree (int)

Return type:

StudentTeacher

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.

All count controls are exact positive non-Boolean integers. Candidate and model matrices must be finite nonempty two-dimensional arrays with exactly one model row per candidate; ref is a finite nonempty model matrix with the same feature width. Raises ValueError if those contracts fail or if n is below the number of model parameters (the information matrix would be singular), and InfeasibleDesignError (a ValueError subclass) if n is otherwise sufficient by count but the candidate pool’s model matrix is rank-deficient – e.g. a constant or otherwise degenerate pool, where n >= p rows can never make up for the pool itself spanning fewer than p independent directions – or if, defensively, every exchange restart still fails to find a finite-merit design.

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.

Raises ValueError if info is not a finite, symmetric, PSD matrix. log det is computed via a Cholesky factorization (PSD-aware: it fails cleanly on a singular M rather than a raw determinant, which would happily return a number for any square matrix regardless of definiteness).

Parameters:
Return type:

float

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

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

Raises ValueError if info is not a finite, symmetric, PSD matrix.

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.

Raises ValueError if info is not a finite, symmetric, PSD matrix, or if ref is given and its column count does not match info’s dimension.

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).

Raises ValueError if info is not a finite, symmetric, PSD matrix, or if ref is given and its column count does not match info’s dimension.

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.

Raises ValueError if info is not a finite, symmetric, PSD matrix.

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().

Raises ValueError (eagerly, at construction) if c is not 1-D, and (per call) if info is not a finite, symmetric, PSD matrix or if c’s length does not match info’s dimension.

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, n_evaluations, failed_evaluations, stopped_reason, 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:
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: row i holds the K constraint moments at point i. A scalar or 1-D input of length n is treated as n points under a single constraint and reshaped to (n, 1) – it is never reinterpreted as one point under n constraints. 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).

Raises:

ValueError – mean and std do not have identical shape once normalized, either is not finite, std is negative anywhere, or there are zero constraint columns.

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(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).

Raises ValueError if y is empty, wrong-width, or contains a non-finite (NaN/Inf) value.

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.

ask() tracks every point it returns as pending (asked, not yet told) until a matching tell() call resolves it. A second, overlapping ask() – issued before an earlier one’s points are told, as in a parallel/async evaluation campaign – folds those still-pending points into the surrogate fit as fantasized (kriging-believer) observations, the same mechanism propose_batch() already uses to keep a single batch’s own picks apart, so acquisition is steered away from proposing them again (MXR-080-0188).

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 pending: ndarray

Return currently pending (asked, not yet told) points as a (P, d) array.

property n_pending: int

Return the number of currently pending (asked, not yet told) points.

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). Every returned point is tracked as pending until tell() resolves it; a second, overlapping ask() before any tell() folds still-pending points into the GP fit as fantasized observations so it is steered away from proposing them again (see the class docstring, MXR-080-0188).

This call is atomic: if it raises (e.g. a GP proposal attempted with zero observations), no internal state changes, so a failed call is always safe to retry.

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. Every point must be finite and within bounds, and every outcome must be finite – a NaN/Inf in either would otherwise poison every downstream GP fit and acquisition score silently. Every point must also match a currently pending proposal, i.e. one this same optimizer’s own ask() actually returned and that has not already been told:

  • a point that was never returned by ask() (an unsolicited observation) is rejected – this is almost always a caller bug (wrong array, a stale point from a different optimizer), so it is rejected by default rather than silently accepted;

  • a point that was already told (whether repeated within this same call or told in a previous one) is rejected as a duplicate – tell() is write-once per point, there is no “update an existing observation” path.

(MXR-080-0188.) Raises ValueError for any of the above. Validation covers the whole batch before anything is applied, so a batch that is only partly invalid never half-applies.

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.

Raises:

ValueError – if mean is empty, cov’s shape does not match mean, either contains a non-finite value, cov is not a valid (finite, PSD-within-tolerance) covariance, or samples is not a positive integer. See _safe_cholesky() for the covariance validation/healing policy – it never silently downgrades a true joint request to an independent one.

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.

Raises:

ValueError – if q, n_candidates, or mc_samples is not a positive integer, if there are zero observations, or if every candidate scores a non-finite q-EI merit at some batch step (rather than silently leaving that step’s pick undefined).

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.

Raises:

ValueError – if q or n_candidates is not a positive integer, if q > n_candidates, if there are zero observations, or if every candidate scores a non-finite acquisition merit (rather than letting argmax – which treats NaN as the maximum – silently select one).

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. A candidate with sd == 0 is deterministic and scores exactly 0: its observation is a constant, which carries no mutual information with y* (MXR-080-1901).

Raises:

ValueError – if mean/std are empty, mismatched in shape, non-finite, or std is negative; if max_samples is empty or non-finite; or if the resulting per-candidate information is non-finite (e.g. an extreme-but-individually-finite mean/std/ max_samples combination can overflow gamma = (y* - mu)/sd to +/-inf, and inf * 0 in the gamma * pdf term then silently yields NaN).

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).

Raises:

ValueError – if mean/std are empty, mismatched in shape, non-finite, or std is negative; or if n_samples is not a positive integer.

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.

n_init (default 2*d), max_evals, batch_size, and n_candidates (default min(2000, 100*d)) must each be exact positive integers – none is silently truncated or coerced. max_evals >= n_init and 1 <= batch_size <= n_candidates are both required up front: the first because the initial design alone spends n_init real evaluations before any GP step can begin, the second because a batch can never select more distinct candidates than the trust region generates for it, no matter how many times it retries (MXR-080-0196).

A well-defined numerical failure while fitting the local GP or sampling its joint posterior (a singular or indefinite covariance – numpy/torch linalg.LinAlgError) shrinks the trust region and retries next iteration, the same policy as an ordinary non-improving batch. Anything else (a configuration or programming bug – a bad fit_kwargs optimizer name, a missing torch install, …) propagates instead of being silently caught and relabeled “model failure” (MXR-080-0197).

Returns {'x', 'y', 'X', 'Y', 'n_restarts', 'n_bayes_batches', 'n_fit_failures', 'last_fit_error', 'degraded_to_random_search', 'stopped_reason'}: the best point/value and the full evaluation history; how many trust-region collapses restarted from a fresh Latin-hypercube design (n_restarts); how many batches were genuine GP-guided Thompson proposals (n_bayes_batches) versus caught numerical failures (n_fit_failures, with the most recent one’s diagnostic string in last_fit_error, else None); degraded_to_random_search is True iff every batch attempt failed numerically, meaning the entire result came from Latin-hypercube restarts and the GP-guided search never once ran; and stopped_reason is always "budget_exhausted" (this loop’s only normal termination) for any call that returns rather than raises.

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)[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).

__post_init__ validates the complete region invariant (MXR-080-0196): dim is an exact positive integer; length, length_min, and length_max are finite with 0 < length_min <= length <= length_max (so the region can never start collapsed, inverted, or unboundedly sized); success_tol is an exact positive integer and failure_tol an exact non-negative one (0 is the “default to dim” sentinel resolved right below). Anything outside that invariant raises immediately instead of accepting geometry with no sensible search meaning – e.g. a length_max < length_min would let a “successful” expansion (min(length * 2, length_max)) silently shrink the region, and a non-finite/non-positive length would make every candidate box _tr_candidates() draws undefined or empty.

Parameters:
property collapsed: bool

Return whether the trust-region length has shrunk below its usable minimum.

update(improved)[source]

Return the next validated state after recording one batch outcome.

States are immutable: callers must retain the returned instance. This prevents a region that passed construction validation from later being assigned invalid geometry.

Parameters:

improved (bool)

Return type:

TrustRegion

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).

n_candidates (and n_reference, for method='alc') must be exact positive integers (MXR-080-0158): a fractional count used to be silently truncated by a bare int() cast, and an n_reference=0 reduced ALC to an empty reference set (all-zero merits, an arbitrary argmax).

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.

n_init (default 2 * len(bounds)) and max_evals must both be exact positive integers with 1 <= n_init <= max_evals (MXR-080-0158): an initial design bigger than the whole budget used to be evaluated in full and returned over budget instead of rejected, a fractional count was silently truncated, and an explicit n_init=0 was silently replaced by the default rather than treated as the caller error it almost certainly is.

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_p + 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.

By Sylvester’s determinant identity, det(I_p + Sigma0 F^T F / noise^2) == det(I_n + F Sigma0 F^T / noise^2) exactly – the two sides differ only in which of the (p, p) or (n, n) matrix gets formed and factorized. This function forms whichever is smaller. That matters beyond speed: scoring a handful of candidate observations (n small, e.g. a few new monitoring sites) against a model with many parameters (p large, e.g. hundreds of grid cells) is exactly this function’s typical caller (mixle_pde.monitoring_design, mixle_pde.geophysics), and forming the full dense (p, p) matrix in that regime is reliably ill-conditioned for slogdet – verified to throw spurious divide-by-zero/overflow RuntimeWarning``s on a real ``p=720 case that the (n, n) formulation computes cleanly (identical value, to float64 precision, with zero warnings).

noise must be finite and strictly positive, and prior_cov (when given) must be a finite, symmetric, positive-semidefinite (p, p) matrix – a real covariance, not merely a matrix that happens to make a determinant come out positive (MXR-080-0159: invalid inputs used to reach slogdet unchecked and could return -inf or a plausible-looking scalar for a model that is not actually a probability distribution). Given valid inputs, the information matrix is guaranteed symmetric positive-definite by construction (an identity plus a PSD term), so its log-determinant is computed from a Cholesky factorization – 2 * sum(log(diag(L))) – rather than slogdet: more numerically stable, and it doubles as a final validity check. F Sigma0 F^T (the n <= p fast path) is symmetric for any symmetric Sigma0, since it is sandwiched by F/F^T, but Sigma0 F^T F (the other order, needed when p < n) is generally not, even though both factors are – the product of two symmetric matrices need not commute – so that branch instead sandwiches F^T F between a symmetric square root of Sigma0, which keeps the determinant exactly (Sylvester’s identity again) while keeping the matrix handed to Cholesky genuinely symmetric.

Parameters:
Return type:

float

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

The nested-Monte-Carlo EIG point estimate.

A thin wrapper around expected_information_gain_nmc_estimate() for callers that only want the point value; see its docstring for the full draw-contract validation (MXR-080-0160) and for how to also get the Monte Carlo standard error / confidence interval.

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. target, if given, must itself be one of fidelities: this function never evaluates a fidelity outside that set, so an out-of-set target could otherwise only ever be “satisfied” by silently substituting a lower-fidelity response for it. costs is the per-fidelity evaluation cost (default: the fidelity value itself) and, if given, must supply exactly one entry per fidelity. 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. Every evaluation – initial per-fidelity seeding and the sequential loop alike – reserves its cost against max_cost before calling objective, so cumulative cost never overshoots max_cost; a max_cost too small to afford even one evaluation is rejected up front.

Returns {'x', 'y', 'X', 'Y', 'cost', 'target_evaluated', 'stopped_reason', 'error', 'failed_evaluations'}. X/Y are the usable finite augmented evaluation history (fidelity in the last column of X) and cost is the total actually spent, including failed calls. failed_evaluations separately preserves each unusable attempt’s point, fidelity, cost, status, and error. x/y are the best target-fidelity point and response, but only when target_evaluated is true: if the budget ran out before any target-fidelity evaluation was affordable, x/y are None rather than silently standing in a lower-fidelity result for the target-fidelity answer the caller asked for. stopped_reason is "budget_exhausted" (the loop ran until no further evaluation fit in max_cost – the normal, successful termination for this optimizer, which has no other stopping criterion) or "surrogate_fit_failed" (a well-defined numerical failure – e.g. numpy/torch linalg.LinAlgError from a singular covariance – aborted the loop early; error then holds a diagnostic string, else None). Any other exception out of surrogate fitting (a bad fit_kwargs optimizer name, a missing torch install, …) propagates rather than being reported as a quiet early stop.

Parameters:
Return type:

dict[str, Any]

sobol_indices(func, bounds, n=4096, *, seed=0, names=None, n_bootstrap=200, confidence=0.95)[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) (bootstrap resampling below is free of additional func calls, so this count is exact).

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

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

  • n_bootstrap (int) – number of bootstrap resamples used for the standard error / confidence interval of each index (see “Returns” below).

  • confidence (float) – confidence level for the bootstrap interval, in (0, 1) – e.g. 0.95 for a two-sided 95% interval.

Returns:

  • 'S1' / 'ST': the raw, unclipped point estimate. A finite-sample Sobol estimator can legitimately land slightly negative, or (for S1) slightly above 1, purely from Monte Carlo noise – most often when the true index is near 0 or 1. That is not a bug in the estimate; it is the honest signal that sampling noise is comparable to the index’s true magnitude, and this raw value is always what is returned here (never silently clamped).

  • 'S1_clipped' / 'ST_clipped': the same estimate clipped into its theoretically valid range (S1 into [0, 1], ST below at 0) – a separate, explicitly-named convenience for e.g. plotting or reporting a single non-negative “share of variance”; never the only value returned, and never used in place of the raw estimate above.

  • 'S1_standard_error' / 'ST_standard_error': a bootstrap standard error, from n_bootstrap resamples (with replacement) of the already-evaluated A/B/AB_i rows, each re-run through the same S1/ST formula – the standard uncertainty-quantification approach for Sobol’ estimators (Archer, Saltelli & Sobol’ 1997), since no exact closed-form sampling distribution exists for them in general.

  • 'S1_ci_low' / 'S1_ci_high' and 'ST_ci_low' / 'ST_ci_high': the corresponding percentile bootstrap confidence interval at the confidence level.

Plus 'names', 'var' (the total output variance), and 'sampling_method' (normally 'scrambled_sobol'; 'pseudorandom_fallback' only when SciPy’s QMC capability is absent). S1[i] is the fraction of output variance from input i alone; ST[i] includes all interactions involving i (so ST[i] - S1[i] measures i’s interaction strength, and ST[i] ~ 0 – relative to its standard error, not just numerically close – means input i can be fixed).

Return type:

A dict with, for both first- and total-order

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

Morris elementary-effects screening – a low-cost first factor ranking.

Walks trajectories randomized-direction one-factor-at-a-time paths on a levels-point grid per dimension (Morris 1991; Campolongo, Cariboni & Saltelli 2007); 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.

Each trajectory walks the grid {0, 1/(p-1), 2/(p-1), ..., 1} (p = levels) and takes d one-at-a-time steps, one per input in a random order, each changing exactly one coordinate by the fixed step delta = p / (2*(p-1)) = (p/2) / (p-1) – i.e. exactly p/2 grid points, chosen (the standard choice) so the step always lands exactly on another grid point. Both the visit order and each dimension’s step direction (up by delta or down by delta) are independently randomized per trajectory. The starting coordinate for a dimension stepping up is drawn only from the grid’s lower half (indices 0..p/2-1, which by construction have room for a +delta move without leaving [0, 1]); for a dimension stepping down, only from the upper half (indices p/2..p-1, symmetric room for -delta). Every step therefore lands on the grid and inside bounds by construction – never clipped or silently shrunk to a smaller step – and, because the two halves are equal-sized and a dimension’s pre-step value is confined to its assigned half while its post-step value always lands in the other, every grid level is visited with equal (1/p) probability, marginally, by any one coordinate at any position along a trajectory: the design is unbiased across the input box, not concentrated away from either boundary.

levels must be even, so the grid-index step levels // 2 is exact (an odd level count has no step that both equals the standard delta formula and lands on a grid point).

Parameters:
Return type:

dict[str, Any]

fast_indices(func, bounds, n=600, *, harmonics=6, seed=0, names=None, n_bootstrap=200, confidence=0.95)[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.

S1 is the raw Tarantola-corrected finite-sample estimate. S1_clipped is a separate presentation convenience. Circular moving-block bootstrap fields report the raw estimator’s standard error and percentile interval without additional model evaluations.

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 low-cost, 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)}.

Raises:

ValueError – if mean is non-finite, cov’s shape doesn’t match mean, cov is not a valid (finite, PSD-within-tolerance) covariance, n is not a positive integer, any quantiles entry is outside [0, 1], or method is not registered.

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.

(MXR-080-1901) Registering a name that is already taken raises instead of replacing it. The old _PROPAGATORS[name] = fn silently overwrote, so re-registering "montecarlo" (by import order, a name collision between two plugins, or a module imported twice) changed what EVERY existing propagate(..., method="montecarlo") call in the process returns – the same call site, the same documented method name, different outputs, and nothing anywhere recording the substitution. This is the registry’s only mutation point, so refusing the collision here is enough to make a resolved method name mean one fixed thing for the life of the process. Adding a NEW name is unaffected: that is the extension point and it stays open. To deliberately replace a built-in, delete the existing entry from _PROPAGATORS first – an explicit, greppable act rather than an invisible one.

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.

Raises:

ValueError – if mean is non-finite, cov’s shape doesn’t match mean, cov is not a valid covariance (see _safe_cholesky()), alpha/beta/kappa are not finite non-Boolean scalars, alpha <= 0, beta < 0, d + lambda <= 0, func’s return doesn’t have a leading dimension of 2d+1 (rejected explicitly, not silently reinterpreted via a total-size-only reshape), or func’s return or derived covariance is non-finite/invalid.

Parameters:
Return type:

tuple[ndarray, ndarray]

noisy_minimize(objective, bounds, n_init=5, n_iter=15, seed=None, *, maximize=False, gp=None, fit_kwargs=None, **minimize_kwargs)[source]

Bayesian optimization of a NOISY objective, reporting the posterior-mean incumbent.

Runs mixle.doe.bayesopt.minimize() (same exploration), then replaces the raw argmin-of- observations incumbent with the posterior-mean incumbent (posterior_incumbent()), so best_x is the model’s believed optimum rather than the luckiest noisy draw and best_y is its denoised estimate. The full (x, y) evaluation history is preserved on the result.

Parameters:
Return type:

BayesOptResult

posterior_incumbent(x, y, *, maximize=False, gp=None, fit_kwargs=None)[source]

Among the evaluated points x, return the one with the best GP POSTERIOR MEAN (min by default).

Fits a GP surrogate on (x, y) – the same surrogate the BO loop uses – and predicts the denoised mean at every evaluated point, so the reported optimum reflects the model’s belief, not a single lucky noisy observation. For a deterministic objective the posterior mean at the observed points is ~the observations, so this reduces to the ordinary argmin/argmax.

Raises ValueError if x contains non-finite values, or if the surrogate’s predicted mean doesn’t match the candidate contract (wrong length, non-finite, or entirely non-finite – see mixle.doe.bayesopt._validate_prediction() / _select_index), the same shared surrogate boundary the sequential/batch/knowledge-gradient proposal paths use (MXR-080-0170).

Parameters:
Return type:

IncumbentResult

class IncumbentResult(best_x, best_mean, best_index)[source]

Bases: object

The posterior-mean incumbent among evaluated points: its location, posterior mean, and index.

Parameters:
calibrate(simulator, x, y, theta0, *, discrepancy=True, discrepancy_lengthscale=None, seed=0, max_iter=1000)[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. x and y must have the same number of rows, and x, y, theta0, and simulator(x, theta0) must all be finite.

  • y (ndarray) – field inputs and observations. x and y must have the same number of rows, and x, y, theta0, and simulator(x, theta0) must all be finite.

  • theta0 (Sequence[float]) – initial calibration parameters (its length sets the parameter count). There must be more observations than free theta parameters, and the fitted simulator sensitivity matrix must have full parameter-column rank.

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

  • discrepancy_lengthscale (float | None) – fixed GP correlation length; None (default) uses 10% of the input domain. If given explicitly it must be positive – unlike None, 0 is not silently treated as “use the default”: a caller-supplied non-positive value is a contract violation, not a request for the default.

  • seed (int) – seeds the small, reproducible random perturbation used to initialize the optimizer, so the same seed value reproduces the same fit and a different seed can escape an exact tie or a degenerate flat ridge at the unperturbed starting point.

  • max_iter (int) – positive optimizer iteration budget. Reused, unmodified, for the internal no-discrepancy warm-start fit when discrepancy=True.

Raises:

ValueError – on an x/y row-count mismatch, non-finite x/y/theta0/ simulator(x, theta0), a simulator(x, theta0) shape that does not match y, a non-positive discrepancy_lengthscale or max_iter, fewer observations than theta parameters, rank-deficient simulator sensitivity, or if the optimizer fails to converge to a finite objective.

Return type:

KOCalibration

class KOCalibration(theta, ls, amp, noise, simulator, x, y, theta_standard_error=None, *, sensitivity_rank=None, sensitivity_singular_values=None, sensitivity_condition_number=None)[source]

Bases: object

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

theta is the maximizer of the (marginal) likelihood, not a full Bayesian inference of it – there is no prior and no posterior distribution over theta here. theta_standard_error is an asymptotic (Laplace / observed-Fisher-information) standard error for each theta component, from the inverse Hessian of the negative log-likelihood at the MLE; theta_ci_low/ci_high are the corresponding two-sided 95% Wald intervals (theta +- 1.96 * theta_standard_error). This is a local, asymptotic approximation (it can be optimistic for small n, strong parameter correlation, or a non-quadratic likelihood near the optimum) – not a substitute for a real posterior when that matters. All three are arrays of nan when the Hessian at the optimum is not numerically invertible (or not positive-definite), which the point estimate itself does not depend on. sensitivity_rank, sensitivity_singular_values, and sensitivity_condition_number certify local structural identifiability; a rank-deficient fit raises CalibrationIdentifiabilityError instead of returning an arbitrary point.

All fitted arrays are privately copied and frozen. Array properties return detached copies, so neither caller-owned inputs nor a returned view can retroactively change the fitted residual, interval, or predictor.

Parameters:

  • sensitivity_rank (int | None)

  • sensitivity_singular_values (np.ndarray | None)

  • sensitivity_condition_number (float | None)

property theta: ndarray

Detached fitted parameter vector.

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

exception CalibrationIdentifiabilityError(*, rank, n_parameters, singular_values, non_identifiable_directions)[source]

Bases: ValueError

Raised when simulator sensitivity cannot identify every calibration-parameter direction.

Parameters:

  • rank (int)

  • n_parameters (int)

  • singular_values (np.ndarray)

  • non_identifiable_directions (np.ndarray)

Return type:

None

Submodules