From fc0a9022145fe64e86f427552496af7e36bb711b Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 24 Jul 2026 22:44:18 +0000 Subject: [PATCH 1/5] Share one desirability implementation between optimizer and plots The Derringer-Suich desirability helpers existed twice, in experiments/optimization.py and in the optimisation plots module, and the two copies had drifted apart. The plotting copy supported an asymmetric weight on a target goal; the optimizer copy did not. The optimizer copy raised on an unrecognised goal and guarded a zero importance sum; the plotting copy silently returned zero and divided by zero respectively. A reader could not tell which behaviour they would get. Move them to experiments/_desirability.py, holding the union: asymmetric target weights, plus the strict argument checking. A typo in a goal name now raises instead of masquerading as an infeasible response. The module docstring separates the two quantities that were easy to confuse: weight shapes an individual ramp, importance sets how much a response counts when the composite is formed. Also fixes the single-response fallback in the desirability contour plot, which set low and high to negative and positive infinity. That made the ramp compute inf/inf, so the whole grid came out NaN and the plot was blank. Bounds are now derived from the range the model spans over the plotted region. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TVbQmQNRwrSerdh7WQX2Aw --- .../experiments/_desirability.py | 214 ++++++++++++++++++ .../experiments/optimization.py | 101 +-------- .../visualization/plots/optimization_plots.py | 164 ++++++-------- tests/test_optimization.py | 58 ++--- tests/test_visualization.py | 65 +++--- 5 files changed, 359 insertions(+), 243 deletions(-) create mode 100644 src/process_improve/experiments/_desirability.py 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/optimization.py b/src/process_improve/experiments/optimization.py index f60d58b4..2245e9fc 100644 --- a/src/process_improve/experiments/optimization.py +++ b/src/process_improve/experiments/optimization.py @@ -32,6 +32,8 @@ import numpy as np from scipy import optimize +from process_improve.experiments._desirability import composite_desirability, individual_desirability + logger = logging.getLogger(__name__) # --------------------------------------------------------------------------- @@ -436,99 +438,6 @@ def _steepest_path( # noqa: PLR0913 } -# --------------------------------------------------------------------------- -# Derringer-Suich desirability functions -# --------------------------------------------------------------------------- - - -def _desirability_maximize(y: float, low: float, high: float, weight: float = 1.0) -> float: - """One-sided desirability for maximisation. - - d = 0 if y <= low, d = ((y - low) / (high - low))^weight if - low < y < high, d = 1 if y >= high. - """ - 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. - - 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. - - 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 - - -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) - - 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 individual desirability values. - - D = (d1^w1 * d2^w2 * ... * dk^wk) ^ (1 / sum(wi)) - - If any d_i is zero, the composite is zero. - """ - if not d_values: - return 0.0 - - weights = importances or [1.0] * len(d_values) - - # If any desirability is zero, composite is zero - 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)) - - def _optimize_desirability( # noqa: PLR0913 fitted_models: list[dict[str, Any]], goals: list[dict[str, Any]], @@ -565,9 +474,9 @@ 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 @@ -604,7 +513,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)}, diff --git a/src/process_improve/experiments/visualization/plots/optimization_plots.py b/src/process_improve/experiments/visualization/plots/optimization_plots.py index c60a2b5d..887540d4 100644 --- a/src/process_improve/experiments/visualization/plots/optimization_plots.py +++ b/src/process_improve/experiments/visualization/plots/optimization_plots.py @@ -11,6 +11,7 @@ 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 ( @@ -28,77 +29,60 @@ from process_improve.visualization.types import AnnotationType, 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 _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 +129,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 +194,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 [] diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 85ba791f..737b1a25 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 # --------------------------------------------------------------------------- diff --git a/tests/test_visualization.py b/tests/test_visualization.py index 10f3980b..d302dda9 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) From 4c8f7842da221f466d1018e5d76697f8e70327c4 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 24 Jul 2026 22:49:08 +0000 Subject: [PATCH 2/5] Compute and shade the sweet spot on the overlay plot The overlay plot's docstring promised "constraint regions shaded" and a "feasible operating region that satisfies all specifications", but nothing computed one. Each response's low/high bounds became an AnnotationType.label, and no adapter has a branch for that annotation type, so the limits were dropped silently at render time. What reached the reader was N filled contour surfaces stacked in a single colorscale. Now the plot computes the region where every bounded response is simultaneously within specification, and shades it. Each bounded response draws its contours at exactly its two specification limits rather than at eight arbitrary levels, so the boundary that matters is the one visible. The region can be empty, which is a real answer to a real question: these specifications cannot be met together. metadata carries sweet_spot_fraction and sweet_spot_empty so a caller can tell without inspecting pixels, and distinguishes "empty" from "no bounds were given" (None). The contour adapter now honours the style keys it was already being handed: 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 no existing plot changes. Without this the sweet-spot layer could not be drawn distinctly from the response contours. _get_overlay_responses now also accepts the output of optimize_responses directly. It previously read only analysis_results["optimization"]["responses"], a shape no library function produces, so every caller had to hand-assemble it. Verified by rendering a two-response case and cross-checking the shaded mask against an independently recomputed feasibility set: no disagreements over 2500 grid cells. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TVbQmQNRwrSerdh7WQX2Aw --- .../visualization/plots/optimization_plots.py | 190 ++++++++++++++---- .../visualization/adapters/plotly_adapter.py | 43 +++- src/process_improve/visualization/colors.py | 11 + tests/test_visualization.py | 108 +++++++++- 4 files changed, 293 insertions(+), 59 deletions(-) diff --git a/src/process_improve/experiments/visualization/plots/optimization_plots.py b/src/process_improve/experiments/visualization/plots/optimization_plots.py index 887540d4..23715b01 100644 --- a/src/process_improve/experiments/visualization/plots/optimization_plots.py +++ b/src/process_improve/experiments/visualization/plots/optimization_plots.py @@ -18,15 +18,15 @@ 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 # --------------------------------------------------------------------------- # Grid helpers shared by the desirability and overlay plots @@ -45,6 +45,50 @@ def _response_surface( 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. @@ -210,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: @@ -244,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, @@ -304,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_visualization.py b/tests/test_visualization.py index d302dda9..5e5506da 100644 --- a/tests/test_visualization.py +++ b/tests/test_visualization.py @@ -1183,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( @@ -1217,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={ @@ -1231,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: From 0d3d387cc8eec687b7cbceac3c9046ada595f87a Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 24 Jul 2026 22:57:46 +0000 Subject: [PATCH 3/5] Report intervals at the optimum, and name importance for what it is Four changes to the desirability path, all user-facing. The optimizer reported a point and nothing about how well that point was known. Coefficients alone cannot say: the residual variance and the design's leverage live on the fitted model object, not in its coefficients. Passing fitted_results alongside fitted_models now yields a confidence interval and a prediction interval per response at the optimum. Omitting it changes nothing. This matters in practice because a predicted optimum can sit inside its required range while the interval around it does not. Goals were consumed in list order while goal["response"] was documented as the key that ties a goal to its model, and was never read. Passing the two lists in different orders optimised the wrong thing and said nothing. They are now matched by name when both sides supply one, falling back to position with a warning when they do not. desirability_weights is renamed to response_importance. The old name described the wrong quantity: these are importances, setting how much each response counts in the composite, not the per-goal weight that shapes an individual ramp. The old name still works and warns. The desirability result also carries a "responses" list pairing coefficients with specification limits, so it can go straight to the overlay plot instead of the caller assembling a dict by hand. New tests assert where the optimum lands rather than only that a number came back: a plane's optimum at its known corner, two conflicting responses settling strictly between their individual optima, and raising one response's importance moving the compromise its way. The existing end-to-end tests asserted only that the composite exceeded zero, which a badly wrong optimum would also satisfy. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TVbQmQNRwrSerdh7WQX2Aw --- CHANGELOG.md | 50 +++- CITATION.cff | 4 +- pyproject.toml | 2 +- .../experiments/_tools/optimize_responses.py | 36 ++- .../experiments/optimization.py | 254 ++++++++++++++++-- tests/test_optimization.py | 238 ++++++++++++++++ tests/test_visualization_spec.py | 68 +++++ 7 files changed, 629 insertions(+), 23 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e561028..e6df8708 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,53 @@ those changes. ## [Unreleased] +## [1.61.0] - 2026-07-24 + +### 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. + +### 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.60.0] - 2026-07-23 ### Added @@ -2652,7 +2699,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.60.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.61.0...HEAD +[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 [1.58.0]: https://github.com/kgdunn/process-improve/compare/v1.57.0...v1.58.0 diff --git a/CITATION.cff b/CITATION.cff index 3b94d0e9..08ba5d92 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.60.0 -date-released: "2026-07-23" +version: 1.61.0 +date-released: "2026-07-24" keywords: - chemometrics - multivariate analysis diff --git a/pyproject.toml b/pyproject.toml index deed236b..cdc0efac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.60.0" +version = "1.61.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/_tools/optimize_responses.py b/src/process_improve/experiments/_tools/optimize_responses.py index 71a03682..0d96f5cb 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,9 +62,23 @@ 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).", + ) 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.", ) @@ -75,7 +92,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 +114,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 +139,8 @@ 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, 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 2245e9fc..db199613 100644 --- a/src/process_improve/experiments/optimization.py +++ b/src/process_improve/experiments/optimization.py @@ -26,10 +26,13 @@ import logging import re +import warnings from collections.abc import Callable 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 @@ -438,6 +441,64 @@ def _steepest_path( # noqa: PLR0913 } +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*. + + 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. + + Parameters + ---------- + fitted_models : list[dict] + Each optionally has ``"response_name"``. + goals : list[dict] + Each optionally has ``"response"``. + + Returns + ------- + list[dict] + Goals in the same order as *fitted_models*. + + Raises + ------ + ValueError + If the two lists differ in length. + """ + 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) + + model_names = [m.get("response_name") for m in fitted_models] + goal_names = [g.get("response") for g in goals] + + 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 + + return [by_name[str(n)] for n in model_names] + + def _optimize_desirability( # noqa: PLR0913 fitted_models: list[dict[str, Any]], goals: list[dict[str, Any]], @@ -453,13 +514,15 @@ 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. Returns ------- @@ -467,6 +530,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: @@ -613,6 +677,118 @@ 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, +) -> 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) + + 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, @@ -620,6 +796,9 @@ 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, desirability_weights: list[float] | None = None, ) -> dict[str, Any]: """Find optimal factor settings for one or multiple responses. @@ -640,16 +819,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"``, @@ -663,9 +850,22 @@ 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. 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 ------- @@ -673,6 +873,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 @@ -705,6 +912,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"] @@ -729,13 +952,14 @@ 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, ) elif method == "ridge_analysis": diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 737b1a25..39452fb2 100644 --- a/tests/test_optimization.py +++ b/tests/test_optimization.py @@ -585,6 +585,244 @@ 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] + ) + + # --------------------------------------------------------------------------- # Stubs # --------------------------------------------------------------------------- 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({})) From 7531a039eb459ede1bfebce365d2e83fc228c260 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 05:41:59 +0000 Subject: [PATCH 4/5] Let the optimizer search the region the design actually covers The desirability search was hard-coded to the coded cube, and a stationary point was judged inside or outside the design space by the same test. That is right for a two-level design, whose runs sit on the cube. It is wrong for a central composite design, whose axial runs sit at plus or minus alpha: the optimizer would refuse to consider settings the experiment had actually covered, and inside_design_space could report a legitimately explored point as outside. search_bounds now sets that region. It takes one (low, high) pair applied to every factor, or a mapping for per-factor control, and defaults to (-1, 1), so nothing changes for callers who do not pass it. The multi-start sampling follows the bounds too, rather than always seeding inside the cube, so widening the region actually widens where the search looks. Malformed bounds are rejected rather than silently coerced: a reversed or degenerate pair, a non-finite limit, or a mapping naming a factor the model does not have. That last one matters because a mistyped factor name would otherwise be ignored without complaint. Checked against the bioreactor example used in the book: its optimum is interior to the cube, so the answer there is unchanged, which is the reassuring case. The tests cover the cases where it is not: a plane whose optimum moves out to the widened bound, and a quadratic whose stationary point at 1.2 is outside the cube but inside a rotatable two-factor composite design's reach. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TVbQmQNRwrSerdh7WQX2Aw --- CHANGELOG.md | 8 ++ .../experiments/_tools/optimize_responses.py | 22 ++++ .../experiments/optimization.py | 116 ++++++++++++++++-- tests/test_optimization.py | 112 +++++++++++++++++ 4 files changed, 249 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e6df8708..b5c33882 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,14 @@ those changes. - `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 diff --git a/src/process_improve/experiments/_tools/optimize_responses.py b/src/process_improve/experiments/_tools/optimize_responses.py index 0d96f5cb..745ee0c3 100644 --- a/src/process_improve/experiments/_tools/optimize_responses.py +++ b/src/process_improve/experiments/_tools/optimize_responses.py @@ -76,12 +76,33 @@ class OptimizeResponsesInput(BaseModel): 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="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=( @@ -141,6 +162,7 @@ def optimize_responses_tool(spec: OptimizeResponsesInput) -> dict[str, Any]: 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 db199613..5767d2fc 100644 --- a/src/process_improve/experiments/optimization.py +++ b/src/process_improve/experiments/optimization.py @@ -27,7 +27,7 @@ import logging import re import warnings -from collections.abc import Callable +from collections.abc import Callable, Sequence from typing import Any import numpy as np @@ -222,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. @@ -267,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)}, @@ -441,6 +447,70 @@ def _steepest_path( # noqa: PLR0913 } +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. + + 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. + + 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. + """ + 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 + + if search_bounds is None: + return [default] * len(factor_names) + + 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 _align_goals_to_models( fitted_models: list[dict[str, Any]], goals: list[dict[str, Any]], @@ -506,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. @@ -523,6 +594,8 @@ def _optimize_desirability( # noqa: PLR0913 importances : list[float] or None 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 ------- @@ -542,8 +615,9 @@ def neg_composite(x: np.ndarray) -> float: d_vals.append(d) 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 @@ -555,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) @@ -753,6 +830,7 @@ def _desirability_result( # noqa: PLR0913 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. @@ -770,7 +848,9 @@ def _desirability_result( # noqa: PLR0913 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) + 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( @@ -799,6 +879,7 @@ def optimize_responses( # noqa: PLR0913, C901 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. @@ -862,6 +943,18 @@ def optimize_responses( # noqa: PLR0913, C901 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 Deprecated alias for *response_importance*. The name was misleading: these values are importances, not the ``weight`` that shapes an @@ -935,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" @@ -960,6 +1057,7 @@ def optimize_responses( # noqa: PLR0913, C901 response_importance=response_importance, fitted_results=fitted_results, significance_level=significance_level, + search_bounds=search_bounds, ) elif method == "ridge_analysis": diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 39452fb2..5c41b723 100644 --- a/tests/test_optimization.py +++ b/tests/test_optimization.py @@ -823,6 +823,118 @@ def test_mismatched_fitted_results_length_is_rejected(self) -> None: ) +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 # --------------------------------------------------------------------------- From b5cca8b69141571175e1ce0e3aafa978d657ca96 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 05:54:24 +0000 Subject: [PATCH 5/5] Cover the tool-level search_bounds path The region arrives over the MCP boundary as a JSON list rather than a tuple, so it is converted before reaching the library. That conversion had no test, including the per-factor mapping form and the error path for a reversed pair. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TVbQmQNRwrSerdh7WQX2Aw --- tests/test_experiments_tools.py | 44 +++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) 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