# (c) Kevin Dunn, 2010-2026. MIT License.
"""Optimal designs (D, I, A and E) over a constrained factor region, by candidate exchange.
A ``Constraint`` such as ``"3*T + 5*D <= 600"`` cuts a corner off the factor box.
Classical designs cannot follow that cut, so the design is chosen from a set of
feasible *candidate* points instead, in three steps:
1. **Candidate set.** A grid over the box, plus the points where each constraint
boundary crosses a grid line (found by bisection) and the vertices where linear
constraints meet each other or the box, minus every infeasible point. The boundary
points matter: a D-optimal design pushes its runs to the edge of the region, and
without them the nearest grid point can sit well inside it.
2. **Model matrix.** Each candidate is expanded into the columns of the model the
design is for (intercept, main effects, two-factor interactions, squares).
3. **Exchange.** A Fedorov exchange swaps design runs for candidates while the
criterion improves (see :class:`Criterion`: the determinant of the information
matrix ``X'X`` for D, a weighted trace of its inverse for A and I, its smallest
eigenvalue for E), with several random starts. For D, A and I each run in turn
takes its best swap, with ``(X'X)^-1`` kept current by rank-one updates; when the
grid was coarse or sampled, each start's design is then polished by moving single
coordinates on a 5-level lattice.
Constraint expressions are parsed into a small arithmetic tree and evaluated with
numpy. Nothing is passed to ``eval``, so an expression from an untrusted caller can
only compute a number.
"""
from __future__ import annotations
import ast
import keyword
import logging
import math
from collections.abc import Callable, Iterable, Iterator, Mapping
from dataclasses import dataclass, field, replace
from typing import TYPE_CHECKING, Any, cast
import numpy as np
import pandas as pd
from scipy.stats import qmc
from process_improve._random import check_random_state
from process_improve.experiments._expression_text import MAX_EXPRESSION_LENGTH, nfkc
from process_improve.experiments._lm import suggest_identifier
from process_improve.experiments._region_geometry import box_vertices
from process_improve.experiments._rules import Rules, compile_rules
from process_improve.experiments._uniform_sampling import (
UniformSampler,
_affine_form,
_compares_numbers_only,
_spread,
)
from process_improve.experiments.factor import FactorType
if TYPE_CHECKING:
from process_improve.experiments.factor import Constraint, Factor
from process_improve.experiments.region import DesignRegion
logger = logging.getLogger(__name__)
#: Model types the candidate exchange builds a model matrix for.
_MODEL_TYPES = ("main_effects", "interactions", "quadratic")
#: Upper limit on the candidate grid; a larger grid is refused before allocation.
MAX_CANDIDATES = 100_000
#: Points taken from a grid of requested levels too large to list (a power of two, for the Sobol sequence): the
#: exchange's cost grows with this. An automatic grid is sampled down to :func:`_candidate_cap` instead.
_SAMPLED_CANDIDATES = 2**15
#: Slack when testing ``g(x) <= 0``, so points found on a boundary are kept.
_FEASIBILITY_TOL = 1e-9
_BISECTION_STEPS = 50
_MAX_EXCHANGES = 500
#: Rank-one updates of ``M^-1`` between full re-factorisations, which bound the drift of the updates.
_REFACTOR_EVERY = 50
#: ``M^-1`` is also rebuilt when the largest candidate variance falls below this fraction of its value at the last
#: rebuild, as it does while a singular design is filled in: the rounding error scales with the old, larger values.
_COLLAPSE = 1e-3
#: Most passes over the design rows the row-wise exchange makes.
_MAX_PASSES = 100
_N_STARTS = 5
#: Most random starts of the D, A and I exchange, and the work (candidates x coefficients x free runs, summed over
#: the starts) that sets how many a problem gets between ``_N_STARTS`` and this (see :func:`_n_starts`).
_MAX_STARTS = 50
_START_WORK = 400_000_000
#: A greedy start step picks at random among candidates whose variance is within this fraction of the largest.
_GREEDY_SLACK = 0.05
# ---------------------------------------------------------------------------
# 1. Constraint expressions -> vectorised inequality functions g(x) <= 0
# ---------------------------------------------------------------------------
_Env = Mapping[str, np.ndarray]
_Evaluator = Callable[[_Env], Any]
_BINARY_OPS: dict[type[ast.operator], Callable[[Any, Any], Any]] = {
ast.Add: np.add,
ast.Sub: np.subtract,
ast.Mult: np.multiply,
ast.Div: np.divide,
ast.Pow: np.power,
}
_UNARY_OPS: dict[type[ast.unaryop], Callable[[Any], Any]] = {ast.USub: np.negative, ast.UAdd: np.positive}
_FUNCTIONS: dict[str, Callable[[Any], Any]] = {
"abs": np.abs,
"sqrt": np.sqrt,
"exp": np.exp,
"log": np.log,
"log10": np.log10,
}
#: For each allowed comparison, whether ``g`` is ``rhs - lhs`` (True) or ``lhs - rhs``.
_FLIP: dict[type[ast.cmpop], bool] = {ast.Lt: False, ast.LtE: False, ast.Gt: True, ast.GtE: True}
def _shown(name: str) -> str:
"""Return ``name`` quoted, with its code points when it is not ASCII, so two names that print alike read apart."""
return repr(name) if name.isascii() else f"{name!r} ({' '.join(f'U+{ord(c):04X}' for c in name)})"
[docs]
def factor_lookup(names: Iterable[str]) -> dict[str, str]:
"""Map each factor name's NFKC form, the name a parsed limit holds, to the name as given.
Raises
------
ValueError
If two names share an NFKC form, so a limit could not tell which factor it names.
"""
lookup: dict[str, str] = {}
for name in sorted(names): # sorted, so a refusal names the two in one order whatever order they came in
other = lookup.setdefault(nfkc(name), name)
if other != name:
raise ValueError(
f"The factors {_shown(other)} and {_shown(name)} read as the same name, {nfkc(name)!r}, in a "
"limit: Python's parser writes every name in its Unicode NFKC form. Rename one of them."
)
return lookup
def _compile(node: ast.expr, names: Mapping[str, str]) -> _Evaluator:
"""Turn one arithmetic node into a function of the factor values.
Only numbers, factor names, ``+ - * / **``, unary signs and the functions in
``_FUNCTIONS`` are accepted; anything else raises ``ValueError``. ``names`` is
:func:`factor_lookup`'s map, from the parsed spelling of a factor's name to its name.
"""
if isinstance(node, ast.Constant) and isinstance(node.value, int | float) and not isinstance(node.value, bool):
number = float(node.value)
return lambda _env: number
if isinstance(node, ast.Name):
name = names.get(node.id)
if name is None:
known = ", ".join(_shown(n) for n in sorted(names.values()))
raise ValueError(f"Unknown name {_shown(node.id)}; constraints may use the continuous factors [{known}].")
return lambda env: env[name]
if isinstance(node, ast.BinOp) and type(node.op) in _BINARY_OPS:
fn, lhs, rhs = _BINARY_OPS[type(node.op)], _compile(node.left, names), _compile(node.right, names)
return lambda env: fn(lhs(env), rhs(env))
if isinstance(node, ast.UnaryOp) and type(node.op) in _UNARY_OPS:
unary, operand = _UNARY_OPS[type(node.op)], _compile(node.operand, names)
return lambda env: unary(operand(env))
if (
isinstance(node, ast.Call)
and isinstance(node.func, ast.Name)
and node.func.id in _FUNCTIONS
and len(node.args) == 1
and not node.keywords
):
func, argument = _FUNCTIONS[node.func.id], _compile(node.args[0], names)
return lambda env: func(argument(env))
raise ValueError(f"Unsupported syntax in constraint: {ast.unparse(node)!r}.")
def _difference(lhs: _Evaluator, rhs: _Evaluator) -> _Evaluator:
"""Return ``g = lhs - rhs``, which is feasible where ``g <= 0``."""
return lambda env: np.asarray(lhs(env), dtype=float) - np.asarray(rhs(env), dtype=float)
[docs]
def parse_constraint(expression: str, names: set[str]) -> list[_Evaluator]:
"""Parse a constraint into inequality functions, each feasible where ``g(x) <= 0``.
Parameters
----------
expression : str
An inequality in actual units, e.g. ``"3*T + 5*D <= 600"``. A chained
comparison such as ``"400 <= 3*T + 5*D <= 600"`` gives two inequalities.
Strict and non-strict comparisons are treated alike, since the boundary
of a continuous region has no volume.
names : set[str]
Names that may appear in the expression (the continuous factors). A name is
matched through its Unicode NFKC form, as Python's parser writes it, so a factor
named with a decomposed accent or a full-width letter can appear in a limit.
Returns
-------
list[Callable]
One function per comparison. Each takes a mapping from factor name to an
array of actual values and returns an array ``g``.
Raises
------
ValueError
For an expression that is too long, is not a comparison, uses ``==`` or
``!=``, or uses a name or syntax outside the allowed set; or for two names
that share an NFKC form.
"""
if len(expression) > MAX_EXPRESSION_LENGTH:
raise ValueError(f"Constraint expression is longer than {MAX_EXPRESSION_LENGTH} characters.")
try:
tree = ast.parse(expression, mode="eval").body
except SyntaxError as err:
raise ValueError(f"Constraint {expression!r} is not a valid expression: {err.msg}.") from err
compare = tree if isinstance(tree, ast.Compare) else None
if compare is None:
raise ValueError(f"Constraint {expression!r} must be an inequality, e.g. 'A + B <= 10'.")
lookup = factor_lookup(names)
terms = [_compile(compare.left, lookup)] + [_compile(c, lookup) for c in compare.comparators]
sides = [compare.left, *compare.comparators]
inequalities = []
for op, lhs, rhs, left, right in zip(compare.ops, terms[:-1], terms[1:], sides[:-1], sides[1:], strict=True):
flip = _FLIP.get(type(op))
if flip is None:
raise ValueError(
f"Constraint {expression!r}: only <, <=, >, >= are supported. An equality leaves "
"a region with no volume to place runs in; for proportions summing to 1, use "
"mixture factors."
)
constant = _cancelled(right, left, names) if flip else _cancelled(left, right, names)
inequalities.append(constant or (_difference(rhs, lhs) if flip else _difference(lhs, rhs)))
return inequalities
def _cancelled(first: ast.expr, second: ast.expr, names: set[str]) -> _Evaluator | None:
"""Return ``g = first - second`` as its exact value when its factors cancel, as in ``0.1*A + 0.2*A <= 0.3*A``.
Evaluated term by term in floating point, that difference is rounding noise in
``A``, which the tests of the region and the design engines would read as a limit
with a gradient. Read in exact arithmetic it compares numbers only. None unless both
sides are affine by their syntax and every factor's coefficient is exactly zero.
"""
from process_improve.experiments.region import _exact_difference # noqa: PLC0415 region imports this module
exact = _exact_difference(first, second, sorted(names))
if exact is None or exact[0].any():
return None
value = -exact[1]
return lambda _env: np.float64(value)
def _factor_names(tree: ast.AST) -> list[str]:
"""Names in ``tree`` that stand for factors (every ``ast.Name`` except a called function), in order of use."""
functions = {id(node.func) for node in ast.walk(tree) if isinstance(node, ast.Call)}
seen: dict[str, None] = {}
for node in ast.walk(tree):
if isinstance(node, ast.Name) and id(node) not in functions:
seen.setdefault(node.id, None)
return list(seen)
class _Renamer(ast.NodeTransformer):
"""Rename factor names in one pass, leaving called function names alone."""
def __init__(self, names: Mapping[str, str]) -> None:
self.names = names
def visit_Call(self, node: ast.Call) -> ast.Call:
node.args = [self.visit(arg) for arg in node.args]
return node
def visit_Name(self, node: ast.Name) -> ast.Name:
return ast.copy_location(ast.Name(id=self.names.get(node.id, node.id), ctx=node.ctx), node)
[docs]
def rename_factors(expression: str, names: Mapping[str, str]) -> str:
"""Rename factors in a constraint expression, all at once, and return it in canonical form.
Every name is replaced in a single pass over the parsed expression, so a swap such
as ``{"A": "B", "B": "A"}`` works. Function names (``sqrt``, ``log``, ...) are left
alone even when a factor shares them. The result is ``ast.unparse`` of the renamed
expression, so ``rename_factors(expression, {})`` gives the canonical text of an
expression, and renaming back with the inverse mapping gives that canonical text.
Parameters
----------
expression : str
A constraint, e.g. ``"3*T + 5*D <= 600"``, as :func:`parse_constraint` accepts.
names : Mapping[str, str]
Old name to new name. Names the expression does not use are ignored. An old
name is matched through its Unicode NFKC form, as :func:`parse_constraint` does.
Returns
-------
str
The renamed expression, e.g. ``"3 * Temp + 5 * D <= 600"``.
Raises
------
ValueError
If ``expression`` is not a constraint :func:`parse_constraint` accepts (before
or after renaming, which also bounds the length of the result), a new name is
not a valid identifier, or two different factors of the expression would end
up with the same name.
Examples
--------
>>> rename_factors("3*T + 5*D <= 600", {"T": "D", "D": "T"})
'3 * D + 5 * T <= 600'
>>> rename_factors("sqrt(log) <= 2", {"log": "L"})
'sqrt(L) <= 2'
"""
if len(expression) > MAX_EXPRESSION_LENGTH:
raise ValueError(f"Constraint expression is longer than {MAX_EXPRESSION_LENGTH} characters.")
try:
tree = ast.parse(expression, mode="eval")
except SyntaxError as err:
raise ValueError(f"Constraint {expression!r} is not a valid expression: {err.msg}.") from err
used = _factor_names(tree)
parse_constraint(expression, set(used))
for new in names.values():
if not (isinstance(new, str) and new.isidentifier() and not keyword.iskeyword(new)):
raise ValueError(
f"{new!r} cannot be a factor name in a constraint: it is not a valid identifier. "
f"Try {suggest_identifier(new)!r}."
)
parsed = {nfkc(old): new for old, new in names.items()} # the tree holds the NFKC spelling of each name
renamed = [parsed.get(name, name) for name in used]
if len({nfkc(name) for name in renamed}) < len(renamed):
merged = sorted(name for name in used if [nfkc(r) for r in renamed].count(nfkc(parsed.get(name, name))) > 1)
raise ValueError(f"Renaming {expression!r} with {dict(names)} would merge the factors {merged} into one.")
text = ast.unparse(_Renamer(parsed).visit(tree))
parse_constraint(text, set(renamed))
return text
# ---------------------------------------------------------------------------
# 2. Candidate set: grid + boundary crossings, filtered to the feasible region
# ---------------------------------------------------------------------------
@dataclass
class _Region:
"""The factor box, split into continuous and categorical parts.
With ``rules`` (constraints that name categorical factors, see
:mod:`~process_improve.experiments._rules`), some combinations of levels are not
allowed, and the inequalities in ``conditional[j]`` apply only on the runs where
``rules.limits[j]`` says they do. ``inequalities`` apply everywhere. ``source`` is
the same region as a :class:`~process_improve.experiments.region.DesignRegion`,
which weights each allowed combination by its volume in :meth:`draw_levels`.
"""
continuous: list[Factor]
categorical: list[Factor]
inequalities: list[_Evaluator]
rules: Rules | None = None
conditional: list[list[_Evaluator]] = field(default_factory=list)
#: Model columns kept when the rules leave some coefficients inestimable; None keeps every column.
kept_columns: np.ndarray | None = None
source: DesignRegion | None = field(default=None, repr=False, compare=False)
_lengths: np.ndarray | None = field(default=None, init=False, repr=False, compare=False)
def _limits(self) -> list[tuple[_Evaluator, int | None]]:
"""Every inequality, with the clause of ``rules.limits`` it applies under (None where it applies everywhere)."""
return [(g, None) for g in self.inequalities] + [
(g, j) for j, group in enumerate(self.conditional) for g in group
]
def labels(self, cats: np.ndarray) -> dict[str, np.ndarray]:
"""Return the level of each run, as text, for each categorical factor the rules name."""
if self.rules is None:
return {}
return {
f.name: np.asarray(self.rules.levels[f.name], dtype=str)[cats[:, j]]
for j, f in enumerate(self.categorical)
if f.name in self.rules.levels
}
def allowed(self, cats: np.ndarray) -> np.ndarray:
"""Return whether each run's combination of levels is allowed."""
return np.ones(len(cats), dtype=bool) if self.rules is None else self.rules.allowed(self.labels(cats))
def applies(self, cats: np.ndarray) -> np.ndarray:
"""Return ``(runs, clauses)``: whether each clause's inequalities apply to each run."""
if self.rules is None:
return np.zeros((len(cats), 0), dtype=bool)
return self.rules.applies(self.labels(cats))
def given(self, applies: np.ndarray) -> _Region:
"""Return the region, without rules, of the runs on which the clauses ``applies`` marks apply."""
chosen = [g for j, group in enumerate(self.conditional) if applies[j] for g in group]
return _Region(self.continuous, [], self.inequalities + chosen)
def draw_levels(self, n: int, rng: np.random.Generator) -> np.ndarray:
"""Draw ``n`` combinations of levels, each allowed combination in proportion to its volume.
The factors the rules name take one of their allowed combinations, each in
proportion to the volume of the region its limits leave (see
:meth:`~process_improve.experiments.region.DesignRegion.sample_frame`), so every
allowed setting is equally likely; every other categorical factor takes each of
its levels with the same probability.
"""
cats = np.zeros((n, len(self.categorical)), dtype=int)
for j, f in enumerate(self.categorical):
cats[:, j] = rng.integers(len(f.levels or []), size=n)
if self.rules is not None:
table, _ = self.rules.allowed_table()
pick = cast("DesignRegion", self.source)._draw_combinations(n, rng)
for j, f in enumerate(self.categorical):
if f.name in self.rules.levels:
position = {level: i for i, level in enumerate(self.rules.levels[f.name])}
cats[:, j] = [position[level] for level in table[f.name].to_numpy()[pick]]
return cats
def actual(self, coded: np.ndarray) -> dict[str, np.ndarray]:
"""Map coded continuous columns to ``{name: actual values}``."""
return {
f.name: f.center + coded[:, j] * (f.high - f.low) / 2.0 # type: ignore[operator]
for j, f in enumerate(self.continuous)
}
def offsets(self) -> np.ndarray:
"""Return ``centre / half-range`` per continuous factor: where its actual values sit next to its coded ones."""
low = np.array([f.low for f in self.continuous], dtype=float)
high = np.array([f.high for f in self.continuous], dtype=float)
return (low + high) / (high - low)
def slack(self, coded: np.ndarray, cats: np.ndarray | None = None) -> np.ndarray:
"""Return the largest ``g`` per row: the point is feasible where this is ``<= 0``.
With rules, ``cats`` (each run's categorical level indices) is needed: a run whose
combination of levels is not allowed has an infinite slack, and an inequality that
does not apply to a run is left out of its maximum.
A straight constraint has its value divided by the length of its gradient in coded
units, which makes it a distance in coded units. So a point within
``_FEASIBILITY_TOL`` of a straight boundary passes, and one further out fails,
however the constraint is scaled and whatever its units, as in
:meth:`~process_improve.experiments.region.DesignRegion.feasible`. A curved one has
its value divided by the largest size it takes over the box. In the constraint's own
units, ``1e-9`` dropped vertices that lie ``1e-12`` (coded) from a steep boundary
such as ``120*T + 33*D + 450*F <= 17000``, and let every setting pass a gentle one,
such as ``X <= 0.5e-9`` on a factor whose range is ``1e-9``.
"""
env = self.actual(coded)
rows = coded.shape[0]
if self.rules is not None and cats is None:
raise ValueError("The constraints name categorical factors, so the slack needs each run's levels.")
applies = self.applies(cats) if cats is not None else np.zeros((rows, 0), dtype=bool)
g = []
for (ineq, clause), size in zip(self._limits(), self._value_scales(), strict=True):
value = np.broadcast_to(ineq(env), (rows,)) / size
g.append(value if clause is None else np.where(applies[:, clause], value, -np.inf))
slack = np.max(g, axis=0) if g else np.full(rows, -np.inf)
return slack if cats is None else np.where(self.allowed(cats), slack, np.inf)
def _value_scales(self) -> np.ndarray:
"""Return per constraint what :meth:`slack` divides its value by (cached).
The length of a straight constraint's gradient in coded units; the largest size a
curved one takes over the box (see :func:`_spread`); and 1 for one that compares
numbers only, as :meth:`~process_improve.experiments.region.DesignRegion.feasible`
widens its tolerance.
"""
if self._lengths is None:
k = len(self.continuous)
ones = np.ones(k)
lengths = []
for g, _clause in self._limits():
coded = _coded_inequality(self, g)
form = _affine_form(coded, -ones, ones, self.offsets()) if k else None
if form is None:
lengths.append(_spread(coded, -ones, ones) if k else 1.0)
else:
lengths.append(1.0 if _compares_numbers_only(*form) else math.hypot(*form[0]))
self._lengths = np.array(lengths, dtype=float)
return self._lengths
#: The automatic grid holds at most ``max(_MIN_CANDIDATE_CAP, _CANDIDATES_PER_COEFFICIENT * p)`` candidates for a
#: model with ``p`` coefficients (and never more than ``MAX_CANDIDATES``); a larger grid is sampled down to that.
_MIN_CANDIDATE_CAP = 5000
_CANDIDATES_PER_COEFFICIENT = 200
#: Criteria on the prediction variance over the region, which place runs inside it as well as on its edges.
_VARIANCE_CRITERIA = ("i_optimal", "g_optimal")
def _candidate_cap(region: _Region, model_type: str) -> int:
"""Most candidates the automatic grid may hold, which grows with the number of model coefficients."""
n_parameters = model_matrix(
region, np.zeros((1, len(region.continuous))), np.zeros((1, len(region.categorical)), dtype=int), model_type
).shape[1]
return min(MAX_CANDIDATES, max(_MIN_CANDIDATE_CAP, _CANDIDATES_PER_COEFFICIENT * n_parameters))
def _level_choices(model_type: str, criterion: str) -> list[int]:
"""Grid resolutions to try, finest first.
A quadratic model needs 3 levels; for criteria that push runs to the edges of the
region (D, A, E, K) an even level count adds no centre and only displaces the
middle level, so 4 is skipped. Models without squared terms may use 2 levels.
"""
if model_type != "quadratic":
return [5, 4, 3, 2]
return [5, 4, 3] if criterion in _VARIANCE_CRITERIA else [5, 3]
def _grid_levels(
region: _Region, n_levels: int | None, model_type: str, criterion: str = "d_optimal"
) -> tuple[int, bool]:
"""Pick the continuous grid resolution, and whether the grid is too large to list in full.
With ``n_levels`` given, that grid is used, and sampled when it has more than
``MAX_CANDIDATES`` points. Otherwise the finest of :func:`_level_choices` whose grid
stays within :func:`_candidate_cap` is used, and the coarsest is sampled when none
does (see :func:`_lattice_points`). The cap grows with the model, so the exchange's
cost does too, and not with the number of grid levels: 5 levels on 7 factors are
78,125 points, and 3 levels lose nothing for a D-optimal quadratic design on a box.
"""
n_cat = int(np.prod([len(f.levels or []) for f in region.categorical]))
if n_levels is not None:
if n_levels < 2:
raise ValueError("n_levels must be at least 2.")
return n_levels, n_levels ** len(region.continuous) * n_cat > MAX_CANDIDATES
cap = _candidate_cap(region, model_type)
choices = _level_choices(model_type, criterion)
for n in choices:
if n ** len(region.continuous) * n_cat <= cap:
return n, False
return choices[-1], True
def _lattice_points(shape: list[int], k_cont: int, n_sampled: int | None, edges: bool = False) -> np.ndarray:
"""Index vectors of the candidate grid: all of it (``n_sampled`` None), or ``n_sampled`` well-spread points of it.
The sample is the unscrambled Sobol sequence rounded down onto the grid: evenly
spread over the grid, and deterministic without drawing on any random generator.
With ``edges`` (for a 3-level grid), a continuous coordinate takes the middle level
with probability ``1 / k_cont`` rather than 1/3, so a typical point has one factor at
its middle level and the rest at the extremes: those are the points a D-optimal
quadratic design is built from, and a uniform sample of a large grid holds few of
them (2% of the points of the 10-factor grid have no factor at its middle level).
A sampled grid on an odd number of levels also holds the centre and the face
centres (one factor at an extreme, the rest at the middle), the points a quadratic
model leans on and a random sample would rarely contain.
"""
if n_sampled is None:
return np.indices(shape).reshape(len(shape), -1).T
unit = qmc.Sobol(d=len(shape), scramble=False).random_base2(math.ceil(math.log2(n_sampled)))[:n_sampled]
idx = np.minimum((unit * shape).astype(int), np.array(shape) - 1)
if edges and k_cont:
p_middle = min(1.0 / 3.0, 1.0 / k_cont)
cont = unit[:, :k_cont]
idx[:, :k_cont] = np.where(cont < (1.0 - p_middle) / 2, 0, np.where(cont < (1.0 + p_middle) / 2, 1, 2))
levels = shape[0] if k_cont else 0
if levels % 2:
middle = np.full((2 * k_cont + 1, len(shape)), levels // 2)
middle[:, k_cont:] = 0
for axis in range(k_cont):
middle[1 + 2 * axis, axis], middle[2 + 2 * axis, axis] = 0, levels - 1
idx = np.vstack([middle, idx])
_, first = np.unique(idx, axis=0, return_index=True)
return idx[np.sort(first)]
def _boundary_points(region: _Region, grid: np.ndarray, cat_idx: np.ndarray, step: float) -> tuple:
"""Find where each constraint boundary crosses a grid edge, by vectorised bisection.
For every edge between neighbouring grid points along one axis where one end is
feasible for a constraint and the other is not, the crossing point is located on
that edge. The feasible side of the final bracket is returned, so the point
satisfies the constraint it came from; other constraints are checked afterwards.
A limit that applies on some levels only is crossed only on those.
"""
points, cats = [np.empty((0, grid.shape[1]))], [np.empty((0, cat_idx.shape[1]), dtype=int)]
for axis in range(grid.shape[1]):
lower = grid[:, axis] < 1.0 - step / 2
start, cat = grid[lower], cat_idx[lower]
applies = region.applies(cat)
direction = np.zeros(grid.shape[1])
direction[axis] = step
for ineq, clause in region._limits():
g_lo = ineq(region.actual(start)) <= 0
g_hi = ineq(region.actual(start + direction)) <= 0
cross = np.broadcast_to(g_lo != g_hi, (start.shape[0],))
if clause is not None:
cross = cross & applies[:, clause]
if not cross.any():
continue
base, lo_feasible = start[cross], np.broadcast_to(g_lo, cross.shape)[cross]
t_feas, t_infeas = np.where(lo_feasible, 0.0, 1.0), np.where(lo_feasible, 1.0, 0.0)
for _ in range(_BISECTION_STEPS):
t_mid = (t_feas + t_infeas) / 2
ok = ineq(region.actual(base + t_mid[:, None] * direction)) <= 0
t_feas, t_infeas = np.where(ok, t_mid, t_feas), np.where(ok, t_infeas, t_mid)
points.append(base + t_feas[:, None] * direction)
cats.append(cat[cross])
return np.vstack(points), np.vstack(cats)
#: Most vertices of a region the candidate set takes; a region with more gives none. Listing them is quick up to
#: here (0.2 s in 10 factors with 8 limits), while the exchange's cost grows with every candidate.
_MAX_VERTEX_CANDIDATES = 2**13
def _coded_inequality(region: _Region, g: _Evaluator) -> Callable[[np.ndarray], np.ndarray]:
"""``g`` as a function of coded points ``(n, k)``."""
return lambda x: np.broadcast_to(np.asarray(g(region.actual(np.atleast_2d(x))), dtype=float), (len(x),))
def _polytope_vertices(region: _Region) -> np.ndarray:
"""Vertices of the box cut by the region's linear constraints, in coded units.
A vertex is a point where ``k`` of the hyperplanes (box faces and linear constraint
boundaries) meet. These are where D-optimal runs go, and two constraints usually
meet away from every grid line, where neither the grid nor the boundary crossings
reach (Atkinson, Donev and Tobias 2007, ch. 12). They are listed by cutting the box
with one constraint at a time
(:func:`~process_improve.experiments._region_geometry.box_vertices`, which
:class:`~process_improve.experiments.region.DesignRegion` uses too), so the work
grows with the number of vertices rather than with the ways of choosing ``k``
boundaries. Non-linear constraints only filter the result. Returns no points when
there is no linear constraint, or when the region has more than
``_MAX_VERTEX_CANDIDATES`` vertices.
"""
k = len(region.continuous)
ones = np.ones(k)
forms = [_affine_form(_coded_inequality(region, g), -ones, ones, region.offsets()) for g in region.inequalities]
rows = [f for f in forms if f is not None]
if not rows or k == 0:
return np.empty((0, k))
a_rows, b_rows = [], []
for a, b in rows:
# Over the box |a @ x| <= sum |a|, so a row with less than |b| there (a limit on numbers only, whose
# fit leaves a gradient of rounding noise) holds everywhere or nowhere, by the sign of b.
if not a.any() or np.abs(a).sum() < abs(b):
if b < 0: # no point meets it
a_rows.append(np.zeros(k))
b_rows.append(-1.0)
continue
size = math.hypot(*a)
a_rows.append(a / size)
b_rows.append(b / size)
try:
return box_vertices(np.array(a_rows).reshape(len(a_rows), k), np.array(b_rows), _MAX_VERTEX_CANDIDATES)
except ValueError: # more vertices than the candidate set takes
return np.empty((0, k))
def _vertex_candidates(region: _Region, cat_shape: list[int]) -> tuple[np.ndarray, np.ndarray]:
"""Return the vertices of the region at every allowed combination of levels.
Without limits that apply on some levels only, the vertices are the same at every
combination. With them, the combinations are grouped by which limits apply, and each
group's vertices are found once: ``m`` such limits give at most ``2**m`` groups.
"""
if not region.conditional:
return _with_every_level(_polytope_vertices(region), cat_shape)
combos = np.indices(cat_shape).reshape(len(cat_shape), -1).T
k = len(region.continuous)
points, cats = [np.empty((0, k))], [np.empty((0, len(cat_shape)), dtype=int)]
if len(combos) > MAX_CANDIDATES:
return points[0], cats[0]
combos = combos[region.allowed(combos)]
applies = region.applies(combos)
for signature in np.unique(applies, axis=0):
members = combos[np.all(applies == signature, axis=1)]
vertices = _polytope_vertices(region.given(signature))
if len(vertices) * len(members) <= MAX_CANDIDATES:
points.append(np.repeat(vertices, len(members), axis=0))
cats.append(np.tile(members, (len(vertices), 1)))
return np.vstack(points), np.vstack(cats)
def _with_every_level(points: np.ndarray, cat_shape: list[int]) -> tuple[np.ndarray, np.ndarray]:
"""Pair every point with every combination of categorical levels, unless that is more than ``MAX_CANDIDATES``."""
combos = np.indices(cat_shape).reshape(len(cat_shape), -1).T if cat_shape else np.empty((1, 0), dtype=int)
if len(points) * len(combos) > MAX_CANDIDATES:
return np.empty((0, points.shape[1])), np.empty((0, len(cat_shape)), dtype=int)
return np.repeat(points, len(combos), axis=0), np.tile(combos, (len(points), 1))
[docs]
def build_candidates(
region: _Region, n_levels: int | None = None, model_type: str = "quadratic", criterion: str = "d_optimal"
) -> tuple[np.ndarray, np.ndarray, dict]:
"""Return the feasible candidate points in coded units.
The candidates are the grid (see :func:`_grid_levels`), the points where each
constraint boundary crosses a grid edge (:func:`_boundary_points`), and the
vertices of the region cut by its linear constraints (:func:`_polytope_vertices`),
each with every categorical level combination.
Parameters
----------
region : _Region
Factors and constraints.
n_levels : int or None
Grid levels per continuous factor; ``None`` picks them (see :func:`_grid_levels`).
model_type : str
The model the design is for: without squared terms, 2 levels may be used.
criterion : str
The optimality criterion the design is for: I- and G-optimality may use 4 levels
on a quadratic model, the others go from 5 levels to 3 (see :func:`_level_choices`).
Returns
-------
tuple[np.ndarray, np.ndarray, dict]
Coded continuous values ``(n, k_cont)``, categorical level indices
``(n, k_cat)``, and counts for the metadata.
"""
levels, sampled = _grid_levels(region, n_levels, model_type, criterion)
shape = [levels] * len(region.continuous) + [len(f.levels or []) for f in region.categorical]
k_cont = len(region.continuous)
n_sampled = None
if sampled:
n_sampled = _SAMPLED_CANDIDATES if n_levels is not None else _candidate_cap(region, model_type)
edges = n_levels is None and levels == 3 and model_type == "quadratic" and criterion not in _VARIANCE_CRITERIA
idx = _lattice_points(shape, k_cont, n_sampled, edges)
if sampled:
logger.info(
"The %d-level grid on %d factors is too large to list; sampling %d points.", levels, k_cont, len(idx)
)
grid = np.linspace(-1.0, 1.0, levels)[idx[:, :k_cont]]
cat_idx = idx[:, k_cont:]
extra, extra_cat = _boundary_points(region, grid, cat_idx, step=2.0 / (levels - 1))
corners, corner_cat = _vertex_candidates(region, shape[k_cont:])
coded = np.vstack([grid, extra, corners])
cats = np.vstack([cat_idx, extra_cat, corner_cat])
feasible = region.slack(coded, cats) <= _FEASIBILITY_TOL
coded, cats = coded[feasible], cats[feasible]
# Drop duplicates (a boundary point can coincide with a grid point).
_, unique = np.unique(np.hstack([coded.round(9), cats]), axis=0, return_index=True)
unique.sort()
counts = {
"n_levels": levels,
"n_grid_points": grid.shape[0],
"grid_sampled": sampled,
"n_boundary_points": extra.shape[0],
"n_vertex_points": corners.shape[0],
"n_candidates": unique.size,
}
return coded[unique], cats[unique], counts
# ---------------------------------------------------------------------------
# 3. Model matrix and Fedorov exchange
# ---------------------------------------------------------------------------
[docs]
def model_matrix(region: _Region, coded: np.ndarray, cats: np.ndarray, model_type: str) -> np.ndarray:
"""Expand points into model columns: intercept, main effects, interactions, squares.
A categorical factor with ``L`` levels contributes ``L - 1`` effect-coded (sum-to-zero)
columns, the coding ``evaluate_design`` uses (see
:func:`~process_improve.experiments.designs_utils.categorical_effect_columns`): a
two-level factor is -1 and +1, and with more levels the last one, in sorted order, is
-1 in every column. The levels are taken in sorted order, as patsy reads a label
column, so the criterion value reported here is the one ``evaluate_design`` gives
for the returned design.
The D-, I- and G-criteria do not depend on which full-rank coding is used, so they
select the same designs under any coding. The A-, E- and K-criteria do depend on
it, and here optimise the effect-coded model, the usual convention in DoE software.
"""
from process_improve.experiments.designs_utils import ( # noqa: PLC0415
categorical_effect_columns,
sorted_level_ranks,
)
blocks = [coded[:, [j]] for j in range(coded.shape[1])]
for j, f in enumerate(region.categorical):
levels = list(f.levels or [])
blocks.append(categorical_effect_columns(sorted_level_ranks(levels)[cats[:, j]], len(levels)))
columns = [np.ones((coded.shape[0], 1)), *blocks]
if model_type in {"interactions", "quadratic"}:
columns.extend(
np.einsum("ni,nj->nij", blocks[i], blocks[j]).reshape(coded.shape[0], -1)
for i in range(len(blocks))
for j in range(i + 1, len(blocks))
)
if model_type == "quadratic":
columns.append(coded**2)
matrix = np.hstack(columns)
return matrix if region.kept_columns is None else matrix[:, region.kept_columns]
#: Candidate rows handled at a time when a quadratic form is taken over the whole candidate set.
_CHUNK_ROWS = 8192
def _quadratic_forms(rows: np.ndarray, matrix: np.ndarray) -> np.ndarray:
"""``r' matrix r`` for every row ``r`` of ``rows``, a block of rows at a time to bound the temporaries."""
out = np.empty(len(rows))
for start in range(0, len(rows), _CHUNK_ROWS):
block = rows[start : start + _CHUNK_ROWS]
out[start : start + _CHUNK_ROWS] = np.einsum("ij,ij->i", block @ matrix, block)
return out
def _greedy_start(f_cand: np.ndarray, f_fixed: np.ndarray, n_free: int, rng: np.random.Generator) -> np.ndarray:
"""Build a non-singular starting design by adding a high-variance candidate each step.
Each step picks at random among the candidates whose prediction variance is within
``_GREEDY_SLACK`` of the largest. On a grid many candidates tie for the largest
variance; always taking the first of them gave nearly the same start every time,
so the restarts of the exchange explored almost nothing (one or two distinct local
optima from five starts). Adding a run is a rank-one update of ``M^-1`` and of the
variance vector, so a step costs one product of the candidate matrix with a vector.
Both are rebuilt from ``M`` every ``_REFACTOR_EVERY`` steps, and whenever the largest
variance has fallen by ``_COLLAPSE`` since the last rebuild: the ridge makes the first
variances about ``1e6`` times the final ones, and the rounding error they leave behind
would otherwise swamp the final values.
"""
p = f_cand.shape[1]
info = f_fixed.T @ f_fixed + 1e-6 * np.eye(p) # small ridge so the first steps are defined
rows = [int(rng.integers(f_cand.shape[0]))]
info += np.outer(f_cand[rows[0]], f_cand[rows[0]])
m_inv = np.linalg.inv(info)
variance = _quadratic_forms(f_cand, m_inv)
refactored, peak = 0, float(variance.max())
for step in range(1, n_free):
near = np.flatnonzero(variance >= (1.0 - _GREEDY_SLACK) * variance.max())
best = int(rng.choice(near))
rows.append(best)
added = f_cand[best]
info += np.outer(added, added)
if step == n_free - 1:
break
m_inv_added = m_inv @ added
scale = 1.0 / (1.0 + float(added @ m_inv_added))
variance -= scale * (f_cand @ m_inv_added) ** 2
m_inv -= scale * np.outer(m_inv_added, m_inv_added)
if step - refactored >= _REFACTOR_EVERY or variance.max() < _COLLAPSE * peak:
m_inv = np.linalg.inv(info)
variance = _quadratic_forms(f_cand, m_inv)
refactored, peak = step, float(variance.max())
return np.array(rows)
[docs]
@dataclass(frozen=True)
class Criterion:
"""The optimality criterion an exchange maximises.
``"d_optimal"`` maximises ``log det(M)``, with ``M = X'X``. ``"a_optimal"`` and
``"i_optimal"`` minimise ``trace(M^-1 W)``: with ``W = I`` that is the summed
variance of the coefficients (A), and with ``W`` the average of ``f(x) f(x)'``
over the region it is the average prediction variance over the region (I), which
is what ``evaluate_design`` reports as I-efficiency. ``"e_optimal"`` maximises the
smallest eigenvalue of ``M``, which bounds the variance of the worst-estimated
linear combination of the coefficients.
Use :meth:`d`, :meth:`a`, :meth:`i` or :meth:`e` to build one.
"""
name: str
weights: np.ndarray | None = None
region_rows: np.ndarray | None = None
segment: np.ndarray | None = None
[docs]
@classmethod
def d(cls) -> Criterion:
"""D-optimality."""
return cls("d_optimal")
[docs]
@classmethod
def a(cls, n_parameters: int) -> Criterion:
"""A-optimality for a model with ``n_parameters`` coefficients."""
return cls("a_optimal", np.eye(n_parameters))
[docs]
@classmethod
def i(cls, region_rows: np.ndarray) -> Criterion:
"""I-optimality, with the moment matrix estimated from model rows sampled uniformly in the region."""
return cls("i_optimal", region_rows.T @ region_rows / len(region_rows))
[docs]
@classmethod
def e(cls) -> Criterion:
"""E-optimality."""
return cls("e_optimal")
[docs]
@classmethod
def g(cls, region_rows: np.ndarray) -> Criterion:
"""G-optimality: the largest prediction variance over the model rows ``region_rows``."""
return cls("g_optimal", region_rows=region_rows)
[docs]
@classmethod
def k(cls) -> Criterion:
"""K-optimality: the condition number of ``X'X`` (Ye and Zhou 2013)."""
return cls("k_optimal")
[docs]
def blocked(self, n_block_columns: int, segment: np.ndarray) -> Criterion:
"""Return this criterion for a design run in blocks: the factor effects, adjusted for the blocks.
The model rows gain ``n_block_columns`` block columns, and ``segment`` gives the
block of each candidate row; a run is only ever swapped for a candidate of its
own block, so the block sizes stay as they started. With them fixed,
``det(M) = det(Z'Z) det(S)``, where ``S`` is the factor information adjusted for
the blocks (the Schur complement of the block part ``Z'Z`` of ``M``), so D
maximises ``det(S)``. The trace criteria give the block columns no weight, and
``trace(M^-1 W) = trace(S^-1 W)``. Only D, A and I are blocked this way.
"""
weights = None if self.weights is None else np.pad(self.weights, (0, n_block_columns))
return replace(self, weights=weights, segment=segment)
[docs]
def value(self, info: np.ndarray) -> float:
"""Score an information matrix; higher is better.
``log det(M)`` (D), ``-trace(M^-1 W)`` (A, I), ``lambda_min`` (E), minus the largest
prediction variance over the region (G), and ``-log(lambda_max / lambda_min)`` (K).
"""
if self.name in ("e_optimal", "k_optimal", "g_optimal"):
return _spectral_or_minimax_value(self, info)
if self.weights is None:
sign, logdet = np.linalg.slogdet(info)
return float(logdet) if sign > 0 else -np.inf
if np.linalg.matrix_rank(info) < info.shape[0]:
return -np.inf
return -float(np.trace(np.linalg.solve(info, self.weights)))
@property
def n_phases(self) -> int:
"""Exchange phases: E-optimality climbs ``phi_p`` first, then polishes ``lambda_min``."""
return 2 if self.name == "e_optimal" else 1
[docs]
def best_swap(
self, info: np.ndarray, f_design: np.ndarray, f_cand: np.ndarray, phase: int = 0
) -> tuple[int, int, float]:
"""Return ``(i, j, gain)`` for the best single swap: design row ``i`` out, candidate ``j`` in."""
if self.name == "e_optimal":
return _best_e_swap(info, f_design, f_cand, polish=phase > 0)
if self.name == "k_optimal":
return _best_k_swap(info, f_design, f_cand)
gains = self.swap_gains(np.linalg.pinv(info), f_design, f_cand)
i, j = np.unravel_index(np.argmax(gains), gains.shape)
return int(i), int(j), float(gains[i, j])
[docs]
def swap_gains(self, m_inv: np.ndarray, f_design: np.ndarray, f_cand: np.ndarray) -> np.ndarray:
"""Return the gain of every swap (design run ``i`` out, candidate ``j`` in), shape ``(n_design, n_cand)``.
With ``A = M^-1`` and ``d(a, b) = f(a)' A f(b)``, D-optimality multiplies the
determinant by ``1 + d(j) - d(i) - [d(i) d(j) - d(i, j)**2]`` (Fedorov 1972).
For the trace criteria, two Sherman-Morrison updates (add ``j``, then remove
``i``) with ``b(a, b) = f(a)' A W A f(b)`` give the new trace in closed form::
s = 1 + d(j)
trace_new = trace - b(j)/s + [b(i) - 2 d(i,j) b(i,j)/s + d(i,j)**2 b(j)/s**2] / [1 - d(i) + d(i,j)**2/s]
so both kinds are scored for all pairs in a few matrix products.
"""
a_design, a_cand = f_design @ m_inv, f_cand @ m_inv
d_i = np.einsum("ij,ij->i", a_design, f_design)[:, None]
d_j = np.einsum("ij,ij->i", a_cand, f_cand)[None, :]
d_ij = a_design @ f_cand.T
if self.weights is None:
return _d_gain(d_i, d_j, d_ij)
b_design, b_cand = a_design @ self.weights, a_cand @ self.weights
b_i = np.einsum("ij,ij->i", b_design, a_design)[:, None]
b_j = np.einsum("ij,ij->i", b_cand, a_cand)[None, :]
b_ij = b_design @ a_cand.T
return _trace_gain((d_i, d_j, d_ij), (b_i, b_j, b_ij))
def _d_gain(d_i: float | np.ndarray, d_j: np.ndarray, d_ij: np.ndarray) -> np.ndarray:
"""Return the relative change in ``det(M)`` from swapping run ``i`` for candidate ``j``."""
return d_j - d_i - (d_i * d_j - d_ij**2)
def _trace_gain(d: tuple, b: tuple) -> np.ndarray:
"""Return the fall in ``trace(M^-1 W)`` from swapping run ``i`` for candidate ``j`` (:meth:`Criterion.swap_gains`).
``d`` holds ``(d(i), d(j), d(i, j))`` and ``b`` holds ``(b(i), b(j), b(i, j))``, as
scalars or arrays that broadcast against each other.
"""
(d_i, d_j, d_ij), (b_i, b_j, b_ij) = d, b
s = 1.0 + d_j
denominator = 1.0 - d_i + d_ij**2 / s
numerator = b_i - 2.0 * d_ij * b_ij / s + d_ij**2 * b_j / s**2
with np.errstate(divide="ignore", invalid="ignore"):
gain = b_j / s - numerator / denominator
return np.where(denominator > 1e-10, gain, -np.inf)
#: Criteria whose swap gains have a closed form (:meth:`Criterion.swap_gains`), climbed by :func:`_row_exchange`.
_CLOSED_FORM_CRITERIA = ("d_optimal", "a_optimal", "i_optimal")
def _regularised_inverse(info: np.ndarray) -> np.ndarray:
"""``M^-1``, with a small ridge when ``M`` is singular so swaps that fill the missing directions score highest."""
eigenvalues = np.linalg.eigvalsh(info)
largest = max(float(eigenvalues[-1]), 1.0)
ridge = 0.0 if eigenvalues[0] > 1e-10 * largest else 1e-8 * largest
return np.linalg.inv(info + ridge * np.eye(len(info)))
#: Entries of the (design rows x candidates) block of swap-gain terms the row-wise exchange keeps.
_BLOCK_ENTRIES = 2**20
#: Most design rows in such a block: each swap updates the terms of the rows after it.
_BLOCK_ROWS = 32
@dataclass
class _RowTerms:
"""Swap-gain terms of a block of design rows: ``d(i)``, ``d(i, j)``, and ``b(i)``, ``b(i, j)`` for a trace."""
d_i: np.ndarray
d_ij: np.ndarray
b_i: np.ndarray | None = None
b_ij: np.ndarray | None = None
def follow(self, x: np.ndarray, step: tuple, start: int) -> None:
"""Bring rows ``start:`` (model rows ``x``) up to date after one rank-one update ``step`` of ``A`` and ``B``.
With ``A <- A - s u u'``, ``d(i, j)`` falls by ``s (u'x)(F u)`` and ``d(i)`` by
``s (u'x)^2``; ``B`` changes by two cross terms and a square term, and so do
``b(i, j)`` and ``b(i)``. Each costs a few passes over the block, where a fresh
score would multiply the block by the whole candidate matrix again.
"""
u, g, s, vbv, f_u, f_g = step
ux = x @ u
self.d_i[start:] -= s * ux**2
self.d_ij[start:] -= s * np.multiply.outer(ux, f_u)
if self.b_i is not None and self.b_ij is not None:
gx = x @ g
self.b_i[start:] += s * (s * vbv * ux**2 - 2.0 * ux * gx)
self.b_ij[start:] += s * (
np.multiply.outer(s * vbv * ux - gx, f_u) - np.multiply.outer(ux, f_g) # type: ignore[arg-type]
)
class _ExchangeState:
"""``M^-1`` of a design and the candidates' variance terms, kept current through each swap.
For D-optimality the state is ``A = M^-1`` and ``d(j) = f(j)' A f(j)`` for every
candidate; the trace criteria also keep ``B = A W A`` and ``b(j) = f(j)' B f(j)``.
A swap is two rank-one (Sherman-Morrison) updates, adding the candidate and then
removing the run. Each costs one product of the candidate matrix with a vector (two
for a trace) instead of a new inverse and a new pass over every (run, candidate)
pair; the removal reuses the products the swap was scored with. The state is rebuilt
from ``X'X`` every ``_REFACTOR_EVERY`` swaps, and when the largest variance collapses
(see ``_COLLAPSE``), to bound rounding drift.
"""
def __init__(self, f_cand: np.ndarray, f_fixed: np.ndarray, rows: np.ndarray, criterion: Criterion) -> None:
self.f_cand, self.f_fixed, self.rows = f_cand, f_fixed, rows.copy()
self.weights, self.segment = criterion.weights, criterion.segment
self.refactor()
def refactor(self) -> None:
"""Rebuild ``A``, ``d`` (and ``B``, ``b``) from the design's information matrix."""
x = np.vstack([self.f_fixed, self.f_cand[self.rows]])
self.m_inv = _regularised_inverse(x.T @ x)
self.variance = _quadratic_forms(self.f_cand, self.m_inv)
if self.weights is not None:
self.b_mat = self.m_inv @ self.weights @ self.m_inv
self.b_variance = _quadratic_forms(self.f_cand, self.b_mat)
self.n_updates, self.peak = 0, float(self.variance.max())
def terms(self, block: np.ndarray) -> _RowTerms:
"""Return the swap-gain terms of design rows ``block`` against every candidate, one matrix product each."""
x = self.f_cand[self.rows[block]]
x_a = x @ self.m_inv
terms = _RowTerms(np.einsum("ij,ij->i", x_a, x), x_a @ self.f_cand.T)
if self.weights is not None:
x_b = x @ self.b_mat
terms.b_i, terms.b_ij = np.einsum("ij,ij->i", x_b, x), x_b @ self.f_cand.T
return terms
def best_swaps(self, terms: _RowTerms, rows: slice, design_rows: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Return the best candidate for block ``rows`` and the gain of each swap (:meth:`Criterion.swap_gains`).
``design_rows`` are those rows' positions in the design. For D the gain
``d(j) (1 - d(i)) + d(i, j)^2 - d(i)`` is ranked without its last term, in one
temporary, where the general form takes several. In a blocked design a run only
competes with the candidates of its own block.
"""
d_i = terms.d_i[rows]
index = np.arange(len(d_i))
if self.weights is None:
score = np.square(terms.d_ij[rows])
score += np.multiply.outer(1.0 - d_i, self.variance)
self._keep_in_block(score, design_rows)
best = np.argmax(score, axis=1)
return best, score[index, best] - d_i
gains = _trace_gain(
(d_i[:, None], self.variance[None, :], terms.d_ij[rows]),
(terms.b_i[rows, None], self.b_variance[None, :], terms.b_ij[rows]), # type: ignore[index]
)
self._keep_in_block(gains, design_rows)
best = np.argmax(gains, axis=1)
return best, gains[index, best]
def _keep_in_block(self, scores: np.ndarray, design_rows: np.ndarray) -> None:
"""Rule out, in place, every swap of a run for a candidate of another block (when the design is blocked)."""
if self.segment is not None:
own_block = self.segment[self.rows[design_rows]]
scores[own_block[:, None] != self.segment[None, :]] = -np.inf
def _rank_one(
self, v: np.ndarray, sign: float, f_u: np.ndarray | None = None, f_g: np.ndarray | None = None
) -> tuple:
"""Add (``sign = 1``) or remove (``sign = -1``) the model row ``v``: ``A <- A - s (A v)(A v)'``.
``f_u = F A v`` and ``f_g = F B v`` are computed unless given. Returns the step
``(A v, B v, s, v' B v, f_u, f_g)``, from which the products of another row with
the updated ``A`` and ``B`` follow without touching ``F`` again (:meth:`_RowTerms.follow`).
"""
u = self.m_inv @ v
s = sign / (1.0 + sign * float(v @ u))
f_u = self.f_cand @ u if f_u is None else f_u
g, vbv = None, 0.0
if self.weights is not None: # B = A W A takes two cross terms and one square term
g = self.b_mat @ v
vbv = float(v @ g)
f_g = self.f_cand @ g if f_g is None else f_g
self.b_variance += s * (s * vbv * f_u**2 - 2.0 * f_u * f_g)
cross = np.outer(u, g)
self.b_mat += s * (s * vbv * np.outer(u, u) - cross - cross.T)
self.variance -= s * f_u**2
self.m_inv -= s * np.outer(u, u)
return u, g, s, vbv, f_u, f_g
def swap(self, i: int, j: int, terms: _RowTerms, r: int) -> list[tuple] | None:
"""Replace design row ``i`` by candidate ``j``; row ``r`` of ``terms`` holds its scored terms.
Returns the two rank-one steps made, or None when the state was rebuilt instead
(and other rows' terms must be scored afresh).
"""
removed = self.f_cand[self.rows[i]]
self.rows[i] = j
self.n_updates += 1
if self.n_updates >= _REFACTOR_EVERY:
self.refactor()
return None
added = self._rank_one(self.f_cand[j], 1.0)
u, g, s, vbv, f_u, f_g = added
# The removed row's products with the updated A and B, from the ones it was scored with.
u_a = float(u @ removed)
f_a, f_ba = terms.d_ij[r] - s * u_a * f_u, None
if terms.b_ij is not None and g is not None:
f_ba = terms.b_ij[r] - s * (f_u * float(g @ removed) + f_g * u_a) + s * s * vbv * u_a * f_u
steps = [added, self._rank_one(removed, -1.0, f_a, f_ba)]
if self.variance.max() < _COLLAPSE * self.peak:
self.refactor()
return None
return steps
def _climb_block(state: _ExchangeState, block: np.ndarray) -> bool:
"""Swap each row of ``block`` in turn for its best candidate, when that helps; return whether any was swapped.
The block's terms come from one matrix product. A swap changes ``A``, and the terms
of the rows after it follow by rank-one updates (:meth:`_RowTerms.follow`), so every
row is scored against the current design without multiplying by the candidates
again. Rows are scored a few at a time, one after a swap and twice as many after each
group with none, so frequent swaps do not rescore rows that a later swap changes.
"""
terms, r, size, swapped = state.terms(block), 0, 1, False
while r < len(block):
best, gain = state.best_swaps(terms, slice(r, r + size), block[r : r + size])
better = np.flatnonzero(gain > 1e-9) # also refuses NaN
if not len(better):
r, size = r + size, 2 * size
continue
r += int(better[0])
steps = state.swap(int(block[r]), int(best[better[0]]), terms, r)
r, size, swapped = r + 1, 1, True
if r == len(block):
break
if steps is None: # the state was rebuilt: score the rest afresh
terms, block, r = state.terms(block[r:]), block[r:], 0
continue
rest = state.f_cand[state.rows[block[r:]]]
for step in steps:
terms.follow(rest, step, r)
return swapped
def _row_exchange(f_cand: np.ndarray, f_fixed: np.ndarray, rows: np.ndarray, criterion: Criterion) -> np.ndarray:
"""Swap each design row for its best candidate in turn: the modified Fedorov exchange (Cook and Nachtsheim 1980).
Each row is scored against every candidate and its best improving swap is made at
once, where the full Fedorov step scored all ``n x N`` (run, candidate) pairs, and
rebuilt ``X'X`` and its inverse, to make one swap. Passes repeat until none improves
the criterion. Rows are scored in blocks of up to ``_BLOCK_ROWS`` rows and
``_BLOCK_ENTRIES`` terms (:func:`_climb_block`).
"""
state = _ExchangeState(f_cand, f_fixed, rows, criterion)
size = int(np.clip(_BLOCK_ENTRIES // max(len(f_cand), 1), 1, _BLOCK_ROWS))
for _ in range(_MAX_PASSES):
improved = False
for start in range(0, len(rows), size):
improved |= _climb_block(state, np.arange(start, min(len(rows), start + size)))
if not improved:
break
return state.rows
#: Exponent of Kiefer's ``phi_p``, the smooth stand-in for ``lambda_min`` that ranks E-optimal swaps.
_E_P = 8
#: Work allowed per E-optimal iteration for exact scoring, in (swaps x coefficients^2): all swaps when they fit.
_E_EXACT_WORK = 4_000_000
def _phi_p(eigenvalues: np.ndarray) -> np.ndarray:
"""Kiefer's ``phi_p = (sum lambda_k^-p)^(-1/p)`` along the last axis of ascending eigenvalues.
It sits just below ``lambda_min`` and rises whenever any small eigenvalue rises,
so it separates designs that tie on ``lambda_min``. Written through the ratios
``lambda_min / lambda_k <= 1`` so the powers cannot overflow.
"""
floor = 1e-12 * max(float(np.max(eigenvalues)), 1.0)
lam = np.maximum(eigenvalues, floor)
ratio = lam[..., :1] / lam
return lam[..., 0] * np.sum(ratio**_E_P, axis=-1) ** (-1.0 / _E_P)
def _exact_e_scores(
info: np.ndarray, f_design: np.ndarray, f_cand: np.ndarray, pairs: tuple[np.ndarray, np.ndarray]
) -> tuple[np.ndarray, np.ndarray]:
"""``lambda_min`` and ``phi_p`` after each swap in ``pairs`` (design rows, candidate rows), in chunks."""
p = info.shape[0]
chunk = max(1, 2_000_000 // p**2) # about 16 MB of matrices at a time
lam_min, phi = [], []
for start in range(0, len(pairs[0]), chunk):
a = f_design[pairs[0][start : start + chunk]]
b = f_cand[pairs[1][start : start + chunk]]
swapped = info[None] - a[:, :, None] * a[:, None, :] + b[:, :, None] * b[:, None, :]
eigenvalues = np.linalg.eigvalsh(swapped)
lam_min.append(eigenvalues[:, 0])
phi.append(_phi_p(eigenvalues))
return np.concatenate(lam_min), np.concatenate(phi)
def _best_e_swap(
info: np.ndarray, f_design: np.ndarray, f_cand: np.ndarray, *, polish: bool = False
) -> tuple[int, int, float]:
"""Best swap for E-optimality: rank every swap by the gradient of ``phi_p``, then score the best exactly.
The exchange runs in two phases. The first climbs ``phi_p`` itself, a smooth
stand-in for ``lambda_min`` that can trade a little of the smallest eigenvalue for
a lot of the next ones, and so escapes designs where no swap raises
``lambda_min`` alone. The second (``polish=True``) then takes only swaps that
raise ``lambda_min``, or hold it and raise ``phi_p``. Each phase climbs one
objective, so neither can cycle. Two things make ``lambda_min`` hard to climb one
swap at a time:
- **Ties.** On a +/-1 grid thousands of swaps share a first-order estimate, so a
short list of "most promising" swaps can miss the one that helps. The list
scored exactly is therefore large (every swap when ``_E_EXACT_WORK`` allows),
and batched eigenvalue solves keep that cheap.
- **Repeated eigenvalues.** If ``lambda_min`` has multiplicity 2 or more, no single
swap can raise it: adding ``b b'`` gives ``lambda_1(M + b b') <= lambda_2(M)``
(Weyl interlacing). A swap that keeps ``lambda_min`` and raises ``phi_p``
splits the repeated eigenvalue, so a later swap can lift it.
Swaps are ranked by the first-order change in ``phi_p``, which weights each
eigen-direction by ``(lambda_min / lambda_k)^(p+1)``: the directions at or near the
minimum count fully, the rest hardly at all.
Returns
-------
tuple[int, int, float]
Design row out, candidate row in, and the gain: in ``phi_p`` in the first phase;
when polishing, in ``lambda_min`` when it rises, otherwise in ``phi_p`` with
``lambda_min`` held. A gain of 0 means no swap helps.
"""
eigenvalues, eigenvectors = np.linalg.eigh(info)
current, current_phi = float(eigenvalues[0]), float(_phi_p(eigenvalues))
weights = (max(current, 1e-12) / np.maximum(eigenvalues, 1e-12)) ** (_E_P + 1)
added = (f_cand @ eigenvectors) ** 2 @ weights
removed = (f_design @ eigenvectors) ** 2 @ weights
estimate = added[None, :] - removed[:, None]
n_exact = max(2000, _E_EXACT_WORK // info.shape[0] ** 2)
order = np.argsort(estimate, axis=None)[::-1][:n_exact]
rows_out, rows_in = np.unravel_index(order, estimate.shape)
pairs = (rows_out, rows_in)
lam_min, phi = _exact_e_scores(info, f_design, f_cand, pairs)
phi_tol = 1e-9 * max(1.0, abs(current_phi))
if not polish:
best = int(np.argmax(phi))
gain = float(phi[best] - current_phi)
return (int(pairs[0][best]), int(pairs[1][best]), gain) if gain > phi_tol else (0, 0, 0.0)
tol = 1e-9 * max(1.0, abs(current))
if (lam_min > current + tol).any():
best = int(np.argmax(np.where(lam_min > current + tol, lam_min + 1e-12 * phi, -np.inf)))
return int(pairs[0][best]), int(pairs[1][best]), float(lam_min[best] - current)
holds = lam_min >= current - tol
if holds.any():
best = int(np.argmax(np.where(holds, phi, -np.inf)))
gain = float(phi[best] - current_phi)
if gain > phi_tol:
return int(pairs[0][best]), int(pairs[1][best]), gain
return 0, 0, 0.0
def _spectral_or_minimax_value(criterion: Criterion, info: np.ndarray) -> float:
"""``lambda_min`` (E), ``-log(lambda_max / lambda_min)`` (K), or minus the largest prediction variance (G)."""
if criterion.name == "g_optimal":
if np.linalg.matrix_rank(info) < info.shape[0]:
return -np.inf
return -float(_prediction_variance(criterion.region_rows, np.linalg.inv(info)).max()) # type: ignore[arg-type]
eigenvalues = np.linalg.eigvalsh(info)
if criterion.name == "e_optimal":
return float(eigenvalues[0])
if eigenvalues[0] <= 1e-12 * max(float(eigenvalues[-1]), 1e-300):
return -np.inf
return -float(np.log(eigenvalues[-1] / eigenvalues[0]))
def _prediction_variance(rows: np.ndarray, m_inv: np.ndarray) -> np.ndarray:
"""``f(x)' M^-1 f(x)`` for every model row ``f(x)`` in ``rows``."""
return np.einsum("ij,jk,ik->i", rows, m_inv, rows)
def _swap_eigenvalues(info: np.ndarray, f_design: np.ndarray, f_cand: np.ndarray, pairs: tuple) -> np.ndarray:
"""Ascending eigenvalues of ``M - a a' + b b'`` for each swap in ``pairs`` (design rows, candidate rows)."""
p = info.shape[0]
chunk = max(1, 2_000_000 // p**2)
out = []
for start in range(0, len(pairs[0]), chunk):
a = f_design[pairs[0][start : start + chunk]]
b = f_cand[pairs[1][start : start + chunk]]
out.append(np.linalg.eigvalsh(info[None] - a[:, :, None] * a[:, None, :] + b[:, :, None] * b[:, None, :]))
return np.concatenate(out)
#: Swaps whose K-criterion is within this fraction of the best swap's tie with it (:func:`_best_k_swap`).
_K_TIE = 1e-9
def _best_k_swap(info: np.ndarray, f_design: np.ndarray, f_cand: np.ndarray) -> tuple[int, int, float]:
"""Best swap for K-optimality (smallest condition number of ``X'X``).
Every swap is ranked by the first-order change in ``log(lambda_min) - log(lambda_max)``,
``((v' b)^2 - (v' a)^2) / lambda`` along the extreme eigenvectors ``v``, and the most
promising ones are scored exactly with batched eigenvalue solves, as for E-optimality.
Of the swaps that tie for the best score, to within ``_K_TIE``, the one with the first design
row, then the first candidate, is made. On a grid many swaps tie exactly, and an extreme
eigenvalue is usually repeated, so its eigenvector is any vector of its eigenspace, chosen by
the BLAS kernel's rounding. Taking the first maximum in the order that vector ranks the swaps
in gave a different search path, and designs of different quality, on different processors
(#1129).
"""
eigenvalues, eigenvectors = np.linalg.eigh(info)
lam_min, lam_max = max(float(eigenvalues[0]), 1e-12), float(eigenvalues[-1])
current = -np.log(lam_max / lam_min)
def pull(rows: np.ndarray) -> np.ndarray:
return (rows @ eigenvectors[:, 0]) ** 2 / lam_min - (rows @ eigenvectors[:, -1]) ** 2 / lam_max
estimate = pull(f_cand)[None, :] - pull(f_design)[:, None]
n_exact = max(2000, _E_EXACT_WORK // info.shape[0] ** 2)
order = np.argsort(estimate, axis=None)[::-1][:n_exact]
pairs = np.unravel_index(order, estimate.shape)
eig = _swap_eigenvalues(info, f_design, f_cand, pairs)
with np.errstate(divide="ignore", invalid="ignore"):
score = np.where(eig[:, 0] > 1e-12 * eig[:, -1], -np.log(eig[:, -1] / eig[:, 0]), -np.inf)
top = float(np.max(score))
tied = np.flatnonzero(score >= top - _K_TIE * max(1.0, abs(top))) # every swap when all are singular
best = int(tied[np.argmin(order[tied])]) # order holds each swap's (design row, candidate) as a flat index
gain = float(score[best] - current)
return (int(pairs[0][best]), int(pairs[1][best]), gain) if gain > 1e-12 else (0, 0, 0.0)
#: Reweighting rounds of the I-lambda search for G-optimality, and the size of its final polish.
_G_ROUNDS = 12
_G_POLISH_CANDIDATES = 400
_G_POLISH_PASSES = 20
def _g_polish(f_cand: np.ndarray, f_fixed: np.ndarray, rows: np.ndarray, region: np.ndarray) -> np.ndarray:
"""Exact exchange on the largest prediction variance over ``region``: take the best single swap while it helps.
With ``A = (M - a a')^-1`` after removing design row ``a``, adding candidate ``b`` gives
``d(r) = r' A r - (r' A b)^2 / (1 + b' A b)`` at every region row ``r`` at once.
"""
rows = rows.copy()
for _ in range(_G_POLISH_PASSES):
x = np.vstack([f_fixed, f_cand[rows]])
m_inv = np.linalg.inv(x.T @ x)
current = float(_prediction_variance(region, m_inv).max())
pool = np.arange(len(f_cand))
if len(pool) > _G_POLISH_CANDIDATES: # the runs that help most sit where the variance is largest
pool = np.argsort(_prediction_variance(f_cand, m_inv))[::-1][:_G_POLISH_CANDIDATES]
best = (current - 1e-9 * current, -1, -1)
for i, row in enumerate(rows):
a = f_cand[row]
denominator = 1.0 - a @ m_inv @ a
if denominator <= 1e-10:
continue
m_inv_a = m_inv @ a
removed = m_inv + np.outer(m_inv_a, m_inv_a) / denominator
region_removed = region @ removed
base = np.einsum("ij,ij->i", region_removed, region)
b = f_cand[pool]
cross = region_removed @ b.T
scale = 1.0 + np.einsum("ij,jk,ik->i", b, removed, b)
worst = (base[:, None] - cross**2 / scale[None, :]).max(axis=0)
j = int(np.argmin(worst))
if worst[j] < best[0]:
best = (float(worst[j]), i, int(pool[j]))
if best[1] < 0:
break
rows[best[1]] = best[2]
return rows
def _g_exchange(
f_cand: np.ndarray, n_free: int, f_fixed: np.ndarray, rng: np.random.Generator, criterion: Criterion
) -> tuple[np.ndarray, float]:
"""G-optimal design by I-lambda optimality (Hernandez and Nachtsheim 2018), then an exact polish.
The largest prediction variance over the region is hard to lower one swap at a time,
so each round runs the closed-form trace exchange for a weighted average of the
variance (I-lambda), and Lawson's rule then moves weight towards the region points
where the variance is largest (``w <- w * d``). The rounds converge on the minimax
design; the best design met is kept and polished by exact swaps on the maximum.
"""
region = criterion.region_rows
if region is None:
raise ValueError("G-optimality needs the region's model rows.")
best_rows, best_value = np.empty(0, dtype=int), -np.inf
for _ in range(_N_STARTS):
rows = _greedy_start(f_cand, f_fixed, n_free, rng)
weights = np.full(len(region), 1.0 / len(region))
start_rows, start_value = rows.copy(), -np.inf
for _ in range(_G_ROUNDS):
rows = _climb(f_cand, f_fixed, rows, Criterion("i_optimal", (region * weights[:, None]).T @ region))
x = np.vstack([f_fixed, f_cand[rows]])
if np.linalg.matrix_rank(x) < x.shape[1]:
break
variance = _prediction_variance(region, np.linalg.inv(x.T @ x))
if -variance.max() > start_value:
start_rows, start_value = rows.copy(), -float(variance.max())
weights = weights * variance
weights /= weights.sum()
if np.isfinite(start_value):
start_rows = _g_polish(f_cand, f_fixed, start_rows, region)
x = np.vstack([f_fixed, f_cand[start_rows]])
value = criterion.value(x.T @ x)
if value > best_value or len(best_rows) == 0:
best_rows, best_value = start_rows, value
return best_rows, best_value
def _n_starts(criterion: Criterion, f_cand: np.ndarray, n_free: int) -> int:
"""Random starts for the exchange: more for a small problem, so long as their total work stays near ``_START_WORK``.
The work of one start grows with candidates x coefficients x free runs. E and K
score every swap at each step and keep ``_N_STARTS``.
"""
if criterion.name not in _CLOSED_FORM_CRITERIA:
return _N_STARTS
work = f_cand.shape[0] * f_cand.shape[1] * max(n_free, 1)
return int(np.clip(_START_WORK // work, _N_STARTS, _MAX_STARTS))
[docs]
def fedorov_exchange(
f_cand: np.ndarray,
n_free: int,
f_fixed: np.ndarray,
rng: np.random.Generator,
criterion: Criterion | None = None,
) -> tuple[np.ndarray, float]:
"""Choose ``n_free`` candidate rows that maximise ``criterion`` (D-optimality by default).
For D, A and I each design run in turn is swapped for the candidate that improves
the criterion most (:func:`_row_exchange`, the modified Fedorov exchange of Cook and
Nachtsheim 1980), with the gains in the closed form of :meth:`Criterion.swap_gains`
(Fedorov 1972). E and K make the best single swap over all pairs at each step. Rows in
``f_fixed`` stay in the design and are never swapped out. Candidates may repeat,
which gives replicated runs where the criterion wants them. The best of several
random starts is kept: ``_N_STARTS`` for E and K, and up to ``_MAX_STARTS`` for a
small D, A or I problem (see :func:`_n_starts`).
Returns
-------
tuple[np.ndarray, float]
Selected candidate indices and the criterion value of the full design (see
:meth:`Criterion.value`). G-optimality runs its own search, :func:`_g_exchange`.
"""
criterion = criterion if criterion is not None else Criterion.d()
if criterion.name == "g_optimal":
return _g_exchange(f_cand, n_free, f_fixed, rng, criterion)
best_rows, best_value = np.empty(0, dtype=int), -np.inf
for rows in _climbed_starts(f_cand, n_free, f_fixed, rng, criterion):
x = np.vstack([f_fixed, f_cand[rows]])
value = criterion.value(x.T @ x)
if value > best_value or len(best_rows) == 0: # keep a design even if every start is singular
best_rows, best_value = rows.copy(), value
return best_rows, best_value
def _climbed_starts(
f_cand: np.ndarray, n_free: int, f_fixed: np.ndarray, rng: np.random.Generator, criterion: Criterion
) -> Iterator[np.ndarray]:
"""Each random start of the exchange (see :func:`_n_starts`), climbed to a local optimum of ``criterion``."""
for _ in range(_n_starts(criterion, f_cand, n_free)):
start = _greedy_start(f_cand, f_fixed, n_free, rng)
if criterion.segment is not None:
start = _spread_over_blocks(start, criterion.segment)
yield _climb(f_cand, f_fixed, start, criterion)
def _spread_over_blocks(rows: np.ndarray, segment: np.ndarray) -> np.ndarray:
"""Deal a start's runs out to the blocks in turn, so the block sizes differ by at most one.
A blocked exchange's candidates are one copy of the points per block, laid end to
end in block order (``segment`` holds each row's block), so a point moves to the
``b``-th block by moving to the same row of the ``b``-th copy.
"""
n_blocks = len(np.unique(segment))
n_points = len(segment) // n_blocks
return rows % n_points + n_points * (np.arange(len(rows)) % n_blocks)
def _climb(f_cand: np.ndarray, f_fixed: np.ndarray, rows: np.ndarray, criterion: Criterion) -> np.ndarray:
"""Improve ``rows`` by swaps while ``criterion`` improves; return the rows.
D, A and I (whose swaps have a closed form) use the row-wise exchange,
:func:`_row_exchange`. E and K make the best single swap at each step, in each of
their phases.
"""
if criterion.name in _CLOSED_FORM_CRITERIA:
return _row_exchange(f_cand, f_fixed, rows, criterion)
rows = rows.copy()
for phase in range(criterion.n_phases):
for _ in range(_MAX_EXCHANGES):
x = np.vstack([f_fixed, f_cand[rows]])
i, j, gain = criterion.best_swap(x.T @ x, f_cand[rows], f_cand, phase)
if not gain > 1e-9: # also stops on NaN
break
rows[i] = j
return rows
[docs]
@dataclass
class PolishLattice:
"""Where :func:`polish_runs` may move a run.
Parameters
----------
n_levels : int
Each continuous coordinate (coded) may move to one of ``n_levels`` evenly spaced
levels from -1 to +1.
model_rows : Callable[[np.ndarray, np.ndarray], np.ndarray]
Model rows of points given as coded continuous values and categorical level indices.
feasible : Callable[[np.ndarray, np.ndarray], np.ndarray] or None
True for each point (coded continuous values, categorical level indices) that
satisfies the constraints; None for the whole box.
"""
n_levels: int
model_rows: Callable[[np.ndarray, np.ndarray], np.ndarray]
feasible: Callable[[np.ndarray, np.ndarray], np.ndarray] | None = None
#: Levels per factor of the lattice a design from a coarser or sampled grid is polished on.
_POLISH_LEVELS = 5
[docs]
def polish_levels(counts: dict) -> int | None:
"""Levels of the lattice to polish on (see :func:`polish_runs`), or None when the candidate grid already holds it.
A grid coarser than ``_POLISH_LEVELS`` levels, or a sampled one, is polished on the
``_POLISH_LEVELS``-level lattice: runs may then use settings between the grid levels,
which can raise the criterion even for D-optimality (a few runs at +/-0.5 when some
runs are fixed, for instance).
"""
return _POLISH_LEVELS if counts["grid_sampled"] or counts["n_levels"] < _POLISH_LEVELS else None
def _coordinate_moves(coded: np.ndarray, levels: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Every run with one continuous coordinate moved to another level: the points, and the run each came from."""
n_runs, n_levels = len(coded), len(levels)
owner = np.repeat(np.arange(n_runs), n_levels)
points, owners = [np.empty((0, coded.shape[1]))], [np.empty(0, dtype=int)]
for axis in range(coded.shape[1]):
moved = coded[owner].copy()
moved[:, axis] = np.tile(levels, n_runs)
changed = ~np.isclose(moved[:, axis], coded[owner, axis])
points.append(moved[changed])
owners.append(owner[changed])
return np.vstack(points), np.concatenate(owners)
[docs]
def polish_runs(
coded: np.ndarray, cats: np.ndarray, f_fixed: np.ndarray, criterion: Criterion, lattice: PolishLattice
) -> tuple[np.ndarray, np.ndarray, int]:
"""Coordinate exchange on a finer or complete grid, after the candidate exchange on a coarse or sampled one.
Each run in turn is moved to the best of its neighbours, the points that differ
from it in one continuous coordinate (set to another level of ``lattice``), when
that improves ``criterion``; passes repeat until no move helps. The gains are those
of :meth:`Criterion.swap_gains`, so only D, A and I are polished. This reaches grid
points the candidate set left out, at a cost that grows with the number of runs
and factors and not with the size of the grid.
Parameters
----------
coded : np.ndarray
The free runs, continuous factors in coded units, shape ``(n, k_cont)``.
cats : np.ndarray
Their categorical level indices, shape ``(n, k_cat)``; these are not moved.
f_fixed : np.ndarray
Model rows of the runs that stay in the design.
criterion : Criterion
D, A or I.
lattice : PolishLattice
The levels, constraints and model rows.
Returns
-------
tuple[np.ndarray, np.ndarray, int]
The polished runs (coded, categorical) and the number of moves made.
"""
coded, cats = coded.copy(), cats.copy()
levels = np.linspace(-1.0, 1.0, lattice.n_levels)
x = np.array(lattice.model_rows(coded, cats), dtype=float) # a copy: patsy's rows can be read-only under pandas 3
n_moves = 0
for _ in range(_MAX_PASSES):
points, owner = _coordinate_moves(coded, levels)
if lattice.feasible is not None:
keep = lattice.feasible(points, cats[owner])
points, owner = points[keep], owner[keep]
rows = lattice.model_rows(points, cats[owner])
moved, start, size = False, 0, 1
m_inv = _regularised_inverse(f_fixed.T @ f_fixed + x.T @ x)
while start < len(coded): # runs a few at a time, as in _climb_block: one after a move, doubling after none
group = np.flatnonzero((owner >= start) & (owner < start + size))
gains = _paired_gains(criterion, m_inv, x[owner[group]], rows[group])
improving = gains > 1e-9 # also refuses NaN
if not improving.any():
start, size = start + size, 2 * size
continue
run = int(owner[group[improving]].min()) # the first run with a move that helps takes its best one
best = int(group[np.argmax(np.where(owner[group] == run, gains, -np.inf))])
m_inv = _swapped_inverse(m_inv, x[run], rows[best])
coded[run], x[run] = points[best], rows[best]
moved, start, size, n_moves = True, run + 1, 1, n_moves + 1
if not moved:
break
return coded, cats, n_moves
def _swapped_inverse(m_inv: np.ndarray, removed: np.ndarray, added: np.ndarray) -> np.ndarray:
"""``M^-1`` after the model row ``removed`` is replaced by ``added``: two Sherman-Morrison updates."""
for v, sign in ((added, 1.0), (removed, -1.0)):
u = m_inv @ v
m_inv = m_inv - sign / (1.0 + sign * float(v @ u)) * np.outer(u, u)
return m_inv
def _paired_gains(criterion: Criterion, m_inv: np.ndarray, x_out: np.ndarray, x_in: np.ndarray) -> np.ndarray:
"""Gain of swapping each row of ``x_out`` for the matching row of ``x_in`` (see :meth:`Criterion.swap_gains`)."""
a_out, a_in = x_out @ m_inv, x_in @ m_inv
d = (np.einsum("ij,ij->i", a_out, x_out), np.einsum("ij,ij->i", a_in, x_in), np.einsum("ij,ij->i", a_out, x_in))
if criterion.weights is None:
return _d_gain(*d)
b_out = a_out @ criterion.weights
b = (
np.einsum("ij,ij->i", b_out, a_out),
np.einsum("ij,ij->i", a_in @ criterion.weights, a_in),
np.einsum("ij,ij->i", b_out, a_in),
)
return _trace_gain(d, b)
[docs]
@dataclass
class CandidatePool:
"""Candidate points, their model rows, and the lattice their runs may be polished on.
Parameters
----------
coded : np.ndarray
Continuous factor settings in coded units, shape ``(N, k_cont)``.
cats : np.ndarray
Categorical level indices, shape ``(N, k_cat)``.
rows : np.ndarray
Model rows of the candidates, shape ``(N, p)``.
lattice : PolishLattice or None
Where each start's design is polished (:func:`polish_runs`); None for no polish.
"""
coded: np.ndarray
cats: np.ndarray
rows: np.ndarray
lattice: PolishLattice | None = None
[docs]
def select_runs(
pool: CandidatePool, n_free: int, f_fixed: np.ndarray, rng: np.random.Generator, criterion: Criterion
) -> tuple[np.ndarray | None, np.ndarray, np.ndarray, float]:
"""Choose ``n_free`` runs by :func:`fedorov_exchange`, polishing each start's design when the pool says so.
The polish can move a run off the candidate set, and the start whose design is best
*after* polishing is kept: on the 10-factor augmentation of the 1.97.0 review, the
best design before polishing was not the best one after it in three seeds out of four.
Returns
-------
tuple
The candidate indices chosen (None once polished, as runs may then lie between
candidates), the runs' coded settings and categorical levels, and the criterion
value of the full design.
"""
lattice = pool.lattice
if lattice is None or criterion.name not in _CLOSED_FORM_CRITERIA:
rows, value = fedorov_exchange(pool.rows, n_free, f_fixed, rng, criterion)
return rows, pool.coded[rows], pool.cats[rows], value
best: tuple[np.ndarray | None, np.ndarray, np.ndarray, float] = (None, pool.coded[:0], pool.cats[:0], -np.inf)
for rows in _climbed_starts(pool.rows, n_free, f_fixed, rng, criterion):
coded, cats, _moves = polish_runs(pool.coded[rows], pool.cats[rows], f_fixed, criterion, lattice)
x = np.vstack([f_fixed, lattice.model_rows(coded, cats)])
value = criterion.value(x.T @ x)
if value > best[3] or not len(best[1]): # keep a design even if every start is singular
best = (None, coded, cats, value)
return best
# ---------------------------------------------------------------------------
# Public entry point
# ---------------------------------------------------------------------------
[docs]
@dataclass
class ConstrainedOptions:
"""Optional settings for :func:`constrained_optimal_design` and the constrained mixture design.
Parameters
----------
model_type : str
``"main_effects"``, ``"interactions"`` or ``"quadratic"`` (a mixture design maps
these to Scheffé models).
criterion : str
``"d_optimal"`` (default), ``"i_optimal"``, ``"a_optimal"``, ``"e_optimal"``,
``"g_optimal"`` or ``"k_optimal"``.
fixed_runs : pandas.DataFrame or None
Runs kept in the design (continuous in coded units, categorical as labels),
already validated by the caller. They count towards the budget.
n_levels : int or None
Grid levels per continuous factor. ``None`` picks the finest of 5, 4 and 3
(and 2, for a model without squared terms; 4 is skipped for a quadratic model
unless the criterion is I or G) that keeps the grid within
``max(5000, 200 p)`` points for ``p`` model coefficients, and samples the
coarsest grid down to that when none does. A sampled or coarse grid is then
polished (see :func:`polish_runs`).
candidates : pandas.DataFrame or None
Settings the runs must be chosen from, in actual units (proportions for a
mixture), one column per factor, instead of a generated grid. Rows that
break a constraint are dropped, and the I-optimality average is taken over
the remaining rows.
n_blocks : int or None
Blocks the design is run in. For D, A and I the blocks enter the exchange as
fixed effects (see :meth:`Criterion.blocked`), and ``metadata["block_labels"]``
gives each run's block; the fixed runs, already made, are a block of their
own. Other criteria leave the blocks to be assigned afterwards.
"""
model_type: str = "interactions"
criterion: str = "d_optimal"
fixed_runs: pd.DataFrame | None = None
n_levels: int | None = None
candidates: pd.DataFrame | None = None
n_blocks: int | None = None
[docs]
def block_columns(labels: np.ndarray, n_blocks: int) -> np.ndarray:
"""Return sum-to-zero columns for the 0-based block ``labels``: ``n_blocks - 1`` of them, the last block -1 in each.
``analyze_experiment`` codes blocks this way, so a criterion over these columns is
the one the analysis estimates with, and the intercept is the average block.
"""
coding = np.vstack([np.eye(n_blocks - 1), -np.ones((1, n_blocks - 1))])
return coding[np.asarray(labels, dtype=int)]
def _blocks_in_exchange(opts: ConstrainedOptions, criterion: Criterion) -> bool:
"""Whether the blocks enter the exchange: there are two or more, and the criterion is D, A or I."""
return opts.n_blocks is not None and opts.n_blocks > 1 and criterion.name in _CLOSED_FORM_CRITERIA
def _blocked_search(
pool: CandidatePool, f_fixed: np.ndarray, criterion: Criterion, n_blocks: int, budget: int
) -> tuple[CandidatePool, np.ndarray, Criterion, int]:
"""Set the exchange up for a design run in ``n_blocks`` blocks: return its pool, fixed rows, criterion and budget.
The fixed runs, already made, are a block of their own, the first; the free runs
fill the others in near-equal shares. Each block the free runs fill gets a copy of
the candidates, with that block's columns (:func:`block_columns`) and its index as
one more categorical column, which the polish never moves.
Raises
------
ValueError
If the block effects leave the budget too few runs to estimate the model, or the
budget leaves fewer than two new runs for each block they fill.
"""
n_fixed = len(f_fixed)
free_blocks = np.arange(1 if n_fixed else 0, n_blocks)
f_fixed = np.hstack([f_fixed, block_columns(np.zeros(n_fixed, dtype=int), n_blocks)])
n_parameters = pool.rows.shape[1] + n_blocks - 1
needed = n_fixed + n_parameters - (int(np.linalg.matrix_rank(f_fixed)) if n_fixed else 0)
if budget < needed:
raise ValueError(
f"Running the design in {n_blocks} blocks adds {n_blocks - 1} block effect(s) to the model's "
f"{pool.rows.shape[1]} coefficients, so it needs at least {needed} runs, more than budget={budget}. "
f"Use budget >= {needed} or fewer blocks."
)
if budget - n_fixed < 2 * len(free_blocks):
raise ValueError(
f"{n_blocks} blocks need at least two new runs in each of the {len(free_blocks)} blocks the new runs "
f"fill, but the budget leaves {budget - n_fixed}. Use fewer blocks or a larger budget."
)
segment = np.repeat(free_blocks, len(pool.rows))
copies = len(free_blocks)
lattice = pool.lattice
if lattice is not None:
factor_rows = lattice.model_rows
lattice = replace(
lattice,
model_rows=lambda coded, cats: np.hstack(
[factor_rows(coded, cats[:, :-1]), block_columns(cats[:, -1], n_blocks)]
),
)
blocked_pool = CandidatePool(
np.tile(pool.coded, (copies, 1)),
np.column_stack([np.tile(pool.cats, (copies, 1)), segment]),
np.hstack([np.tile(pool.rows, (copies, 1)), block_columns(segment, n_blocks)]),
lattice,
)
return blocked_pool, f_fixed, criterion.blocked(n_blocks - 1, segment), budget
def _adjusted_for_blocks(criterion: Criterion, value: float, labels: np.ndarray, n_blocks: int) -> float:
"""Return ``value`` for the factor effects adjusted for the blocks (see :meth:`Criterion.blocked`).
For D, ``log det(M)`` less ``log det(Z'Z)`` of the block columns ``Z`` is
``log det(S)``; the trace criteria already score ``S``.
"""
if criterion.name != "d_optimal":
return value
z = block_columns(labels, n_blocks)
return value - float(np.linalg.slogdet(z.T @ z)[1])
def _coded_runs(region: _Region, runs: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
"""Return runs given continuous in coded units and categorical as labels, as coded values and level indices."""
coded = runs[[f.name for f in region.continuous]].to_numpy(dtype=float)
labels = {f.name: [str(lv) for lv in f.levels or []] for f in region.categorical}
cats = np.array(
[[labels[f.name].index(str(v)) for v in runs[f.name]] for f in region.categorical], dtype=int
).T.reshape(len(runs), len(region.categorical))
return coded, cats
def _fixed_rows(region: _Region, fixed_runs: pd.DataFrame | None, model_type: str, n_columns: int) -> np.ndarray:
"""Model rows for the fixed runs (none when ``fixed_runs`` is None), warning about any outside the region."""
if fixed_runs is None:
return np.empty((0, n_columns))
fixed_coded, fixed_cats = _coded_runs(region, fixed_runs)
n_outside = int((region.slack(fixed_coded, fixed_cats) > _FEASIBILITY_TOL).sum())
if n_outside:
logger.warning("%d fixed run(s) lie outside the constrained region; they are kept as given.", n_outside)
return model_matrix(region, fixed_coded, fixed_cats, model_type)
def _fixed_run_metadata(region: _Region, fixed_runs: pd.DataFrame) -> dict:
"""Metadata on the fixed runs: how many, and how many lie outside the region."""
meta = {"n_fixed_runs": len(fixed_runs)}
coded, cats = _coded_runs(region, fixed_runs)
n_outside = int((region.slack(coded, cats) > _FEASIBILITY_TOL).sum())
if n_outside:
meta["n_fixed_runs_outside_region"] = n_outside
return meta
def _user_candidates(region: _Region, candidates: pd.DataFrame) -> tuple[np.ndarray, np.ndarray, dict, list]:
"""Code a user's candidate settings, drop the infeasible and repeated ones.
Returns
-------
tuple
Coded continuous values, categorical level indices, counts for the
metadata, and the ``candidates`` index label of each kept row.
"""
names = [f.name for f in region.continuous] + [f.name for f in region.categorical]
missing = [n for n in names if n not in candidates.columns]
if missing:
raise ValueError(f"candidates is missing columns for factors: {missing}.")
actual = candidates[[f.name for f in region.continuous]].apply(pd.to_numeric, errors="coerce")
if actual.isna().any().any():
raise ValueError("candidates has missing or non-numeric values in a continuous factor column.")
low = np.array([f.low for f in region.continuous], dtype=float)
high = np.array([f.high for f in region.continuous], dtype=float)
coded = (actual.to_numpy(dtype=float) - (low + high) / 2.0) / ((high - low) / 2.0)
cats = np.zeros((len(candidates), len(region.categorical)), dtype=int)
for j, f in enumerate(region.categorical):
labels = [str(lv) for lv in f.levels or []]
values = candidates[f.name].astype(str)
unknown = sorted(set(values) - set(labels))
if unknown:
raise ValueError(f"candidates has unknown levels {unknown} for categorical factor {f.name!r}.")
cats[:, j] = [labels.index(v) for v in values]
if (np.abs(coded) > 1.0 + 1e-9).any():
logger.warning("Some candidates lie outside the factors' low/high range; they are kept (coded beyond +/-1).")
feasible = region.slack(coded, cats) <= _FEASIBILITY_TOL
_, first = np.unique(np.hstack([coded.round(9), cats]), axis=0, return_index=True)
keep = np.zeros(len(coded), dtype=bool)
keep[first] = True
keep &= feasible
counts = {
"n_candidates_supplied": len(candidates),
"n_candidates_infeasible": int((~feasible).sum()),
"n_candidates": int(keep.sum()),
}
return coded[keep], cats[keep], counts, list(candidates.index[keep])
def _corners_suffice(region: _Region, opts: ConstrainedOptions) -> bool:
"""Whether some D-optimal design uses only the corners of the box: no squared terms and no constraints."""
return (
opts.criterion == "d_optimal"
and opts.model_type != "quadratic"
and not region.inequalities
and not region.conditional
)
def _polish_lattice(region: _Region, opts: ConstrainedOptions, counts: dict) -> PolishLattice | None:
"""Return the lattice to polish a design from a generated grid on, or None when there is nothing to add.
See :func:`polish_levels`. Supplied candidates and an explicit ``n_levels`` are kept
to, and a design that only needs the box's corners is not polished.
"""
if opts.candidates is not None or opts.n_levels is not None or opts.criterion not in _CLOSED_FORM_CRITERIA:
return None
n_levels = polish_levels(counts)
if n_levels is None or _corners_suffice(region, opts):
return None
return PolishLattice(
n_levels,
lambda coded, cats: model_matrix(region, coded, cats, opts.model_type),
lambda coded, cats: region.slack(coded, cats) <= _FEASIBILITY_TOL,
)
def _candidate_pool(
region: _Region, opts: ConstrainedOptions, rng: np.random.Generator
) -> tuple[np.ndarray, np.ndarray, dict, list | None, Callable[[], np.ndarray]]:
"""Return the candidates (a generated grid, or the user's), and the rows I-optimality averages over.
Raises
------
ValueError
If no candidate point satisfies the constraints.
"""
coded, cats, counts, labels, region_rows = _candidates_and_region(region, opts, rng)
if counts["n_candidates"] == 0:
raise ValueError("No candidate point satisfies all the constraints; check them for conflicts.")
return coded, cats, counts, labels, region_rows
def _candidates_and_region(
region: _Region, opts: ConstrainedOptions, rng: np.random.Generator
) -> tuple[np.ndarray, np.ndarray, dict, list | None, Callable[[], np.ndarray]]:
"""See :func:`_candidate_pool`, which checks that some candidate is feasible."""
if opts.candidates is not None:
coded, cats, counts, labels = _user_candidates(region, opts.candidates)
return coded, cats, counts, labels, lambda: model_matrix(region, coded, cats, opts.model_type)
n_levels = opts.n_levels
if n_levels is None and _corners_suffice(region, opts):
# Without squared terms each model row is affine in every single coordinate, so det(X'X) is a convex
# quadratic in it and is largest at -1 or +1: some exact D-optimal design uses only the box's corners.
n_levels = 2
coded, cats, counts = build_candidates(region, n_levels, opts.model_type, opts.criterion)
return coded, cats, counts, None, lambda: _uniform_rows(region, opts.model_type, rng, coded, cats)
[docs]
def selection_counts(labels: list, rows: np.ndarray) -> dict[str, int]:
"""Count how often each supplied candidate (by index label) was chosen."""
chosen = pd.Series([labels[r] for r in rows]).value_counts(sort=False)
return {str(label): int(n) for label, n in chosen.items()}
def _missing_rank(f_fixed: np.ndarray, n_parameters: int) -> int:
"""Count the model coefficients the fixed runs leave unestimated: the rank the free runs must supply."""
return n_parameters - (int(np.linalg.matrix_rank(f_fixed)) if len(f_fixed) else 0)
[docs]
def refuse_budget_short_of_fixed_runs(f_fixed: np.ndarray, n_parameters: int, budget: int) -> None:
"""Raise when the budget leaves too few new runs to supply the rank the fixed runs lack.
Fixed runs that repeat a point (three centre runs, say) carry less information than
their count, so a budget above the number of coefficients can still be too small.
Raises
------
ValueError
If ``budget`` minus the fixed runs is below the rank they lack.
"""
n_fixed = len(f_fixed)
missing_rank = _missing_rank(f_fixed, n_parameters)
if budget - n_fixed < missing_rank:
raise ValueError(
f"The {n_fixed} fixed run(s) span only {n_parameters - missing_rank} of the {n_parameters} model "
f"coefficients, so at least {missing_rank} more run(s) are needed: budget >= {n_fixed + missing_rank}, "
f"not {budget}."
)
#: Uniform draws used to estimate the region's moment matrix for I-optimality.
_N_MOMENT_SAMPLES = 20_000
def _uniform_rows(
region: _Region,
model_type: str,
rng: np.random.Generator,
seeds: np.ndarray,
seed_cats: np.ndarray | None = None,
) -> np.ndarray:
"""Model rows at points drawn uniformly from the feasible region.
Thin regions are handled by hit-and-run started from the feasible candidates
(``seeds``, whose level indices are ``seed_cats``), so the moment matrix covers the
whole region, not the few points a rejection pass happens to hit.
With rules, each allowed combination of levels carries the weight of its volume (see
:meth:`_Region.draw_levels`), and the continuous settings of a run are drawn
uniformly from the region its levels leave: the box cut by the inequalities that
apply there. Without rules, every combination of levels is allowed and the region
is the same for each.
"""
if region.rules is None:
coded = _uniform_points(region, _N_MOMENT_SAMPLES, rng, seeds)
cats = np.empty((len(coded), 0), dtype=int)
if region.categorical:
cats = np.column_stack([rng.integers(len(f.levels or []), size=len(coded)) for f in region.categorical])
return model_matrix(region, coded, cats, model_type)
cats = region.draw_levels(_N_MOMENT_SAMPLES, rng)
applies = region.applies(cats)
seed_applies = region.applies(seed_cats) if seed_cats is not None else None
coded = np.empty((len(cats), len(region.continuous)))
for signature in np.unique(applies, axis=0):
rows = np.all(applies == signature, axis=1)
mine = seeds if seed_applies is None else seeds[np.all(seed_applies == signature, axis=1)]
coded[rows] = _uniform_points(region.given(signature), int(rows.sum()), rng, mine)
return model_matrix(region, coded, cats, model_type)
def _uniform_points(region: _Region, n: int, rng: np.random.Generator, seeds: np.ndarray) -> np.ndarray:
"""Draw ``n`` coded points uniformly from the box cut by the region's inequalities (its rules are not read)."""
k = len(region.continuous)
def coded_inequality(g: Callable) -> Callable[[np.ndarray], np.ndarray]:
return lambda x: np.broadcast_to(np.asarray(g(region.actual(x)), dtype=float), (len(x),))
sampler = UniformSampler(
[coded_inequality(g) for g in region.inequalities],
(np.full(k, -1.0), np.full(k, 1.0)),
lambda m: rng.uniform(-1.0, 1.0, size=(m, k)),
seeds=lambda: seeds,
tolerances=region._value_scales()[: len(region.inequalities)], # so a drawn point passes slack()
)
return sampler.draw(n, rng)
#: Uniform draws added to the candidates for the G-criterion's maximum (the candidates hold the boundary points).
_N_G_REGION_SAMPLES = 3000
[docs]
def make_criterion(
name: str,
n_parameters: int,
region_rows: Callable[[], np.ndarray],
candidate_rows: np.ndarray | None = None,
) -> Criterion:
"""Build the :class:`Criterion` called ``name``.
``region_rows`` (model rows drawn uniformly in the region) is only evaluated for
I- and G-optimality. G takes its maximum over those rows plus ``candidate_rows``,
which hold the vertices and boundary points where the worst case usually sits.
"""
if name == "d_optimal":
return Criterion.d()
if name == "a_optimal":
return Criterion.a(n_parameters)
if name == "i_optimal":
return Criterion.i(region_rows())
if name == "e_optimal":
return Criterion.e()
if name == "k_optimal":
return Criterion.k()
if name == "g_optimal":
sample = region_rows()[:_N_G_REGION_SAMPLES]
rows = sample if candidate_rows is None else np.vstack([candidate_rows, sample])
return Criterion.g(np.unique(rows.round(12), axis=0))
raise ValueError(
f"Unknown criterion {name!r}; choose 'd_optimal', 'i_optimal', 'a_optimal', 'e_optimal', 'g_optimal' or "
"'k_optimal'."
)
def _check_factors_and_model(factors: list[Factor], model_type: str) -> None:
"""Refuse mixture components and a model type the model matrix does not build."""
if model_type not in _MODEL_TYPES:
raise ValueError(f"model_type={model_type!r} is not supported; choose from {', '.join(_MODEL_TYPES)}.")
if any(f.type == FactorType.mixture for f in factors):
raise ValueError(
"Mixture components need the mixture engine: give only mixture factors, and generate_design routes "
"them there. Mixture-process designs are not supported."
)
[docs]
def constrained_optimal_design(
factors: list[Factor],
budget: int,
constraints: list[Constraint],
options: ConstrainedOptions | None = None,
random_state: int | np.random.Generator | None = None,
) -> tuple[np.ndarray, dict]:
"""Generate a D-, I-, A- or E-optimal design from a candidate set, every chosen run satisfying ``constraints``.
This is also the optimal-design backend without constraints, unless the design is a
split-plot design or pyoptex is asked for. Runs in ``options.fixed_runs`` are kept as given, even
outside the region (with a warning); ``metadata["n_fixed_runs_outside_region"]``
counts them.
Parameters
----------
factors : list[Factor]
Continuous and categorical factors. Mixture factors are not supported here.
budget : int
Total number of runs, including ``fixed_runs``. A budget that leaves too few free
runs for the rank the fixed runs lack, or for the block effects, raises a ``ValueError``.
constraints : list[Constraint]
Inequalities in actual units over the continuous factors, e.g.
``Constraint(expression="3*T + 5*D <= 600")``. May be empty.
options : ConstrainedOptions or None
Model type, criterion, fixed runs and grid resolution; defaults to a
D-optimal design for an interactions model, no fixed runs and an automatic grid.
random_state : int, numpy.random.Generator or None
Seed for the random starts of the exchange.
Returns
-------
tuple[np.ndarray, dict]
The design in coded units (categorical factors as labels) and metadata.
Raises
------
ValueError
If ``model_type`` is unknown, a constraint cannot be parsed, no candidate point
is feasible, or the feasible region cannot support the requested model.
"""
rng = check_random_state(random_state)
opts = options if options is not None else ConstrainedOptions()
model_type, fixed_runs = opts.model_type, opts.fixed_runs
_check_factors_and_model(factors, model_type)
region = region_with_rules(factors, constraints)
coded, cats, counts, labels, region_rows = _candidate_pool(region, opts, rng)
f_cand, inestimable = _candidate_rows(region, coded, cats, model_type)
f_fixed = _fixed_rows(region, fixed_runs, model_type, f_cand.shape[1])
n_fixed = len(f_fixed)
_check_candidates_support_model(counts, f_fixed, f_cand, model_type)
p = f_cand.shape[1]
refuse_budget_short_of_fixed_runs(f_fixed, p, budget)
criterion = make_criterion(opts.criterion, p, region_rows, f_cand)
pool = CandidatePool(coded, cats, f_cand, _polish_lattice(region, opts, counts))
block_labels = None
if _blocks_in_exchange(opts, criterion):
n_blocks = int(opts.n_blocks) # type: ignore[arg-type]
search, search_fixed, criterion, budget = _blocked_search(pool, f_fixed, criterion, n_blocks, budget)
rows, chosen, chosen_cats, value = select_runs(search, budget - n_fixed, search_fixed, rng, criterion)
block_labels = np.concatenate([np.zeros(n_fixed, dtype=int), chosen_cats[:, -1]])
rows, chosen_cats = (None if rows is None else rows % len(f_cand)), chosen_cats[:, :-1]
value = _adjusted_for_blocks(criterion, value, block_labels, n_blocks)
else:
rows, chosen, chosen_cats, value = select_runs(pool, budget - n_fixed, f_fixed, rng, criterion)
chosen_rows = model_matrix(region, chosen, chosen_cats, model_type)
if len(chosen) != budget - n_fixed or np.linalg.matrix_rank(np.vstack([f_fixed, chosen_rows])) < p:
raise ValueError(
f"No design of {budget} runs from these candidates can estimate the '{model_type}' model "
f"({p} coefficients). Use a simpler model or more distinct candidate points."
)
if pool.lattice is not None:
counts["polish_levels"] = pool.lattice.n_levels
meta = {
"backend": "candidate_exchange",
"optimality_criterion": opts.criterion,
"model_type": model_type,
**counts,
}
if constraints:
meta["constraints"] = [c.expression for c in constraints]
meta["constraints_enforced"] = True
if inestimable:
meta["inestimable_columns"] = inestimable
meta.update(criterion_metadata(criterion, value))
if block_labels is not None:
meta["block_labels"] = (block_labels + 1).tolist()
if fixed_runs is not None:
meta.update(_fixed_run_metadata(region, fixed_runs))
if labels is not None and rows is not None:
meta["candidate_source"] = "user"
meta["selected_candidates"] = selection_counts(labels, rows)
return _assemble_design(region, factors, chosen, chosen_cats, fixed_runs), meta
[docs]
def region_with_rules(factors: list[Factor], constraints: list[Constraint]) -> _Region:
"""Return the region of continuous and categorical ``factors`` cut by ``constraints``, their rules included.
A constraint that names a categorical factor, or joins comparisons with ``and``,
``or`` or ``not``, is read by :func:`~process_improve.experiments._rules.compile_rules`;
every other one by :func:`parse_constraint`, as before.
"""
continuous = [f for f in factors if f.type != FactorType.categorical]
categorical = [f for f in factors if f.type == FactorType.categorical]
names = {f.name for f in continuous}
plain, rules = compile_rules(constraints, names, {f.name: list(f.levels or []) for f in categorical})
inequalities = [g for _source, text in plain for g in parse_constraint(text, names)]
if rules is None:
return _Region(continuous, categorical, inequalities)
from process_improve.experiments.region import DesignRegion # noqa: PLC0415 region imports this module
conditional = [parse_constraint(clause.limit or "", names) for clause in rules.limits]
return _Region(continuous, categorical, inequalities, rules, conditional, source=DesignRegion(factors, constraints))
def _candidate_rows(
region: _Region, coded: np.ndarray, cats: np.ndarray, model_type: str
) -> tuple[np.ndarray, dict[str, int]]:
"""Return the candidates' model rows, and per term the columns the rules leave inestimable (none without rules).
With rules, the columns that cannot be estimated are dropped from the model, here and
in every later call of :func:`model_matrix` on ``region``.
"""
f_cand = model_matrix(region, coded, cats, model_type)
if region.rules is None:
return f_cand, {}
_check_every_level_runs(region, cats)
region.kept_columns, inestimable = _estimable_columns(region, f_cand, model_type)
return f_cand[:, region.kept_columns], inestimable
def _check_every_level_runs(region: _Region, cats: np.ndarray) -> None:
"""Refuse candidates that leave a level of a categorical factor with no setting the constraints allow."""
for j, f in enumerate(region.categorical):
levels = list(f.levels or [])
empty = [levels[i] for i in sorted(set(range(len(levels))) - set(cats[:, j].tolist()))]
if empty:
shown = ", ".join(f"{f.name} = {level!r}" for level in empty)
raise ValueError(
f"No candidate the constraints allow has {shown}: the limits that apply there (or the candidates "
"supplied) leave that level without a run. Loosen them, or remove that level from the factor."
)
#: A model column is estimable from the candidates when its part not explained by the columns before it is larger
#: than this share of its length.
_ESTIMABLE_TOL = 1e-8
def _estimable_columns(region: _Region, f_cand: np.ndarray, model_type: str) -> tuple[np.ndarray, dict[str, int]]:
"""Return the model columns the candidates can estimate, and per term the number of columns lost, with a warning.
A combination of levels that is not allowed is a missing cell: with ``catalyst`` at 3
levels and ``solvent`` at 2, one disallowed pair leaves 5 cells for the 6
coefficients of the main effects and their interaction, so one interaction column
cannot be estimated. The columns are kept in model order (intercept, main effects,
interactions, squares) while each adds rank, so the lost columns are the highest-order
ones. Every set of columns kept this way spans the same model space over the
candidates, so the D- and I-optimal designs do not depend on which column of a term
is dropped.
"""
keep = independent_columns(f_cand)
kept = keep.tolist()
terms = _column_terms(region, model_type)
lost: dict[str, int] = {}
for j in sorted(set(range(len(terms))) - set(kept)):
lost[terms[j]] = lost.get(terms[j], 0) + 1
for term, n_lost in lost.items():
width = terms.count(term)
logger.warning(
"The constraints leave the %s term with %d of its %d column(s) estimable%s; the design is built for "
"the columns that can be estimated, and metadata['inestimable_columns'] records the %d lost.",
term,
width - n_lost,
width,
_missing_cells(region, term),
n_lost,
)
return keep, lost
[docs]
def independent_columns(matrix: np.ndarray) -> np.ndarray:
"""Return the columns of ``matrix`` kept, in order, while each adds rank: a basis of its column space.
Worked on the triangular factor of a QR decomposition, which has the column space's
geometry at the size of the number of columns. A column is kept when its part not
explained by the columns kept before it is longer than ``_ESTIMABLE_TOL`` of its own
length.
"""
r = np.linalg.qr(matrix, mode="r")
basis = np.empty((r.shape[0], 0))
kept = []
for j in range(r.shape[1]):
column = r[:, j]
residual = column - basis @ (basis.T @ column)
size = float(np.linalg.norm(residual))
if size > _ESTIMABLE_TOL * max(float(np.linalg.norm(column)), 1.0):
basis = np.column_stack([basis, residual / size])
kept.append(j)
return np.array(kept, dtype=int)
def _column_terms(region: _Region, model_type: str) -> list[str]:
"""Return the term each column of :func:`model_matrix` belongs to, before any column is dropped."""
blocks = [(f.name, 1) for f in region.continuous] + [(f.name, len(f.levels or []) - 1) for f in region.categorical]
terms = ["Intercept"] + [name for name, width in blocks for _ in range(width)]
if model_type in {"interactions", "quadratic"}:
for i, (first, width_i) in enumerate(blocks):
for second, width_j in blocks[i + 1 :]:
terms += [f"{first}:{second}"] * (width_i * width_j)
if model_type == "quadratic":
terms += [f"{f.name}^2" for f in region.continuous]
return terms
def _missing_cells(region: _Region, term: str) -> str:
"""Name the pairs of levels no allowed run has, for an interaction of two factors the rules name."""
if region.rules is None or term.count(":") != 1:
return ""
first, second = term.split(":")
if first not in region.rules.levels or second not in region.rules.levels:
return ""
allowed, _ = region.rules.allowed_table()
seen = set(zip(allowed[first], allowed[second], strict=True))
missing = [
f"{first} = {a!r} with {second} = {b!r}"
for a in region.rules.levels[first]
for b in region.rules.levels[second]
if (a, b) not in seen
]
return f" (no allowed run has {'; '.join(missing)})" if missing else ""
def _check_candidates_support_model(counts: dict, f_fixed: np.ndarray, f_cand: np.ndarray, model_type: str) -> None:
"""Raise unless the candidates, with the fixed runs, can estimate the model."""
if np.linalg.matrix_rank(np.vstack([f_fixed, f_cand])) < f_cand.shape[1]:
raise ValueError(
f"The feasible region ({counts['n_candidates']} candidate points) cannot support a "
f"'{model_type}' model with {f_cand.shape[1]} coefficients. Use a simpler model, supply "
"candidates, or loosen the constraints."
)
def _assemble_design(
region: _Region, factors: list[Factor], chosen: np.ndarray, chosen_cats: np.ndarray, fixed_runs: pd.DataFrame | None
) -> np.ndarray:
"""Return the design: the fixed runs, then the chosen ones, continuous factors coded and categorical as labels."""
design = pd.DataFrame(chosen, columns=[f.name for f in region.continuous])
for j, f in enumerate(region.categorical):
design[f.name] = np.asarray(f.levels, dtype=object)[chosen_cats[:, j]]
if fixed_runs is not None:
design = pd.concat([fixed_runs[design.columns].reset_index(drop=True), design], ignore_index=True)
design = design[[f.name for f in factors]]
return design.to_numpy() if region.categorical else design.to_numpy(dtype=float)
