# (c) Kevin Dunn, 2010-2026. MIT License.
"""Recommend the next round of a designed experiment from the results of the rounds so far.
A study is a sequence of rounds: screen the factors, resolve what the screen left open (an
effect aliased with an interaction, curvature it could not model), move towards the optimum
when it lies past the ranges studied, and confirm the predicted optimum. After every round the
question is the same: given every run so far, what should the next round be?
:func:`recommend_next_step` answers it with one action, the evidence for it, the runs to make,
and how to judge them once they are in.
The decision is deterministic and follows the evidence in this order:
1. **Confirmation runs** at the predicted optimum are judged by their mean against the
prediction interval for a mean of that many runs. Confirmed, the study **stops**, unless
the goal still improves past an edge of a range that may move. Missed, the next round adds
what the design could not show (axial runs, a fold-over); else, when the model chosen
again with the missed runs as data puts its optimum elsewhere, that optimum is confirmed;
else a miss that may be chance is confirmed again, judged on the new runs alone, and a
gross or repeated miss calls for a look at the process (**investigate**).
2. A model that leaves fewer than two **degrees of freedom for error** cannot be judged: add
centre runs to a two-level design of numeric factors, else replicate the runs.
3. **Nothing stands out from the noise**: with fewer than four degrees of freedom for error,
add runs first; else widen the ranges, or replicate when they cannot move; and when that
was done before (``history``), stop: no factor moves the response by more than the noise.
4. An active effect is **aliased with a plausible interaction**, so the data cannot say which
of the two is real: fold the design over, or add runs chosen for the model that holds both.
5. **Curvature** (centre runs off the plane of the corners, lack of fit, or a square in the
model) that the design cannot model: add axial runs, which make a face-centred central
composite design, or runs chosen for the full quadratic model. Curvature the centre runs
show could be any factor's, so the axial runs go in every numeric factor (up to six).
6. A **lack of fit** that the full quadratic model in the active factors does not remove: the
model misses effects. Every factor with some evidence of an effect joins the active ones,
and runs are chosen for the full quadratic model in all of them.
7. The best settings sit at an **edge of a range** and the goal improves past it by more than
twice the noise: when the ranges may move, run the path of steepest ascent beyond it.
8. Otherwise, **confirm** at the predicted optimum.
How the terms are judged matters as much as the rules. Each term is tested in the largest
standard model the runs estimate with two degrees of freedom for error (one column per chain of
exact aliases), or by Lenth's method, with simulated critical values, on every contrast of a
saturated two-level design, or by the staged analysis of a definitive screening design. Whether
anything stands out at all is one test at ``alpha`` whose critical value is simulated for the
design (:func:`_stands_out`), so pure noise raises a false alarm in one study in twenty, not in
most of them. A model that does not fit the repeated runs has its terms tested against the
repeats' own variation instead of its residual, which holds what it misses. Interactions of the
active factors (strong heredity) that the judging could not test are then tested in the model of
the active terms.
A factor the selected model leaves out does not change the prediction, so any setting in its
range will do. The runs recommended set it at the end its small, not significant, effect
favours: a setting the study has run, and never a guess.
"""
from __future__ import annotations
import functools
import itertools
import math
import warnings
from dataclasses import dataclass, field
from typing import Any, Literal
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from patsy import dmatrix
from scipy import optimize, stats
from statsmodels.regression.linear_model import RegressionResultsWrapper
from statsmodels.stats.multitest import multipletests
from process_improve._random import check_random_state
from process_improve.experiments._analyses.aliasing import alias_chains_of
from process_improve.experiments._analyses.curvature import _run_curvature_test
from process_improve.experiments._analyses.lack_of_fit import _run_lack_of_fit
from process_improve.experiments._analyses.model_selection import _proposed_terms
from process_improve.experiments._analyses.prediction import _mean_tests
from process_improve.experiments.augment import augment_design
from process_improve.experiments.factor import Factor, FactorType
from process_improve.experiments.optimization import _extract_b_and_B
from process_improve.experiments.true_process import _canonical
#: The actions :func:`recommend_next_step` can recommend.
ACTIONS = (
"stop",
"confirm",
"add_center_points",
"replicate",
"widen_ranges",
"de_alias",
"augment_to_rsm",
"add_runs_for_model",
"move_region",
"investigate",
)
Goal = Literal["maximize", "minimize", "target"]
#: A coded setting this close to +/-1 is at the end of its range.
_AT_EDGE = 1e-6
#: An interaction whose column the model's columns reproduce this well (R squared; a variance inflation
#: factor of 10) cannot be told apart from the model's terms by these runs.
_ALIASED = 0.9
#: An alias weight at least this large names a model term the interaction is confused with.
_ALIAS_WEIGHT = 0.25
#: A confirmation setting within this much of the optimum, in coded units, was made at the optimum.
_SAME_SETTING = 0.1
#: The goal must improve past an edge by this many noise standard deviations, per coded unit, to move.
_WORTH_MOVING = 2.0
#: A confirmation mean more than this many half-widths of its interval from the prediction is no chance miss.
_GROSS_MISS = 3.0
#: Steps along the path of steepest ascent beyond the best settings, in coded units of the factor moving most.
_PATH_STEPS = (0.5, 1.0, 1.5, 2.0)
#: The search for the optimum starts from every corner of up to this many factors, and from this many random points.
_ALL_CORNERS = 6
_RANDOM_STARTS = 8
#: Centre runs added to a two-level design that has none.
_CENTER_RUNS = 3
#: Runs beyond a model's coefficients that runs chosen for it leave, to estimate the noise.
_SPARE_RUNS = 3
#: The fewest degrees of freedom for error with which the model's terms are judged.
_MIN_ERROR_DF = 2
#: Below this many degrees of freedom for error a real effect can fail to stand out for want of them: the
#: t quantile at 0.05 is 2.78 with four, 4.30 with two, against 2.0 with many.
_FEW_ERROR_DF = 4
#: A range widened for lack of any visible effect grows by this factor about its middle.
_WIDEN = 1.5
# --------------------------------------------------------------------------- the study, coded
@dataclass(frozen=True)
class _Study:
"""The factors, and the runs coded with internal names ``x0, x1, ...`` and the response as ``y``."""
factors: list[Factor]
data: pd.DataFrame
response: str
@property
def names(self) -> list[str]:
return [f"x{i}" for i in range(len(self.factors))]
@property
def numeric(self) -> list[bool]:
return [f.type != FactorType.categorical for f in self.factors]
def real(self, text: str) -> str:
"""Write a term or formula with the factors' own names in place of the internal ones."""
for i in reversed(range(len(self.factors))): # x10 before x1
text = text.replace(f"x{i}", self.factors[i].name)
return text.replace("y ~", f"{self.response} ~", 1)
def actual(self, coded: np.ndarray) -> dict[str, Any]:
"""Return coded settings in the factors' own units, or a categorical factor's level."""
out: dict[str, Any] = {}
for f, value in zip(self.factors, coded, strict=True):
if f.type == FactorType.categorical:
out[f.name] = f.levels[1] if value >= 0 else f.levels[0] # type: ignore[index]
else:
mid, half = (f.high + f.low) / 2, (f.high - f.low) / 2 # type: ignore[operator]
out[f.name] = float(mid + value * half)
return out
def _check_factors(factors: list[Factor]) -> None:
if not factors:
raise ValueError("At least one factor is required.")
for f in factors:
if f.type == FactorType.mixture:
raise ValueError(
f"Factor {f.name!r} is a mixture component; recommend_next_step does not take mixtures, whose "
"runs are blends."
)
if f.type == FactorType.categorical and len(f.levels or []) != 2:
raise ValueError(
f"Categorical factor {f.name!r} has {len(f.levels or [])} levels; recommend_next_step takes "
"categorical factors at two levels."
)
def _code(factors: list[Factor], frame: pd.DataFrame, response: str, block: str | None = None) -> pd.DataFrame:
"""Return the runs coded to -1/+1 under the internal names, the response as ``y`` and the block as ``block``.
A run with no block (a confirmation run, say) has None as its block.
"""
missing = [f.name for f in factors if f.name not in frame.columns]
if missing:
raise ValueError(f"The runs have no column for the factor(s) {missing}.")
coded: dict[str, pd.Series] = {}
for i, f in enumerate(factors):
column = frame[f.name]
if f.type == FactorType.categorical:
levels = list(f.levels or [])
unknown = sorted({str(v) for v in column if v not in levels})
if unknown:
raise ValueError(f"Factor {f.name!r} has levels {unknown} in the runs, not among {levels}.")
coded[f"x{i}"] = column.map({levels[0]: -1.0, levels[1]: 1.0}).astype(float)
else:
mid, half = (f.high + f.low) / 2, (f.high - f.low) / 2 # type: ignore[operator]
coded[f"x{i}"] = (pd.to_numeric(column, errors="raise") - mid) / half
out = pd.DataFrame(coded, index=frame.index)
if response not in frame.columns:
raise ValueError(f"The runs have no response column {response!r}.")
out["y"] = pd.to_numeric(frame[response], errors="raise").astype(float)
if block is not None:
labels = frame[block] if block in frame.columns else pd.Series(None, index=frame.index, dtype=object)
out[_BLOCK_LABEL] = labels.map(lambda v: None if pd.isna(v) else str(v))
return out
#: The column holding each run's block, as text, and the prefix of the columns that code the blocks.
_BLOCK_LABEL = "block"
_BLOCK = "blk"
def _with_blocks(data: pd.DataFrame) -> pd.DataFrame:
"""Return the runs with the blocks coded as contrasts that sum to zero over the blocks, for every fit to hold.
The blocks are taken in the order they first appear in the runs, the latest last. Block ``j``
of ``m`` is 1 in contrast ``j``, and the latest block is -1 in every contrast: the process as
it was when the last runs were made, at which the confirmation runs are judged (see
:func:`_points`). One block needs no contrast.
"""
data = data.drop(columns=_nuisance(data))
if _BLOCK_LABEL not in data.columns:
return data
levels = list(dict.fromkeys(v for v in data[_BLOCK_LABEL] if v is not None))
for j, level in enumerate(levels[:-1]):
data[f"{_BLOCK}{j + 1}"] = [1.0 if v == level else -1.0 if v == levels[-1] else 0.0 for v in data[_BLOCK_LABEL]]
return data
def _nuisance(data: pd.DataFrame) -> list[str]:
"""Return the block contrasts among the runs' columns: every fit holds them, and no decision reads them."""
return [str(c) for c in data.columns if str(c).startswith(_BLOCK)]
# --------------------------------------------------------------------------- what the design can show
@dataclass(frozen=True)
class _Design:
"""What the runs, as a design, can and cannot show."""
n_runs: int
center: np.ndarray # runs with every numeric factor at its middle (and at least one numeric factor)
two_level: bool # every run that is not a centre run has every factor at an end
three_level: bool # not two-level, and every numeric factor at three levels or more: squares can be fitted
pure_error_df: int
def _design(study: _Study) -> _Design:
x = study.data[study.names].to_numpy(dtype=float)
numeric = np.array(study.numeric)
center = np.all(np.abs(x[:, numeric]) < _AT_EDGE, axis=1) if numeric.any() else np.zeros(len(x), dtype=bool)
away = x[~center]
two_level = bool(len(away)) and bool(np.all(np.abs(np.abs(away) - 1) < _AT_EDGE))
three_level = (
not two_level
and bool(numeric.any())
and all(len(np.unique(np.round(x[:, j], 9))) >= 3 for j in np.flatnonzero(numeric))
)
settings = pd.DataFrame(np.round(x, 9))
if _BLOCK_LABEL in study.data.columns: # repeats in different blocks differ by the blocks too
settings["block"] = study.data[_BLOCK_LABEL].fillna("").to_numpy()
groups = settings.value_counts()
return _Design(
n_runs=len(x),
center=center,
two_level=two_level,
three_level=three_level,
pure_error_df=int((groups - 1).sum()),
)
# --------------------------------------------------------------------------- the model
#: The column marking the centre runs of a two-level design, so the curvature they show is fitted
#: rather than left in the error, where it would hide the effects.
_CURVATURE = "curvature"
#: The fewest contrasts with which Lenth's method judges a saturated two-level design.
_LENTH_MIN_CONTRASTS = 7
@dataclass
class _Model:
"""The model the decision reads, and how its terms were judged.
``judged`` holds every term the runs could test, with its p-value; Lenth's method gives
none, and records 0 for an effect beyond its margin of error and 1 for one within it.
``signal`` says whether anything stands out from the noise once the number of terms
judged is allowed for. ``terms`` is the model for prediction: the active terms and the
main effects their interactions and squares are built on. ``fit`` is that model's fit,
with the curvature column when the design has centre runs among two-level runs.
``misses`` is the lack-of-fit test of the model the terms were judged in, when that model
did not fit the repeated runs; empty otherwise.
"""
terms: list[str]
fit: RegressionResultsWrapper
method: str
judged: dict[str, float]
signal: bool
noise_sd: float
noise_df: int
misses: dict[str, Any] = field(default_factory=dict)
b0: float = 0.0
b: np.ndarray = field(default_factory=lambda: np.zeros(0))
B: np.ndarray = field(default_factory=lambda: np.zeros((0, 0)))
@property
def error_df(self) -> int:
return self.noise_df
@property
def sigma(self) -> float:
return self.noise_sd
def p_value(self, term: str) -> float:
return float(self.judged.get(term, math.nan))
def active(self, alpha: float) -> list[str]:
"""Return the judged terms, other than the curvature column, that are significant; none without a signal."""
if not self.signal:
return []
return [t for t, p in self.judged.items() if t != _CURVATURE and p < alpha]
def value(self, x: np.ndarray) -> np.ndarray:
x = np.atleast_2d(x)
return self.b0 + x @ self.b + np.einsum("ij,jk,ik->i", x, self.B, x)
def gradient(self, x: np.ndarray) -> np.ndarray:
return self.b + 2 * self.B @ x
def _term_factors(term: str) -> set[str]:
return set(_canonical(term)[1]) if term != _CURVATURE else set()
def _with_parents(terms: list[str], names: list[str]) -> list[str]:
"""Return the terms after the main effects of every factor their interactions and squares use (heredity)."""
used = set().union(*(_term_factors(t) for t in terms)) if terms else set()
mains = [n for n in names if n in used]
return list(dict.fromkeys([*mains, *terms]))
def _rhs(names: list[str], numeric: list[bool], model: str) -> str:
"""Return the right-hand side of a standard model in the internal names; squares only of numeric factors."""
return " + ".join(_rhs_terms(names, numeric, model))
def _rhs_terms(names: list[str], numeric: list[bool], model: str) -> list[str]:
"""Return the terms of a standard model, one per column, in the internal names."""
terms = list(names)
if model != "main_effects":
terms += [f"{a}:{b}" for a, b in itertools.combinations(names, 2)]
if model == "quadratic":
terms += [f"I({n} ** 2)" for n, num in zip(names, numeric, strict=True) if num]
return terms
def _rhs_of(data: pd.DataFrame, terms: list[str]) -> str:
"""Return the right-hand side of a fit: the block contrasts first, so they lead their chains, then the terms."""
every = [*[c for c in _nuisance(data) if c not in terms], *terms]
return " + ".join(every) if every else "1"
def _ols(data: pd.DataFrame, terms: list[str]) -> RegressionResultsWrapper:
return smf.ols("y ~ " + _rhs_of(data, terms), data=data).fit()
def _matrix(data: pd.DataFrame, terms: list[str]) -> tuple[np.ndarray, list[str]]:
"""Return the model matrix of the terms (and the block contrasts), the intercept first, and its column names.
The matrix is read-only: it is built once for each set of runs and kept (see :func:`_built`).
"""
settings = data.drop(columns="y", errors="ignore")
key = (tuple(map(str, settings.columns)), pd.util.hash_pandas_object(settings, index=False).to_numpy().tobytes())
matrix, names = _built(_rhs_of(data, terms), _Settings(key, settings))
return matrix, list(names)
@dataclass(frozen=True)
class _Settings:
"""The runs' settings, without the response, compared and hashed by their columns and values alone."""
key: tuple[tuple[str, ...], bytes]
frame: pd.DataFrame = field(compare=False, hash=False)
@functools.lru_cache(maxsize=64)
def _built(rhs: str, settings: _Settings) -> tuple[np.ndarray, tuple[str, ...]]:
"""Build a model matrix once for each set of runs.
Judging one study builds the same matrices several times, and so do the simulated studies
of one design (see :mod:`process_improve.experiments.simulate`): only the response differs.
"""
matrix = dmatrix(rhs, settings.frame, return_type="dataframe")
values = matrix.to_numpy(dtype=float)
values.flags.writeable = False
return values, tuple(str(c) for c in matrix.columns)
def _chain_leaders(data: pd.DataFrame, rhs: str) -> list[str] | None:
"""Return one term per chain of exact aliases in a model; None when even those cannot all be estimated.
A block contrast leads its chain, so a term confounded with the blocks is not among the terms returned.
"""
exog, names = _matrix(data, [rhs])
leaders = [chain[0][0] for chain in alias_chains_of(exog, names)]
kept = exog[:, [names.index(t) for t in leaders]]
if int(np.linalg.matrix_rank(kept)) < kept.shape[1]:
return None
nuisance = _nuisance(data)
return [t for t in leaders if t != "Intercept" and t not in nuisance]
def _t_tests(data: pd.DataFrame, terms: list[str], extra: list[str], method: str, alpha: float) -> dict[str, Any]:
fit = _ols(data, [*terms, *extra])
judged = {t: float(fit.pvalues[t]) for t in [*terms, *extra]}
return {
"method": method,
"judged": judged,
"signal": _signal(fit, [*terms, *extra], alpha),
"considered": set(terms),
"noise_sd": float(np.sqrt(fit.scale)),
"noise_df": int(fit.df_resid),
"fit": fit,
}
def _against_pure_error(study: _Study, judged: dict[str, Any], alpha: float) -> dict[str, Any]:
"""Judge the terms against the repeats' own variation when the judging model does not fit the repeats.
Effects the model leaves out land in its residual, which then overstates the noise: every term
looks smaller than it is, and the model keeps too few terms to find what it misses. The
lack-of-fit test against the repeated runs shows it, and the terms are then tested against the
repeats' variation (pure error) instead. Only when something stands out (else there are no
terms to judge), and only with :data:`_FEW_ERROR_DF` degrees of freedom of pure error: the
test finds a lack of fit most often when the repeats happen to agree closely, and with fewer
such a chance spread would make many terms look active.
"""
fit = judged.get("fit")
if fit is None or not judged["signal"]: # nothing stands out: no terms to judge, whatever the noise
return judged
test = _run_lack_of_fit(fit, study.data, "y", _repeats(study), alpha)["lack_of_fit"]
if "error" in test or not test["significant"] or test["df_pure_error"] < _FEW_ERROR_DF:
return judged
t = fit.params / fit.bse * math.sqrt(fit.scale / test["ms_pure_error"])
p = pd.Series(2 * stats.t.sf(np.abs(t), test["df_pure_error"]), index=t.index)
return {
**judged,
"method": (
f"{judged['method']}, against the variation of the repeated runs: the model does not fit them "
f"(lack of fit p = {_g(test['p_value'])})"
),
"judged": {term: float(p[term]) for term in judged["judged"]},
"noise_sd": float(np.sqrt(test["ms_pure_error"])),
"noise_df": int(test["df_pure_error"]),
"misses": {"tested": True, "significant": True, "p_value": float(test["p_value"])},
}
#: Draws of pure noise behind the critical values of the tests for a signal, and the p-values of Lenth's method.
_NULL_DRAWS = 4_000
_LENTH_DRAWS = 20_000
def _signal(fit: RegressionResultsWrapper, tested: list[str], alpha: float) -> bool:
"""Say whether any of the ``tested`` terms of a fit stands out from the noise (see :func:`_stands_out`)."""
x = np.asarray(fit.model.exog, dtype=float)
columns = list(fit.model.exog_names)
return _stands_out(x, np.asarray(fit.model.endog, dtype=float), [columns.index(t) for t in tested], alpha)
def _stands_out(
x: np.ndarray, y: np.ndarray, tested: list[int], alpha: float, together: list[int] | None = None
) -> bool:
"""Say whether anything stands out from the noise in the regression of ``y`` on ``x``, at level ``alpha``.
Two tests see different signals: the smallest p-value of the ``tested`` columns times their
number (Bonferroni) sees one large effect, and the F-test of the columns ``together`` (by
default the tested ones), against the model without them, sees several moderate ones,
which Bonferroni misses with few degrees of freedom for error. The statistic is the smaller
of the two p-values; its critical value is its ``alpha`` quantile over responses of pure
noise on these runs, so the false alarm rate is ``alpha`` for this design, not up to twice
it. Columns in neither (the intercept, the block contrasts) are fitted and never tested.
"""
n, p = x.shape
df = n - p
together = tested if together is None else together
if not tested or df < 1:
return False
pinv = np.linalg.pinv(x)
scale = np.diag(pinv @ pinv.T)[tested]
rest = x[:, [j for j in range(p) if j not in together]]
rest_pinv = np.linalg.pinv(rest)
def statistic(responses: np.ndarray) -> np.ndarray:
beta = pinv @ responses
rss = ((responses - x @ beta) ** 2).sum(axis=0)
s2 = rss / df
t = np.abs(beta[tested]) / np.sqrt(np.outer(scale, s2))
bonferroni = 2 * stats.t.sf(t.max(axis=0), df) * len(tested)
without = ((responses - rest @ (rest_pinv @ responses)) ** 2).sum(axis=0)
f = (without - rss) / len(together) / s2
return np.minimum(bonferroni, stats.f.sf(f, len(together), df))
null = statistic(np.random.default_rng(n).standard_normal((n, _NULL_DRAWS)))
return bool(statistic(y[:, None])[0] < np.quantile(null, alpha))
def _pse(contrasts: np.ndarray) -> np.ndarray:
"""Return Lenth's pseudo standard error of each row of absolute contrasts."""
s0 = 1.5 * np.median(contrasts, axis=-1, keepdims=True)
return 1.5 * np.nanmedian(np.where(contrasts < 2.5 * s0, contrasts, np.nan), axis=-1)
@functools.lru_cache(maxsize=16)
def _lenth_null(m: int) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Return Lenth's ratios ``|contrast| / PSE`` over ``m`` contrasts of pure noise, each sample sorted.
The samples are one contrast's ratio, the largest, the second largest, and the smaller
of the simulated p-values of the largest and the second largest. Lenth's own margins, a t
quantile on ``m / 3`` degrees of freedom, are approximate (with seven contrasts the
individual one is 3.76 where the exact value is 2.31), and the simultaneous one is so wide
that two clear effects in an 8-run design rarely pass it. Simulated quantiles, as Ye and
Hamada (2000) tabulate them, are exact for any ``m``. The evidence that something is
active comes from the largest ratio, which sees one large effect, or the second largest,
which sees two moderate ones, with the critical value set so that the false alarm rate of
the pair is the significance level.
"""
z = np.abs(np.random.default_rng(m).standard_normal((_LENTH_DRAWS, m)))
ratios = z / _pse(z)[:, None]
ordered = np.sort(ratios, axis=1)
largest, second = np.sort(ordered[:, -1]), np.sort(ordered[:, -2])
both = np.minimum(_beyond(largest, ordered[:, -1]), _beyond(second, ordered[:, -2]))
return np.sort(ratios[:, 0]), largest, second, np.sort(both) # the contrasts are exchangeable
def _beyond(null: np.ndarray, values: np.ndarray) -> np.ndarray:
"""Return the share of a sorted null sample at or beyond each value: its simulated p-value."""
return 1 - np.searchsorted(null, values) / len(null)
def _lenth(data: pd.DataFrame, terms: list[str], alpha: float) -> dict[str, Any] | None:
"""Judge the terms of a saturated two-level design by Lenth's method, on every contrast the runs hold.
The contrasts are the terms' coefficients and, to make up ``n - 1``, the projections of the
response on an orthonormal basis of what the model leaves (scaled to the coefficients'
variance): a Plackett-Burman design's unused columns, or a full factorial's three-factor
interaction. The block contrasts are fitted and are no contrasts of their own. None when the
terms are not balanced, orthogonal to each other and to the blocks, or the contrasts are too few.
"""
x, every = _matrix(data, terms)
nuisance = _nuisance(data)
effect = [j for j, name in enumerate(every) if j and name not in nuisance] # patsy orders terms by degree
names = [every[j] for j in effect]
n, columns = len(x), x[:, effect]
blocks = x[:, [j for j, name in enumerate(every) if name in nuisance]]
if (
not np.allclose(columns.T @ columns, n * np.eye(columns.shape[1]))
or not np.allclose(columns.sum(axis=0), 0)
or not np.allclose(blocks.T @ columns, 0)
):
return None
y = data["y"].to_numpy(dtype=float)
effects = columns.T @ y / n
rest = np.linalg.qr(x, mode="complete")[0][:, x.shape[1] :].T @ y / np.sqrt(n)
m = len(effects) + len(rest)
pse = float(_pse(np.abs(np.concatenate([effects, rest]))))
if m < _LENTH_MIN_CONTRASTS or not pse > 0:
return None
one, largest, second, both = _lenth_null(m)
ratios = np.abs(effects) / pse
p_values = _beyond(one, ratios)
top = np.sort(ratios)[::-1] # the terms' ratios only: a signal must be in a term the decision can act on
observed = min(_beyond(largest, top[:1])[0], _beyond(second, top[1:2])[0] if len(top) > 1 else 1.0)
signal = bool(observed < np.quantile(both, alpha))
return {
"method": f"Lenth's method on the {m} contrasts of the saturated two-level design, with simulated p-values",
"judged": dict(zip(names, (float(p) for p in p_values), strict=True)),
"signal": signal,
"considered": set(terms),
"noise_sd": pse,
"noise_df": max(1, round(m / 3)),
"contrasts": m,
}
def _lenth_on_the_cube(
study: _Study, design: _Design, terms: list[str], alpha: float, t_df: int
) -> dict[str, Any] | None:
"""Judge a two-level design with centre runs by Lenth's method on its two-level runs, when t-tests see less.
Centre runs added to a saturated two-level design leave a t-test few degrees of freedom for
error: two, with three centre runs in a 16-run fraction, where a t of 4.3 is needed to pass.
The pseudo standard error of the two-level runs' m contrasts is worth about m / 3 degrees of
freedom (Lenth, 1989), so it judges the terms instead when that is more, and the centre runs
test curvature against their own repeats. Each of the two tests is made at ``alpha / 2`` for
the signal, so together they raise a false alarm at most ``alpha`` of the time. None when
t-tests have :data:`_FEW_ERROR_DF` or more, or Lenth's method does not apply to the two-level runs.
"""
if t_df >= _FEW_ERROR_DF or not (design.two_level and design.center.any()):
return None
judged = _lenth(study.data.loc[~design.center], terms, alpha / 2)
if judged is None or judged["noise_df"] <= t_df:
return None
test = _centre_runs_test(study, design, alpha)
p = math.nan if "error" in test else float(test["p_value"])
judged["judged"][_CURVATURE] = p
judged["signal"] = judged["signal"] or bool(p < alpha / 2)
judged["method"] = (
f"Lenth's method on the {judged['contrasts']} contrasts of the two-level runs, with simulated p-values"
)
if "error" not in test:
judged["method"] += ", and a curvature test of the centre runs against their repeats"
return judged
def _omars(data: pd.DataFrame, names: list[str], alpha: float) -> dict[str, Any] | None:
from process_improve.experiments.omars import _has_omars_structure, analyze_omars # noqa: PLC0415 - one caller
if _nuisance(data):
return None # the staged analysis has no place for blocks
if not _has_omars_structure(data[names].to_numpy(dtype=float)):
return None # its tests hold only for that structure; checked first, as its search can take long elsewhere
with warnings.catch_warnings():
warnings.simplefilter("ignore") # a structure warning cannot arise after the check above
result = analyze_omars(
data[names],
data["y"],
alpha_main=alpha,
alpha_second_overall=alpha,
alpha_second_subset=alpha,
interaction_heredity="strong", # a DSD's A:B is partly aliased with C:D and E:F; heredity picks A:B
)
if not result.success:
return None
judged = {name: float(p) for name, p in result.main_effect_p_values.items()}
judged.update(dict.fromkeys(result.active_interactions, 0.0))
judged.update({f"I({q[:-2]} ** 2)": 0.0 for q in result.active_quadratics})
return {
"method": "the staged analysis of a definitive screening design (analyze_omars)",
"judged": judged,
"signal": _omars_signal(data, names, alpha),
"considered": set(_rhs_terms(names, [True] * len(names), "quadratic")),
"noise_sd": float(result.updated_rmse),
"noise_df": int(result.updated_error_df),
}
def _omars_signal(data: pd.DataFrame, names: list[str], alpha: float) -> bool:
"""Say whether a definitive screening design shows anything: its main effects, or its second-order terms.
The main effects are tested by Bonferroni, and the second-order terms together by the
F-test, against the error of the full second-order model, as the staged analysis does;
the pair's critical value is simulated (:func:`_stands_out`). The staged analysis's own
second-order test pools the main effects it found inactive into the error, which makes
that error too small and the test fire on noise.
"""
mains = np.column_stack([np.ones(len(data)), data[names].to_numpy(dtype=float)])
second = _matrix(data, _rhs_terms(names, [True] * len(names), "quadratic"))[0]
second = second - mains @ np.linalg.lstsq(mains, second, rcond=None)[0]
u, sv, _ = np.linalg.svd(second, full_matrices=False)
basis = u[:, sv > 1e-8 * sv.max()]
x = np.hstack([mains, basis])
y = data["y"].to_numpy(dtype=float)
return _stands_out(x, y, list(range(1, len(names) + 1)), alpha, together=list(range(1, x.shape[1])))
def _judge_terms(study: _Study, design: _Design, alpha: float) -> dict[str, Any] | None:
"""Judge every term the runs can test, in the largest standard model they estimate.
Testing the terms a stepwise search picked, on the data that picked them, makes them look
stronger than they are: in simulated noise it found an effect in most studies. So each
term is tested in the largest model the runs estimate with at least two degrees of freedom
for error (one column per chain of exact aliases): the full quadratic, the interactions, or
the main effects. A two-level design too small for that is judged by Lenth's method on all
its contrasts. A design whose second-order terms are partly aliased with each other, which
no standard model separates, is judged by its own staged analysis when it is a definitive
screening design; with centre runs among two-level runs (a Plackett-Burman design with
centre runs), a stepwise search (AICc) proposes its interactions, which are then tested
beside every main effect. Testing the main effects alone would leave an active interaction
in the error, where it hides the effects. Returns None when the runs are too few to judge
anything.
"""
data, names = study.data, study.names
extra = _extra(data, design, study.numeric)
error_runs = len(data) - len(_nuisance(data)) - 1 - len(extra) # each block contrast takes a run
saturated = design.two_level and not extra # judged by Lenth's method when no model leaves error
judged, aliased = _judge_in_standard_models(study, design, alpha, extra, error_runs)
if not aliased:
return judged
judged = _judge_proposed_terms(study, design, alpha, extra, error_runs)
if judged is None and saturated: # the main effects, and the columns the design leaves for error
judged = _lenth(data, names, alpha)
if judged is not None:
return judged
mains = _chain_leaders(data, _rhs(names, study.numeric, "main_effects"))
if mains is not None and error_runs - len(mains) >= _MIN_ERROR_DF:
return _t_tests(data, mains, extra, "t-tests in the main effects model", alpha)
return None
def _judge_in_standard_models(
study: _Study, design: _Design, alpha: float, extra: list[str], error_runs: int
) -> tuple[dict[str, Any] | None, bool]:
"""Judge the terms in the largest standard model the runs estimate, or by Lenth's method on all the contrasts.
Returns the judgement, None when even the main effects leave too few runs, and whether the
terms are partly aliased with each other, which no standard model separates (the judgement
is then None, and the caller judges them another way). ``error_runs`` is the number of runs
left for error once the intercept, the block contrasts and ``extra`` are fitted.
"""
data, names, numeric = study.data, study.names, study.numeric
models = (["quadratic"] if design.three_level else []) + ["interactions", "main_effects"]
for model in models:
leaders = _chain_leaders(data, _rhs(names, numeric, model))
if leaders is None:
return None, True
if error_runs - len(leaders) >= _MIN_ERROR_DF:
method = f"t-tests in the {model.replace('_', ' ')} model"
judged = _lenth_on_the_cube(study, design, leaders, alpha, error_runs - len(leaders))
return judged or _t_tests(data, leaders, extra, method, alpha), False
if design.two_level and not extra and model == "interactions":
judged = _lenth(data, leaders, alpha)
if judged is not None:
return judged, False
return None, False
def _judge_proposed_terms(
study: _Study, design: _Design, alpha: float, extra: list[str], error_runs: int
) -> dict[str, Any] | None:
"""Judge partly aliased terms: by the staged analysis of a definitive screening design, or the AICc proposal.
The search proposes second-order terms, which are then tested beside every main effect.
None when neither applies, or the proposal leaves too few runs for error.
"""
data, names = study.data, study.names
if design.three_level:
judged = _omars(data, names, alpha)
if judged is not None:
return judged
searched = "quadratic" if design.three_level else "interactions"
proposed = _proposed_terms(_less_blocks(data), "y", names, searched)
second = _estimable(data, names, [t for t in proposed if t not in names], extra)
if error_runs - len(names) - len(second) < _MIN_ERROR_DF:
return None
method = f"t-tests of every main effect and the {searched} terms a stepwise search (AICc) proposed"
judged = _t_tests(data, [*names, *second], extra, method, alpha)
judged["considered"] = set(_rhs_terms(names, study.numeric, searched)) # the search weighed every term
return judged
def _estimable(data: pd.DataFrame, names: list[str], second: list[str], extra: list[str]) -> list[str]:
"""Return the proposed second-order terms, in order, that keep the model with every main effect full rank.
The search proposed them beside the main effects it kept; beside every main effect, a
Plackett-Burman design's partial aliasing can leave one of them a combination of the others.
"""
kept: list[str] = []
for term in second:
x = _matrix(data, [*names, *kept, term, *extra])[0]
if int(np.linalg.matrix_rank(x)) == x.shape[1]:
kept.append(term)
return kept
def _extra(data: pd.DataFrame, design: _Design, numeric: list[bool]) -> list[str]:
"""Return the curvature column when the runs are two-level with centre runs, and the blocks do not hold it.
Centre runs made in a block of their own would show the shift between the blocks as curvature: their
column is then a combination of the block contrasts, and is left out.
"""
if not (design.two_level and design.center.any() and any(numeric)):
return []
blocks = _matrix(data, [])[0] # the intercept and the block contrasts
return [_CURVATURE] if _explained(blocks, data[_CURVATURE].to_numpy(dtype=float)) < 1 - 1e-9 else []
def _less_blocks(data: pd.DataFrame) -> pd.DataFrame:
"""Return the runs with the shifts between the blocks taken out of the response, for a search that has no blocks."""
if not _nuisance(data):
return data
blocks = _matrix(data, [])[0]
y = data["y"].to_numpy(dtype=float)
return data.assign(y=y - blocks[:, 1:] @ np.linalg.lstsq(blocks, y, rcond=None)[0][1:])
def _interaction_candidates(study: _Study, active: list[str], *, weak: bool = False) -> list[str]:
"""Return the two-factor interactions heredity makes plausible, not yet active.
Strong heredity needs both factors in active terms; ``weak`` heredity needs one.
"""
factors = set().union(*(_term_factors(t) for t in active)) if active else set()
present = {frozenset(_term_factors(t)) for t in active if ":" in t}
return [
f"{a}:{b}"
for a, b in itertools.combinations(study.names, 2)
if ({a, b} & factors if weak else {a, b} <= factors) and frozenset((a, b)) not in present
]
def _separable(study: _Study, terms: list[str], candidates: list[str]) -> list[str]:
"""Return the candidate interactions the runs can tell apart from the model's terms and from each other.
Each is kept when the model's columns, with those of the candidates kept before it, do not
reproduce its column (see :func:`_aliasing`): in a resolution IV fraction B:D and C:G can
each be told from the main effects, but not from each other.
"""
x = _matrix(study.data, terms)[0]
kept = []
for candidate in candidates:
column = _product(study, candidate)
if _explained(x, column) < _ALIASED:
kept.append(candidate)
x = np.column_stack([x, column])
return kept
def _product(study: _Study, interaction: str) -> np.ndarray:
a, b = interaction.split(":")
return (study.data[a] * study.data[b]).to_numpy(dtype=float)
def _explained(x: np.ndarray, column: np.ndarray) -> float:
"""Return the share of a column's variation about its mean that the columns of ``x`` reproduce (R squared)."""
centred = column - column.mean()
total = float(centred @ centred)
if total <= 0:
return 1.0 # a constant column: the intercept holds it
residual = column - x @ np.linalg.lstsq(x, column, rcond=None)[0]
return 1.0 - float(residual @ residual) / total
def _second_look(study: _Study, judged: dict[str, Any], active: list[str], extra: list[str], alpha: float) -> None:
"""Test the plausible interactions of the active factors that the judging never tested, and add the active ones.
A Plackett-Burman design judged on its main effects never tests A:C, though A and C are
active and A:C is nearly orthogonal to every main effect. In the model of the active terms
it has the degrees of freedom to be tested, and strong heredity keeps the family small.
When the model the terms were judged in misses the repeated runs (``judged["misses"]``), the
effects it leaves out are likeliest to be interactions of an active factor with one judged
inactive: such an effect changes sign across the active factor's range and averages out of
the other factor's main effect. Those are tested too (weak heredity), even when a stepwise
search weighed them, since it weighed them against a residual that holds what the model
misses; Holm's method keeps the larger family's false alarms at ``alpha``. Once one is found,
the model holds a term the judging model lacked, and its own lack-of-fit test decides.
"""
misses = bool(judged.get("misses"))
terms = _with_parents(active, study.names)
tested = judged["judged"] if misses else judged["considered"]
untested = [c for c in _interaction_candidates(study, active, weak=misses) if c not in tested]
untested = _separable(study, [*terms, *extra], untested)
spare = len(study.data) - len(_nuisance(study.data)) - len(terms) - len(untested) - 1 - len(extra)
if not untested or spare < _MIN_ERROR_DF:
return
fit = _ols(study.data, [*terms, *untested, *extra])
p_values = [float(fit.pvalues[c]) for c in untested]
if misses:
p_values = [float(p) for p in multipletests(p_values, method="holm")[1]]
found = dict(zip(untested, p_values, strict=True))
judged["judged"].update(found)
judged["considered"] |= set(untested)
judged["method"] += (
f"; then the {len(untested)} interaction(s) of active factors it could not test, by t-tests in the model "
"of the active terms" + (" (every factor with an active one, Holm's method)" if misses else "")
)
new = [c for c, p in found.items() if p < alpha]
active.extend(new)
if new:
judged.pop("misses", None)
def _fit(study: _Study, terms: list[str], judged: dict[str, Any], design: _Design) -> _Model:
"""Fit the model for prediction, and write it as the polynomial ``b0 + x'b + x'Bx`` in coded units.
With blocks, ``b0`` is the latest block's level (every contrast -1), where the next runs are made.
"""
extra = _extra(study.data, design, study.numeric)
fit = _ols(study.data, [*terms, *extra])
nuisance = _nuisance(study.data)
params = [
{"term": str(name), "coefficient": float(value)}
for name, value in fit.params.items()
if str(name) != _CURVATURE and str(name) not in nuisance
]
b0, b, B = _extract_b_and_B(params, study.names)
b0 -= sum(float(fit.params[c]) for c in nuisance) # the latest block is -1 in every contrast
return _Model(
terms=terms,
fit=fit,
method=judged["method"],
judged=judged["judged"],
signal=judged["signal"],
noise_sd=judged["noise_sd"],
noise_df=judged["noise_df"],
misses=judged.get("misses", {}),
b0=b0,
b=b,
B=B,
)
def _select(
study: _Study, design: _Design, model_terms: list[str] | None, alpha: float, *, against_repeats: bool = True
) -> _Model:
"""Judge the terms (see :func:`_judge_terms`), or the caller's model's, and fit the model for prediction.
``against_repeats`` allows the terms to be judged against the repeated runs when the model does not fit them
(:func:`_against_pure_error`).
"""
extra = _extra(study.data, design, study.numeric)
if model_terms is not None:
terms = _internal_terms(study, model_terms)
return _fit(study, terms, _t_tests(study.data, terms, extra, "t-tests in the model given", alpha), design)
judged = _judge_terms(study, design, alpha)
if judged is None:
empty = {
"method": "none: too few runs to judge any term",
"judged": {},
"signal": False,
"noise_sd": math.nan,
"noise_df": 0,
}
return _fit(study, [], empty, design)
if against_repeats:
judged = _against_pure_error(study, judged, alpha)
active = [t for t, p in judged["judged"].items() if t != _CURVATURE and p < alpha] if judged["signal"] else []
if active:
_second_look(study, judged, active, extra, alpha)
return _fit(study, _with_parents(active, study.names), judged, design)
def _internal_terms(study: _Study, model_terms: list[str]) -> list[str]:
"""Return the caller's model terms in the internal names; the intercept is implied."""
internal = {f.name: f"x{i}" for i, f in enumerate(study.factors)}
terms = []
for term in model_terms:
canonical, names = _canonical(term)
unknown = [n for n in names if n not in internal]
if unknown:
raise ValueError(f"The model term {term!r} names {unknown}, which are not factors.")
if canonical == "Intercept":
continue
parts = [internal[n] for n in names]
if "**" in canonical:
terms.append(f"I({parts[0]} ** 2)")
else:
terms.append(":".join(parts))
return terms
# --------------------------------------------------------------------------- the evidence
def _favoured_ends(study: _Study, sense: float) -> np.ndarray:
"""Each factor's end that its marginal effect favours for the goal (+1 or -1), from the runs at its ends."""
ends = np.ones(len(study.factors))
y = _less_blocks(study.data)["y"].to_numpy(dtype=float)
for j, name in enumerate(study.names):
x = study.data[name].to_numpy(dtype=float)
high, low = y[x > _AT_EDGE], y[x < -_AT_EDGE]
contrast = (high.mean() - low.mean()) if len(high) and len(low) else 0.0
ends[j] = 1.0 if contrast * (sense or 1.0) >= 0 else -1.0
return ends
def _sense(goal: Goal, target: float | None, predicted: float) -> float:
"""Which way the response should move: +1 up, -1 down, 0 for a target already met."""
if goal == "maximize":
return 1.0
if goal == "minimize":
return -1.0
assert target is not None
return float(np.sign(target - predicted))
def _optimum(study: _Study, model: _Model, goal: Goal, target: float | None, seed: int) -> np.ndarray:
"""Return the coded settings inside the ranges at which the model best meets the goal.
A polynomial of at most second order is searched from the centre, the corner its linear
terms point to, every corner up to six factors, and eight random starts. A categorical
factor is searched at each of its two levels; a factor the model leaves out is set at
the end it favours (:func:`_favoured_ends`). Many settings can meet a target; the search
then prefers the one nearest the middle of the ranges.
"""
in_model = set().union(*(_term_factors(t) for t in model.terms)) if model.terms else set()
free = [j for j, name in enumerate(study.names) if name in in_model and study.numeric[j]]
categorical = [j for j, name in enumerate(study.names) if name in in_model and not study.numeric[j]]
nearest = 1e-6 * float(study.data["y"].var()) # the weight that breaks ties among settings that meet a target
def loss(x: np.ndarray) -> tuple[float, np.ndarray]:
"""Return the loss at ``x`` and its gradient, exact for a polynomial of second order."""
value, slope = float(model.value(x)[0]), model.gradient(x)
if goal == "target":
miss = value - float(target) # type: ignore[arg-type]
return miss**2 + nearest * float(x @ x), 2 * miss * slope + 2 * nearest * x
return (-value, -slope) if goal == "maximize" else (value, slope)
rng = check_random_state(seed)
base = _favoured_ends(study, _sense(goal, target, float(model.b0)))
best_x, best_loss = base.copy(), math.inf
for levels in itertools.product((-1.0, 1.0), repeat=len(categorical)):
start = base.copy()
start[categorical] = levels
for s in _starts(model, free, goal, rng) if free else [np.zeros(0)]:
def partial(z: np.ndarray, start: np.ndarray = start) -> tuple[float, np.ndarray]:
x = start.copy()
x[free] = z
value, slope = loss(x)
return value, slope[free]
if free:
found = optimize.minimize(partial, s, jac=True, bounds=[(-1.0, 1.0)] * len(free), method="L-BFGS-B")
z, value = np.clip(found.x, -1.0, 1.0), float(found.fun)
else: # the model holds no numeric factor: nothing to search
z, value = s, partial(s)[0]
if value < best_loss - 1e-12:
best_loss, best_x = value, start.copy()
best_x[free] = z
return best_x
def _starts(model: _Model, free: list[int], goal: Goal, rng: np.random.Generator) -> list[np.ndarray]:
"""Return the starting points of the search for the optimum over the ``free`` factors."""
starts = [np.zeros(len(free)), np.sign(model.b[free] * (-1 if goal == "minimize" else 1))]
if len(free) <= _ALL_CORNERS:
starts += [np.array(v, dtype=float) for v in itertools.product((-1.0, 1.0), repeat=len(free))]
return starts + list(rng.uniform(-1, 1, size=(_RANDOM_STARTS, len(free))))
def _curvature(study: _Study, design: _Design, model: _Model, alpha: float) -> dict[str, Any]:
"""Say whether the response curves, and how the runs show it: centre runs, or squares in the model."""
squares = [t for t in model.terms if "**" in t]
significant_squares = [t for t in squares if model.p_value(t) < alpha]
if design.three_level:
return {
"source": "squares",
"significant": bool(significant_squares),
"terms": significant_squares,
}
if design.two_level and design.center.any():
p = model.p_value(_CURVATURE)
if np.isfinite(p): # the curvature column's test in the judging fit, against its residual error
return {"source": "center_runs", "significant": bool(p < alpha), "p_value": p}
if design.pure_error_df > 0: # judged by Lenth's method: the centre runs' own repeats instead
test = _centre_runs_test(study, design, alpha)
if "error" not in test:
return {"source": "center_runs", "significant": bool(test["significant"]), "p_value": test["p_value"]}
return {"source": None, "significant": False}
def _centre_runs_test(study: _Study, design: _Design, alpha: float) -> dict[str, Any]:
"""Test the centre runs' mean against the two-level runs' mean, against the repeats (:func:`_run_curvature_test`).
The test ignores the blocks, so it is not made when the blocks hold the centre runs
(:func:`_blocks_hold_the_centre_runs`): the result then holds an ``error``, as when there are no repeats.
"""
if _blocks_hold_the_centre_runs(study, design):
return {"error": "The centre runs share no block with a two-level run: the blocks hide any curvature."}
return _run_curvature_test(study.data, "y", study.names, alpha, _repeats(study))["curvature_test"]
def _blocks_hold_the_centre_runs(study: _Study, design: _Design) -> bool:
"""Say whether the centre runs of a two-level design share no block with a two-level run: a round of their own.
They then differ from the two-level runs by the shift between the blocks too, and their column is a
combination of the block contrasts (see :func:`_extra`): no test can tell that shift from curvature.
"""
return design.two_level and bool(design.center.any()) and not _extra(study.data, design, study.numeric)
def _repeats(study: _Study) -> list[str]:
"""Return the columns whose values mark repeated runs: the settings, and the block when there are blocks."""
return [*study.names, _BLOCK_LABEL] if _nuisance(study.data) else study.names
def _lack_of_fit(study: _Study, model: _Model, alpha: float) -> dict[str, Any]:
"""Test the model against the repeated runs, or carry over the judging model's lack of fit.
The model the terms were judged in holds every term of this one; when it misses the repeats
(``model.misses``), so does this one, whatever its own test says. That test pools the effects
of the terms the model drops into its lack of fit, and so has less power.
"""
own: dict[str, Any] = {"tested": False}
if model.error_df > 0:
result = _run_lack_of_fit(model.fit, study.data, "y", _repeats(study), alpha)["lack_of_fit"]
if result.get("p_value") is not None: # no p-value: no replicates, no lack-of-fit df, or zero pure error
own = {"tested": True, "significant": bool(result["significant"]), "p_value": result["p_value"]}
return model.misses if model.misses and not own.get("significant") else own
def _aliasing(study: _Study, model: _Model, active: list[str]) -> list[dict[str, Any]]:
"""List each plausible interaction outside the model that the runs cannot tell apart from its active terms.
A two-factor interaction is plausible when both its factors are in active terms (strong
heredity): in a resolution III fraction, main effect C aliased with A:B is in doubt when A
and B are active, not when they are not. The runs cannot tell it apart from the model when
the model's columns reproduce its column (R squared of :data:`_ALIASED` or more): wholly
in a regular fraction, where it is aliased, and not at all in a Plackett-Burman or
definitive screening design, whose partial aliasing the analysis separates. The active terms
it is confused with are those its column loads on, with the weight that is each term's
bias per unit of the interaction's coefficient, ``(X1'X1)^-1 X1' x2``.
"""
candidates = _interaction_candidates(study, [t for t in model.terms if t in active or ":" in t])
if not active or not candidates:
return []
exog = np.asarray(model.fit.model.exog, dtype=float)
columns = list(model.fit.model.exog_names)
found = []
for candidate in candidates:
column = _product(study, candidate)
if _explained(exog, column) < _ALIASED:
continue
weights = np.linalg.pinv(exog) @ column
loads = [j for j, name in enumerate(columns) if name in active and abs(weights[j]) >= _ALIAS_WEIGHT]
if not loads: # confused with the model as a whole, not one term: name the term it loads on most
loads = [max((j for j, name in enumerate(columns) if name in active), key=lambda j: abs(weights[j]))]
found += [{"term": columns[j], "aliased_with": candidate, "weight": float(weights[j])} for j in loads]
return found
def _edges(study: _Study, model: _Model, x: np.ndarray, sense: float) -> tuple[list[dict[str, Any]], float]:
"""List the numeric factors at an end of their range at ``x``, and the gain per coded unit past them.
A factor improves past its end when the gradient there points out of the range, the way
the goal wants the response to move. The gain is the directional derivative along the
outward directions of those factors, scaled to one coded unit in the factor that moves
most: the change the first step of the path of steepest ascent should bring. A factor
whose range is a hard limit (``Factor.hard_limits``) cannot go past it, so it is no
part of that direction.
"""
gradient = model.gradient(x)
in_model = set().union(*(_term_factors(t) for t in model.terms)) if model.terms else set()
edges, outward = [], np.zeros(len(x))
for j, name in enumerate(study.names):
if not study.numeric[j] or name not in in_model or abs(x[j]) < 1 - _AT_EDGE:
continue
pushes = sense * gradient[j] * np.sign(x[j])
improves = bool(pushes > 1e-9)
limited = study.factors[j].hard_limits
edges.append(
{
"factor": study.factors[j].name,
"end": "high" if x[j] > 0 else "low",
"improves_beyond": improves,
"hard_limit": limited,
}
)
if improves and not limited:
outward[j] = sense * gradient[j]
if not outward.any():
return edges, 0.0
direction = outward / np.abs(outward).max()
return edges, float(abs(gradient @ direction))
# --------------------------------------------------------------------------- the runs of each action
def _points(study: _Study, coded: np.ndarray) -> pd.DataFrame:
"""Return coded settings as the model's fit reads them: the curvature column 1 at the centre, the latest block.
Runs made next are made after every run so far, so they are predicted as the process was in
the latest block (every block contrast -1), not at the average of the blocks.
"""
frame = pd.DataFrame(np.atleast_2d(coded), columns=study.names)
numeric = [n for n, num in zip(study.names, study.numeric, strict=True) if num]
frame[_CURVATURE] = (frame[numeric].abs() < _AT_EDGE).all(axis=1).astype(float) if numeric else 0.0
for contrast in _nuisance(study.data):
frame[contrast] = -1.0
return frame
def _rows(study: _Study, coded: np.ndarray) -> list[dict[str, Any]]:
return [study.actual(row) for row in np.atleast_2d(coded)]
def _two_level_part(study: _Study, design: _Design) -> pd.DataFrame:
return study.data.loc[~design.center, study.names]
def _with_others(base: np.ndarray, runs: pd.DataFrame, columns: list[int], study: _Study) -> np.ndarray:
"""Put augment_design's runs over ``columns`` into full coded rows, the other factors at ``base``."""
out = np.tile(base, (len(runs), 1))
for j in columns:
out[:, j] = runs[study.names[j]].to_numpy(dtype=float)
return out
# --------------------------------------------------------------------------- the context of one decision
@dataclass
class _Context:
"""The study, its model, the evidence and the caller's choices: what every action builder reads."""
study: _Study
design: _Design
model: _Model
goal: Goal
target: float | None
alpha: float
n_confirmation_runs: int
ranges_can_move: bool
seed: int
history: tuple[str, ...] = ()
optimum: np.ndarray = field(default_factory=lambda: np.zeros(0))
sense: float = 1.0
active: list[str] = field(default_factory=list)
active_factors: list[int] = field(default_factory=list)
curvature: dict[str, Any] = field(default_factory=dict)
lack_of_fit: dict[str, Any] = field(default_factory=dict)
aliasing: list[dict[str, Any]] = field(default_factory=list)
edges: list[dict[str, Any]] = field(default_factory=list)
gain_beyond: float = 0.0
@property
def k(self) -> int:
return len(self.study.factors)
def name(self, j: int) -> str:
return self.study.factors[j].name
def predicted(self) -> float:
return float(self.model.value(self.optimum)[0])
def mean_interval(self, n: int) -> tuple[float, float]:
"""Return the prediction interval at the optimum for the mean of ``n`` new runs."""
frame = self.model.fit.get_prediction(_points(self.study, self.optimum)).summary_frame(alpha=self.alpha).iloc[0]
if self.model.error_df <= 0:
return math.nan, math.nan
t = float(stats.t.ppf(1 - self.alpha / 2, self.model.error_df))
half = t * float(np.sqrt(frame["mean_se"] ** 2 + self.model.fit.scale / n))
return float(frame["mean"]) - half, float(frame["mean"]) + half
def _gather(ctx: _Context) -> _Context:
"""Fill in the evidence: the optimum, the active terms and factors, curvature, lack of fit, aliasing, edges."""
study, model = ctx.study, ctx.model
ctx.optimum = _optimum(study, model, ctx.goal, ctx.target, ctx.seed)
ctx.sense = _sense(ctx.goal, ctx.target, ctx.predicted()) or 1.0
ctx.active = model.active(ctx.alpha)
factors = set().union(*(_term_factors(t) for t in ctx.active)) if ctx.active else set()
ctx.active_factors = [j for j, name in enumerate(study.names) if name in factors]
ctx.curvature = _curvature(study, ctx.design, model, ctx.alpha)
ctx.lack_of_fit = _lack_of_fit(study, model, ctx.alpha)
ctx.aliasing = _aliasing(study, model, ctx.active)
ctx.edges, ctx.gain_beyond = _edges(study, model, ctx.optimum, ctx.sense)
return ctx
# --------------------------------------------------------------------------- the actions
@dataclass
class _Step:
"""One recommended action: what to do, why, the runs to make (coded), and how to judge them."""
action: str
headline: str
reasons: list[str]
runs: np.ndarray
after: str
extra: dict[str, Any] = field(default_factory=dict)
def _g(value: float) -> str:
"""Write a number to four significant figures."""
return f"{value:.4g}"
def _names(ctx: _Context, indices: list[int]) -> str:
return ", ".join(ctx.name(j) for j in indices)
def _confirm(ctx: _Context) -> _Step:
n = ctx.n_confirmation_runs
low, high = ctx.mean_interval(n)
left_out = [j for j in range(ctx.k) if j not in _model_factor_indices(ctx)]
reasons = [f"The model predicts {_g(ctx.predicted())} at these settings."]
if np.isfinite(low):
reasons.append(
f"The mean of {n} runs there should fall between {_g(low)} and {_g(high)} "
f"({100 * (1 - ctx.alpha):.0f}% prediction interval for a mean of {n} runs)."
)
if left_out:
reasons.append(
f"{_names(ctx, left_out)} {'is' if len(left_out) == 1 else 'are'} not in the model, so set at the end "
"the small, not significant, effect favours; any setting in the range should do as well."
)
reasons.extend(_held_at_limits(ctx))
after = (
f"If the mean of the {n} runs lies between {_g(low)} and {_g(high)}, accept these settings. If not, the "
"model has missed something: call recommend_next_step again with these runs as the confirmation runs."
)
return _Step(
"confirm",
f"Confirm the predicted optimum with {n} runs at these settings.",
reasons,
np.tile(ctx.optimum, (n, 1)),
after,
)
def _model_factor_indices(ctx: _Context) -> list[int]:
factors = set().union(*(_term_factors(t) for t in ctx.model.terms)) if ctx.model.terms else set()
return [j for j, name in enumerate(ctx.study.names) if name in factors]
def _add_center_points(ctx: _Context) -> _Step | None:
"""Three centre runs, for a two-level design of numeric factors that has none."""
if not all(ctx.study.numeric) or not ctx.design.two_level or ctx.design.center.any():
return None
return _Step(
"add_center_points",
f"Add {_CENTER_RUNS} runs at the centre of every range.",
[
(
"Every run so far sets each factor at an end of its range, so nothing shows whether the response "
"curves between the ends, and no run is repeated to estimate the noise."
),
],
np.zeros((_CENTER_RUNS, ctx.k)),
"Analyse again with the centre runs: if their mean is off the plane of the corners "
"(analysis_type='curvature_test'), the next round adds axial runs; their spread estimates the noise.",
)
def _replicate(ctx: _Context) -> _Step:
runs = ctx.study.data[ctx.study.names].to_numpy(dtype=float)
return _Step(
"replicate",
f"Make the {len(runs)} runs again.",
[
(
f"The model leaves {ctx.model.error_df} degree(s) of freedom to estimate the noise, too few to judge "
"which effects are real."
),
],
runs,
"Analyse every run together: the repeats estimate the noise, and every coefficient's variance halves.",
)
def _more_error_df(ctx: _Context) -> _Step:
"""Centre runs, or the runs again, when nothing stands out but the noise is estimated from too few runs."""
step = _add_center_points(ctx) or _replicate(ctx)
step.reasons = [
(
f"No effect stands out from the noise ({ctx.model.method}), but the noise is estimated with only "
f"{ctx.model.error_df} degree(s) of freedom: a real effect needs to be several times larger to stand out "
"than it would with more. More runs make the noise estimate firm before the ranges are questioned."
),
*step.reasons,
]
return step
def _widen_ranges(ctx: _Context) -> _Step:
"""Widen the numeric ranges about their middles and run the design again; replicate when none can move.
A range that is a hard limit (``Factor.hard_limits``) stays as it is.
"""
spread = float(ctx.study.data["y"].std(ddof=1)) if len(ctx.study.data) > 1 else math.nan
smallest = min(ctx.model.judged.values(), default=math.nan)
reasons = [
(
f"No effect stands out from the noise ({ctx.model.method}): with {len(ctx.model.judged)} terms judged, the "
f"smallest p-value, {_g(smallest)}, is what chance alone gives at {ctx.alpha:g} (the response's standard "
f"deviation over the runs is {_g(spread)})."
),
(
"Either the ranges are too narrow for the factors to move the response, or the measurement is too "
"noisy to see them move it: check its repeatability before spending more runs."
),
]
soft = [j for j, f in enumerate(ctx.study.factors) if f.type != FactorType.categorical and not f.hard_limits]
if not ctx.ranges_can_move or not soft:
step = _replicate(ctx)
why = "cannot move" if not ctx.ranges_can_move else "are hard limits"
step.reasons = [*reasons, f"The ranges {why}, so more runs at the same settings average the noise down."]
return step
ranges = {}
for j in soft:
f = ctx.study.factors[j]
mid, half = (f.high + f.low) / 2, (f.high - f.low) / 2 # type: ignore[operator]
ranges[f.name] = [mid - _WIDEN * half, mid + _WIDEN * half]
runs = ctx.study.data[ctx.study.names].to_numpy(dtype=float, copy=True)
runs[:, soft] *= _WIDEN # the same coded design, read on the old scale: the new ends are 1.5 old half-ranges out
numeric = sum(f.type != FactorType.categorical for f in ctx.study.factors)
which = "every numeric range" if len(soft) == numeric else f"the ranges of {_names(ctx, soft)}"
kept = [j for j in range(ctx.k) if ctx.study.numeric[j] and j not in soft]
reasons.append("Wider ranges make each factor's effect larger against the same noise.")
if kept:
one = len(kept) == 1
held = (
"range of {} is a hard limit, and stays as it is"
if one
else "ranges of {} are hard limits, and stay as they are"
)
reasons.append(f"The {held.format(_names(ctx, kept))}.")
return _Step(
"widen_ranges",
f"Widen {which} by half again about the middle, and run the {len(runs)}-run design again.",
reasons,
runs,
"Analyse the new runs on their own ranges; if still nothing stands out, check the measurement system.",
{"new_ranges": ranges},
)
def _de_alias(ctx: _Context) -> _Step | None:
"""Return the runs that separate each ambiguous active term from the interaction it is aliased with."""
if not ctx.aliasing:
return None
pairs = ", ".join(f"{ctx.study.real(a['term'])} with {ctx.study.real(a['aliased_with'])}" for a in ctx.aliasing)
reasons = [
(
f"Active effects are aliased with interactions that heredity makes plausible ({pairs}): the data cannot "
"say which of each pair is real, and the optimum depends on it."
),
]
two_level = _two_level_part(ctx.study, ctx.design)
mains = any(":" not in a["term"] for a in ctx.aliasing)
if ctx.design.two_level and mains:
new = augment_design(two_level, "foldover", factor_names=ctx.study.names)["new_runs"]
headline = f"Fold the design over: {len(new)} runs, each a run made before with every setting switched."
reasons.append("The fold-over separates every main effect from every two-factor interaction.")
else:
partners = sorted({a["aliased_with"] for a in ctx.aliasing})
rhs = " + ".join(dict.fromkeys([*ctx.model.terms, *partners]))
n_terms = len(ctx.model.terms) + len(partners) + 1
n_new = max(len(partners), n_terms - len(ctx.study.data)) + _SPARE_RUNS
frame = ctx.study.data[ctx.study.names]
new = augment_design(
frame,
"add_runs_optimal",
target_model=rhs,
n_additional_runs=n_new,
factor_names=ctx.study.names,
random_state=ctx.seed,
)["new_runs"]
headline = (
f"Add {len(new)} runs chosen to estimate each active term apart from the interaction it is aliased with."
)
reasons.append(
"Runs chosen for the model that holds both terms of each pair (D-optimal) let it tell them apart."
)
runs = pd.DataFrame(new)[ctx.study.names].to_numpy(dtype=float)
return _Step(
"de_alias",
headline,
reasons,
runs,
"Analyse every run together: the model can then tell each pair apart, and the next step follows from it.",
)
def _quadratic_runs(ctx: _Context, factors: list[int], suspects: list[int], *, search: bool = True) -> _Step | None:
"""Return runs that let the full quadratic model in ``factors`` be fitted: axial runs, or runs chosen for it.
None when there is no numeric factor among them, or the runs already fit that model (see
:func:`_quadratic_fits`). The other numeric factors are held at the middle of their ranges,
as in a central composite design, so the runs added tell the squares of ``factors`` from
theirs; categorical factors are held at the level the optimum takes. The ``suspects``,
factors whose squares the runs so far cannot tell from the curvature (see
:func:`_square_suspects`), vary too, in runs chosen to separate their squares. Without
``search``, runs chosen D-optimally are left blank (see :func:`_d_optimal_runs`).
"""
numeric = [j for j in factors if ctx.study.numeric[j]]
if not numeric or (not suspects and _quadratic_fits(ctx, numeric)):
return None
added, how = (
_separating_runs(ctx, numeric, suspects, search) if suspects else _quadratic_design(ctx, numeric, search)
)
varied = [*numeric, *suspects]
held = [j for j in range(ctx.k) if ctx.study.numeric[j] and j not in varied]
if held:
how += f" Hold {_names(ctx, held)} at the middle of {'its range' if len(held) == 1 else 'their ranges'}."
runs = _with_others(_held_settings(ctx), added, varied, ctx.study)
return _Step(
"augment_to_rsm",
how,
[],
runs,
"Fit the full quadratic model to every run; its optimum, inside the ranges or on an edge, is "
"where the confirmation runs go next.",
)
def _separating_runs(ctx: _Context, numeric: list[int], suspects: list[int], search: bool) -> tuple[pd.DataFrame, str]:
"""Return runs chosen to fit the full quadratic model in ``numeric`` and each suspect's own square, and why."""
chosen = [ctx.study.names[j] for j in numeric]
others = [ctx.study.names[j] for j in suspects]
rhs = " + ".join(
[
*chosen,
*(f"{a}:{b}" for a, b in itertools.combinations(chosen, 2)),
*(f"I({n} ** 2)" for n in chosen),
*others,
*(f"I({n} ** 2)" for n in others),
]
)
frame = ctx.study.data[[*chosen, *others]]
columns = np.asarray(dmatrix(rhs, frame), dtype=float)
n_new = columns.shape[1] - int(np.linalg.matrix_rank(columns)) + _SPARE_RUNS
added = _d_optimal_runs(ctx, frame, rhs, n_new, search)
how = (
f"Add {len(added)} runs chosen (D-optimal) to fit the full quadratic model in {_names(ctx, numeric)} and to "
f"tell the squares of {_names(ctx, suspects)} apart: so far the runs give "
f"{'its square' if len(suspects) == 1 else 'their squares'} the same column as the curvature they show, "
"which could be any of theirs."
)
return added, how
def _quadratic_design(ctx: _Context, numeric: list[int], search: bool) -> tuple[pd.DataFrame, str]:
"""Return the runs that complete a face-centred central composite design in ``numeric``, or D-optimal runs."""
names = [ctx.study.names[j] for j in numeric]
frame = ctx.study.data[names]
quadratic = _model_columns(frame, "quadratic")
away = ~ctx.design.center
two_level_here = bool(np.all(np.abs(np.abs(frame[away].to_numpy(dtype=float)) - 1) < _AT_EDGE))
interactions = _model_columns(frame, "interactions")
interactions_fit = len(names) >= 2 and int(np.linalg.matrix_rank(interactions)) == interactions.shape[1]
if two_level_here and interactions_fit:
result = augment_design(frame, "upgrade_to_rsm", alpha="face_centered", factor_names=names)
added = pd.DataFrame(result["new_runs"])
how = (
f"Add {len(added)} axial and centre runs in {_names(ctx, numeric)}: with the runs before, a face-centred "
"central composite design, which fits the full quadratic model and keeps every run inside the ranges."
)
else:
n_new = quadratic.shape[1] - int(np.linalg.matrix_rank(quadratic)) + _SPARE_RUNS
added = _d_optimal_runs(ctx, frame, "quadratic", n_new, search)
how = f"Add {len(added)} runs chosen (D-optimal) to fit the full quadratic model in {_names(ctx, numeric)}."
return added, how
def _d_optimal_runs(ctx: _Context, frame: pd.DataFrame, model: str, n_new: int, search: bool) -> pd.DataFrame:
"""Return ``n_new`` runs chosen (D-optimal) to fit ``model`` with the runs in ``frame``, in its columns.
Without ``search`` the runs are blank, for :func:`_alternatives`, which reads only how many
there are: the search always adds ``n_new`` (or refuses before it starts), and it is most of
the time a large study's step takes, 15 s of 16 for 11 factors in a Plackett-Burman design.
"""
if not search:
return pd.DataFrame(np.nan, index=range(n_new), columns=frame.columns)
result = augment_design(
frame,
"add_runs_optimal",
target_model=model,
n_additional_runs=n_new,
factor_names=list(frame.columns),
random_state=ctx.seed,
)
return pd.DataFrame(result["new_runs"])
def _quadratic_fits(ctx: _Context, numeric: list[int]) -> bool:
"""Say whether the runs fit the full quadratic model in these numeric factors, apart from every other square.
Centre runs among two-level runs give every factor three levels, so the quadratic model in
any one factor looks estimable; but its square is then the same column as every other
factor's, and the curvature the centre runs show could be any of theirs. The squares of
these factors must therefore be told apart from the other factors' squares. (A definitive
screening design's squares are; that they are partly aliased with its interactions is
what its own analysis resolves.)
"""
data = ctx.study.data
chosen = data[[ctx.study.names[j] for j in numeric]].to_numpy(dtype=float)
quadratic = _model_columns(pd.DataFrame(chosen), "quadratic")
if int(np.linalg.matrix_rank(quadratic)) < quadratic.shape[1]:
return False
others = [n for j, n in enumerate(ctx.study.names) if ctx.study.numeric[j] and j not in numeric]
if not others:
return True
own = np.column_stack([np.ones(len(chosen)), chosen**2])
squares = data[others].to_numpy(dtype=float) ** 2
squares -= squares.mean(axis=0)
together = int(np.linalg.matrix_rank(np.hstack([own, squares])))
return together == own.shape[1] + int(np.linalg.matrix_rank(squares))
def _model_columns(frame: pd.DataFrame, model: str) -> np.ndarray:
x = frame.to_numpy(dtype=float)
parts = [np.ones((len(x), 1)), x]
if model in {"interactions", "quadratic"}:
parts += [x[:, [i]] * x[:, [j]] for i, j in itertools.combinations(range(x.shape[1]), 2)]
if model == "quadratic":
parts += [x**2]
return np.column_stack(parts)
def _held_settings(ctx: _Context) -> np.ndarray:
"""Return where added runs hold the other factors: numeric at the middle, categorical as at the optimum."""
held = ctx.optimum.copy()
held[np.array(ctx.study.numeric)] = 0.0
return held
def _move_region(ctx: _Context) -> _Step | None:
"""Return runs along the path of steepest ascent, past the edges of the ranges the goal improves beyond."""
if not ctx.ranges_can_move or ctx.gain_beyond <= 0:
return None
noise = ctx.model.sigma
if not np.isfinite(noise) or ctx.gain_beyond < _WORTH_MOVING * noise:
return None
gradient = ctx.model.gradient(ctx.optimum)
outward = np.zeros(ctx.k)
for j, edge in enumerate(_edge_flags(ctx)):
if edge:
outward[j] = ctx.sense * gradient[j]
direction = outward / np.abs(outward).max()
runs = np.array([ctx.optimum + t * direction for t in _PATH_STEPS])
predicted = ctx.model.value(runs)
moving = [j for j in range(ctx.k) if direction[j] != 0]
return _Step(
"move_region",
f"Move past the ranges of {_names(ctx, moving)}: {len(runs)} runs along the path of steepest ascent.",
[
(
f"The best settings in the ranges sit at their edge, and the model's slope there gains "
f"{_g(ctx.gain_beyond)} per coded unit beyond it, more than {_WORTH_MOVING:g} times the noise "
f"(standard deviation {_g(noise)})."
),
(
"Past the ranges the model is extrapolating: its predictions along the path "
f"({', '.join(_g(v) for v in predicted)}) say where to look, not what will be found."
),
*_held_at_limits(ctx),
],
runs,
"Make the runs in order and stop at the first one that does worse than the run before it. Then centre a "
"new study on the best run, with ranges as wide as these, and start again with its screening or "
"response surface design.",
{
"path_predictions": [float(v) for v in predicted],
"direction": {ctx.name(j): float(direction[j]) for j in moving},
},
)
def _edge_flags(ctx: _Context) -> list[bool]:
"""For each factor, whether it is at an edge the goal improves beyond, and may go past it."""
beyond = {e["factor"] for e in ctx.edges if e["improves_beyond"] and not e["hard_limit"]}
return [f.name in beyond for f in ctx.study.factors]
def _held_at_limits(ctx: _Context) -> list[str]:
"""Say which factors the goal would take past their hard limits, which hold them at that end."""
held = [e["factor"] for e in ctx.edges if e["improves_beyond"] and e["hard_limit"]]
if not held:
return []
names = ", ".join(held)
one = len(held) == 1
return [
(
f"The response would improve past the {'end' if one else 'ends'} of {names}, but "
f"{'its range is a hard limit' if one else 'their ranges are hard limits'}: "
f"{'it stays at that end' if one else 'they stay at those ends'}."
)
]
def _no_effect(ctx: _Context) -> _Step:
"""Stop: nothing stood out again after the ranges were widened, or after the runs were made again."""
tried = "widened" if "widen_ranges" in ctx.history else "made again"
return _Step(
"stop",
"Stop: no factor moves the response by more than the noise.",
[
f"Nothing stands out from the noise ({ctx.model.method}), as before the runs were {tried}.",
(
"Over the ranges the process allows, these factors do not change the response by more than its "
f"run-to-run noise (standard deviation {_g(ctx.model.sigma)}). Check the measurement's repeatability, "
"and look for factors this study did not vary."
),
],
np.zeros((0, ctx.k)),
"Set the factors for cost or convenience. A new study needs other factors, or a better measurement.",
)
def _investigate_confirmation(ctx: _Context, why: str) -> _Step:
"""No more runs yet: the confirmation runs miss the prediction by far, though the runs fit the model they need."""
return _Step(
"investigate",
"Look at the process before running more: the confirmation runs miss the prediction by far.",
[
why,
(
"The runs already fit the model the active factors need, so more runs of the same kind would not "
"explain the miss. Check whether the process shifted between the rounds (material, equipment, "
"season), whether a run is an outlier, and whether a factor held constant changed."
),
],
np.zeros((0, ctx.k)),
"Once the cause is found, call recommend_next_step again; a shift between rounds is modelled as a block.",
)
def _repeat_confirmation(ctx: _Context, why: str) -> _Step:
"""More confirmation runs at the same settings, after a miss small enough to be chance."""
step = _confirm(ctx)
n = ctx.n_confirmation_runs
step.headline = f"Repeat the confirmation: {n} more runs at the same settings."
step.reasons = [
why,
(
f"The miss is small enough to be chance: a {100 * (1 - ctx.alpha):.0f}% interval misses about one time in "
f"{round(1 / ctx.alpha)} when the model is right."
),
*step.reasons[:1],
]
step.after = (
"Add the confirmation runs that missed to the runs, and call recommend_next_step with only these new runs as "
"the confirmation runs: the model then knows the response at these settings, and the new runs are judged on "
"their own."
)
return step
def _stop(ctx: _Context, judged: dict[str, Any]) -> _Step:
return _Step(
"stop",
"Stop: the confirmation runs agree with the model at its optimum.",
[
(
f"Their mean, {_g(judged['observed_mean'])}, lies inside the {100 * (1 - ctx.alpha):.0f}% prediction "
f"interval for a mean of {judged['n_runs']} runs, {_g(judged['pi_low'])} to {_g(judged['pi_high'])}."
),
],
np.zeros((0, ctx.k)),
"Use these settings. Run them a few more times in production to see the process hold there.",
)
# --------------------------------------------------------------------------- the decision
def _curved(ctx: _Context) -> bool:
"""Whether the runs show curvature, or a lack of fit, that the selected model does not capture."""
if ctx.curvature.get("source") == "center_runs" and ctx.curvature.get("significant"):
return True
return bool(ctx.lack_of_fit.get("tested") and ctx.lack_of_fit.get("significant"))
#: Curvature that the runs cannot attribute gets axial runs in every numeric factor up to this many. Near an
#: optimum a factor's linear effect vanishes while its square does not, so the active factors are a poor guide.
_AXIAL_FOR_ALL = 6
def _rsm_step(ctx: _Context, why: str) -> _Step | None:
numeric = [j for j in range(ctx.k) if ctx.study.numeric[j]]
attributed = ctx.curvature.get("source") == "squares"
if len(numeric) <= _AXIAL_FOR_ALL and not attributed:
factors = numeric
else:
factors = ctx.active_factors or _model_factor_indices(ctx) or numeric
step = _quadratic_runs(ctx, factors, _square_suspects(ctx, factors))
if step is not None:
step.reasons = [why, *step.reasons]
return step
def _square_suspects(ctx: _Context, factors: list[int]) -> list[int]:
"""Return the numeric factors outside ``factors`` whose square the runs cannot tell from the curvature they show.
Centre runs among two-level runs give every factor's square the same column, so curvature
at the centre could be any factor's; and a square judged active stands for every square
with the same column, of which the analysis names only the first. Each such factor needs
runs that separate its square, or the optimum is sought in the wrong factor.
"""
data = ctx.study.data
references = [data[_square_factor(t)].to_numpy(dtype=float) ** 2 for t in ctx.curvature.get("terms") or []]
if ctx.curvature.get("source") == "center_runs":
references.append(data[_CURVATURE].to_numpy(dtype=float))
return [
j
for j, name in enumerate(ctx.study.names)
if ctx.study.numeric[j]
and j not in factors
and any(_identical(data[name].to_numpy(dtype=float) ** 2, r) for r in references)
]
def _square_factor(term: str) -> str:
"""Return the factor a square term squares."""
return _canonical(term)[1][0]
def _identical(a: np.ndarray, b: np.ndarray) -> bool:
"""Say whether two columns that vary are the same up to sign, scale and a constant: exactly aliased."""
a, b = a - a.mean(), b - b.mean()
scale = float(np.linalg.norm(a)) * float(np.linalg.norm(b))
return scale > 1e-12 and abs(float(a @ b)) > (1 - 1e-9) * scale
def _full_quadratic_fits(ctx: _Context) -> bool:
"""Say whether the runs already fit the full quadratic model in the active numeric factors."""
numeric = [j for j in ctx.active_factors if ctx.study.numeric[j]]
return bool(numeric) and _quadratic_fits(ctx, numeric)
def _decide(ctx: _Context) -> _Step:
"""Apply the rules of the module docstring, in order, to a study without confirmation runs to judge."""
if ctx.model.error_df < _MIN_ERROR_DF and ctx.design.pure_error_df < _MIN_ERROR_DF:
return _add_center_points(ctx) or _replicate(ctx)
if not ctx.model.signal:
return _nothing_stands_out(ctx)
step = _de_alias(ctx) or _curvature_step(ctx) or _lack_of_fit_step(ctx)
if step is not None:
return step
step = _move_region(ctx) or _confirm(ctx)
if ctx.lack_of_fit.get("significant"):
step.reasons.append(
f"The lack-of-fit test is significant (p = {_g(ctx.lack_of_fit['p_value'])}), though the runs fit the full "
"quadratic model in every factor with some evidence of an effect: the model is off somewhere in the "
"region, or the repeats agree more closely than the noise usually allows. Runs at the optimum test it "
"where it matters; look for an outlying run if they miss."
)
if _blocks_hold_the_centre_runs(ctx.study, ctx.design):
step.reasons.append(
"The centre runs share no block with a corner run, so they differ from the corners by the shift between "
"the blocks as well: the runs cannot show curvature. Centre runs made in a block with corner runs test it."
)
return step
#: When the model does not fit the repeats, a factor with a term this significant joins the active ones: something is
#: missing, and a missed effect then costs more than a spare factor, as in screening.
_SOME_EVIDENCE = 0.2
def _lack_of_fit_step(ctx: _Context) -> _Step | None:
"""Return runs for the full quadratic model in every factor with some evidence of an effect, when the model misses.
The curvature step has found the active factors' full quadratic model fitted. A lack of fit then
points at factors judged inactive: the effects they have went into the residual. When no factor
beyond the active ones shows some evidence (p below :data:`_SOME_EVIDENCE`), or the runs already
fit that larger model, the runs for the interactions the second look could not test follow
(:func:`_interaction_runs`). None when the model fits.
"""
if not ctx.lack_of_fit.get("significant"):
return None
evidence = {
f for t, p in ctx.model.judged.items() if t != _CURVATURE and p < _SOME_EVIDENCE for f in _term_factors(t)
}
joining = [
j
for j, name in enumerate(ctx.study.names)
if name in evidence and j not in ctx.active_factors and ctx.study.numeric[j]
]
step = _quadratic_runs(ctx, sorted([*ctx.active_factors, *joining]), []) if joining else None
if step is None:
return _interaction_runs(ctx)
step.action = "add_runs_for_model"
step.reasons = [
(
f"The model does not fit the repeats (lack of fit p = {_g(ctx.lack_of_fit['p_value'])}): it misses "
f"effects. {_names(ctx, joining)} {'shows' if len(joining) == 1 else 'show'} some evidence of an effect "
f"(p below {_SOME_EVIDENCE:g}) and {'joins' if len(joining) == 1 else 'join'} the active factors."
),
*step.reasons,
]
return step
def _interaction_runs(ctx: _Context) -> _Step | None:
"""Return runs that let the interactions of the active factors with the others be tested, when the model misses.
No factor beyond the active ones shows some evidence of an effect, or the runs already fit their
full quadratic model, so the effects the model misses are likeliest to be interactions of an
active factor with one judged inactive (see :func:`_second_look`). Some of those the runs could
not test: aliased with each other, as in a definitive screening design, or short of degrees of
freedom. None when the model the terms were judged in fits the repeats, since the second look
then did not look for them, or when every one was tested.
"""
if not ctx.model.misses:
return None
names, numeric = ctx.study.names, ctx.study.numeric
numeric_names = {n for n, is_numeric in zip(names, numeric, strict=True) if is_numeric}
untested = [
c
for c in _interaction_candidates(ctx.study, ctx.active, weak=True)
if set(c.split(":")) <= numeric_names and c not in ctx.model.judged
]
if not untested:
return None
rhs = " + ".join(dict.fromkeys([*names, *ctx.model.terms, *untested]))
frame = ctx.study.data[names]
columns = np.asarray(dmatrix(rhs, frame), dtype=float)
rank = int(np.linalg.matrix_rank(columns))
blocks = len(_nuisance(ctx.study.data)) + 1 # the new runs are a block of their own
n_new = max(columns.shape[1] - rank, columns.shape[1] + blocks + _MIN_ERROR_DF - len(frame), 0) + _SPARE_RUNS
new = augment_design(
frame, "add_runs_optimal", target_model=rhs, n_additional_runs=n_new, factor_names=names, random_state=ctx.seed
)["new_runs"]
active = _names(ctx, [j for j in ctx.active_factors if numeric[j]])
return _Step(
"add_runs_for_model",
f"Add {len(new)} runs chosen (D-optimal) to test the interactions of {active} with the other factors.",
[
(
f"The model does not fit the repeats (lack of fit p = {_g(ctx.lack_of_fit['p_value'])}): it misses "
f"effects, and no other factor shows them through its main effect. The likeliest are interactions "
f"of {active} with the others, and the runs so far cannot test {len(untested)} of them."
),
],
pd.DataFrame(new)[names].to_numpy(dtype=float),
"Analyse every run together: the interactions of the active factors are then tested, and the next step "
"follows from them.",
)
def _nothing_stands_out(ctx: _Context) -> _Step:
"""More error degrees of freedom, wider ranges, or a stop, when no effect stands out from the noise."""
if ctx.model.error_df < _FEW_ERROR_DF:
return _more_error_df(ctx)
if "widen_ranges" in ctx.history or (not ctx.ranges_can_move and "replicate" in ctx.history):
return _no_effect(ctx)
return _widen_ranges(ctx)
def _curvature_step(ctx: _Context) -> _Step | None:
"""Return the runs a curved response needs before its optimum can be located, or None when the runs suffice."""
if _curved(ctx):
why = (
f"The centre runs are off the plane of the corners (curvature test p = {_g(ctx.curvature['p_value'])})."
if ctx.curvature.get("source") == "center_runs"
else f"The model does not fit the repeats (lack of fit p = {_g(ctx.lack_of_fit['p_value'])})."
)
step = _rsm_step(ctx, why) # None when the runs already fit the full quadratic model
if step is not None:
return step
squares = ctx.curvature.get("terms") or []
if squares and not _full_quadratic_fits(ctx):
why = (
f"Squares are active ({', '.join(ctx.study.real(t) for t in squares)}), but the runs cannot fit the full "
"quadratic model in the active factors, so the optimum cannot be located."
)
step = _rsm_step(ctx, why)
if step is not None:
step.action = "add_runs_for_model"
return step
return None
def _learn_more(ctx: _Context, why: str) -> _Step | None:
"""After a failed confirmation: the runs that add the information the model lacked, or None when none would."""
for build in (
_de_alias,
_add_center_points,
lambda c: _rsm_step(c, why) if not _full_quadratic_fits(c) else None,
):
step = build(ctx)
if step is not None:
step.reasons = [why, *step.reasons]
return step
return None
def _judge(ctx: _Context, confirmation: pd.DataFrame, raw: pd.Series) -> dict[str, Any]:
"""Judge the confirmation runs: the mean at each setting against the interval for such a mean."""
points = _points(ctx.study, confirmation[ctx.study.names].to_numpy(dtype=float))
summary = ctx.model.fit.get_prediction(points).summary_frame(alpha=ctx.alpha)
observed = raw.to_numpy(dtype=float)
tests = _mean_tests(ctx.model.fit, points[ctx.study.names], summary, (observed, observed), ctx.alpha)
in_model = _model_factor_indices(ctx)
for test in tests:
at = np.array([test["settings"][name] for name in ctx.study.names])
test["at_optimum"] = bool(np.all(np.abs(at[in_model] - ctx.optimum[in_model]) <= _SAME_SETTING))
at_optimum = [t for t in tests if t["at_optimum"]]
return {"tests": tests, "at_optimum": at_optimum[0] if at_optimum else None}
# --------------------------------------------------------------------------- the public function
[docs]
def recommend_next_step( # noqa: PLR0913 - each setting is a separate input to the decision
factors: list[Factor],
runs: pd.DataFrame,
response: str,
*,
goal: Goal = "maximize",
target: float | None = None,
confirmation_runs: pd.DataFrame | None = None,
model_terms: list[str] | None = None,
budget: int | None = None,
ranges_can_move: bool = True,
n_confirmation_runs: int = 3,
alpha: float = 0.05,
random_state: int = 42,
history: list[str] | None = None,
block: str | None = None,
) -> dict[str, Any]:
"""Recommend the next round of a designed experiment from every run so far.
Parameters
----------
factors : list[Factor]
The study's factors: continuous, or categorical at two levels. Their ``low`` and
``high`` define the coded scale (-1 to +1) and the region the optimum is sought in.
runs : pandas.DataFrame
Every run of the rounds so far, one row each, with a column per factor in its own
units (a categorical factor's level) and the response. Runs with no response yet are
left out.
response : str
The response column.
goal : {"maximize", "minimize", "target"}
What the study wants of the response.
target : float or None
The value aimed at, for ``goal="target"``.
confirmation_runs : pandas.DataFrame or None
Runs made to confirm a predicted optimum, judged against the model the other runs
give, with the same columns as ``runs``.
model_terms : list[str] or None
The model to use, as terms (``"A"``, ``"A:B"``, ``"I(A ** 2)"``), its terms tested in it;
``None`` judges every term the runs can test and keeps the active ones (see the module
documentation).
budget : int or None
Runs left. An action that needs more falls back to confirmation runs, or to stopping.
ranges_can_move : bool
Whether the next round may set factors past their ranges: the path of steepest ascent,
or wider ranges. ``False`` keeps every recommended run inside them.
n_confirmation_runs : int
Runs made to confirm an optimum.
alpha : float
Significance level of every test the decision reads.
random_state : int
Seed for the search for the optimum and for runs chosen by exchange.
history : list[str] or None
The actions of the earlier rounds of this study, oldest first, as this function
returned them. With ``"widen_ranges"`` among them (or ``"replicate"``, when the ranges
cannot move), a study in which nothing stands out stops rather than widen the ranges
again.
block : str or None
The column of ``runs`` that holds each run's block: the round it was made in, say. Every
model then holds a shift for each block, so a process that drifted between the rounds
does not bias the effects or swell the noise. The blocks are taken in the order they
first appear, and runs made next are predicted as the process was in the latest one.
Confirmation runs need no block; when they miss and join the runs, they count in the
latest block.
Returns
-------
dict[str, Any]
``action`` (one of :data:`ACTIONS`), ``headline``, ``reasons`` (the evidence, with its
numbers), ``runs`` (the next round, in the factors' own units, one dict per run),
``n_runs``, ``after`` (how to judge those runs), ``model`` (the terms, coefficients and
p-values of the model the decision read, in the factors' names), ``evidence`` (active
factors, curvature, lack of fit, aliasing, the optimum with its prediction and the edges
it sits on, the confirmation runs' verdict), ``alternatives`` (other actions that apply,
with their run counts) and, for some actions, ``new_ranges`` (``widen_ranges``), or
``path_predictions`` and ``direction``, each moving factor's step in coded units per unit
of the factor that moves most (``move_region``).
Raises
------
ValueError
For a mixture factor, a categorical factor not at two levels, a factor or response
missing from the runs, a level not among a factor's levels, a target goal without a
target, an action in ``history`` not among :data:`ACTIONS`, a ``block`` that is not a
column of the runs, or fewer than two runs with a response.
Examples
--------
>>> import pandas as pd
>>> from process_improve.experiments import Factor, TrueProcess, generate_design
>>> factors = [Factor(name=n, low=-1, high=1) for n in "ABCDEFGH"]
>>> design = generate_design(factors, design_type="dsd", n_center_points=5).design_actual
>>> truth = TrueProcess({"Intercept": 100, "A": 10, "C": 7, "G": -6, "A:C": 8}, sigma=3)
>>> runs = design.assign(y=truth.simulate(design.drop(columns="RunOrder"), random_state=1))
>>> step = recommend_next_step(factors, runs, "y", ranges_can_move=False)
>>> step["action"], step["evidence"]["active_terms"]
('confirm', ['A', 'C', 'G', 'A:C'])
"""
choices = _Choices(goal, target, alpha, n_confirmation_runs, ranges_can_move, random_state, tuple(history or ()))
_check_choices(factors, runs, choices, block)
data = _with_blocks(_code(factors, runs, response, block).dropna(subset=["y"]).reset_index(drop=True))
if len(data) < 2:
raise ValueError("At least two runs with a response are needed to recommend the next step.")
ctx = _context(_Study(factors, data, response), model_terms, choices)
judged = None
if confirmation_runs is not None and len(confirmation_runs):
confirming = _code(factors, confirmation_runs, response, block).dropna(subset=["y"]).reset_index(drop=True)
judged = _judge(ctx, confirming, confirming["y"])
ctx, step = _after_confirmation(ctx, confirming, judged["at_optimum"], model_terms)
else:
step = _decide(ctx)
step = _within_budget(ctx, step, budget)
return _reply(ctx, step, judged)
def _check_choices(factors: list[Factor], runs: pd.DataFrame, choices: _Choices, block: str | None) -> None:
"""Refuse choices :func:`recommend_next_step` cannot act on, with what is wrong."""
_check_factors(factors)
if choices.goal not in ("maximize", "minimize", "target"):
raise ValueError(f"goal must be 'maximize', 'minimize' or 'target', not {choices.goal!r}.")
if choices.goal == "target" and choices.target is None:
raise ValueError("A target goal needs its target.")
if choices.n_confirmation_runs < 1:
raise ValueError(f"n_confirmation_runs must be at least 1, not {choices.n_confirmation_runs}.")
unknown = sorted(set(choices.history) - set(ACTIONS))
if unknown:
raise ValueError(f"history holds {unknown}, which are not actions recommend_next_step returns: {ACTIONS}.")
if block is not None and block not in runs.columns:
raise ValueError(f"The runs have no block column {block!r}.")
def _after_confirmation(
ctx: _Context, confirming: pd.DataFrame, verdict: dict[str, Any] | None, model_terms: list[str] | None
) -> tuple[_Context, _Step]:
"""Return the context the next step reads and the step, once confirmation runs are judged (``verdict``).
Runs inside the interval stop the study, or move it on when the ranges may move. Otherwise
they are data, made after the other runs (in the last block, when there are blocks), and the
model is chosen again with them in.
"""
if verdict is not None and verdict["within_PI"]:
return ctx, _move_region(ctx) or _stop(ctx, verdict)
data = ctx.study.data
if _BLOCK_LABEL in data.columns:
last = next((v for v in reversed(data[_BLOCK_LABEL].tolist()) if v is not None), None)
confirming[_BLOCK_LABEL] = confirming[_BLOCK_LABEL].map(lambda v: last if v is None else v)
together = pd.concat([data.drop(columns=[*_nuisance(data), _CURVATURE]), confirming], ignore_index=True)
choices = _Choices(
ctx.goal, ctx.target, ctx.alpha, ctx.n_confirmation_runs, ctx.ranges_can_move, ctx.seed, ctx.history
)
judged_ctx = ctx
# After a miss the missed runs are repeats far from the model: judged against them, terms would chase the miss
ctx = _context(
_Study(ctx.study.factors, _with_blocks(together), ctx.study.response),
model_terms,
choices,
against_repeats=verdict is None,
)
if verdict is not None:
return ctx, _after_a_miss(ctx, judged_ctx, verdict)
step = _decide(ctx)
step.reasons.append(
"The confirmation runs were not at the optimum the model now predicts, so they join the other runs."
)
return ctx, step
def _after_a_miss(ctx: _Context, judged_ctx: _Context, verdict: dict[str, Any]) -> _Step:
"""Return the step after confirmation runs that miss the prediction.
Runs that add what the design behind the prediction could not show come first. Then, if
the model chosen again with the missed runs among its data puts its optimum elsewhere (a
term the first model held by chance, say), that optimum is confirmed. Otherwise a miss
within :data:`_GROSS_MISS` half-widths of the interval may be chance (the noise estimate
behind the interval can itself be low), and the confirmation is made again; ``ctx`` holds
the missed runs as data, so the model already knows the response there, and the new runs
are judged on their own. Pooling them with the runs that missed would carry the first miss
into the second test. A larger miss, or a miss after a repeat (two confirmations in
``history``), calls for a look at the process instead.
"""
why = (
f"The confirmation runs' mean, {_g(verdict['observed_mean'])}, lies outside the prediction interval "
f"{_g(verdict['pi_low'])} to {_g(verdict['pi_high'])} for a mean of {verdict['n_runs']} runs: the model "
"missed something, or the miss is chance."
)
step = _learn_more(judged_ctx, why) # what the design the prediction came from could not show
if step is not None:
return step
in_model = _model_factor_indices(ctx)
if np.any(np.abs(ctx.optimum[in_model] - judged_ctx.optimum[in_model]) > _SAME_SETTING):
step = _confirm(ctx)
step.headline = f"Confirm the optimum the model now predicts, with {ctx.n_confirmation_runs} runs there."
step.reasons = [
why,
"With the confirmation runs among its data, the model is chosen again, and its optimum has moved.",
*step.reasons,
]
return step
half = (verdict["pi_high"] - verdict["pi_low"]) / 2
gross = abs(verdict["observed_mean"] - verdict["predicted"]) > _GROSS_MISS * half
if gross or ctx.history.count("confirm") >= 2:
return _investigate_confirmation(ctx, why)
return _repeat_confirmation(ctx, why)
@dataclass(frozen=True)
class _Choices:
"""The caller's choices that every decision reads."""
goal: Goal
target: float | None
alpha: float
n_confirmation_runs: int
ranges_can_move: bool
seed: int
history: tuple[str, ...]
def _context(
study: _Study, model_terms: list[str] | None, choices: _Choices, *, against_repeats: bool = True
) -> _Context:
design = _design(study)
study.data[_CURVATURE] = design.center.astype(float)
model = _select(study, design, model_terms, choices.alpha, against_repeats=against_repeats)
ctx = _Context(study, design, model, **vars(choices))
return _gather(ctx)
def _within_budget(ctx: _Context, step: _Step, budget: int | None) -> _Step:
"""Fall back to confirmation runs, then to stopping, when the step needs more runs than are left."""
if budget is None or len(step.runs) <= budget:
return step
if step.action != "confirm" and ctx.n_confirmation_runs <= budget and ctx.active:
fallback = _confirm(ctx)
fallback.reasons.insert(0, f"{step.headline} That needs {len(step.runs)} runs; {budget} are left.")
return fallback
return _Step(
"stop",
f"Stop: {budget} run(s) are left, fewer than the next step needs.",
[
f"{step.headline} That needs {len(step.runs)} runs.",
f"The best settings the model knows predict {_g(ctx.predicted())}.",
],
np.zeros((0, ctx.k)),
"Use the best settings the model knows, and plan the next study when runs are available.",
)
def _alternatives(ctx: _Context, chosen: str) -> list[dict[str, Any]]:
"""Return the other actions that apply to this study, each with its run count and headline."""
options = []
for build in (
_confirm,
_add_center_points,
_de_alias,
lambda c: _quadratic_runs(c, c.active_factors, _square_suspects(c, c.active_factors), search=False),
_move_region,
):
try:
step = build(ctx)
except ValueError: # an augmentation the runs cannot take
continue
if step is not None and step.action != chosen:
options.append({"action": step.action, "n_runs": len(step.runs), "headline": step.headline})
return options
def _reply(ctx: _Context, step: _Step, judged: dict[str, Any] | None) -> dict[str, Any]:
study, model = ctx.study, ctx.model
low, high = ctx.mean_interval(ctx.n_confirmation_runs)
nuisance = _nuisance(study.data)
coefficients = [
{"term": study.real(str(term)), "coefficient": float(value), "p_value": float(model.fit.pvalues[term])}
for term, value in model.fit.params.items()
if str(term) not in nuisance
]
evidence: dict[str, Any] = {
"n_runs": ctx.design.n_runs,
"n_blocks": len(nuisance) + 1,
"active_terms": [study.real(t) for t in ctx.active],
"active_factors": [ctx.name(j) for j in ctx.active_factors],
"noise_sd": model.sigma,
"error_df": model.error_df,
"curvature": {**ctx.curvature, "terms": [study.real(t) for t in ctx.curvature.get("terms", [])]},
"lack_of_fit": ctx.lack_of_fit,
"aliasing": [
{**a, "term": study.real(a["term"]), "aliased_with": study.real(a["aliased_with"])} for a in ctx.aliasing
],
"optimum": {
"settings": study.actual(ctx.optimum),
"coded": dict(zip([f.name for f in study.factors], (float(v) for v in ctx.optimum), strict=True)),
"predicted": ctx.predicted(),
"interval_for_mean": [low, high],
"n_runs_in_mean": ctx.n_confirmation_runs,
"edges": ctx.edges,
},
}
if judged is not None:
evidence["confirmation"] = [
{**t, "settings": study.actual(np.array([t["settings"][n] for n in study.names]))} for t in judged["tests"]
]
reply = {
"action": step.action,
"headline": step.headline,
"reasons": step.reasons,
"runs": _rows(study, step.runs) if len(step.runs) else [],
"n_runs": len(step.runs),
"after": step.after,
"model": {
"terms": [study.real(t) for t in model.terms],
"judged_by": model.method,
"judged": {study.real(t): p for t, p in model.judged.items()},
"signal": model.signal,
"coefficients": coefficients,
"r_squared_adj": float(model.fit.rsquared_adj) if model.error_df > 0 else math.nan,
},
"evidence": evidence,
"alternatives": _alternatives(ctx, step.action),
}
reply.update(step.extra)
return reply
