diff --git a/CHANGELOG.md b/CHANGELOG.md index e78e309a..c32cf9a4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,60 @@ those changes. ## [Unreleased] +## [1.62.0] - 2026-07-26 + +### Added + +- `optimize_responses(..., fitted_results=[...], significance_level=0.05)` reports a + confidence interval and a prediction interval for each response at the optimum. + Coefficients alone locate an optimum but cannot say how well it is known, so the + fitted statsmodels results objects are passed alongside the coefficient dicts. + Behaviour is unchanged when `fitted_results` is omitted. +- The `"desirability"` result now carries a `"responses"` list, pairing each model's + coefficients with its specification limits, so it can be handed directly to + `visualize_doe(plot_type="overlay")` without hand-assembling a dict. +- `process_improve.experiments._desirability`, one shared implementation of the + Derringer-Suich functions, replacing two copies that had drifted apart. Target goals + now accept a separate `weight_high` for the falling side. +- `optimize_responses(..., search_bounds=...)` sets the coded region to search, and the + region against which a stationary point is judged inside or outside. It accepts one + `(low, high)` pair for every factor, or a mapping for per-factor control. The previous + behaviour was hard-coded to the factorial cube, which suits a two-level design but + understates a central composite design, whose axial runs sit at plus or minus alpha: + the optimizer could not consider settings the experiment had actually covered, and + `inside_design_space` could report a legitimately explored point as outside. The + default is unchanged at `(-1, 1)`. + +### Changed + +- **The overlay plot now computes and shades the sweet spot**: the region where every + bounded response is simultaneously within specification. Each bounded response draws + its contours at exactly its two specification limits rather than eight arbitrary + levels. New metadata `sweet_spot_fraction` and `sweet_spot_empty` report whether any + compromise exists, distinguishing an empty region from one where no bounds were given. +- The contour renderer honours the style keys it was already being passed: `colorscale`, + `zmin`, `zmax`, `showscale`, `ncontours`, `contours_coloring`, explicit contour + `start`/`end`/`size`, layer opacity, and the layer colour for line-drawn contours. + Defaults match the previous hard-coded values, so existing plots are unchanged. +- `optimize_responses` matches goals to fitted models by response name when both sides + supply one, falling back to list order with a warning. Goals were previously consumed + positionally while `goal["response"]` was documented as the key and never read, so + passing the two lists in different orders optimised the wrong thing without complaint. +- `optimize_responses(..., desirability_weights=)` is renamed to `response_importance`. + The old name is still accepted and now emits a `DeprecationWarning`. The values are + importances, which set how much each response counts in the composite, not the + per-goal `weight` that shapes an individual desirability ramp. + +### Fixed + +- The desirability contour plot's single-response fallback set the bounds to negative and + positive infinity, so the ramp computed `inf/inf` and the whole grid came out NaN. + Bounds are now derived from the range the model spans over the plotted region. +- The overlay plot's `low`/`high` bounds became an `AnnotationType.label`, which no + adapter renders, so the specification limits never reached the reader. +- `OverlayPlot` read only `analysis_results["optimization"]["responses"]`, a shape no + library function produced. It now also accepts the output of `optimize_responses`, and + falls back to a single-response contour when only top-level coefficients are present. ## [1.61.0] - 2026-07-25 ### Fixed @@ -2707,7 +2761,8 @@ this entry records them together. - Reworked the README with a sharper value proposition and a "Why not scikit-learn?" comparison table. -[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.61.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.62.0...HEAD +[1.62.0]: https://github.com/kgdunn/process-improve/compare/v1.61.0...v1.62.0 [1.61.0]: https://github.com/kgdunn/process-improve/compare/v1.60.0...v1.61.0 [1.60.0]: https://github.com/kgdunn/process-improve/compare/v1.59.0...v1.60.0 [1.59.0]: https://github.com/kgdunn/process-improve/compare/v1.58.0...v1.59.0 diff --git a/CITATION.cff b/CITATION.cff index 8990021e..14553329 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,8 +12,8 @@ authors: repository-code: "https://github.com/kgdunn/process-improve" url: "https://kgdunn.github.io/process-improve/" license: MIT -version: 1.61.0 -date-released: "2026-07-25" +version: 1.62.0 +date-released: "2026-07-26" keywords: - chemometrics - multivariate analysis diff --git a/pyproject.toml b/pyproject.toml index cdc0efac..b8b154c9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.61.0" +version = "1.62.0" description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.' readme = "README.md" license = "MIT" diff --git a/src/process_improve/experiments/_desirability.py b/src/process_improve/experiments/_desirability.py new file mode 100644 index 00000000..795ffe4c --- /dev/null +++ b/src/process_improve/experiments/_desirability.py @@ -0,0 +1,214 @@ +# (c) Kevin Dunn, 2010-2026. MIT License. +"""Derringer-Suich desirability functions, shared by the optimizer and the plots. + +Two quantities in this module are easy to confuse, so they are named apart: + +*weight* + The exponent that shapes an *individual* desirability ramp, per response. + A weight of 1 gives a linear ramp from 0 to 1 across the acceptable range. + A weight above 1 concentrates desirability near the good end, so only + values close to the target score well. A weight below 1 flattens the ramp, + so anything inside the range scores nearly as well as the target. + +*importance* + The exponent applied when the individual desirabilities are combined into + the composite. It sets how much one response counts relative to another, + and it has no effect on any individual ramp. + +This module was previously duplicated in ``optimization.py`` and in +``visualization/plots/optimization_plots.py``, and the two copies drifted. +It now holds the union of their behaviour: asymmetric target weights from the +plotting copy, and the strict argument checking from the optimizer copy. + +References +---------- +Derringer, G. and Suich, R. (1980). "Simultaneous Optimization of Several +Response Variables", *Journal of Quality Technology*, **12**, 214-219. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + +__all__ = [ + "composite_desirability", + "desirability_maximize", + "desirability_minimize", + "desirability_target", + "individual_desirability", +] + + +def desirability_maximize(y: float, low: float, high: float, weight: float = 1.0) -> float: + """One-sided desirability for maximisation. + + Parameters + ---------- + y : float + Predicted response value. + low : float + Value at or below which the response is unacceptable, d = 0. + high : float + Value at or above which the response is fully acceptable, d = 1. + weight : float + Exponent shaping the ramp between *low* and *high*. See the module + docstring for how this differs from importance. + + Returns + ------- + float + Desirability in [0, 1]. + """ + if y <= low: + return 0.0 + if y >= high: + return 1.0 + return ((y - low) / (high - low)) ** weight + + +def desirability_minimize(y: float, low: float, high: float, weight: float = 1.0) -> float: + """One-sided desirability for minimisation. + + Parameters + ---------- + y : float + Predicted response value. + low : float + Value at or below which the response is fully acceptable, d = 1. + high : float + Value at or above which the response is unacceptable, d = 0. + weight : float + Exponent shaping the ramp between *low* and *high*. + + Returns + ------- + float + Desirability in [0, 1]. + """ + if y <= low: + return 1.0 + if y >= high: + return 0.0 + return ((high - y) / (high - low)) ** weight + + +def desirability_target( # noqa: PLR0913 + y: float, + low: float, + target: float, + high: float, + weight_low: float = 1.0, + weight_high: float = 1.0, +) -> float: + """Two-sided desirability for a target value. + + Ramps up from *low* to *target*, then back down from *target* to *high*. + Values outside ``[low, high]`` score zero. + + Parameters + ---------- + y : float + Predicted response value. + low, target, high : float + Lower bound, most desirable value, and upper bound. + weight_low : float + Exponent shaping the rising side, below the target. + weight_high : float + Exponent shaping the falling side, above the target. + + Returns + ------- + float + Desirability in [0, 1]. + """ + if y <= low or y >= high: + return 0.0 + if y <= target: + return ((y - low) / (target - low)) ** weight_low + return ((high - y) / (high - target)) ** weight_high + + +def individual_desirability(y: float, goal: dict[str, Any]) -> float: + """Desirability of one predicted response, given its goal. + + Parameters + ---------- + y : float + Predicted response value. + goal : dict + Keys: ``goal`` (one of ``"maximize"``, ``"minimize"``, ``"target"``), + ``low``, ``high``, ``target`` (required for a target goal), ``weight``, + and optionally ``weight_high`` to make a target goal asymmetric. When + ``weight_high`` is absent, ``weight`` applies to both sides. + + Returns + ------- + float + Desirability in [0, 1]. + + Raises + ------ + ValueError + If ``goal`` is missing, unrecognised, or a target goal has no target. + """ + if "goal" not in goal: + msg = "Goal dict is missing the required 'goal' key. Use 'maximize', 'minimize', or 'target'." + raise ValueError(msg) + + goal_type = goal["goal"] + low = goal.get("low", 0.0) + high = goal.get("high", 1.0) + weight = goal.get("weight", 1.0) + + if goal_type == "maximize": + return desirability_maximize(y, low, high, weight) + if goal_type == "minimize": + return desirability_minimize(y, low, high, weight) + if goal_type == "target": + if "target" not in goal: + msg = "A goal of 'target' requires a 'target' value in the goal dict." + raise ValueError(msg) + return desirability_target(y, low, goal["target"], high, weight, goal.get("weight_high", weight)) + + msg = f"Unknown goal type: {goal_type!r}. Use 'maximize', 'minimize', or 'target'." + raise ValueError(msg) + + +def composite_desirability(d_values: list[float], importances: list[float] | None = None) -> float: + """Weighted geometric mean of the individual desirabilities. + + ``D = (d1^w1 * d2^w2 * ... * dk^wk) ^ (1 / sum(wi))`` + + A single zero drives the composite to zero, whatever the other responses + do. That is the defining property of the geometric mean here: a setting + that fails one specification outright cannot be rescued by excelling + elsewhere. + + Parameters + ---------- + d_values : list of float + Individual desirabilities, each in [0, 1]. + importances : list of float, optional + Relative importance per response. Defaults to equal importance. + + Returns + ------- + float + Composite desirability in [0, 1]. Zero for an empty input, for any + zero individual value, or when the importances sum to zero. + """ + if not d_values: + return 0.0 + + weights = importances if importances is not None else [1.0] * len(d_values) + + if any(d == 0.0 for d in d_values): + return 0.0 + + log_d = sum(w * np.log(d) for d, w in zip(d_values, weights, strict=True)) + w_sum = sum(weights) + if w_sum == 0: + return 0.0 + return float(np.exp(log_d / w_sum)) diff --git a/src/process_improve/experiments/_tools/optimize_responses.py b/src/process_improve/experiments/_tools/optimize_responses.py index 71a03682..745ee0c3 100644 --- a/src/process_improve/experiments/_tools/optimize_responses.py +++ b/src/process_improve/experiments/_tools/optimize_responses.py @@ -28,9 +28,12 @@ class OptimizeResponsesInput(BaseModel): goals: list[dict[str, Any]] | None = Field( None, description=( - "Per-response optimisation goals. Each entry: response (str), " + "Per-response optimisation goals. Each entry: response (str, matched against the model's " + "response_name so the two lists need not be in the same order), " "goal ('maximize'|'minimize'|'target'), low, high, optional target, " - "weight, importance. Required for 'desirability' method." + "weight, weight_high, importance. Required for 'desirability' method. " + "'weight' shapes that response's own desirability ramp; 'importance' sets how much the " + "response counts relative to the others in the composite. They are different things." ), ) method: Literal[ @@ -59,12 +62,47 @@ class OptimizeResponsesInput(BaseModel): ge=1, description="Number of steps for steepest ascent/descent (default 10).", ) + response_importance: list[float] | None = Field( + None, + description=( + "Relative importance per response in the composite desirability, aligned with fitted_models. " + "Overrides the per-goal 'importance'. This is not the per-goal 'weight', which shapes an " + "individual response's ramp." + ), + ) + significance_level: float = Field( + 0.05, + gt=0.0, + lt=1.0, + description="Alpha for intervals reported at the optimum (default 0.05, giving 95% intervals).", + ) + search_bounds: list[float] | dict[str, list[float]] | None = Field( + None, + description=( + "Coded region to search, as [low, high] applied to every factor, or a mapping from factor " + "name to its own [low, high]. Defaults to the factorial cube, [-1, 1]. That default suits a " + "two-level design but understates a central composite design, whose axial runs sit at plus " + "or minus alpha: pass [-1.41, 1.41] for a two-factor rotatable central composite design so " + "the search covers the region the experiment actually explored." + ), + ) desirability_weights: list[float] | None = Field( None, - description="Importance weights for composite desirability (overrides per-goal importance).", + description="Deprecated alias for 'response_importance'. Use 'response_importance' instead.", ) +def _as_bounds( + raw: list[float] | dict[str, list[float]] | None, +) -> tuple[float, float] | dict[str, tuple[float, float]] | None: + """Convert the JSON-friendly list form into the tuples the library expects.""" + if raw is None: + return None + if isinstance(raw, dict): + return {name: (float(pair[0]), float(pair[1])) for name, pair in raw.items()} + return (float(raw[0]), float(raw[1])) + + @tool_spec( name="optimize_responses", description=( @@ -75,7 +113,10 @@ class OptimizeResponsesInput(BaseModel): "'canonical_analysis' (eigenvalue decomposition to classify the response surface shape). " "Each fitted_model must include coefficients (as returned by analyze_experiment with " "analysis_type='coefficients'), factor_names, and response_name. " - "For desirability, each goal specifies whether to maximize, minimize, or target a value." + "For desirability, each goal specifies whether to maximize, minimize, or target a value. " + "The desirability result also carries a 'responses' list, pairing each model's coefficients " + "with its specification limits, which can be passed straight to visualize_doe(plot_type='overlay') " + "to see the region where every response is simultaneously within specification." ), input_model=OptimizeResponsesInput, examples=""" @@ -94,6 +135,12 @@ class OptimizeResponsesInput(BaseModel): {"response": "cost", "goal": "minimize", "low": 10, "high": 40}], method="desirability")`` + # "Optimize two responses, counting yield twice as heavily as cost" + -> ``optimize_responses(fitted_models=[model1, model2], + goals=[{"response": "yield", "goal": "maximize", "low": 30, "high": 50}, + {"response": "cost", "goal": "minimize", "low": 10, "high": 40}], + method="desirability", response_importance=[2.0, 1.0])`` + # "Generate a steepest ascent path from a first-order model" -> ``optimize_responses(fitted_models=[model], method="steepest_ascent", step_size=0.5, n_steps=8, @@ -113,6 +160,9 @@ def optimize_responses_tool(spec: OptimizeResponsesInput) -> dict[str, Any]: factor_ranges=spec.factor_ranges, step_size=spec.step_size, n_steps=spec.n_steps, + response_importance=spec.response_importance, + significance_level=spec.significance_level, + search_bounds=_as_bounds(spec.search_bounds), desirability_weights=spec.desirability_weights, ) return clean(result) diff --git a/src/process_improve/experiments/optimization.py b/src/process_improve/experiments/optimization.py index f60d58b4..5767d2fc 100644 --- a/src/process_improve/experiments/optimization.py +++ b/src/process_improve/experiments/optimization.py @@ -26,12 +26,17 @@ import logging import re -from collections.abc import Callable +import warnings +from collections.abc import Callable, Sequence from typing import Any import numpy as np +import pandas as pd +from patsy import PatsyError from scipy import optimize +from process_improve.experiments._desirability import composite_desirability, individual_desirability + logger = logging.getLogger(__name__) # --------------------------------------------------------------------------- @@ -217,6 +222,7 @@ def _find_stationary_point( coefficients: list[dict[str, Any]], factor_names: list[str], factor_ranges: dict[str, dict[str, float]] | None = None, + search_bounds: tuple[float, float] | dict[str, tuple[float, float]] | None = None, ) -> dict[str, Any]: """Find the stationary point of a second-order response surface model. @@ -262,8 +268,13 @@ def _find_stationary_point( else: classification = "saddle_point" - # Check if stationary point is inside the design space (coded [-1, 1]) - inside_design_space = bool(np.all(np.abs(x_s) <= 1.0)) + # Is the stationary point inside the region the experiment covered? The + # default region is the factorial cube; a central composite design reaches + # further, so its axial distance can be supplied via search_bounds. + region = _resolve_search_bounds(search_bounds, factor_names) + inside_design_space = bool( + all(low <= value <= high for value, (low, high) in zip(x_s, region, strict=True)) + ) result: dict[str, Any] = { "stationary_point_coded": {n: float(x_s[i]) for i, n in enumerate(factor_names)}, @@ -436,97 +447,126 @@ def _steepest_path( # noqa: PLR0913 } -# --------------------------------------------------------------------------- -# Derringer-Suich desirability functions -# --------------------------------------------------------------------------- +def _resolve_search_bounds( + search_bounds: tuple[float, float] | dict[str, tuple[float, float]] | None, + factor_names: list[str], +) -> list[tuple[float, float]]: + """Return per-factor coded bounds for the region to search. + The default of (-1, 1) is the factorial cube, which is the right region for + a two-level design. It is not the right region for a central composite + design, whose axial runs sit at plus or minus alpha: restricting the search + to the cube there would refuse to consider settings the experiment actually + covered. Pass the design's axial distance to search the whole region. -def _desirability_maximize(y: float, low: float, high: float, weight: float = 1.0) -> float: - """One-sided desirability for maximisation. + Parameters + ---------- + search_bounds : tuple, dict, or None + A single ``(low, high)`` pair applied to every factor, or a mapping from + factor name to its own pair. Factors absent from the mapping fall back to + (-1, 1). ``None`` means (-1, 1) throughout. + factor_names : list[str] + Ordered factor names. - d = 0 if y <= low, d = ((y - low) / (high - low))^weight if - low < y < high, d = 1 if y >= high. + Returns + ------- + list[tuple[float, float]] + One ``(low, high)`` pair per factor, in *factor_names* order. + + Raises + ------ + ValueError + If a pair is malformed, non-finite, or has low >= high, or if the + mapping names a factor the model does not have. """ - if y <= low: - return 0.0 - if y >= high: - return 1.0 - return ((y - low) / (high - low)) ** weight - - -def _desirability_minimize(y: float, low: float, high: float, weight: float = 1.0) -> float: - """One-sided desirability for minimisation. + default = (-1.0, 1.0) + + def _check(pair: Sequence[float], where: str) -> tuple[float, float]: + try: + low, high = (float(pair[0]), float(pair[1])) + except (TypeError, ValueError, IndexError, KeyError) as exc: + msg = f"search_bounds{where} must be a (low, high) pair of numbers; got {pair!r}." + raise ValueError(msg) from exc + if not (np.isfinite(low) and np.isfinite(high)): + msg = f"search_bounds{where} must be finite; got ({low}, {high})." + raise ValueError(msg) + if low >= high: + msg = f"search_bounds{where} must have low < high; got ({low}, {high})." + raise ValueError(msg) + return low, high - d = 1 if y <= low, d = ((high - y) / (high - low))^weight if - low < y < high, d = 0 if y >= high. - """ - if y <= low: - return 1.0 - if y >= high: - return 0.0 - return ((high - y) / (high - low)) ** weight - - -def _desirability_target( # noqa: PLR0913 - y: float, - low: float, - target: float, - high: float, - weight_low: float = 1.0, - weight_high: float = 1.0, -) -> float: - """Two-sided desirability for a target value. + if search_bounds is None: + return [default] * len(factor_names) - Ramps up from *low* to *target*, then ramps down from *target* to - *high*. Returns 0 outside ``[low, high]``. - """ - if y <= low or y >= high: - return 0.0 - if y <= target: - return ((y - low) / (target - low)) ** weight_low - return ((high - y) / (high - target)) ** weight_high + if isinstance(search_bounds, dict): + unknown = set(search_bounds) - set(factor_names) + if unknown: + msg = f"search_bounds names unknown factor(s) {sorted(unknown)}; the model has {factor_names}." + raise ValueError(msg) + return [ + _check(search_bounds[name], f"[{name!r}]") if name in search_bounds else default + for name in factor_names + ] + return [_check(search_bounds, "")] * len(factor_names) -def _individual_desirability(y: float, goal: dict[str, Any]) -> float: - """Compute individual desirability for a single response value.""" - goal_type = goal["goal"] - low = goal.get("low", 0.0) - high = goal.get("high", 1.0) - weight = goal.get("weight", 1.0) - if goal_type == "maximize": - return _desirability_maximize(y, low, high, weight) - if goal_type == "minimize": - return _desirability_minimize(y, low, high, weight) - if goal_type == "target": - target = goal["target"] - return _desirability_target(y, low, target, high, weight, weight) +def _align_goals_to_models( + fitted_models: list[dict[str, Any]], + goals: list[dict[str, Any]], +) -> list[dict[str, Any]]: + """Return *goals* reordered to match *fitted_models*. - msg = f"Unknown goal type: {goal_type!r}. Use 'maximize', 'minimize', or 'target'." - raise ValueError(msg) + Goals were previously consumed in list order while ``goal["response"]`` was + documented as the key that ties a goal to its model. Passing the two lists + in different orders therefore optimised the wrong thing without complaint. + When every model names its response and every goal names a matching one, the + goals are reordered by name. Otherwise the original positional order is kept, + with a warning, since that is the only interpretation left. -def _composite_desirability(d_values: list[float], importances: list[float] | None = None) -> float: - """Weighted geometric mean of individual desirability values. + Parameters + ---------- + fitted_models : list[dict] + Each optionally has ``"response_name"``. + goals : list[dict] + Each optionally has ``"response"``. - D = (d1^w1 * d2^w2 * ... * dk^wk) ^ (1 / sum(wi)) + Returns + ------- + list[dict] + Goals in the same order as *fitted_models*. - If any d_i is zero, the composite is zero. + Raises + ------ + ValueError + If the two lists differ in length. """ - if not d_values: - return 0.0 + if len(goals) != len(fitted_models): + msg = f"Got {len(fitted_models)} fitted model(s) but {len(goals)} goal(s); they must correspond one to one." + raise ValueError(msg) - weights = importances or [1.0] * len(d_values) + model_names = [m.get("response_name") for m in fitted_models] + goal_names = [g.get("response") for g in goals] - # If any desirability is zero, composite is zero - if any(d == 0.0 for d in d_values): - return 0.0 + if any(n is None for n in model_names) or any(n is None for n in goal_names): + logger.warning( + "Matching goals to fitted models by position: not every model has 'response_name' and not every " + "goal has 'response'. Name both to have them matched by name instead." + ) + return goals + + by_name = {str(g["response"]): g for g in goals} + if len(by_name) != len(goals) or set(by_name) != {str(n) for n in model_names}: + logger.warning( + "Matching goals to fitted models by position: the goal 'response' names %s do not correspond " + "one to one with the model 'response_name' values %s.", + sorted(str(n) for n in goal_names), + sorted(str(n) for n in model_names), + ) + return goals - log_d = sum(w * np.log(d) for d, w in zip(d_values, weights, strict=True)) - w_sum = sum(weights) - if w_sum == 0: - return 0.0 - return float(np.exp(log_d / w_sum)) + return [by_name[str(n)] for n in model_names] def _optimize_desirability( # noqa: PLR0913 @@ -536,6 +576,7 @@ def _optimize_desirability( # noqa: PLR0913 factor_ranges: dict[str, dict[str, float]] | None = None, importances: list[float] | None = None, random_state: int | np.random.Generator | None = 42, + search_bounds: tuple[float, float] | dict[str, tuple[float, float]] | None = None, ) -> dict[str, Any]: """Optimise composite desirability using scipy SLSQP. @@ -544,13 +585,17 @@ def _optimize_desirability( # noqa: PLR0913 fitted_models : list[dict] Each has ``"coefficients"`` and ``"response_name"``. goals : list[dict] - Per-response goals. + Per-response goals. Matched to *fitted_models* by response name when + both sides supply one, otherwise by position. factor_names : list[str] Ordered factor names. factor_ranges : dict or None Factor bounds in actual units. importances : list[float] or None - Relative importance weights for composite desirability. + Relative importance of each response in the composite. This is not the + same as a goal's ``weight``, which shapes that response's own ramp. + search_bounds : tuple, dict, or None + Coded region to search. Defaults to the factorial cube, (-1, 1). Returns ------- @@ -558,6 +603,7 @@ def _optimize_desirability( # noqa: PLR0913 Optimal settings, predicted responses, individual and composite desirability. """ + goals = _align_goals_to_models(fitted_models, goals) evaluators = [_build_model_evaluator(m["coefficients"], factor_names) for m in fitted_models] def neg_composite(x: np.ndarray) -> float: @@ -565,12 +611,13 @@ def neg_composite(x: np.ndarray) -> float: d_vals = [] for evaluator, goal in zip(evaluators, goals, strict=True): y_pred = evaluator(x) - d = _individual_desirability(y_pred, goal) + d = individual_desirability(y_pred, goal) d_vals.append(d) - return -_composite_desirability(d_vals, importances) + return -composite_desirability(d_vals, importances) - k = len(factor_names) - bounds = [(-1.0, 1.0)] * k + bounds = _resolve_search_bounds(search_bounds, factor_names) + lows = np.array([b[0] for b in bounds]) + highs = np.array([b[1] for b in bounds]) # Multi-start: try centre + random points. # SEC-33 (#282): the hard-coded ``42`` moved to the public signature @@ -582,7 +629,10 @@ def neg_composite(x: np.ndarray) -> float: best_result = None best_value = np.inf - starting_points = [np.zeros(k), *[rng.uniform(-1, 1, size=k) for _ in range(9)]] + # Start from the centre of the searched region, then sample across it, so + # that widening the bounds actually widens where the search looks. + centre = (lows + highs) / 2.0 + starting_points = [centre, *[rng.uniform(lows, highs) for _ in range(9)]] for x0 in starting_points: res = optimize.minimize(neg_composite, x0, method="SLSQP", bounds=bounds) @@ -604,7 +654,7 @@ def neg_composite(x: np.ndarray) -> float: resp_name = model_dict.get("response_name", "response") y_pred = float(evaluator(x_opt)) predictions[resp_name] = y_pred - individual_d[resp_name] = _individual_desirability(y_pred, goal) + individual_d[resp_name] = individual_desirability(y_pred, goal) result: dict[str, Any] = { "optimal_coded": {n: float(x_opt[i]) for i, n in enumerate(factor_names)}, @@ -704,6 +754,121 @@ def _coded_to_actual(coded: dict[str, float], factor_ranges: dict[str, dict[str, # --------------------------------------------------------------------------- +def _intervals_at_point( + fitted_results: list[Any], + fitted_models: list[dict[str, Any]], + factor_names: list[str], + point_coded: dict[str, float], + significance_level: float, +) -> dict[str, Any]: + """Confidence and prediction intervals for each response at one point. + + The optimizer works from coefficients alone, which is enough to locate an + optimum but not to say how well it is known. The residual variance and the + design's leverage at that point are needed for that, and both live on the + fitted model object rather than in its coefficients. + + Parameters + ---------- + fitted_results : list + Statsmodels results objects, aligned with *fitted_models*, fitted on + the coded factors. + fitted_models : list[dict] + Used only for the response names. + factor_names : list[str] + Ordered factor names, matching the columns the models were fitted on. + point_coded : dict[str, float] + Coded factor settings at which to report the intervals. + significance_level : float + Alpha. 0.05 gives 95% intervals. + + Returns + ------- + dict + Keyed by response name. Each entry has ``predicted``, + ``confidence_interval``, ``prediction_interval``, and + ``confidence_level``. A response whose model cannot be evaluated + carries an ``error`` string instead, so one failure does not discard + the intervals for the others. + """ + from process_improve.experiments._analyses.prediction import _run_prediction # noqa: PLC0415 + + if len(fitted_results) != len(fitted_models): + msg = ( + f"Got {len(fitted_models)} fitted model(s) but {len(fitted_results)} fitted result(s); " + "they must correspond one to one and be in the same order." + ) + raise ValueError(msg) + + new_point = pd.DataFrame([{name: point_coded[name] for name in factor_names}]) + + intervals: dict[str, Any] = {} + for i, (results_obj, model) in enumerate(zip(fitted_results, fitted_models, strict=True)): + resp_name = model.get("response_name", f"Response {i + 1}") + try: + record = _run_prediction(results_obj, new_point, alpha=significance_level)["predictions"][0] + except (AttributeError, KeyError, TypeError, ValueError, PatsyError) as exc: + logger.warning("Could not compute intervals for response %r: %s", resp_name, exc) + intervals[resp_name] = {"error": str(exc)} + continue + + intervals[resp_name] = { + "predicted": record["predicted"], + "confidence_interval": [record["ci_low"], record["ci_high"]], + "prediction_interval": [record["pi_low"], record["pi_high"]], + "confidence_level": 1.0 - significance_level, + } + return intervals + + +def _desirability_result( # noqa: PLR0913 + *, + fitted_models: list[dict[str, Any]], + goals: list[dict[str, Any]], + factor_names: list[str], + factor_ranges: dict[str, dict[str, float]] | None, + response_importance: list[float] | None, + fitted_results: list[Any] | None, + significance_level: float, + search_bounds: tuple[float, float] | dict[str, tuple[float, float]] | None = None, +) -> dict[str, Any]: + """Assemble the full desirability result: optimum, intervals, and plot input. + + Returns + ------- + dict + The optimum from :func:`_optimize_desirability`, plus + ``"response_intervals"`` when *fitted_results* is supplied, plus + ``"responses"``, which pairs each model's coefficients with its + specification limits so the result can be passed straight to the + overlay plot. + """ + aligned_goals = _align_goals_to_models(fitted_models, goals) + importances = response_importance + if importances is None: + importances = [g.get("importance", 1.0) for g in aligned_goals] + + desirability = _optimize_desirability( + fitted_models, aligned_goals, factor_names, factor_ranges, importances, search_bounds=search_bounds + ) + + if fitted_results is not None: + desirability["response_intervals"] = _intervals_at_point( + fitted_results, fitted_models, factor_names, desirability["optimal_coded"], significance_level + ) + + carried = ("goal", "low", "high", "target", "weight", "weight_high", "importance") + desirability["responses"] = [ + { + "name": model.get("response_name", f"Response {i + 1}"), + "coefficients": model.get("coefficients", []), + **{key: goal[key] for key in carried if key in goal}, + } + for i, (model, goal) in enumerate(zip(fitted_models, aligned_goals, strict=True)) + ] + return desirability + + def optimize_responses( # noqa: PLR0913, C901 fitted_models: list[dict[str, Any]], goals: list[dict[str, Any]] | None = None, @@ -711,6 +876,10 @@ def optimize_responses( # noqa: PLR0913, C901 factor_ranges: dict[str, dict[str, float]] | None = None, step_size: float = 0.5, n_steps: int = 10, + response_importance: list[float] | None = None, + fitted_results: list[Any] | None = None, + significance_level: float = 0.05, + search_bounds: tuple[float, float] | dict[str, tuple[float, float]] | None = None, desirability_weights: list[float] | None = None, ) -> dict[str, Any]: """Find optimal factor settings for one or multiple responses. @@ -731,16 +900,24 @@ def optimize_responses( # noqa: PLR0913, C901 goals : list[dict] or None Per-response optimisation goals. Each dict has keys: - - ``"response"`` (str) - response name (must match a model). + - ``"response"`` (str) - response name. Matched against each model's + ``"response_name"``; when both sides name their responses the goals + are reordered to match, so the two lists need not be in the same + order. When either side omits a name, goals are taken in list order. - ``"goal"`` (str) - ``"maximize"``, ``"minimize"``, or ``"target"``. - ``"target"`` (float, optional) - target value (required when ``goal="target"``). - ``"low"`` (float) - lower acceptable bound. - ``"high"`` (float) - upper acceptable bound. - - ``"weight"`` (float, default 1) - desirability shape parameter. - - ``"importance"`` (float, default 1) - relative importance for - composite desirability. + - ``"weight"`` (float, default 1) - the exponent shaping *this* + response's desirability ramp between ``low`` and ``high``. Above 1 + concentrates desirability near the good end; below 1 flattens it. + - ``"weight_high"`` (float, optional) - a separate exponent for the + falling side of a ``"target"`` goal. Defaults to ``"weight"``. + - ``"importance"`` (float, default 1) - how much this response counts + relative to the others when the composite is formed. Unlike + ``weight``, it has no effect on this response's own ramp. method : str Optimisation method: ``"desirability"``, @@ -754,9 +931,34 @@ def optimize_responses( # noqa: PLR0913, C901 Step magnitude for steepest ascent/descent (coded units). n_steps : int Number of steps for steepest ascent/descent. + response_importance : list[float] or None + Relative importance per response, overriding the per-goal + ``"importance"`` values. Aligned with *fitted_models*. + fitted_results : list or None + Optional statsmodels results objects, one per entry in *fitted_models* + and in the same order, as returned by ``lm()`` or by + ``analyze_experiment``. When supplied, a confidence interval and a + prediction interval for each response are reported at the optimum. + The models must have been fitted on the coded factors, since the + optimum is located in coded units. + significance_level : float + Alpha for those intervals. The default of 0.05 gives 95% intervals. + search_bounds : tuple, dict, or None + The coded region to search, and the region against which a stationary + point is judged inside or outside. Defaults to the factorial cube, + ``(-1, 1)`` on every factor. + + That default suits a two-level design but understates a central + composite design, whose axial runs sit at plus or minus alpha: leaving + it at the cube would refuse to consider settings the experiment + actually covered. Pass ``(-1.41, 1.41)`` for a two-factor rotatable + central composite design, or a mapping such as + ``{"T": (-1.41, 1.41)}`` to widen one factor only. Factors left out of + a mapping keep the (-1, 1) default. desirability_weights : list[float] or None - Importance weights for composite desirability (overrides per-goal - ``"importance"`` values). + Deprecated alias for *response_importance*. The name was misleading: + these values are importances, not the ``weight`` that shapes an + individual ramp. Returns ------- @@ -764,6 +966,13 @@ def optimize_responses( # noqa: PLR0913, C901 Results keyed by method. Always includes ``"method"`` and ``"factor_names"``. + Raises + ------ + ValueError + If *method* is unknown, if *fitted_models* is empty, if a method that + needs goals is called without them, or if both *response_importance* + and *desirability_weights* are given. + Examples -------- >>> from process_improve.experiments.optimization import optimize_responses @@ -796,6 +1005,22 @@ def optimize_responses( # noqa: PLR0913, C901 msg = "At least one fitted model is required." raise ValueError(msg) + if desirability_weights is not None: + if response_importance is not None: + msg = ( + "Pass either 'response_importance' or the deprecated 'desirability_weights', not both. " + "They set the same thing: how much each response counts in the composite." + ) + raise ValueError(msg) + warnings.warn( + "'desirability_weights' is deprecated; use 'response_importance'. The values are importances, " + "which set how much each response counts in the composite, not the per-goal 'weight' that shapes " + "an individual desirability ramp.", + DeprecationWarning, + stacklevel=2, + ) + response_importance = desirability_weights + # Use factor_names from the first model as the canonical ordering factor_names = fitted_models[0]["factor_names"] coefficients = fitted_models[0]["coefficients"] @@ -803,12 +1028,16 @@ def optimize_responses( # noqa: PLR0913, C901 result: dict[str, Any] = {"method": method, "factor_names": factor_names} if method == "stationary_point": - result["stationary_point"] = _find_stationary_point(coefficients, factor_names, factor_ranges) + result["stationary_point"] = _find_stationary_point( + coefficients, factor_names, factor_ranges, search_bounds + ) elif method == "canonical_analysis": result["canonical_analysis"] = _canonical_analysis(coefficients, factor_names) # Also include the stationary point for context - result["stationary_point"] = _find_stationary_point(coefficients, factor_names, factor_ranges) + result["stationary_point"] = _find_stationary_point( + coefficients, factor_names, factor_ranges, search_bounds + ) elif method in ("steepest_ascent", "steepest_descent"): direction = "ascent" if method == "steepest_ascent" else "descent" @@ -820,13 +1049,15 @@ def optimize_responses( # noqa: PLR0913, C901 if goals is None: msg = "Goals are required for desirability optimization." raise ValueError(msg) - - importances = desirability_weights - if importances is None: - importances = [g.get("importance", 1.0) for g in goals] - - result["desirability"] = _optimize_desirability( - fitted_models, goals, factor_names, factor_ranges, importances + result["desirability"] = _desirability_result( + fitted_models=fitted_models, + goals=goals, + factor_names=factor_names, + factor_ranges=factor_ranges, + response_importance=response_importance, + fitted_results=fitted_results, + significance_level=significance_level, + search_bounds=search_bounds, ) elif method == "ridge_analysis": diff --git a/src/process_improve/experiments/visualization/plots/optimization_plots.py b/src/process_improve/experiments/visualization/plots/optimization_plots.py index c60a2b5d..23715b01 100644 --- a/src/process_improve/experiments/visualization/plots/optimization_plots.py +++ b/src/process_improve/experiments/visualization/plots/optimization_plots.py @@ -11,94 +11,122 @@ import numpy as np +from process_improve.experiments._desirability import composite_desirability, individual_desirability from process_improve.experiments.visualization.plots.registry import BasePlot, register_plot from process_improve.experiments.visualization.plots.surfaces import _build_coef_map, _compute_grid, _evaluate_model from process_improve.visualization.colors import ( DESIRABILITY_COLORSCALE, DOE_PALETTE, FACTOR_COLORS, + SWEET_SPOT_COLORSCALE, ) from process_improve.visualization.spec import ( - Annotation, ChartSpec, Encoding, LayerSpec, PanelSpec, ) -from process_improve.visualization.types import AnnotationType, MarkType +from process_improve.visualization.types import MarkType # --------------------------------------------------------------------------- -# Shared desirability helpers +# Grid helpers shared by the desirability and overlay plots # --------------------------------------------------------------------------- -def _desirability_maximize(y: float, low: float, high: float, weight: float = 1.0) -> float: - """Individual desirability for a maximise goal.""" - if y <= low: - return 0.0 - if y >= high: - return 1.0 - return ((y - low) / (high - low)) ** weight - - -def _desirability_minimize(y: float, low: float, high: float, weight: float = 1.0) -> float: - """Individual desirability for a minimise goal.""" - if y <= low: - return 1.0 - if y >= high: - return 0.0 - return ((high - y) / (high - low)) ** weight - - -def _desirability_target( # noqa: PLR0913 - y: float, - low: float, - target: float, - high: float, - weight_low: float = 1.0, - weight_high: float = 1.0, -) -> float: - """Individual desirability for a target goal.""" - if y <= low or y >= high: - return 0.0 - if y <= target: - return ((y - low) / (target - low)) ** weight_low - return ((high - y) / (high - target)) ** weight_high - - -def _individual_desirability(y: float, goal: dict[str, Any]) -> float: - """Compute individual desirability for a single response.""" - goal_type = goal.get("goal", "maximize") - weight = goal.get("weight", 1.0) - - if goal_type == "maximize": - return _desirability_maximize(y, goal.get("low", 0.0), goal.get("high", 1.0), weight) - if goal_type == "minimize": - return _desirability_minimize(y, goal.get("low", 0.0), goal.get("high", 1.0), weight) - if goal_type == "target": - return _desirability_target( - y, - goal.get("low", 0.0), - goal.get("target", 0.5), - goal.get("high", 1.0), - weight, - goal.get("weight_high", weight), - ) - return 0.0 +def _response_surface( + coef_map: dict[str, float], + factor_x: str, + factor_y: str, + hold_values: dict[str, float], + n_grid: int, +) -> list[list[float]]: + """Predict one response over the coded grid, row-major as ``z[y][x]``.""" + _, _, z_matrix = _compute_grid(coef_map, factor_x, factor_y, hold_values, n_grid) + return z_matrix + + +def _is_bounded(low: float | None, high: float | None) -> bool: + """Report whether both specification limits are present, finite, and ordered.""" + if low is None or high is None: + return False + return bool(np.isfinite(low) and np.isfinite(high) and high > low) + + +def _feasibility_mask( + surfaces: list[tuple[dict[str, Any], list[list[float]]]], + n_grid: int, +) -> tuple[list[list[float]] | None, int]: + """Mark the grid cells where every bounded response is within specification. + + Parameters + ---------- + surfaces : list of (goal dict, predicted surface) + One entry per response that could be evaluated. + n_grid : int + Grid resolution, assumed square. + + Returns + ------- + mask : list[list[float]] or None + ``1.0`` inside the sweet spot and ``nan`` outside, so the renderer draws + nothing where the region does not apply. ``None`` when no response + carries usable bounds, in which case there is no sweet spot to speak of + rather than an empty one. + n_bounded : int + How many responses contributed bounds. + """ + bounded = [(resp, z) for resp, z in surfaces if _is_bounded(resp.get("low"), resp.get("high"))] + if not bounded: + return None, 0 + + mask: list[list[float]] = [] + for j in range(n_grid): + row: list[float] = [] + for i in range(n_grid): + inside = all(resp["low"] <= z[j][i] <= resp["high"] for resp, z in bounded) + row.append(1.0 if inside else float("nan")) + mask.append(row) + return mask, len(bounded) + + +def _resolve_goal_bounds(goal: dict[str, Any], surface: list[list[float]]) -> dict[str, Any]: + """Return *goal* with finite ``low`` and ``high``, derived if not supplied. + + A response with no stated specification has nothing to ramp between. Rather + than defaulting to an arbitrary [0, 1] or to infinities (which make the + desirability undefined), the ramp spans the range the fitted model actually + covers over the plotted region. + + Parameters + ---------- + goal : dict + Goal dict, possibly missing ``low`` / ``high`` or carrying non-finite + values for them. + surface : list[list[float]] + Predicted response over the grid. + + Returns + ------- + dict + A copy of *goal* with finite bounds. If the surface is flat, the bounds + are nudged apart so the ramp stays well defined. + """ + resolved = dict(goal) + low = resolved.get("low") + high = resolved.get("high") + if low is not None and high is not None and np.isfinite(low) and np.isfinite(high): + return resolved + flat = [v for row in surface for v in row if np.isfinite(v)] + if not flat: + resolved["low"], resolved["high"] = 0.0, 1.0 + return resolved -def _composite_desirability(d_values: list[float], importances: list[float] | None = None) -> float: - """Weighted geometric mean of individual desirabilities.""" - if not d_values: - return 0.0 - if any(d == 0.0 for d in d_values): - return 0.0 - weights = importances or [1.0] * len(d_values) - total_w = sum(weights) - product = 1.0 - for d, w in zip(d_values, weights): # noqa: B905 - product *= d ** (w / total_w) - return product + z_min, z_max = float(min(flat)), float(max(flat)) + if z_min == z_max: + z_min, z_max = z_min - 0.5, z_max + 0.5 + resolved["low"], resolved["high"] = z_min, z_max + return resolved # --------------------------------------------------------------------------- @@ -145,28 +173,26 @@ def to_spec(self) -> ChartSpec: importances = [r.get("importance", 1.0) for r in responses] - # Build coefficient maps for each response - coef_maps = [] - for resp in responses: - coeffs = resp.get("coefficients", []) - coef_maps.append(_build_coef_map(coeffs)) + # Predict each response over the grid first, so that a response with no + # stated specification can have its ramp derived from the range the + # model actually spans here (see _resolve_goal_bounds). + surfaces = [ + _response_surface( + _build_coef_map(resp.get("coefficients", [])), factor_x, factor_y, self.hold_values, n_grid + ) + for resp in responses + ] + goals = [_resolve_goal_bounds(resp, surface) for resp, surface in zip(responses, surfaces, strict=True)] - # Evaluate composite desirability over grid z_matrix: list[list[float]] = [] - for y_val in y_grid: + for j in range(n_grid): row: list[float] = [] - for x_val in x_grid: - point = dict(self.hold_values) - point[factor_x] = x_val - point[factor_y] = y_val - - d_values = [] - for cm, resp in zip(coef_maps, responses): # noqa: B905 - y_hat = _evaluate_model(cm, point) - d = _individual_desirability(y_hat, resp) - d_values.append(d) - - row.append(_composite_desirability(d_values, importances)) + for i in range(n_grid): + d_values = [ + individual_desirability(surface[j][i], goal) + for surface, goal in zip(surfaces, goals, strict=True) + ] + row.append(composite_desirability(d_values, importances)) z_matrix.append(row) contour_layer = LayerSpec( @@ -212,16 +238,12 @@ def _get_optimization_responses(self) -> list[dict[str, Any]]: if responses: return responses - # Fall back: single-response from top-level coefficients + # Fall back: single-response from top-level coefficients. Bounds are + # left unset deliberately; _resolve_goal_bounds derives them from the + # predicted range over the plotted region. coefficients = self._get_coefficients() if coefficients: - return [{ - "coefficients": coefficients, - "goal": "maximize", - "low": float("-inf"), - "high": float("inf"), - "weight": 1.0, - }] + return [{"coefficients": coefficients, "goal": "maximize", "weight": 1.0}] return [] @@ -232,17 +254,25 @@ def _get_optimization_responses(self) -> list[dict[str, Any]]: @register_plot("overlay") class OverlayPlot(BasePlot): - """Overlay contour plot for multiple responses. + """Overlay contour plot for multiple responses, showing the sweet spot. + + Each response model is drawn as contour lines on shared axes. Where a + response carries ``low`` / ``high`` bounds, its two specification limits are + drawn as its contour levels, so the reader sees the boundary that matters + rather than eight arbitrary levels. - Draws contour lines from each response model on the same axes, - with optional constraint regions shaded. Useful for identifying - a feasible operating region that satisfies all specifications. + The region where *every* bounded response is simultaneously inside its + limits is shaded. That region is the sweet spot: the set of operating points + that meet all the specifications at once. It can be empty, which is a + meaningful answer and is reported as such rather than drawn as blank space. Data sources ------------ Requires ``analysis_results`` with ``"optimization"`` containing - ``"responses"`` - each with ``"coefficients"``, and optionally - ``"low"``/``"high"`` bounds for feasibility shading. + ``"responses"``, each with ``"coefficients"`` and optionally ``"low"`` / + ``"high"``. The output of + :func:`~process_improve.experiments.optimization.optimize_responses` is also + accepted directly. """ def to_spec(self) -> ChartSpec: @@ -266,53 +296,76 @@ def to_spec(self) -> ChartSpec: y_grid = np.linspace(-1, 1, n_grid).tolist() layers: list[LayerSpec] = [] + surfaces: list[tuple[dict[str, Any], list[list[float]]]] = [] for i, resp in enumerate(responses): coeffs = resp.get("coefficients", []) if not coeffs: continue - coef_map = _build_coef_map(coeffs) - _, _, z_matrix = _compute_grid( - coef_map, factor_x, factor_y, self.hold_values, n_grid, - ) + z_matrix = _response_surface(_build_coef_map(coeffs), factor_x, factor_y, self.hold_values, n_grid) + surfaces.append((resp, z_matrix)) resp_name = resp.get("name", f"Response {i + 1}") color = FACTOR_COLORS[i % len(FACTOR_COLORS)] + style: dict[str, Any] = { + "x_grid": x_grid, + "y_grid": y_grid, + "z_matrix": z_matrix, + "contours_coloring": "lines", + "showscale": False, + } + + if _is_bounded(resp.get("low"), resp.get("high")): + # Draw exactly the two specification limits. `size` is the gap + # between levels, so this yields contours at low and at high. + spec_low, spec_high = float(resp["low"]), float(resp["high"]) + style["contours"] = {"start": spec_low, "end": spec_high, "size": spec_high - spec_low} + else: + style["ncontours"] = 8 + + layers.append( + LayerSpec( + mark=MarkType.contour, + data=[], + x=Encoding(field="x", title=factor_x), + y=Encoding(field="y", title=factor_y), + name=resp_name, + color=color, + style=style, + ) + ) - layer = LayerSpec( - mark=MarkType.contour, - data=[], - x=Encoding(field="x", title=factor_x), - y=Encoding(field="y", title=factor_y), - name=resp_name, - color=color, - style={ - "x_grid": x_grid, - "y_grid": y_grid, - "z_matrix": z_matrix, - "contours_coloring": "lines", - "ncontours": 8, - }, + mask, n_bounded = _feasibility_mask(surfaces, n_grid) + sweet_spot_cells = 0 if mask is None else sum(1 for row in mask for v in row if v == 1.0) + sweet_spot_fraction = sweet_spot_cells / float(n_grid * n_grid) if mask is not None else None + + if mask is not None: + # Drawn first so it sits behind the response contours. + layers.insert( + 0, + LayerSpec( + mark=MarkType.contour, + data=[], + x=Encoding(field="x", title=factor_x), + y=Encoding(field="y", title=factor_y), + name="Sweet spot", + opacity=0.35, + style={ + "x_grid": x_grid, + "y_grid": y_grid, + "z_matrix": mask, + "colorscale": SWEET_SPOT_COLORSCALE, + "zmin": 0.0, + "zmax": 1.0, + "showscale": False, + "showlabels": False, + }, + ), ) - layers.append(layer) - - annotations: list[Annotation] = [] - # Add constraint region annotations if bounds are specified - for resp in responses: - low = resp.get("low") - high = resp.get("high") - if low is not None and high is not None: - resp_name = resp.get("name", "Response") - annotations.append(Annotation( - annotation_type=AnnotationType.label, - label=f"{resp_name}: [{low:.2f}, {high:.2f}]", - style={"color": DOE_PALETTE["neutral"]}, - )) panel = PanelSpec( layers=layers, - annotations=annotations, title=f"Overlay: {factor_x} x {factor_y}", x_title=factor_x, y_title=factor_y, @@ -326,15 +379,44 @@ def to_spec(self) -> ChartSpec: "factors": [factor_x, factor_y], "hold_values": self.hold_values, "n_responses": len(responses), + "n_bounded_responses": n_bounded, + "specification_limits": { + resp.get("name", f"Response {i + 1}"): [resp.get("low"), resp.get("high")] + for i, resp in enumerate(responses) + if _is_bounded(resp.get("low"), resp.get("high")) + }, + "sweet_spot_fraction": sweet_spot_fraction, + "sweet_spot_empty": None if mask is None else sweet_spot_cells == 0, }, ) def _get_overlay_responses(self) -> list[dict[str, Any]]: - """Extract multi-response data for overlay.""" + """Extract multi-response data for overlay. + + Accepts either the hand-assembled ``{"optimization": {"responses": ...}}`` + shape or the result of ``optimize_responses``, which nests the same + information under ``"desirability"``. + """ if not self.analysis_results: return [] + opt = self.analysis_results.get("optimization", {}) - return opt.get("responses", []) + responses = opt.get("responses", []) + if responses: + return responses + + # optimize_responses(...) output, passed through directly. + desirability = self.analysis_results.get("desirability", {}) + responses = desirability.get("responses", []) + if responses: + return responses + + # Fall back: a single response built from top-level coefficients, so + # the plot degrades to a plain contour rather than an empty frame. + coefficients = self._get_coefficients() + if coefficients: + return [{"name": "Response 1", "coefficients": coefficients}] + return [] # --------------------------------------------------------------------------- diff --git a/src/process_improve/visualization/adapters/plotly_adapter.py b/src/process_improve/visualization/adapters/plotly_adapter.py index 8512205b..e1b9f92c 100644 --- a/src/process_improve/visualization/adapters/plotly_adapter.py +++ b/src/process_improve/visualization/adapters/plotly_adapter.py @@ -286,14 +286,41 @@ def _contour_trace(self, layer: LayerSpec) -> go.Contour: x_vals = layer.style.get("x_grid", []) y_vals = layer.style.get("y_grid", []) z_matrix = layer.style.get("z_matrix", []) - return go.Contour( - x=x_vals, - y=y_vals, - z=z_matrix, - name=layer.name, - colorscale=_SURFACE_COLORSCALE, - contours=dict(showlabels=True), - ) + + # Contour levels. `contours` may carry explicit start/end/size, which is + # how a caller pins the levels to specification limits rather than + # letting plotly choose them. + contours: dict[str, Any] = dict(layer.style.get("contours", {})) + contours.setdefault("showlabels", layer.style.get("showlabels", True)) + coloring = layer.style.get("contours_coloring") + if coloring is not None: + contours["coloring"] = coloring + + trace_kwargs: dict[str, Any] = { + "x": x_vals, + "y": y_vals, + "z": z_matrix, + "name": layer.name, + "colorscale": layer.style.get("colorscale", _SURFACE_COLORSCALE), + "contours": contours, + "opacity": layer.opacity, + "showscale": layer.style.get("showscale", True), + } + + # Only forward the optional keys when set: passing zmin/zmax of None + # would override plotly's own auto-ranging. + for key in ("zmin", "zmax", "ncontours"): + value = layer.style.get(key) + if value is not None: + trace_kwargs[key] = value + + # With line-drawn contours the layer colour selects the line colour, so + # several responses overlaid on one panel stay distinguishable. + if layer.color and contours.get("coloring") == "lines": + trace_kwargs["line"] = dict(color=layer.color) + trace_kwargs["showscale"] = False + + return go.Contour(**trace_kwargs) def _surface_trace(self, layer: LayerSpec) -> go.Surface: x_vals = layer.style.get("x_grid", []) diff --git a/src/process_improve/visualization/colors.py b/src/process_improve/visualization/colors.py index e8d0f18f..ff8a8959 100644 --- a/src/process_improve/visualization/colors.py +++ b/src/process_improve/visualization/colors.py @@ -91,6 +91,17 @@ [1.0, "#059669"], # d = 1 (best) ] +# --------------------------------------------------------------------------- +# Sweet-spot fill: a single flat green, used behind overlaid response contours +# to shade the region where every response is within specification. It is one +# colour rather than a scale because the mask it renders is binary. +# --------------------------------------------------------------------------- + +SWEET_SPOT_COLORSCALE: list[list[float | str]] = [ + [0.0, "#059669"], + [1.0, "#059669"], +] + # --------------------------------------------------------------------------- # Font stacks (shared by the registered Plotly base themes) # --------------------------------------------------------------------------- diff --git a/tests/test_experiments_tools.py b/tests/test_experiments_tools.py index 88ed5229..bc3ea11a 100644 --- a/tests/test_experiments_tools.py +++ b/tests/test_experiments_tools.py @@ -346,6 +346,50 @@ def test_stationary_point_quadratic(self) -> None: ) assert "error" not in result + @staticmethod + def _rising_plane_call(**extra: object) -> dict: + """Maximize a plane rising with A, so the optimum sits at the upper bound.""" + return execute_tool_call( + "optimize_responses", + { + "fitted_models": [ + { + "response_name": "y", + "factor_names": ["A", "B"], + "coefficients": [ + {"term": "Intercept", "coefficient": 0.0}, + {"term": "A", "coefficient": 1.0}, + ], + } + ], + "goals": [{"response": "y", "goal": "maximize", "low": -5.0, "high": 5.0}], + "method": "desirability", + **extra, + }, + ) + + def test_search_bounds_as_a_single_pair(self) -> None: + """Over the wire the region arrives as a JSON list, not a tuple.""" + result = self._rising_plane_call(search_bounds=[-1.41, 1.41]) + assert "error" not in result + assert result["desirability"]["optimal_coded"]["A"] == pytest.approx(1.41, abs=1e-4) + + def test_search_bounds_per_factor(self) -> None: + """A mapping widens one factor and leaves the other at the cube.""" + result = self._rising_plane_call(search_bounds={"A": [-2.0, 2.0]}) + assert "error" not in result + assert result["desirability"]["optimal_coded"]["A"] == pytest.approx(2.0, abs=1e-4) + + def test_search_bounds_default_is_the_cube(self) -> None: + """Omitting it leaves the previous behaviour in place.""" + result = self._rising_plane_call() + assert result["desirability"]["optimal_coded"]["A"] == pytest.approx(1.0, abs=1e-4) + + def test_malformed_search_bounds_returns_an_error(self) -> None: + """A reversed pair is reported, not silently accepted.""" + result = self._rising_plane_call(search_bounds=[1.0, -1.0]) + assert "low < high" in result["error"] + def test_invalid_method_returns_error(self) -> None: """Unknown method is rejected by the pydantic Literal.""" from process_improve.tool_safety import ToolInputInvalidError diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 85ba791f..5c41b723 100644 --- a/tests/test_optimization.py +++ b/tests/test_optimization.py @@ -5,16 +5,18 @@ import numpy as np import pytest +from process_improve.experiments._desirability import ( + composite_desirability, + desirability_maximize, + desirability_minimize, + desirability_target, + individual_desirability, +) from process_improve.experiments.optimization import ( _build_model_evaluator, _canonical_analysis, - _composite_desirability, - _desirability_maximize, - _desirability_minimize, - _desirability_target, _extract_b_and_B, _find_stationary_point, - _individual_desirability, _parse_term, _steepest_path, evaluate_model, @@ -371,20 +373,20 @@ class TestDesirabilityMaximize: def test_below_low(self) -> None: """Value below low bound gives d=0.""" - assert _desirability_maximize(5.0, 10.0, 20.0) == 0.0 + assert desirability_maximize(5.0, 10.0, 20.0) == 0.0 def test_above_high(self) -> None: """Value above high bound gives d=1.""" - assert _desirability_maximize(25.0, 10.0, 20.0) == 1.0 + assert desirability_maximize(25.0, 10.0, 20.0) == 1.0 def test_at_midpoint(self) -> None: """Midpoint gives d=0.5 with linear weight.""" - assert _desirability_maximize(15.0, 10.0, 20.0) == pytest.approx(0.5) + assert desirability_maximize(15.0, 10.0, 20.0) == pytest.approx(0.5) def test_weight_effect(self) -> None: """Weight < 1 (concave) gives higher d at midpoint than linear.""" - d_linear = _desirability_maximize(15.0, 10.0, 20.0, weight=1.0) - d_concave = _desirability_maximize(15.0, 10.0, 20.0, weight=0.5) + d_linear = desirability_maximize(15.0, 10.0, 20.0, weight=1.0) + d_concave = desirability_maximize(15.0, 10.0, 20.0, weight=0.5) assert d_concave > d_linear @@ -393,15 +395,15 @@ class TestDesirabilityMinimize: def test_below_low(self) -> None: """Value below low bound gives d=1.""" - assert _desirability_minimize(5.0, 10.0, 20.0) == 1.0 + assert desirability_minimize(5.0, 10.0, 20.0) == 1.0 def test_above_high(self) -> None: """Value above high bound gives d=0.""" - assert _desirability_minimize(25.0, 10.0, 20.0) == 0.0 + assert desirability_minimize(25.0, 10.0, 20.0) == 0.0 def test_at_midpoint(self) -> None: """Midpoint gives d=0.5 with linear weight.""" - assert _desirability_minimize(15.0, 10.0, 20.0) == pytest.approx(0.5) + assert desirability_minimize(15.0, 10.0, 20.0) == pytest.approx(0.5) class TestDesirabilityTarget: @@ -409,50 +411,50 @@ class TestDesirabilityTarget: def test_at_target(self) -> None: """At target value, d=1.""" - assert _desirability_target(15.0, 10.0, 15.0, 20.0) == pytest.approx(1.0) + assert desirability_target(15.0, 10.0, 15.0, 20.0) == pytest.approx(1.0) def test_below_low(self) -> None: """Below low bound, d=0.""" - assert _desirability_target(5.0, 10.0, 15.0, 20.0) == 0.0 + assert desirability_target(5.0, 10.0, 15.0, 20.0) == 0.0 def test_above_high(self) -> None: """Above high bound, d=0.""" - assert _desirability_target(25.0, 10.0, 15.0, 20.0) == 0.0 + assert desirability_target(25.0, 10.0, 15.0, 20.0) == 0.0 def test_between_low_and_target(self) -> None: """Between low and target, 0 < d < 1.""" - d = _desirability_target(12.5, 10.0, 15.0, 20.0) + d = desirability_target(12.5, 10.0, 15.0, 20.0) assert 0.0 < d < 1.0 def test_between_target_and_high(self) -> None: """Between target and high, 0 < d < 1.""" - d = _desirability_target(17.5, 10.0, 15.0, 20.0) + d = desirability_target(17.5, 10.0, 15.0, 20.0) assert 0.0 < d < 1.0 class TestIndividualDesirability: - """Verify _individual_desirability dispatch to correct function.""" + """Verify individual_desirability dispatch to correct function.""" def test_maximize_goal(self) -> None: """Maximize goal above high gives d=1.""" goal = {"goal": "maximize", "low": 10.0, "high": 20.0} - assert _individual_desirability(25.0, goal) == 1.0 + assert individual_desirability(25.0, goal) == 1.0 def test_minimize_goal(self) -> None: """Minimize goal below low gives d=1.""" goal = {"goal": "minimize", "low": 10.0, "high": 20.0} - assert _individual_desirability(5.0, goal) == 1.0 + assert individual_desirability(5.0, goal) == 1.0 def test_target_goal(self) -> None: """Target goal at target gives d=1.""" goal = {"goal": "target", "low": 10.0, "high": 20.0, "target": 15.0} - assert _individual_desirability(15.0, goal) == pytest.approx(1.0) + assert individual_desirability(15.0, goal) == pytest.approx(1.0) def test_unknown_goal_raises(self) -> None: """Unknown goal type raises ValueError.""" goal = {"goal": "unknown", "low": 10.0, "high": 20.0} with pytest.raises(ValueError, match="Unknown goal"): - _individual_desirability(15.0, goal) + individual_desirability(15.0, goal) class TestCompositeDesirability: @@ -460,26 +462,26 @@ class TestCompositeDesirability: def test_all_ones(self) -> None: """All d=1 gives composite D=1.""" - assert _composite_desirability([1.0, 1.0, 1.0]) == pytest.approx(1.0) + assert composite_desirability([1.0, 1.0, 1.0]) == pytest.approx(1.0) def test_any_zero_gives_zero(self) -> None: """Any d=0 makes composite D=0.""" - assert _composite_desirability([1.0, 0.0, 1.0]) == 0.0 + assert composite_desirability([1.0, 0.0, 1.0]) == 0.0 def test_geometric_mean(self) -> None: """Unweighted: D = sqrt(0.5 * 0.8) = sqrt(0.4).""" - d = _composite_desirability([0.5, 0.8]) + d = composite_desirability([0.5, 0.8]) assert d == pytest.approx(np.sqrt(0.4)) def test_weighted(self) -> None: """Weighted geometric mean with importances [2, 1].""" - d = _composite_desirability([0.5, 0.8], importances=[2.0, 1.0]) + d = composite_desirability([0.5, 0.8], importances=[2.0, 1.0]) expected = np.exp((2.0 * np.log(0.5) + 1.0 * np.log(0.8)) / 3.0) assert d == pytest.approx(expected) def test_empty_list(self) -> None: """Empty list returns 0.""" - assert _composite_desirability([]) == 0.0 + assert composite_desirability([]) == 0.0 # --------------------------------------------------------------------------- @@ -583,6 +585,356 @@ def test_random_state_is_configurable(self) -> None: assert out_c["composite_desirability"] > 0.5 +# --------------------------------------------------------------------------- +# Where the optimum actually lands +# --------------------------------------------------------------------------- + + +class TestDesirabilityOptimumLocation: + """Pin down where the optimiser lands, not just that it produced a number. + + The other desirability tests assert only that keys are present and that the + composite is above zero, which a badly wrong optimum would also satisfy. + """ + + @staticmethod + def _plane(name: str, intercept: float, slope_a: float, slope_b: float) -> dict: + """Build a plane in A and B, so the optimum is known without solving anything.""" + return { + "response_name": name, + "coefficients": [ + {"term": "Intercept", "coefficient": intercept}, + {"term": "A", "coefficient": slope_a}, + {"term": "B", "coefficient": slope_b}, + ], + "factor_names": ["A", "B"], + } + + def test_single_response_lands_on_the_known_corner(self) -> None: + """Maximising a plane drives both factors to the corner that maximises it.""" + model = self._plane("y", intercept=0.0, slope_a=1.0, slope_b=-1.0) + goals = [{"response": "y", "goal": "maximize", "low": -2.0, "high": 2.0}] + out = optimize_responses([model], goals=goals, method="desirability")["desirability"] + assert out["optimal_coded"]["A"] == pytest.approx(1.0, abs=1e-4) + assert out["optimal_coded"]["B"] == pytest.approx(-1.0, abs=1e-4) + assert out["predicted_responses"]["y"] == pytest.approx(2.0, abs=1e-4) + assert out["composite_desirability"] == pytest.approx(1.0, abs=1e-4) + + def test_two_responses_compromise_between_their_optima(self) -> None: + """Conflicting responses settle strictly between their individual optima. + + y1 wants A at +1, y2 wants A at -1, and both are indifferent to B. The + compromise must therefore sit strictly inside the A range. + """ + y1 = self._plane("y1", intercept=0.0, slope_a=1.0, slope_b=0.0) + y2 = self._plane("y2", intercept=0.0, slope_a=-1.0, slope_b=0.0) + goals = [ + {"response": "y1", "goal": "maximize", "low": -1.0, "high": 1.0}, + {"response": "y2", "goal": "maximize", "low": -1.0, "high": 1.0}, + ] + out = optimize_responses([y1, y2], goals=goals, method="desirability")["desirability"] + assert out["optimal_coded"]["A"] == pytest.approx(0.0, abs=1e-3) + + def test_importance_pulls_the_optimum_toward_the_favoured_response(self) -> None: + """Raising one response's importance moves the compromise its way.""" + y1 = self._plane("y1", intercept=0.0, slope_a=1.0, slope_b=0.0) + y2 = self._plane("y2", intercept=0.0, slope_a=-1.0, slope_b=0.0) + goals = [ + {"response": "y1", "goal": "maximize", "low": -1.0, "high": 1.0}, + {"response": "y2", "goal": "maximize", "low": -1.0, "high": 1.0}, + ] + balanced = optimize_responses([y1, y2], goals=goals, method="desirability")["desirability"] + favoured = optimize_responses( + [y1, y2], goals=goals, method="desirability", response_importance=[5.0, 1.0] + )["desirability"] + assert favoured["optimal_coded"]["A"] > balanced["optimal_coded"]["A"] + + +class TestGoalMatching: + """Goals should follow their response name, not their list position.""" + + @staticmethod + def _models() -> list[dict]: + return [ + { + "response_name": "yield", + "coefficients": [{"term": "Intercept", "coefficient": 0.0}, {"term": "A", "coefficient": 1.0}], + "factor_names": ["A", "B"], + }, + { + "response_name": "cost", + "coefficients": [{"term": "Intercept", "coefficient": 0.0}, {"term": "A", "coefficient": -1.0}], + "factor_names": ["A", "B"], + }, + ] + + def test_goal_order_does_not_change_the_answer(self) -> None: + """Reordering goals relative to models used to silently invert the problem.""" + yield_goal = {"response": "yield", "goal": "maximize", "low": -1.0, "high": 1.0} + cost_goal = {"response": "cost", "goal": "minimize", "low": -1.0, "high": 1.0} + + in_order = optimize_responses( + self._models(), goals=[yield_goal, cost_goal], method="desirability" + )["desirability"] + reversed_order = optimize_responses( + self._models(), goals=[cost_goal, yield_goal], method="desirability" + )["desirability"] + + assert in_order["optimal_coded"]["A"] == pytest.approx(reversed_order["optimal_coded"]["A"], abs=1e-6) + # Both goals push A to +1: yield rises with A, and cost falls with A. + assert in_order["optimal_coded"]["A"] == pytest.approx(1.0, abs=1e-4) + + def test_mismatched_length_is_rejected(self) -> None: + """One goal per model, or the pairing is undefined.""" + goals = [{"response": "yield", "goal": "maximize", "low": -1.0, "high": 1.0}] + with pytest.raises(ValueError, match="correspond one to one"): + optimize_responses(self._models(), goals=goals, method="desirability") + + def test_unnamed_goals_fall_back_to_position(self, caplog: pytest.LogCaptureFixture) -> None: + """Without names on both sides, position is the only reading available.""" + goals = [ + {"goal": "maximize", "low": -1.0, "high": 1.0}, + {"goal": "minimize", "low": -1.0, "high": 1.0}, + ] + with caplog.at_level("WARNING"): + out = optimize_responses(self._models(), goals=goals, method="desirability") + assert "by position" in caplog.text + assert out["desirability"]["optimal_coded"]["A"] == pytest.approx(1.0, abs=1e-4) + + +class TestResponseImportanceNaming: + """The old kwarg name said 'weights' but carried importances.""" + + @staticmethod + def _model() -> dict: + return { + "response_name": "y", + "coefficients": [{"term": "Intercept", "coefficient": 0.0}, {"term": "A", "coefficient": 1.0}], + "factor_names": ["A", "B"], + } + + def test_deprecated_alias_still_works(self) -> None: + """desirability_weights keeps working, with a warning.""" + goals = [{"response": "y", "goal": "maximize", "low": -1.0, "high": 1.0}] + with pytest.warns(DeprecationWarning, match="response_importance"): + out = optimize_responses( + [self._model()], goals=goals, method="desirability", desirability_weights=[1.0] + ) + assert out["desirability"]["composite_desirability"] > 0.0 + + def test_both_names_together_is_an_error(self) -> None: + """Passing both leaves the intent ambiguous.""" + goals = [{"response": "y", "goal": "maximize", "low": -1.0, "high": 1.0}] + with pytest.raises(ValueError, match="not both"): + optimize_responses( + [self._model()], + goals=goals, + method="desirability", + response_importance=[1.0], + desirability_weights=[2.0], + ) + + def test_result_carries_responses_for_the_overlay_plot(self) -> None: + """The desirability result is directly consumable by the overlay plot.""" + goals = [{"response": "y", "goal": "maximize", "low": -1.0, "high": 1.0}] + out = optimize_responses([self._model()], goals=goals, method="desirability")["desirability"] + assert out["responses"][0]["name"] == "y" + assert out["responses"][0]["low"] == -1.0 + assert out["responses"][0]["high"] == 1.0 + assert out["responses"][0]["coefficients"] + + +class TestIntervalsAtOptimum: + """Uncertainty at the optimum, when the fitted model objects are supplied.""" + + @staticmethod + def _fit() -> tuple[dict, object]: + """Fit a small two-factor model on coded factors and return both forms.""" + import pandas as pd + import statsmodels.formula.api as smf + + design = pd.DataFrame({ + "A": [-1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0], + "B": [-1.0, -1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0], + }) + design["y"] = 40.0 + 5.0 * design["A"] - 2.0 * design["B"] + [0.3, -0.2, 0.1, -0.1, 0.2, -0.3, 0.1, 0.0] + ols_result = smf.ols("y ~ A + B", data=design).fit() + model = { + "response_name": "y", + "coefficients": [ + {"term": term, "coefficient": float(value)} for term, value in ols_result.params.items() + ], + "factor_names": ["A", "B"], + } + return model, ols_result + + def test_no_intervals_without_fitted_results(self) -> None: + """Behaviour is unchanged when the fitted objects are not supplied.""" + model, _ = self._fit() + goals = [{"response": "y", "goal": "maximize", "low": 30.0, "high": 50.0}] + out = optimize_responses([model], goals=goals, method="desirability")["desirability"] + assert "response_intervals" not in out + + def test_intervals_are_reported_and_ordered(self) -> None: + """The prediction interval contains the confidence interval, which contains the fit.""" + model, fitted = self._fit() + goals = [{"response": "y", "goal": "maximize", "low": 30.0, "high": 50.0}] + out = optimize_responses( + [model], goals=goals, method="desirability", fitted_results=[fitted] + )["desirability"] + + interval = out["response_intervals"]["y"] + ci_low, ci_high = interval["confidence_interval"] + pi_low, pi_high = interval["prediction_interval"] + predicted = interval["predicted"] + + assert ci_low < predicted < ci_high + assert pi_low < ci_low + assert pi_high > ci_high + assert interval["confidence_level"] == pytest.approx(0.95) + # The optimizer and the fitted model must agree on the predicted value. + assert predicted == pytest.approx(out["predicted_responses"]["y"], abs=1e-6) + + def test_significance_level_widens_the_interval(self) -> None: + """A smaller alpha gives a wider interval.""" + model, fitted = self._fit() + goals = [{"response": "y", "goal": "maximize", "low": 30.0, "high": 50.0}] + + def width(alpha: float) -> float: + out = optimize_responses( + [model], + goals=goals, + method="desirability", + fitted_results=[fitted], + significance_level=alpha, + )["desirability"] + low, high = out["response_intervals"]["y"]["confidence_interval"] + return high - low + + assert width(0.01) > width(0.05) > width(0.20) + + def test_mismatched_fitted_results_length_is_rejected(self) -> None: + """One fitted result per model, in the same order.""" + model, fitted = self._fit() + goals = [{"response": "y", "goal": "maximize", "low": 30.0, "high": 50.0}] + with pytest.raises(ValueError, match="correspond one to one"): + optimize_responses( + [model], goals=goals, method="desirability", fitted_results=[fitted, fitted] + ) + + +class TestSearchBounds: + """The searched region defaults to the cube, but need not be the cube. + + A two-level design covers the factorial cube, so (-1, 1) is right for it. + A central composite design reaches further: its axial runs sit at plus or + minus alpha. Searching only the cube there refuses to consider settings the + experiment actually covered. + """ + + @staticmethod + def _rising_plane() -> dict: + """Return a plane rising with A, so the optimum sits at the upper bound.""" + return { + "response_name": "y", + "coefficients": [{"term": "Intercept", "coefficient": 0.0}, {"term": "A", "coefficient": 1.0}], + "factor_names": ["A", "B"], + } + + @staticmethod + def _goals() -> list[dict]: + return [{"response": "y", "goal": "maximize", "low": -2.0, "high": 2.0}] + + def test_default_is_the_factorial_cube(self) -> None: + """Unchanged behaviour when nothing is passed.""" + out = optimize_responses([self._rising_plane()], goals=self._goals(), method="desirability") + assert out["desirability"]["optimal_coded"]["A"] == pytest.approx(1.0, abs=1e-4) + + def test_widening_the_region_moves_the_optimum_out(self) -> None: + """A central composite design's axial reach is searchable.""" + out = optimize_responses( + [self._rising_plane()], + goals=self._goals(), + method="desirability", + search_bounds=(-1.41, 1.41), + ) + assert out["desirability"]["optimal_coded"]["A"] == pytest.approx(1.41, abs=1e-4) + + def test_per_factor_bounds(self) -> None: + """One factor can be widened without widening the others.""" + model = { + "response_name": "y", + "coefficients": [ + {"term": "Intercept", "coefficient": 0.0}, + {"term": "A", "coefficient": 1.0}, + {"term": "B", "coefficient": 1.0}, + ], + "factor_names": ["A", "B"], + } + # The ramp is deliberately wider than the region can reach, so the + # desirability never saturates and the optimum stays unique. + goals = [{"response": "y", "goal": "maximize", "low": -5.0, "high": 5.0}] + out = optimize_responses( + [model], goals=goals, method="desirability", search_bounds={"A": (-1.41, 1.41)} + ) + coded = out["desirability"]["optimal_coded"] + assert coded["A"] == pytest.approx(1.41, abs=1e-4) + assert coded["B"] == pytest.approx(1.0, abs=1e-4) + + def test_stationary_point_region_test_respects_the_bounds(self) -> None: + """A point outside the cube can still be inside a composite design's region. + + The quadratic below has its maximum at A = 1.2, which is outside the + factorial cube but well within the axial reach of a rotatable + two-factor central composite design. + """ + model = { + "response_name": "y", + "coefficients": [ + {"term": "Intercept", "coefficient": 0.0}, + {"term": "A", "coefficient": 2.4}, + {"term": "B", "coefficient": 0.0}, + {"term": "I(A ** 2)", "coefficient": -1.0}, + {"term": "I(B ** 2)", "coefficient": -1.0}, + ], + "factor_names": ["A", "B"], + } + cube = optimize_responses([model], method="stationary_point")["stationary_point"] + assert cube["stationary_point_coded"]["A"] == pytest.approx(1.2, abs=1e-6) + assert cube["inside_design_space"] is False + + composite = optimize_responses( + [model], method="stationary_point", search_bounds=(-1.41, 1.41) + )["stationary_point"] + assert composite["inside_design_space"] is True + + @pytest.mark.parametrize( + ("bad", "match"), + [ + ((1.0, -1.0), "low < high"), + ((0.0, 0.0), "low < high"), + ((float("-inf"), 1.0), "finite"), + ((1.0,), "pair of numbers"), + ("wide", "pair of numbers"), + ], + ) + def test_malformed_bounds_are_rejected(self, bad: object, match: str) -> None: + with pytest.raises(ValueError, match=match): + optimize_responses( + [self._rising_plane()], goals=self._goals(), method="desirability", search_bounds=bad + ) + + def test_unknown_factor_in_bounds_is_rejected(self) -> None: + """A typo in a factor name would otherwise be silently ignored.""" + with pytest.raises(ValueError, match="unknown factor"): + optimize_responses( + [self._rising_plane()], + goals=self._goals(), + method="desirability", + search_bounds={"Temperature": (-2.0, 2.0)}, + ) + + # --------------------------------------------------------------------------- # Stubs # --------------------------------------------------------------------------- diff --git a/tests/test_visualization.py b/tests/test_visualization.py index 10f3980b..5e5506da 100644 --- a/tests/test_visualization.py +++ b/tests/test_visualization.py @@ -13,14 +13,14 @@ import pytest -from process_improve.experiments.visualization import main_effects_plot, visualize_doe -from process_improve.experiments.visualization.plots.optimization_plots import ( - _composite_desirability, - _desirability_maximize, - _desirability_minimize, - _desirability_target, - _individual_desirability, +from process_improve.experiments._desirability import ( + composite_desirability, + desirability_maximize, + desirability_minimize, + desirability_target, + individual_desirability, ) +from process_improve.experiments.visualization import main_effects_plot, visualize_doe from process_improve.experiments.visualization.plots.registry import ( create_plot, get_available_plot_types, @@ -1096,42 +1096,55 @@ class TestDesirabilityHelpers: """Unit tests for the pure desirability functions in optimization_plots.""" def test_maximize_boundaries(self) -> None: - assert _desirability_maximize(0.0, low=1.0, high=2.0) == 0.0 - assert _desirability_maximize(3.0, low=1.0, high=2.0) == 1.0 - assert _desirability_maximize(1.5, low=1.0, high=2.0) == pytest.approx(0.5) + assert desirability_maximize(0.0, low=1.0, high=2.0) == 0.0 + assert desirability_maximize(3.0, low=1.0, high=2.0) == 1.0 + assert desirability_maximize(1.5, low=1.0, high=2.0) == pytest.approx(0.5) def test_minimize_boundaries(self) -> None: - assert _desirability_minimize(0.5, low=1.0, high=2.0) == 1.0 - assert _desirability_minimize(2.5, low=1.0, high=2.0) == 0.0 - assert _desirability_minimize(1.5, low=1.0, high=2.0) == pytest.approx(0.5) + assert desirability_minimize(0.5, low=1.0, high=2.0) == 1.0 + assert desirability_minimize(2.5, low=1.0, high=2.0) == 0.0 + assert desirability_minimize(1.5, low=1.0, high=2.0) == pytest.approx(0.5) def test_target_all_branches(self) -> None: # Outside [low, high] on either side -> 0 - assert _desirability_target(0.5, low=1.0, target=2.0, high=3.0) == 0.0 - assert _desirability_target(3.5, low=1.0, target=2.0, high=3.0) == 0.0 + assert desirability_target(0.5, low=1.0, target=2.0, high=3.0) == 0.0 + assert desirability_target(3.5, low=1.0, target=2.0, high=3.0) == 0.0 # Below target: rising ramp - assert _desirability_target(1.5, low=1.0, target=2.0, high=3.0) == pytest.approx(0.5) + assert desirability_target(1.5, low=1.0, target=2.0, high=3.0) == pytest.approx(0.5) # Above target: falling ramp - assert _desirability_target(2.5, low=1.0, target=2.0, high=3.0) == pytest.approx(0.5) + assert desirability_target(2.5, low=1.0, target=2.0, high=3.0) == pytest.approx(0.5) def test_individual_desirability_goal_dispatch(self) -> None: - assert _individual_desirability(0.75, {"goal": "maximize", "low": 0.0, "high": 1.0}) == pytest.approx(0.75) - assert _individual_desirability(0.25, {"goal": "minimize", "low": 0.0, "high": 1.0}) == pytest.approx(0.75) + assert individual_desirability(0.75, {"goal": "maximize", "low": 0.0, "high": 1.0}) == pytest.approx(0.75) + assert individual_desirability(0.25, {"goal": "minimize", "low": 0.0, "high": 1.0}) == pytest.approx(0.75) target_goal = {"goal": "target", "low": 0.0, "target": 0.5, "high": 1.0} - assert _individual_desirability(0.25, target_goal) == pytest.approx(0.5) - # Unknown goal type falls through to 0.0 - assert _individual_desirability(0.5, {"goal": "unknown_goal"}) == 0.0 + assert individual_desirability(0.25, target_goal) == pytest.approx(0.5) + + def test_individual_desirability_rejects_unknown_goal(self) -> None: + """An unrecognised goal is an error, not a silent zero. + + The plotting copy of this helper used to return 0.0 here, which made a + typo in a goal name look like an infeasible response. + """ + with pytest.raises(ValueError, match="Unknown goal type"): + individual_desirability(0.5, {"goal": "unknown_goal"}) + + def test_individual_desirability_asymmetric_target_weights(self) -> None: + """A target goal may ramp at different rates on each side.""" + goal = {"goal": "target", "low": 0.0, "target": 0.5, "high": 1.0, "weight": 1.0, "weight_high": 2.0} + assert individual_desirability(0.25, goal) == pytest.approx(0.5) + assert individual_desirability(0.75, goal) == pytest.approx(0.25) def test_composite_empty_list(self) -> None: - assert _composite_desirability([]) == 0.0 + assert composite_desirability([]) == 0.0 def test_composite_any_zero(self) -> None: - assert _composite_desirability([0.9, 0.0, 0.8]) == 0.0 + assert composite_desirability([0.9, 0.0, 0.8]) == 0.0 def test_composite_importance_weighting(self) -> None: """Importance weights should change the geometric mean.""" - unweighted = _composite_desirability([0.25, 1.0]) - weighted = _composite_desirability([0.25, 1.0], importances=[3.0, 1.0]) + unweighted = composite_desirability([0.25, 1.0]) + weighted = composite_desirability([0.25, 1.0], importances=[3.0, 1.0]) assert unweighted == pytest.approx(0.25**0.5) assert weighted == pytest.approx(0.25**0.75) assert weighted != pytest.approx(unweighted) @@ -1170,10 +1183,35 @@ def test_no_analysis_results(self) -> None: spec = plot.to_spec() assert "no response data" in spec.title.lower() - def test_missing_optimization_key(self, coefficients_2f: list) -> None: - plot = create_plot("overlay", analysis_results={"coefficients": coefficients_2f}) + def test_missing_optimization_key_falls_back_to_single_response(self, coefficients_2f: list) -> None: + """Top-level coefficients degrade to a one-response contour. + + This matches the desirability contour plot's behaviour, and is more + useful than an empty frame. There are no bounds, so there is no sweet + spot to shade and the metadata says so. + """ + plot = create_plot( + "overlay", analysis_results={"coefficients": coefficients_2f}, factors_to_plot=["A", "B"] + ) spec = plot.to_spec() - assert "no response data" in spec.title.lower() + assert spec.plot_type == "overlay" + assert len(spec.panels[0].layers) == 1 + assert spec.metadata["n_bounded_responses"] == 0 + assert spec.metadata["sweet_spot_empty"] is None + + def test_accepts_optimize_responses_output(self, coefficients_2f: list) -> None: + """The result of optimize_responses can be passed straight through.""" + plot = create_plot( + "overlay", + analysis_results={ + "desirability": { + "responses": [{"name": "Yield", "coefficients": coefficients_2f, "low": 30.0, "high": 50.0}], + }, + }, + factors_to_plot=["A", "B"], + ) + spec = plot.to_spec() + assert spec.metadata["n_bounded_responses"] == 1 def test_single_factor(self, coefficients_2f: list) -> None: plot = create_plot( @@ -1204,8 +1242,12 @@ def test_response_with_empty_coefficients_skipped(self, coefficients_2f: list) - assert len(spec.panels[0].layers) == 1 assert spec.metadata["n_responses"] == 2 - def test_constraint_labels_from_bounds(self, coefficients_2f: list) -> None: - """Responses carrying low/high bounds produce constraint annotations.""" + def test_bounds_pin_the_contour_levels(self, coefficients_2f: list) -> None: + """A bounded response draws its contours at the specification limits. + + Previously the bounds became a text label that no adapter rendered, so + the limits never reached the reader. + """ plot = create_plot( "overlay", analysis_results={ @@ -1218,9 +1260,72 @@ def test_constraint_labels_from_bounds(self, coefficients_2f: list) -> None: factors_to_plot=["A", "B"], ) spec = plot.to_spec() - anns = spec.panels[0].annotations - assert len(anns) == 1 - assert anns[0].label == "Yield: [30.00, 50.00]" + # Layer 0 is the sweet-spot shading; layer 1 is the response itself. + response_layer = next(layer for layer in spec.panels[0].layers if layer.name == "Yield") + assert response_layer.style["contours"] == {"start": 30.0, "end": 50.0, "size": 20.0} + assert spec.metadata["specification_limits"] == {"Yield": [30.0, 50.0]} + + def test_sweet_spot_matches_a_hand_computed_mask(self) -> None: + """The shaded region is exactly where both responses are in range. + + Two planes are used so the feasible set can be worked out by hand: + y1 = A rises with A alone and is wanted at or above 0, so it is in + specification for A >= 0; y2 = B rises with B alone and is wanted at or + below 0, so it is in specification for B <= 0. The sweet spot is + therefore the quadrant A >= 0 and B <= 0, one quarter of the region. + """ + y1 = [{"term": "Intercept", "coefficient": 0.0}, {"term": "A", "coefficient": 1.0}] + y2 = [{"term": "Intercept", "coefficient": 0.0}, {"term": "B", "coefficient": 1.0}] + plot = create_plot( + "overlay", + analysis_results={ + "optimization": { + "responses": [ + {"name": "y1", "coefficients": y1, "low": 0.0, "high": 10.0}, + {"name": "y2", "coefficients": y2, "low": -10.0, "high": 0.0}, + ], + }, + }, + factors_to_plot=["A", "B"], + ) + spec = plot.to_spec() + mask_layer = spec.panels[0].layers[0] + assert mask_layer.name == "Sweet spot" + + # The 50-point grid spans [-1, 1] inclusive, so 25 of the 50 values are + # <= 0 and 25 are >= 0: a quarter of the cells, plus the shared axes. + mask = mask_layer.style["z_matrix"] + x_grid = mask_layer.style["x_grid"] + y_grid = mask_layer.style["y_grid"] + for j, y_val in enumerate(y_grid): + for i, x_val in enumerate(x_grid): + expected_inside = x_val >= 0.0 and y_val <= 0.0 + is_inside = mask[j][i] == 1.0 + assert is_inside == expected_inside, f"cell ({x_val:.3f}, {y_val:.3f})" + + assert spec.metadata["n_bounded_responses"] == 2 + assert spec.metadata["sweet_spot_empty"] is False + assert spec.metadata["sweet_spot_fraction"] == pytest.approx(0.25, abs=0.02) + + def test_empty_sweet_spot_is_reported(self) -> None: + """Specifications that cannot be met together report an empty region.""" + y1 = [{"term": "Intercept", "coefficient": 0.0}, {"term": "A", "coefficient": 1.0}] + plot = create_plot( + "overlay", + analysis_results={ + "optimization": { + "responses": [ + # A never exceeds +1 over the coded region, so a + # specification of 5 to 10 is unreachable everywhere. + {"name": "y1", "coefficients": y1, "low": 5.0, "high": 10.0}, + ], + }, + }, + factors_to_plot=["A", "B"], + ) + spec = plot.to_spec() + assert spec.metadata["sweet_spot_empty"] is True + assert spec.metadata["sweet_spot_fraction"] == 0.0 class TestRidgeTracePlotEdgeCases: diff --git a/tests/test_visualization_spec.py b/tests/test_visualization_spec.py index 2a78d417..bcc17489 100644 --- a/tests/test_visualization_spec.py +++ b/tests/test_visualization_spec.py @@ -84,3 +84,71 @@ def test_incomplete_reference_band_is_skipped(self) -> None: spec = ChartSpec(panels=[_scatter_panel([band])]) result = PlotlyAdapter().render(spec) assert "data" in result + + +def _contour_panel(style: dict, *, color: str | None = None, opacity: float = 1.0) -> PanelSpec: + layer = LayerSpec( + mark=MarkType.contour, + data=[], + x=Encoding(field="x"), + y=Encoding(field="y"), + name="surface", + color=color, + opacity=opacity, + style={"x_grid": [0.0, 1.0], "y_grid": [0.0, 1.0], "z_matrix": [[0.0, 1.0], [1.0, 2.0]], **style}, + ) + return PanelSpec(layers=[layer]) + + +def _first_trace(panel: PanelSpec) -> dict: + rendered = PlotlyAdapter().render(ChartSpec(panels=[panel])) + return rendered["data"][0] + + +class TestPlotlyContourStyling: + """The contour trace must honour the style keys it is handed. + + These keys were previously accepted into the spec and then discarded by the + adapter, so an overlay of several responses rendered as filled surfaces all + in one colorscale. + """ + + def test_defaults_are_unchanged_when_no_style_given(self) -> None: + """A bare contour keeps the previous appearance.""" + trace = _first_trace(_contour_panel({})) + assert trace["contours"]["showlabels"] is True + assert trace["colorscale"] is not None + assert trace["showscale"] is True + + def test_colorscale_and_z_limits_are_forwarded(self) -> None: + trace = _first_trace( + _contour_panel({"colorscale": [[0.0, "#000000"], [1.0, "#FFFFFF"]], "zmin": 0.0, "zmax": 1.0}) + ) + assert [list(stop) for stop in trace["colorscale"]] == [[0.0, "#000000"], [1.0, "#FFFFFF"]] + assert trace["zmin"] == 0.0 + assert trace["zmax"] == 1.0 + + def test_explicit_contour_levels_are_forwarded(self) -> None: + """Pinning start/end/size is how specification limits become contours.""" + trace = _first_trace(_contour_panel({"contours": {"start": 30.0, "end": 50.0, "size": 20.0}})) + assert trace["contours"]["start"] == 30.0 + assert trace["contours"]["end"] == 50.0 + assert trace["contours"]["size"] == 20.0 + + def test_line_coloring_uses_the_layer_colour(self) -> None: + """Overlaid responses stay distinguishable only if each keeps its colour.""" + trace = _first_trace(_contour_panel({"contours_coloring": "lines"}, color="#123456")) + assert trace["contours"]["coloring"] == "lines" + assert trace["line"]["color"] == "#123456" + # A per-response colour bar would be meaningless, so it is suppressed. + assert trace["showscale"] is False + + def test_showscale_and_opacity_are_forwarded(self) -> None: + trace = _first_trace(_contour_panel({"showscale": False}, opacity=0.35)) + assert trace["showscale"] is False + assert trace["opacity"] == 0.35 + + def test_ncontours_is_forwarded_and_omitted_when_unset(self) -> None: + assert _first_trace(_contour_panel({"ncontours": 8}))["ncontours"] == 8 + # Passing None would override plotly's own auto-ranging, so it is dropped. + assert "ncontours" not in _first_trace(_contour_panel({}))