From d21f2aee2338895eb2454d8aac82774bdfbf75fc Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 06:18:53 +0000 Subject: [PATCH 01/16] Fix DSD conference-matrix construction: prime-power Paley, C16 table, no silent fallback Placeholder commit; work follows. From f06014cb8951846dbc12a9eb5956651b2f8c080c Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 06:30:04 +0000 Subject: [PATCH 02/16] Fix DSD conference-matrix construction so every factor count gives a real DSD The conference-matrix constructor implemented Paley's construction only for a prime q = m - 1. At every other order it fell back to a cyclic approximation that is not a conference matrix, so the definitive screening designs for 9, 10, 15, 16, 21, 22, 25, 26, 27 and 28 factors had main effects correlated with each other at up to r = 0.9. generate_design returned those with a log warning only, and design_type="omars" shares the construction: it already reported omars_verified: False for the same counts. Paley's construction now works over any odd prime-power field, which covers 9, 10, 25, 26, 27 and 28 factors at minimal run size, and a verified order-16 table covers 15 and 16. The fallback is removed: when the minimal order has no conference matrix the constructor steps up to the next order it can build and records that in the metadata, and when nothing can be built it raises. Order 22 is the first order with no conference matrix at all, so 21 and 22 factors now cost 49 runs rather than 45. Every constructed matrix is checked against C.T @ C == (m - 1) * I before use. The upstream definitive_screening_design package shows why: its order-10 table has one mistyped entry, which leaves its 9- and 10-factor designs slightly non-orthogonal with nothing to catch it. Two adjacent errors found while auditing the same path: the planner estimated 2k + 1 runs for every DSD, undercounting every odd factor count, and it advertised a nonexistent fake_factor parameter on the strength of a comment that had the parity argument backwards. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 40 +- CITATION.cff | 4 +- pyproject.toml | 2 +- .../experiments/designs_response_surface.py | 481 ++++++++++++++---- .../experiments/strategy/budget.py | 7 +- .../experiments/strategy/engine.py | 7 +- tests/test_conference_matrices.py | 323 ++++++++++++ tests/test_design_generation.py | 29 +- tests/test_design_properties.py | 67 ++- tests/test_recommend_strategy.py | 19 +- 10 files changed, 869 insertions(+), 110 deletions(-) create mode 100644 tests/test_conference_matrices.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e561028..dc92c2f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,43 @@ those changes. ## [Unreleased] +## [1.60.1] - 2026-07-25 + +### Fixed + +- Definitive screening designs are now genuine DSDs at every factor count. The + conference-matrix constructor only implemented Paley's construction for a prime + `q = m - 1`; at every other order it fell back to a cyclic approximation that is not + a conference matrix, so `generate_design(..., design_type="dsd")` returned designs + whose main effects were correlated with each other at up to r = 0.9 for 9, 10, 15, + 16, 21, 22, 25, 26, 27 and 28 factors. The same applied to `design_type="omars"`, + which shares the construction and already reported `omars_verified: False` for those + counts. Three changes: + - Paley's construction now works over any odd prime-power field (GF(9), GF(25), + GF(27), GF(49), ...), covering 9, 10, 25, 26, 27 and 28 factors at minimal run size. + - A verified order-16 conference matrix table covers 15 and 16 factors. + - The cyclic fallback is gone. When no conference matrix exists at the minimal order, + the constructor steps up to the next order it can build and reports the choice in + the design metadata (`conference_order` and `minimal_conference_order`); when + nothing can be built it raises. Every constructed matrix is checked against + `C.T @ C == (m - 1) * I` before use, so a mistyped table entry cannot slip through. + - Consequence for run counts: 21 and 22 factors now need 49 runs rather than 45, + because no conference matrix of order 22 exists. All other counts are unchanged. +- `estimate_screening_runs(k, "definitive_screening")` returned `2k + 1` for every `k`, + which undercounted every odd factor count (a 7-factor DSD has 17 runs, not 15) and + every count where the minimal conference matrix does not exist. It now asks the + constructor via the new `dsd_run_count`. +- The recommended-strategy plan advertised a `fake_factor` design parameter for + definitive screening designs, justified by a comment saying a DSD "needs odd factor + count". No such parameter exists and the claim is backwards: an odd factor count is + handled by building the next even conference order and dropping a column. + +### Added + +- `dsd_run_count(n_factors)` and `dsd_conference_order(n_factors)` in + `process_improve.experiments.designs_response_surface`, for planning a DSD's size + without building it. + ## [1.60.0] - 2026-07-23 ### Added @@ -2652,7 +2689,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.60.1...HEAD +[1.60.1]: https://github.com/kgdunn/process-improve/compare/v1.60.0...v1.60.1 [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..63cba69e 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.60.1 +date-released: "2026-07-25" keywords: - chemometrics - multivariate analysis diff --git a/pyproject.toml b/pyproject.toml index deed236b..65a6a341 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.60.0" +version = "1.60.1" 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/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index cc0e58bc..9002bd9b 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -293,31 +293,44 @@ def dispatch_box_behnken( def dispatch_dsd(factors: list[Factor]) -> tuple[np.ndarray, dict]: """Generate a Definitive Screening Design (DSD). - Follows the conference-matrix-based construction of Jones & Nachtsheim - (2011). For *k* factors the DSD has ``2k + 1`` runs when *k* is even - (using a conference matrix of order *k*) and ``2k + 3`` runs when *k* - is odd (using a conference matrix of order ``k + 1`` and dropping the - last column; Xiao, Lin & Bai 2012). The design can estimate all main - effects and quadratic effects and detect two-factor interactions with - minimal confounding, provided the underlying conference matrix is - genuine (``C.T @ C == (m-1) * I``). - - The Paley construction used by :func:`_conference_matrix` produces a - genuine conference matrix whenever ``m - 1`` is an odd prime. For other - *m* the function falls back to a cyclic approximation and logs a - warning; the resulting DSD will still run but its main-effects - orthogonality may be degraded. + Follows the conference-matrix foldover of Jones & Nachtsheim (2011): for a + conference matrix ``C`` of order *m* the design is ``[C; -C; 0]``, giving + ``2m + 1`` runs. For *k* factors the smallest usable order is ``m = k`` + when *k* is even and ``m = k + 1`` when *k* is odd (the surplus column is + dropped; Xiao, Lin & Bai 2012), so a DSD normally has ``2k + 1`` runs for + even *k* and ``2k + 3`` runs for odd *k*. + + A conference matrix does not exist for every even order. Order 22 is the + first exception: a conference matrix of order ``m ≡ 2 (mod 4)`` exists only + if ``m - 1`` is a sum of two squares (Belevitch 1950; van Lint & Seidel + 1966), and ``21 = 3 x 7`` is not. When the minimal order is unavailable + this function steps up to the next order it can construct and drops the + surplus columns, which costs runs but keeps the design a genuine DSD. The + order actually used is reported in the metadata as ``"conference_order"``. + + Every constructed matrix is checked against the defining property + ``C.T @ C == (m - 1) * I`` before it is used, so a degraded design can + never reach the caller. Parameters ---------- factors : list[Factor] - Continuous factors. + Continuous factors. At least three are required. Returns ------- tuple[np.ndarray, dict] - Coded design matrix and metadata. Metadata includes the name of - the conference-matrix construction that was used. + Coded design matrix and metadata. Metadata keys are + ``"construction"`` (which conference-matrix construction was used), + ``"conference_order"`` (the order *m* of that matrix) and, when the + minimal order was not constructible, ``"minimal_conference_order"`` + (the order that would have been used had it existed). + + Raises + ------ + ValueError + If fewer than three factors are supplied, or if no conference matrix + of usable order can be constructed. References ---------- @@ -332,19 +345,29 @@ def dispatch_dsd(factors: list[Factor]) -> tuple[np.ndarray, dict]: if k < 3: raise ValueError("Definitive Screening Designs require at least 3 factors.") - # For even k, use a conference matrix of order k (-> 2k + 1 runs). - # For odd k, use a conference matrix of order k + 1 and drop the last - # column so the design has k factors and 2(k+1) + 1 = 2k + 3 runs. - m = k if k % 2 == 0 else k + 1 - C, construction = _conference_matrix(m) - - zero_row = np.zeros((1, m)) - coded_matrix = np.vstack([C, -C, zero_row]) + minimal_order = k if k % 2 == 0 else k + 1 + order = _smallest_constructible_conference_order(minimal_order) + conference, construction = _conference_matrix(order) - if k % 2 == 1: + zero_row = np.zeros((1, order)) + coded_matrix = np.vstack([conference, -conference, zero_row]) + if order > k: coded_matrix = coded_matrix[:, :k] - return coded_matrix, {"construction": construction} + meta = {"construction": construction, "conference_order": order} + if order > minimal_order: + meta["minimal_conference_order"] = minimal_order + logger.info( + "No conference matrix of order %d exists, so the definitive screening design for %d factors " + "uses order %d instead: %d runs rather than the %d a minimal design would need.", + minimal_order, + k, + order, + 2 * order + 1, + 2 * minimal_order + 1, + ) + + return coded_matrix, meta def _is_prime(n: int) -> bool: @@ -363,18 +386,173 @@ def _is_prime(n: int) -> bool: return True +def _prime_power_factorization(q: int) -> tuple[int, int] | None: + """Factor *q* as ``p ** n`` for a prime *p*, or return None. + + Parameters + ---------- + q : int + The integer to factor. + + Returns + ------- + tuple[int, int] or None + ``(p, n)`` with ``p ** n == q``, or None when *q* is not a prime power + (including ``q < 2``). + """ + if q < 2: + return None + smallest_divisor = next(d for d in range(2, q + 1) if q % d == 0) + exponent, remaining = 0, q + while remaining % smallest_divisor == 0: + remaining //= smallest_divisor + exponent += 1 + return (smallest_divisor, exponent) if remaining == 1 else None + + +def _monic_polynomials(degree: int, p: int) -> list[list[int]]: + """Enumerate every monic polynomial of the given *degree* over GF(*p*). + + Polynomials are coefficient lists in ascending degree order, so + ``[2, 0, 1]`` is ``x**2 + 2``. + """ + polynomials = [] + for index in range(p**degree): + coefficients, remaining = [], index + for _ in range(degree): + coefficients.append(remaining % p) + remaining //= p + polynomials.append([*coefficients, 1]) + return polynomials + + +def _polynomial_remainder(dividend: list[int], divisor: list[int], p: int) -> list[int]: + """Remainder of *dividend* divided by the monic *divisor*, over GF(*p*). + + Both arguments are coefficient lists in ascending degree order. The result + is padded to exactly ``len(divisor) - 1`` coefficients. + """ + remainder = [c % p for c in dividend] + degree = len(divisor) - 1 + for i in range(len(remainder) - 1, degree - 1, -1): + coefficient = remainder[i] + if coefficient: + for j in range(degree + 1): + remainder[i - degree + j] = (remainder[i - degree + j] - coefficient * divisor[j]) % p + remainder = remainder[:degree] + return remainder + [0] * (degree - len(remainder)) + + +def _is_irreducible(polynomial: list[int], p: int) -> bool: + """Return True iff the monic *polynomial* is irreducible over GF(*p*). + + Uses trial division by every monic polynomial of degree up to half the + degree of *polynomial*, which is inexpensive at the field sizes needed + here (``p ** n`` of a few hundred at most). + """ + degree = len(polynomial) - 1 + for divisor_degree in range(1, degree // 2 + 1): + for divisor in _monic_polynomials(divisor_degree, p): + if not any(_polynomial_remainder(polynomial, divisor, p)): + return False + return True + + +def _irreducible_polynomial(p: int, n: int) -> list[int]: + """Return the first monic irreducible polynomial of degree *n* over GF(*p*).""" + for candidate in _monic_polynomials(n, p): + if _is_irreducible(candidate, p): + return candidate + raise ValueError(f"No irreducible polynomial of degree {n} found over GF({p}).") + + +def _decode_field_element(element: int, p: int, n: int) -> list[int]: + """Expand the integer label of a GF(``p ** n``) element into its coefficients.""" + coefficients, remaining = [], element + for _ in range(n): + coefficients.append(remaining % p) + remaining //= p + return coefficients + + +def _encode_field_element(coefficients: list[int], p: int) -> int: + """Collapse GF(``p ** n``) coefficients back into an integer label.""" + return sum(c * p**i for i, c in enumerate(coefficients)) + + +def _gf_multiplication_table(p: int, n: int) -> np.ndarray: + """Multiplication table of GF(``p ** n``), indexed by integer element labels. + + Element *e* stands for the polynomial whose base-*p* digits are the + coefficients of *e*, taken modulo a fixed irreducible polynomial of degree + *n*. For ``n == 1`` this is ordinary multiplication modulo *p*. + """ + q = p**n + if n == 1: + return np.outer(np.arange(p), np.arange(p)) % p + + modulus = _irreducible_polynomial(p, n) + table = np.zeros((q, q), dtype=int) + for a in range(q): + left = _decode_field_element(a, p, n) + for b in range(a, q): + right = _decode_field_element(b, p, n) + product = [0] * (2 * n - 1) + for i, x in enumerate(left): + if x: + for j, y in enumerate(right): + product[i + j] = (product[i + j] + x * y) % p + table[a, b] = table[b, a] = _encode_field_element(_polynomial_remainder(product, modulus, p), p) + return table + + +def _gf_subtraction_table(p: int, n: int) -> np.ndarray: + """Subtraction table of GF(``p ** n``): entry ``[a, b]`` is ``a - b``. + + Addition in GF(``p ** n``) is coefficient-wise modulo *p*, independent of + the choice of irreducible polynomial. + """ + q = p**n + table = np.zeros((q, q), dtype=int) + for a in range(q): + left = _decode_field_element(a, p, n) + for b in range(q): + right = _decode_field_element(b, p, n) + table[a, b] = _encode_field_element([(x - y) % p for x, y in zip(left, right, strict=True)], p) + return table + + +def _quadratic_character(p: int, n: int) -> np.ndarray: + """Legendre symbol on GF(``p ** n``): 0 at zero, +1 on squares, -1 otherwise.""" + q = p**n + multiplication = _gf_multiplication_table(p, n) + squares = {int(multiplication[x, x]) for x in range(1, q)} + character = np.full(q, -1, dtype=int) + character[0] = 0 + for square in squares: + character[square] = 1 + return character + + def _paley_conference_matrix(q: int) -> np.ndarray: """Build a conference matrix of order ``q + 1`` via Paley's construction. - Requires *q* to be an odd prime. Works for both ``q ≡ 1 (mod 4)`` + Requires *q* to be an odd prime power, which covers both ``q ≡ 1 (mod 4)`` (symmetric conference matrix, "Paley type II") and ``q ≡ 3 (mod 4)`` - (skew-symmetric conference matrix, "Paley type I"). In both cases - ``C.T @ C == q * I``. + (skew-symmetric conference matrix, "Paley type I"). In both cases the + result satisfies ``C.T @ C == q * I``. + + Prime powers with exponent above one (9, 25, 27, 49, ...) need arithmetic + in GF(``p ** n``) rather than integers modulo *q*; the field is built from + the first monic irreducible polynomial of the required degree. These are + the orders that make the difference for a DSD: without GF(9), GF(25) and + GF(27) there is no minimal construction for 9, 10, 25, 26, 27 or 28 + factors. Parameters ---------- q : int - An odd prime. + An odd prime power. Returns ------- @@ -382,81 +560,212 @@ def _paley_conference_matrix(q: int) -> np.ndarray: ``(q + 1) x (q + 1)`` matrix with 0s on the diagonal and ±1 off-diagonal. """ - # Legendre symbol χ : GF(q) -> {-1, 0, 1} - quadratic_residues = {(x * x) % q for x in range(1, q)} - chi = np.zeros(q, dtype=int) - for x in range(1, q): - chi[x] = 1 if x in quadratic_residues else -1 - - # Jacobsthal matrix Q[a, b] = χ(b - a) - q_matrix = np.zeros((q, q), dtype=int) - for a in range(q): - for b in range(q): - q_matrix[a, b] = chi[(b - a) % q] + factorization = _prime_power_factorization(q) + if factorization is None or q < 3 or q % 2 == 0: + raise ValueError(f"Paley's construction needs an odd prime power; got q={q}.") + p, n = factorization + + character = _quadratic_character(p, n) + subtraction = _gf_subtraction_table(p, n) + + # Jacobsthal matrix Q[a, b] = chi(b - a), where the subtraction is in the field. + jacobsthal = character[subtraction.T] - n = q + 1 - c_matrix = np.zeros((n, n), dtype=int) - c_matrix[0, 1:] = 1 + size = q + 1 + matrix = np.zeros((size, size), dtype=int) + matrix[0, 1:] = 1 if q % 4 == 1: # Symmetric Paley conference matrix. - c_matrix[1:, 0] = 1 + matrix[1:, 0] = 1 else: # Skew-symmetric Paley conference matrix (q ≡ 3 mod 4). - c_matrix[1:, 0] = -1 - c_matrix[1:, 1:] = q_matrix - return c_matrix - - -def _cyclic_conference_matrix(k: int) -> np.ndarray: - """Legacy cyclic approximation of a conference matrix. - - Does **not** satisfy ``C.T @ C == (k - 1) * I`` in general; used only as - a fallback when no Paley construction is available for the requested - order. + matrix[1:, 0] = -1 + matrix[1:, 1:] = jacobsthal + return matrix + + +# Conference matrix of order 16. 15 is not a prime power, so Paley's +# construction does not reach this order, and it is the only order below 22 +# that needs a table. Taken from the conference-matrix catalogue of Xiao, Lin +# & Bai (2012), by way of Jacob Albrecht's BSD-3-licensed MATLAB port of the +# JMP add-in (Bristol-Myers Squibb, 2015) and its Python translation by Daniele +# Ongari in the `definitive_screening_design` package (MIT-spirited, same +# BSD-3 header retained). Verified here against C.T @ C == 15 * I before use. +_CONFERENCE_MATRIX_16 = np.array( + [ + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [-1, 0, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1], + [-1, -1, 0, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1], + [-1, -1, -1, 0, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1], + [-1, 1, -1, -1, 0, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1], + [-1, -1, 1, -1, -1, 0, 1, 1, 1, -1, 1, -1, -1, -1, 1, 1], + [-1, 1, -1, 1, -1, -1, 0, 1, 1, 1, -1, 1, -1, -1, -1, 1], + [-1, 1, 1, -1, 1, -1, -1, 0, 1, 1, 1, -1, 1, -1, -1, -1], + [-1, -1, -1, -1, -1, -1, -1, -1, 0, 1, 1, 1, 1, 1, 1, 1], + [-1, 1, 1, 1, -1, 1, -1, -1, -1, 0, -1, -1, 1, -1, 1, 1], + [-1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 0, -1, -1, 1, -1, 1], + [-1, -1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 0, -1, -1, 1, -1], + [-1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 0, -1, -1, 1], + [-1, -1, 1, -1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 0, -1, -1], + [-1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 0, -1], + [-1, 1, 1, -1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 0], + ] +) + +_TABULATED_CONFERENCE_MATRICES: dict[int, np.ndarray] = {16: _CONFERENCE_MATRIX_16} + +# Upper bound on the conference-matrix order the search will consider. A DSD +# of order 200 already needs 401 runs, well past any practical screening study. +_MAX_CONFERENCE_ORDER = 200 + + +def _validate_conference_matrix(matrix: np.ndarray, order: int) -> None: + """Check the defining property ``C.T @ C == (m - 1) * I`` of a conference matrix. + + Raises + ------ + ValueError + If the matrix is the wrong shape, has a non-zero diagonal, holds values + other than 0 and ±1, or fails the orthogonality identity. A tabulated + matrix with a single mistyped entry fails here rather than silently + producing a design with correlated main effects. """ - c_matrix = np.zeros((k, k)) - half = (k - 1) // 2 - sequence = [1] * half + [-1] * (k - 1 - half) - for i in range(k): - for j in range(k): - if i != j: - idx = (j - i - 1) % (k - 1) if j > i else (j - i) % (k - 1) - c_matrix[i, j] = sequence[idx] - return c_matrix + if matrix.shape != (order, order): + raise ValueError(f"Conference matrix of order {order} has shape {matrix.shape}.") + if np.any(np.diag(matrix) != 0): + raise ValueError(f"Conference matrix of order {order} has a non-zero diagonal.") + if not np.all(np.isin(matrix, (-1, 0, 1))): + raise ValueError(f"Conference matrix of order {order} holds values outside {{-1, 0, +1}}.") + expected = (order - 1) * np.eye(order) + if not np.allclose(matrix.T @ matrix, expected, atol=1e-9): + raise ValueError( + f"Conference matrix of order {order} fails C.T @ C == {order - 1} * I; " + "the design built from it would not be a definitive screening design." + ) def _conference_matrix(m: int) -> tuple[np.ndarray, str]: """Construct an ``m x m`` conference matrix. - Uses Paley's construction when ``m - 1`` is an odd prime (covering - ``m ∈ {4, 6, 8, 12, 14, 18, 20, 24, 30, 32, 38, 42, 44, 48, 54, 60, 62, - 68, 72, 74, 80, 84, 90, 98, ...}``), which returns a genuine conference - matrix with ``C.T @ C == (m - 1) * I``. For other orders (e.g. - ``m ∈ {10, 16, 22, 26, 28, 34, 36, 40, ...}``) no Paley construction - with a prime *q* is available; the function falls back to a cyclic - approximation and logs a warning. + Two constructions are available. Paley's covers every even order *m* whose + predecessor ``m - 1`` is an odd prime power, which is most of them. A small + table covers order 16, the one gap below order 22. Parameters ---------- m : int - Desired order of the conference matrix. + Desired (even) order of the conference matrix. Returns ------- tuple[np.ndarray, str] The matrix and a short string identifying the construction used - (e.g. ``"paley_q=13"`` or ``"cyclic_fallback"``). + (``"paley_q=13"``, ``"tabulated_order_16"``). + + Raises + ------ + ValueError + If *m* is odd, or if no construction is available at that order. No + approximate matrix is ever returned: a design built from one would not + have orthogonal main effects, which is the defining property of a DSD. """ + if m < 2 or m % 2 != 0: + raise ValueError(f"A conference matrix has even order of at least 2; got m={m}.") + q = m - 1 - if q >= 3 and q % 2 == 1 and _is_prime(q): - return _paley_conference_matrix(q).astype(float), f"paley_q={q}" - - logger.warning( - "No Paley conference-matrix construction known for order m=%d " - "(q = m - 1 = %d is not an odd prime); falling back to a cyclic " - "approximation. The resulting DSD's main-effects orthogonality " - "may be degraded.", - m, - q, + factorization = _prime_power_factorization(q) + if q >= 3 and factorization is not None: + matrix = _paley_conference_matrix(q).astype(float) + construction = f"paley_q={q}" + elif m in _TABULATED_CONFERENCE_MATRICES: + matrix = _TABULATED_CONFERENCE_MATRICES[m].astype(float) + construction = f"tabulated_order_{m}" + else: + raise ValueError(f"No conference-matrix construction is available for order m={m}.") + + _validate_conference_matrix(matrix, m) + return matrix, construction + + +def _is_constructible_conference_order(m: int) -> bool: + """Return True iff :func:`_conference_matrix` can build order *m*.""" + if m < 2 or m % 2 != 0: + return False + q = m - 1 + return (q >= 3 and _prime_power_factorization(q) is not None) or m in _TABULATED_CONFERENCE_MATRICES + + +def _smallest_constructible_conference_order(minimum_order: int) -> int: + """Smallest order at least *minimum_order* for which a conference matrix can be built. + + Parameters + ---------- + minimum_order : int + Smallest acceptable (even) order. + + Returns + ------- + int + The order that will be used. + + Raises + ------ + ValueError + If no order up to ``_MAX_CONFERENCE_ORDER`` can be constructed. + """ + for order in range(minimum_order, _MAX_CONFERENCE_ORDER + 1, 2): + if _is_constructible_conference_order(order): + return order + + raise ValueError( + f"No conference matrix could be constructed with order between {minimum_order} and " + f"{_MAX_CONFERENCE_ORDER}; a definitive screening design is not available at this size." ) - return _cyclic_conference_matrix(m), "cyclic_fallback" + + +def dsd_conference_order(n_factors: int) -> int: + """Return the order of the conference matrix a DSD for *n_factors* factors is built from. + + Normally ``n_factors`` for an even count and ``n_factors + 1`` for an odd + count, but larger when no conference matrix exists at that order (21 and 22 + factors, for instance, need order 24 rather than 22). + + Parameters + ---------- + n_factors : int + Number of factors, at least three. + + Returns + ------- + int + The conference-matrix order. + """ + if n_factors < 3: + raise ValueError("Definitive Screening Designs require at least 3 factors.") + return _smallest_constructible_conference_order(n_factors if n_factors % 2 == 0 else n_factors + 1) + + +def dsd_run_count(n_factors: int) -> int: + """Return the number of runs in the definitive screening design for *n_factors* factors. + + This is ``2m + 1`` for the conference-matrix order *m* returned by + :func:`dsd_conference_order`, which equals the familiar ``2k + 1`` (even + *k*) or ``2k + 3`` (odd *k*) except at the factor counts where the minimal + conference matrix does not exist. + + Parameters + ---------- + n_factors : int + Number of factors, at least three. + + Returns + ------- + int + Run count, including the design's single centre run. + + Examples + -------- + >>> dsd_run_count(6), dsd_run_count(7), dsd_run_count(22) + (13, 17, 49) + """ + return 2 * dsd_conference_order(n_factors) + 1 diff --git a/src/process_improve/experiments/strategy/budget.py b/src/process_improve/experiments/strategy/budget.py index 80ffa333..97151ee0 100644 --- a/src/process_improve/experiments/strategy/budget.py +++ b/src/process_improve/experiments/strategy/budget.py @@ -69,7 +69,12 @@ def estimate_screening_runs(n_factors: int, design_type: str) -> int: # noqa: P return int(math.ceil(n / 4) * 4) if design_type == "definitive_screening": - return 2 * n_factors + 1 + # 2k + 1 runs for an even factor count and 2k + 3 for an odd one, except + # at the counts where the minimal conference matrix does not exist and a + # larger one has to be used. Ask the constructor rather than guessing. + from process_improve.experiments.designs_response_surface import dsd_run_count # noqa: PLC0415 + + return dsd_run_count(n_factors) if n_factors >= 3 else 2 * n_factors + 1 if design_type == "fractional_factorial": # Smallest 2^(k-p) with resolution >= IV diff --git a/src/process_improve/experiments/strategy/engine.py b/src/process_improve/experiments/strategy/engine.py index c3cba685..b7d818b0 100644 --- a/src/process_improve/experiments/strategy/engine.py +++ b/src/process_improve/experiments/strategy/engine.py @@ -226,7 +226,7 @@ def _screening_design_params(design_type: str, n: int) -> dict[str, Any]: if design_type == "plackett_burman": return {"center_points": 0} # PB typically without center points if design_type == "definitive_screening": - return {"fake_factor": n % 2 == 0} # DSD needs odd factor count + return {"center_points": 0} # the DSD's foldover already carries a centre run return {} @@ -554,9 +554,8 @@ def _build_alternatives(spec: DOEProblemSpec, classification: dict[str, Any]) -> n = classification["n_factors"] if n >= 6: - alternatives.append( - f"Definitive Screening Design ({2 * n + 1} runs) to combine screening and curvature detection." - ) + runs = estimate_screening_runs(n, "definitive_screening") + alternatives.append(f"Definitive Screening Design ({runs} runs) to combine screening and curvature detection.") if n <= 5: alternatives.append(f"Full factorial 2^{n} ({2**n} runs) if budget allows complete information.") if n >= 4: diff --git a/tests/test_conference_matrices.py b/tests/test_conference_matrices.py new file mode 100644 index 00000000..1bc6fd7d --- /dev/null +++ b/tests/test_conference_matrices.py @@ -0,0 +1,323 @@ +"""Tests for the conference-matrix constructions behind the definitive screening design. + +A DSD is only a DSD if the matrix it folds is a genuine conference matrix: +``C.T @ C == (m - 1) * I``. If that fails, main effects come out correlated +with each other and the design silently stops doing the one thing it exists to +do. These tests check the constructions, the finite-field arithmetic they rest +on, and the guard that refuses anything that does not satisfy the identity. +""" + +from __future__ import annotations + +import re + +import numpy as np +import pytest + +from process_improve.experiments.designs_response_surface import ( + _conference_matrix, + _gf_multiplication_table, + _irreducible_polynomial, + _is_constructible_conference_order, + _polynomial_remainder, + _prime_power_factorization, + _quadratic_character, + _smallest_constructible_conference_order, + _validate_conference_matrix, + dsd_conference_order, + dsd_run_count, +) + +# Odd prime powers that matter for DSDs at practical factor counts. +_PRIME_POWERS = [3, 5, 7, 9, 11, 13, 17, 19, 23, 25, 27, 29, 31, 37, 41, 43, 47, 49, 81, 121, 125] + + +class TestPrimePowerFactorization: + """Recognising ``q = p ** n``.""" + + @pytest.mark.parametrize( + ("q", "expected"), + [ + (3, (3, 1)), + (9, (3, 2)), + (25, (5, 2)), + (27, (3, 3)), + (49, (7, 2)), + (81, (3, 4)), + (125, (5, 3)), + (169, (13, 2)), + ], + ) + def test_prime_powers_recognised(self, q: int, expected: tuple[int, int]) -> None: + assert _prime_power_factorization(q) == expected + + @pytest.mark.parametrize("q", [1, 0, -4, 15, 21, 33, 35, 45, 51, 55, 57]) + def test_non_prime_powers_rejected(self, q: int) -> None: + """15, 21, 33, ... are the orders that force a table or a step up.""" + assert _prime_power_factorization(q) is None + + def test_agrees_with_brute_force(self) -> None: + """Cross-check against an independent definition over a wide range.""" + primes = [p for p in range(2, 500) if all(p % d for d in range(2, int(p**0.5) + 1))] + for q in range(2, 500): + brute = None + for p in primes: + n, remaining = 0, q + while remaining % p == 0: + remaining //= p + n += 1 + if remaining == 1 and n > 0: + brute = (p, n) + break + assert _prime_power_factorization(q) == brute, f"disagreement at q={q}" + + +class TestFiniteFieldArithmetic: + """GF(p**n) built from the first monic irreducible polynomial of degree n.""" + + @pytest.mark.parametrize("q", _PRIME_POWERS) + def test_multiplication_is_a_group_on_nonzero_elements(self, q: int) -> None: + """Every non-zero row of the table permutes the non-zero elements. + + This is the property that fails if the modulus polynomial is reducible: + the ring would then have zero divisors and rows would repeat values. + """ + p, n = _prime_power_factorization(q) + table = _gf_multiplication_table(p, n) + assert table.shape == (q, q) + assert np.all(table[0, :] == 0) + assert np.all(table[:, 0] == 0) + for row in range(1, q): + assert sorted(table[row, 1:]) == list(range(1, q)) + + @pytest.mark.parametrize("q", _PRIME_POWERS) + def test_multiplication_is_commutative_and_associative(self, q: int) -> None: + p, n = _prime_power_factorization(q) + table = _gf_multiplication_table(p, n) + assert np.array_equal(table, table.T) + rng = np.random.default_rng(seed=42) + for _ in range(50): + a, b, c = rng.integers(0, q, size=3) + assert table[table[a, b], c] == table[a, table[b, c]] + + @pytest.mark.parametrize("q", _PRIME_POWERS) + def test_quadratic_character_is_multiplicative(self, q: int) -> None: + """chi(ab) == chi(a) chi(b), and exactly half the non-zero elements are squares.""" + p, n = _prime_power_factorization(q) + table = _gf_multiplication_table(p, n) + chi = _quadratic_character(p, n) + + assert chi[0] == 0 + assert set(np.unique(chi[1:])) == {-1, 1} + assert (chi[1:] == 1).sum() == (q - 1) // 2 + + for a in range(1, q): + for b in range(1, q): + assert chi[table[a, b]] == chi[a] * chi[b] + + def test_prime_field_matches_integer_arithmetic(self) -> None: + """For n == 1 the table must be plain multiplication modulo p.""" + for p in (3, 5, 7, 11, 13): + table = _gf_multiplication_table(p, 1) + assert np.array_equal(table, np.outer(np.arange(p), np.arange(p)) % p) + + def test_legendre_symbol_matches_euler_criterion(self) -> None: + """On a prime field, chi(a) must equal a**((p-1)/2) mod p, mapped to ±1.""" + for p in (3, 5, 7, 11, 13, 17, 19, 23): + chi = _quadratic_character(p, 1) + for a in range(1, p): + euler = pow(a, (p - 1) // 2, p) + assert chi[a] == (1 if euler == 1 else -1), f"p={p}, a={a}" + + +class TestPolynomialHelpers: + """Polynomial remainder and irreducibility over GF(p).""" + + def test_remainder_of_known_division(self) -> None: + """(x**2 + 1) mod (x + 1) over GF(3) is 2, since (-1)**2 + 1 = 2.""" + assert _polynomial_remainder([1, 0, 1], [1, 1], 3) == [2] + + def test_remainder_degree_is_bounded(self) -> None: + rng = np.random.default_rng(seed=7) + for p in (3, 5): + for _ in range(50): + dividend = [int(c) for c in rng.integers(0, p, size=6)] + divisor = [*[int(c) for c in rng.integers(0, p, size=3)], 1] + remainder = _polynomial_remainder(dividend, divisor, p) + assert len(remainder) == len(divisor) - 1 + assert all(0 <= c < p for c in remainder) + + def test_remainder_of_a_multiple_is_zero(self) -> None: + """Dividing the product f * g by g must leave no remainder.""" + p = 5 + g = [2, 3, 1] # x**2 + 3x + 2 + f = [4, 1] # x + 4 + product = [0] * (len(f) + len(g) - 1) + for i, x in enumerate(f): + for j, y in enumerate(g): + product[i + j] = (product[i + j] + x * y) % p + assert not any(_polynomial_remainder(product, g, p)) + + @pytest.mark.parametrize(("p", "n"), [(3, 2), (3, 3), (3, 4), (5, 2), (5, 3), (7, 2), (11, 2)]) + def test_irreducible_polynomial_has_no_roots(self, p: int, n: int) -> None: + """A necessary condition, and sufficient for degree 2 and 3.""" + polynomial = _irreducible_polynomial(p, n) + assert len(polynomial) == n + 1 + assert polynomial[-1] == 1 + for x in range(p): + value = sum(c * pow(x, i, p) for i, c in enumerate(polynomial)) % p + assert value != 0, f"GF({p}**{n}) modulus has root x={x}" + + +class TestConferenceMatrices: + """The constructions themselves.""" + + @pytest.mark.parametrize("m", [4, 6, 8, 10, 12, 14, 16, 18, 20, 24, 26, 28, 30, 32, 38, 42, 44, 48, 50]) + def test_defining_property(self, m: int) -> None: + """C.T @ C == (m - 1) I, zero diagonal, entries in {-1, 0, +1}.""" + matrix, _construction = _conference_matrix(m) + assert matrix.shape == (m, m) + assert np.all(np.diag(matrix) == 0) + assert np.all(np.isin(matrix, (-1.0, 0.0, 1.0))) + assert np.allclose(matrix.T @ matrix, (m - 1) * np.eye(m)) + assert np.allclose(matrix @ matrix.T, (m - 1) * np.eye(m)) + assert np.all((matrix == 0).sum(axis=0) == 1) + assert np.all((matrix == 0).sum(axis=1) == 1) + + @pytest.mark.parametrize(("m", "expected"), [(10, "paley_q=9"), (16, "tabulated_order_16"), (26, "paley_q=25")]) + def test_construction_reported(self, m: int, expected: str) -> None: + """The orders that the prime-only Paley construction used to miss.""" + _matrix, construction = _conference_matrix(m) + assert construction == expected + + def test_order_16_table_is_exact(self) -> None: + """The tabulated matrix is checked, not trusted. + + The upstream `definitive_screening_design` package carries an order-10 + table with a single mistyped entry, which leaves its 9- and 10-factor + designs slightly non-orthogonal. A table is only as good as its + verification, so verify it. + """ + matrix, _construction = _conference_matrix(16) + assert np.array_equal(matrix.T @ matrix, 15 * np.eye(16)) + # Order 16 is skew-symmetric: C.T == -C away from the border convention. + assert np.array_equal(matrix[1:, 1:], -matrix[1:, 1:].T) + + @pytest.mark.parametrize("m", [22, 34, 36, 40, 46]) + def test_unconstructible_orders_raise(self, m: int) -> None: + """No approximation is ever returned. + + Order 22 is the first even order with no conference matrix at all: one + exists for ``m = 2 (mod 4)`` only when ``m - 1`` is a sum of two + squares, and 21 is not. The others are orders this module has no + construction for. Both cases must raise so the caller steps up. + """ + with pytest.raises(ValueError, match="No conference-matrix construction"): + _conference_matrix(m) + + @pytest.mark.parametrize("m", [3, 5, 7, 15, 0, -2]) + def test_odd_orders_rejected(self, m: int) -> None: + with pytest.raises(ValueError, match="even order"): + _conference_matrix(m) + + def test_constructibility_predicate_agrees_with_the_constructor(self) -> None: + """The cheap predicate and the actual construction must never disagree.""" + for m in range(2, 121): + predicted = _is_constructible_conference_order(m) + try: + _conference_matrix(m) + actual = True + except ValueError: + actual = False + assert predicted == actual, f"disagreement at m={m}" + + def test_smallest_order_steps_up_only_when_it_must(self) -> None: + assert _smallest_constructible_conference_order(6) == 6 + assert _smallest_constructible_conference_order(10) == 10 + assert _smallest_constructible_conference_order(16) == 16 + assert _smallest_constructible_conference_order(22) == 24 + assert _smallest_constructible_conference_order(28) == 28 + assert _smallest_constructible_conference_order(34) == 38 + + +class TestConferenceMatrixValidation: + """The guard that stops a bad matrix reaching a design.""" + + def test_accepts_a_valid_matrix(self) -> None: + matrix, _construction = _conference_matrix(6) + _validate_conference_matrix(matrix, 6) # must not raise + + def test_rejects_a_single_flipped_entry(self) -> None: + """Exactly the upstream order-10 defect: one sign wrong out of 100.""" + matrix, _construction = _conference_matrix(10) + matrix[7, 3] = -matrix[7, 3] + with pytest.raises(ValueError, match=re.escape("fails C.T @ C")): + _validate_conference_matrix(matrix, 10) + + def test_rejects_a_non_zero_diagonal(self) -> None: + matrix, _construction = _conference_matrix(6) + matrix[2, 2] = 1 + with pytest.raises(ValueError, match="non-zero diagonal"): + _validate_conference_matrix(matrix, 6) + + def test_rejects_out_of_range_values(self) -> None: + matrix, _construction = _conference_matrix(6) + matrix[1, 2] = 2 + with pytest.raises(ValueError, match="outside"): + _validate_conference_matrix(matrix, 6) + + def test_rejects_the_wrong_shape(self) -> None: + matrix, _construction = _conference_matrix(6) + with pytest.raises(ValueError, match="has shape"): + _validate_conference_matrix(matrix[:, :4], 6) + + +class TestRunCountHelpers: + """The public run-count helpers used by the strategy layer.""" + + @pytest.mark.parametrize( + ("k", "order", "runs"), + [ + (3, 4, 9), + (4, 4, 9), + (5, 6, 13), + (6, 6, 13), + (9, 10, 21), + (10, 10, 21), + (15, 16, 33), + (16, 16, 33), + (21, 24, 49), # order 22 does not exist, so 49 runs rather than 45 + (22, 24, 49), + (25, 26, 53), + (26, 26, 53), + (27, 28, 57), + (28, 28, 57), + ], + ) + def test_known_values(self, k: int, order: int, runs: int) -> None: + assert dsd_conference_order(k) == order + assert dsd_run_count(k) == runs + + def test_matches_the_generated_design(self) -> None: + """The prediction must equal what the generator actually produces.""" + from process_improve.experiments.designs_response_surface import dispatch_dsd + from process_improve.experiments.factor import Factor + + for k in range(3, 41): + matrix, meta = dispatch_dsd([Factor(name=f"X{i}", low=-1, high=1) for i in range(k)]) + assert matrix.shape[0] == dsd_run_count(k), f"run count mismatch at k={k}" + assert meta["conference_order"] == dsd_conference_order(k), f"order mismatch at k={k}" + + def test_usual_formula_holds_away_from_the_exceptions(self) -> None: + """2k + 1 for even k and 2k + 3 for odd k, except where no matrix exists.""" + exceptions = {21, 22, 33, 34, 35, 36, 39, 40} + for k in range(3, 41): + if k in exceptions: + assert dsd_run_count(k) > (2 * k + 1 if k % 2 == 0 else 2 * k + 3) + else: + assert dsd_run_count(k) == (2 * k + 1 if k % 2 == 0 else 2 * k + 3) + + @pytest.mark.parametrize("k", [0, 1, 2]) + def test_too_few_factors(self, k: int) -> None: + with pytest.raises(ValueError, match="at least 3 factors"): + dsd_run_count(k) diff --git a/tests/test_design_generation.py b/tests/test_design_generation.py index 2188009f..41e63e4b 100644 --- a/tests/test_design_generation.py +++ b/tests/test_design_generation.py @@ -508,16 +508,41 @@ def test_4_factors(self) -> None: assert result.n_runs == 9 def test_main_effects_orthogonal(self) -> None: - """Paley-constructed DSD should have mutually orthogonal main effects.""" + """A DSD should have mutually orthogonal main effects.""" factors = _continuous_factors(6) result = generate_design(factors, design_type="dsd", center_points=0) x = result.design[result.factor_names].values.astype(float) gram = x.T @ x off_diag = gram - np.diag(np.diag(gram)) assert np.abs(off_diag).max() < 1e-9 - # Paley construction should have been used (no cyclic fallback warning). assert result.metadata["construction"].startswith("paley") + def test_ten_factors_uses_the_prime_power_field(self) -> None: + """15 factors and 10 factors have no prime-order Paley construction. + + 10 factors needs GF(9) and 15 needs the tabulated order-16 matrix. + Both used to fall through to an approximation that left main effects + correlated. + """ + result = generate_design(_continuous_factors(10), design_type="dsd", center_points=0) + assert result.metadata["construction"] == "paley_q=9" + assert result.n_runs == 21 + + result = generate_design(_continuous_factors(15), design_type="dsd", center_points=0) + assert result.metadata["construction"] == "tabulated_order_16" + assert result.n_runs == 33 + + def test_twenty_two_factors_steps_up_to_a_larger_conference_matrix(self) -> None: + """No conference matrix of order 22 exists, so the design costs 4 extra runs.""" + result = generate_design(_continuous_factors(22), design_type="dsd", center_points=0) + assert result.n_runs == 49 + assert result.metadata["conference_order"] == 24 + assert result.metadata["minimal_conference_order"] == 22 + + x = result.design[result.factor_names].values.astype(float) + gram = x.T @ x + assert np.abs(gram - np.diag(np.diag(gram))).max() < 1e-9 + def test_requires_3_factors(self) -> None: """DSD should require at least 3 factors.""" factors = _continuous_factors(2, "AB") diff --git a/tests/test_design_properties.py b/tests/test_design_properties.py index 00ec8cce..9b0ec28f 100644 --- a/tests/test_design_properties.py +++ b/tests/test_design_properties.py @@ -231,22 +231,33 @@ class TestDSDProperties: """Jones-Nachtsheim structural invariants of the DSD.""" @_settings - @given(k=st.integers(min_value=3, max_value=14)) + @given(k=st.integers(min_value=3, max_value=30)) def test_run_count_matches_jones_nachtsheim(self, k: int) -> None: - """Even k -> 2k+1 runs; odd k -> 2k+3 runs (Jones-Nachtsheim 2011).""" + """Even k -> 2k+1 runs; odd k -> 2k+3 runs (Jones-Nachtsheim 2011). + + The exception is a factor count whose minimal conference matrix does not + exist, where a larger one is used and the design costs more runs. 21 + and 22 factors are the first such counts: no conference matrix of order + 22 exists, so order 24 is used and the design has 49 runs, not 45. + """ result = generate_design(_factors(k), design_type="dsd", center_points=0) expected = 2 * k + 1 if k % 2 == 0 else 2 * k + 3 - assert result.n_runs == expected + if k in (21, 22): + assert result.n_runs == 49 + assert result.metadata["conference_order"] == 24 + assert result.metadata["minimal_conference_order"] == 22 + else: + assert result.n_runs == expected @_settings - @given(k=st.integers(min_value=3, max_value=14)) + @given(k=st.integers(min_value=3, max_value=30)) def test_columns_balanced(self, k: int) -> None: """Every DSD factor column sums to 0 (foldover structure).""" x = _coded(generate_design(_factors(k), design_type="dsd", center_points=0)) assert np.allclose(x.sum(axis=0), 0.0, atol=1e-9) @_settings - @given(k=st.integers(min_value=3, max_value=14)) + @given(k=st.integers(min_value=3, max_value=30)) def test_exactly_one_zero_pair_per_factor(self, k: int) -> None: """Each factor takes value 0 in exactly 2 rows for odd k and 1 row for even k (main rows).""" result = generate_design(_factors(k), design_type="dsd", center_points=0) @@ -259,18 +270,52 @@ def test_exactly_one_zero_pair_per_factor(self, k: int) -> None: assert np.all(zeros_per_column == 3) @_settings - @given(k=st.sampled_from([3, 4, 5, 6, 7, 8, 12, 13, 14, 18, 19, 20])) - def test_paley_construction_is_orthogonal(self, k: int) -> None: - """When a Paley conference matrix is used, main effects are exactly orthogonal.""" + @given(k=st.integers(min_value=3, max_value=30)) + def test_main_effects_are_exactly_orthogonal(self, k: int) -> None: + """Main effects are mutually orthogonal at every factor count, with no exceptions. + + This invariant used to be conditional on the Paley construction having + been reached, because the cyclic fallback returned a matrix that is not + a conference matrix at all. It now holds unconditionally: the + construction either succeeds or raises. + """ result = generate_design(_factors(k), design_type="dsd", center_points=0) - if not result.metadata.get("construction", "").startswith("paley"): - # Cyclic fallback is known-approximate, not covered by this invariant. - return x = _coded(result) gram = x.T @ x off_diag = gram - np.diag(np.diag(gram)) assert np.abs(off_diag).max() < 1e-9 + @_settings + @given(k=st.integers(min_value=3, max_value=30)) + def test_main_effects_clear_of_second_order_terms(self, k: int) -> None: + """The defining DSD property: no main effect is aliased with any quadratic or two-factor interaction.""" + x = _coded(generate_design(_factors(k), design_type="dsd", center_points=0)) + quadratics = x**2 + quadratics = quadratics - quadratics.mean(axis=0) + assert np.abs(x.T @ quadratics).max() < 1e-9 + + interactions = np.array([x[:, i] * x[:, j] for i in range(k) for j in range(i + 1, k)]).T + interactions = interactions - interactions.mean(axis=0) + assert np.abs(x.T @ interactions).max() < 1e-9 + + @pytest.mark.parametrize("k", [9, 10, 15, 16, 21, 22, 25, 26, 27, 28]) + def test_previously_degraded_factor_counts(self, k: int) -> None: + """Regression test for the factor counts the cyclic fallback used to mishandle. + + Before the prime-power Paley construction and the order-16 table, these + counts fell through to an approximation whose main effects correlated + with each other at up to r = 0.9, and `generate_design` returned it with + only a log warning. + """ + result = generate_design(_factors(k), design_type="dsd", center_points=0) + assert result.metadata["construction"] != "cyclic_fallback" + + x = _coded(result) + gram = x.T @ x + off_diag = gram - np.diag(np.diag(gram)) + assert np.abs(off_diag).max() < 1e-9 + assert np.all((x == 0).sum(axis=0) == 3) + # --------------------------------------------------------------------------- # evaluate_design metrics diff --git a/tests/test_recommend_strategy.py b/tests/test_recommend_strategy.py index dccc786a..2743113e 100644 --- a/tests/test_recommend_strategy.py +++ b/tests/test_recommend_strategy.py @@ -247,8 +247,23 @@ def test_pb_12_runs_for_11_factors(self): assert runs == 12 # next mult of 4 >= 12 def test_dsd_runs(self): - runs = estimate_screening_runs(7, "definitive_screening") - assert runs == 15 # 2*7 + 1 + # An odd factor count needs a conference matrix of order k + 1, so a + # 7-factor DSD has 2*8 + 1 = 17 runs, not 2*7 + 1. The estimate now + # comes from the constructor rather than from the even-k formula. + assert estimate_screening_runs(7, "definitive_screening") == 17 + assert estimate_screening_runs(8, "definitive_screening") == 17 # 2*8 + 1 + # No conference matrix of order 22 exists, so 21 factors need order 24. + assert estimate_screening_runs(21, "definitive_screening") == 49 + + def test_dsd_run_estimate_matches_the_generated_design(self): + """The planner's estimate must not undercount what generate_design produces.""" + from process_improve.experiments.designs import generate_design + from process_improve.experiments.factor import Factor + + for k in (5, 6, 7, 10, 15, 21): + factors = [Factor(name=f"X{i}", low=0, high=1) for i in range(k)] + design = generate_design(factors, design_type="dsd", center_points=0) + assert estimate_screening_runs(k, "definitive_screening") == design.n_runs def test_full_factorial_runs(self): runs = estimate_screening_runs(3, "full_factorial") From cb16cf76e30718a840f0a6bc797c4767f938772c Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 06:35:32 +0000 Subject: [PATCH 03/16] Update DSD documentation to match the corrected construction The coverage note described the cyclic fallback as a known limitation, the OMARS docstring pointed at it as the reason for its verification step, and the knowledge base advertised 2k + 1 runs for every factor count. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- docs/doe/coverage.md | 2 +- src/process_improve/experiments/designs_omars.py | 14 +++++++++----- .../experiments/knowledge/data/design_types.yaml | 4 ++-- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/docs/doe/coverage.md b/docs/doe/coverage.md index fec9d443..c8fd0748 100644 --- a/docs/doe/coverage.md +++ b/docs/doe/coverage.md @@ -54,7 +54,7 @@ All 14 metrics live in `experiments/evaluate.py` behind the `_METRIC_REGISTRY`: ## Caveats -- **DSD conference matrix.** `_conference_matrix` in `designs_response_surface.py` uses Paley's construction when `m − 1` is an odd prime (covers `m ∈ {4, 6, 8, 12, 14, 18, 20, 24, 30, 32, 38, 42, 44, 48, 54, 60, 62, 68, 72, 74, 80, 84, 90, 98, …}`). For other *m* (including `m ∈ {10, 16, 22, 26, 28, 34, 36, 40, 46, 50, 52, 56, …}`) the function falls back to a cyclic approximation that does **not** satisfy `Cᵀ C = (m − 1) I`, and logs a warning. Main-effects orthogonality of the resulting DSD may be degraded in those sizes. +- **DSD conference matrix.** `_conference_matrix` in `designs_response_surface.py` builds order *m* from Paley's construction whenever `m − 1` is an odd prime power (`m ∈ {4, 6, 8, 10, 12, 14, 18, 20, 24, 26, 28, 30, 32, 38, 42, 44, 48, 50, 54, 60, 62, …}`), and from a tabulated matrix at `m = 16`. Every matrix is verified against `Cᵀ C = (m − 1) I` before use. At the remaining orders (`m ∈ {22, 34, 36, 40, 46, 52, 56, 58, …}`) the constructor raises, and `dispatch_dsd` steps up to the next order it can build, dropping the surplus columns. That costs runs but never orthogonality: 21 and 22 factors use order 24, so 49 runs rather than the 45 a minimal design would need. Order 22 is not a gap in this module but a fact about conference matrices: one of order `m ≡ 2 (mod 4)` exists only when `m − 1` is a sum of two squares, and 21 is not. Orders 34 and 58 fail the same test; 36, 40, 46 and 52 are orders this module has no construction for, though a matrix may exist. Use `dsd_run_count(k)` to find the size a given factor count will produce. - **Optimal designs without `pyoptex`.** D-optimal has a point-exchange fallback; I/A-optimal raise `ImportError` instead. Hard-to-change factors (split-plot structure) are currently ignored with a `logger.warning` when `pyoptex` is not available. - **Taguchi OA auto-selection** picks the smallest standard array that covers the requested factor count and level counts, but the underlying `pyDOE3.taguchi_design` requires the number of `levels_per_factor` entries to match the OA column count exactly, so requesting a Taguchi design with fewer factors than any available OA will fail. Users typically pick *k* to match one of the standard arrays (3, 4, 7, 11, 15, …). diff --git a/src/process_improve/experiments/designs_omars.py b/src/process_improve/experiments/designs_omars.py index e6d06ebb..1d4e2f57 100644 --- a/src/process_improve/experiments/designs_omars.py +++ b/src/process_improve/experiments/designs_omars.py @@ -236,9 +236,12 @@ def dispatch_omars(factors: list[Factor], *, verify: bool = True) -> tuple[np.nd folds a conference matrix ``C`` of order ``m`` over its negative and adds a center run, giving ``[C; -C; 0]``. ``m = k`` for even ``k`` (``2k + 1`` runs) and ``m = k + 1`` with the last column dropped for odd ``k`` - (``2k + 3`` runs). This yields the minimal OMARS design for ``k`` factors; - a future enumerator will expand coverage to the larger, non-foldover - members of the OMARS family. + (``2k + 3`` runs), except at the factor counts where no conference matrix + exists at that order and a larger one has to be used; see + :func:`~process_improve.experiments.designs_response_surface.dsd_run_count`. + This yields the minimal OMARS design for ``k`` factors; a future enumerator + will expand coverage to the larger, non-foldover members of the OMARS + family. Parameters ---------- @@ -247,8 +250,9 @@ def dispatch_omars(factors: list[Factor], *, verify: bool = True) -> tuple[np.nd verify : bool When ``True`` (default) the generated matrix is checked with :func:`is_omars` and the result is recorded in the metadata under - ``"omars_verified"``. The check is cheap and guards against the - degraded orthogonality of the cyclic conference-matrix fallback. + ``"omars_verified"``. The check is cheap and is a second line of + defence behind the conference-matrix verification in + :func:`~process_improve.experiments.designs_response_surface.dispatch_dsd`. Returns ------- diff --git a/src/process_improve/experiments/knowledge/data/design_types.yaml b/src/process_improve/experiments/knowledge/data/design_types.yaml index b7de854e..8f36919c 100644 --- a/src/process_improve/experiments/knowledge/data/design_types.yaml +++ b/src/process_improve/experiments/knowledge/data/design_types.yaml @@ -303,8 +303,8 @@ supports_quadratic: "Partial (if <= 3 active factors)" supports_interaction: "Partial (if <= 3 active factors)" min_runs: - formula: "2k + 1" - examples: {"6": 13, "8": 17, "10": 21, "12": 25} + formula: "2k + 1 (even k) or 2k + 3 (odd k); more where no conference matrix exists at that order" + examples: {"6": 13, "7": 17, "8": 17, "10": 21, "12": 25, "21": 49, "22": 49} supports_models: - main_effects - "quadratic (limited active factors)" From 36941cf92e217e2b5d63c4e782c080ebc8663ca3 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 06:37:25 +0000 Subject: [PATCH 04/16] Cover the conference-matrix guard paths Tests for the two reachable guards (Paley on a non-prime-power, and a request past the largest order the search considers). The third, a missing irreducible polynomial, cannot happen over a finite field and is marked as such. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- .../experiments/designs_response_surface.py | 4 +++- tests/test_conference_matrices.py | 12 ++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index 9002bd9b..cc4535c0 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -463,7 +463,9 @@ def _irreducible_polynomial(p: int, n: int) -> list[int]: for candidate in _monic_polynomials(n, p): if _is_irreducible(candidate, p): return candidate - raise ValueError(f"No irreducible polynomial of degree {n} found over GF({p}).") + # Unreachable: an irreducible polynomial of every degree exists over every + # finite field. Kept so a bug in the search cannot return None silently. + raise ValueError(f"No irreducible polynomial of degree {n} found over GF({p}).") # pragma: no cover def _decode_field_element(element: int, p: int, n: int) -> list[int]: diff --git a/tests/test_conference_matrices.py b/tests/test_conference_matrices.py index 1bc6fd7d..23a93463 100644 --- a/tests/test_conference_matrices.py +++ b/tests/test_conference_matrices.py @@ -15,10 +15,12 @@ import pytest from process_improve.experiments.designs_response_surface import ( + _MAX_CONFERENCE_ORDER, _conference_matrix, _gf_multiplication_table, _irreducible_polynomial, _is_constructible_conference_order, + _paley_conference_matrix, _polynomial_remainder, _prime_power_factorization, _quadratic_character, @@ -239,6 +241,16 @@ def test_smallest_order_steps_up_only_when_it_must(self) -> None: assert _smallest_constructible_conference_order(28) == 28 assert _smallest_constructible_conference_order(34) == 38 + def test_beyond_the_search_bound_raises(self) -> None: + """A request past the largest order considered must fail loudly.""" + with pytest.raises(ValueError, match="No conference matrix could be constructed"): + _smallest_constructible_conference_order(_MAX_CONFERENCE_ORDER + 2) + + @pytest.mark.parametrize("q", [15, 21, 2, 4, 1, 0]) + def test_paley_rejects_anything_but_an_odd_prime_power(self, q: int) -> None: + with pytest.raises(ValueError, match="odd prime power"): + _paley_conference_matrix(q) + class TestConferenceMatrixValidation: """The guard that stops a bad matrix reaching a design.""" From d9d6c3fd8ecb7d401a20255f6d6927cf271fee09 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 06:38:20 +0000 Subject: [PATCH 05/16] Test the irreducible-polynomial guard rather than excluding it The repo's .coveragerc replaces coverage's default exclude_lines, so a bare "pragma: no cover" is not honoured here. Reach the guard with a stubbed irreducibility test instead. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- .../experiments/designs_response_surface.py | 7 ++++--- tests/test_conference_matrices.py | 14 ++++++++++++++ 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index cc4535c0..51d03a07 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -463,9 +463,10 @@ def _irreducible_polynomial(p: int, n: int) -> list[int]: for candidate in _monic_polynomials(n, p): if _is_irreducible(candidate, p): return candidate - # Unreachable: an irreducible polynomial of every degree exists over every - # finite field. Kept so a bug in the search cannot return None silently. - raise ValueError(f"No irreducible polynomial of degree {n} found over GF({p}).") # pragma: no cover + # Not reachable through the search above: an irreducible polynomial of every + # degree exists over every finite field. Kept so that a bug in the + # irreducibility test surfaces here instead of silently returning None. + raise ValueError(f"No irreducible polynomial of degree {n} found over GF({p}).") def _decode_field_element(element: int, p: int, n: int) -> list[int]: diff --git a/tests/test_conference_matrices.py b/tests/test_conference_matrices.py index 23a93463..16ca4ba0 100644 --- a/tests/test_conference_matrices.py +++ b/tests/test_conference_matrices.py @@ -160,6 +160,20 @@ def test_remainder_of_a_multiple_is_zero(self) -> None: product[i + j] = (product[i + j] + x * y) % p assert not any(_polynomial_remainder(product, g, p)) + def test_irreducible_polynomial_search_failure_is_loud(self, monkeypatch: pytest.MonkeyPatch) -> None: + """The guard behind the search must raise, not return None. + + An irreducible polynomial of every degree exists over every finite + field, so this path needs a broken irreducibility test to reach. A + silent None here would surface much later as a wrong design. + """ + monkeypatch.setattr( + "process_improve.experiments.designs_response_surface._is_irreducible", + lambda *_args, **_kwargs: False, + ) + with pytest.raises(ValueError, match="No irreducible polynomial"): + _irreducible_polynomial(3, 2) + @pytest.mark.parametrize(("p", "n"), [(3, 2), (3, 3), (3, 4), (5, 2), (5, 3), (7, 2), (11, 2)]) def test_irreducible_polynomial_has_no_roots(self, p: int, n: int) -> None: """A necessary condition, and sufficient for degree 2 and 3.""" From daa6f11f7a944b4f109c889e94c298fa65c89886 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 06:59:55 +0000 Subject: [PATCH 06/16] Correct the attribution for the order-16 conference matrix table Two errors in the previous comment. The 2015 MATLAB port is Copyright (c) 2015 Jacob Albrecht, not Bristol-Myers Squibb; BMS is his affiliation and the organisation named in the BSD-3 no-endorsement clause. The Python package the table was read from is BSD 3-Clause, Copyright (c) 2022 Daniele Ongari, not MIT. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- .../experiments/designs_response_surface.py | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index 51d03a07..c4943681 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -589,11 +589,21 @@ def _paley_conference_matrix(q: int) -> np.ndarray: # Conference matrix of order 16. 15 is not a prime power, so Paley's # construction does not reach this order, and it is the only order below 22 -# that needs a table. Taken from the conference-matrix catalogue of Xiao, Lin -# & Bai (2012), by way of Jacob Albrecht's BSD-3-licensed MATLAB port of the -# JMP add-in (Bristol-Myers Squibb, 2015) and its Python translation by Daniele -# Ongari in the `definitive_screening_design` package (MIT-spirited, same -# BSD-3 header retained). Verified here against C.T @ C == 15 * I before use. +# that needs a table. +# +# Provenance: the conference-matrix catalogue of Xiao, Lin & Bai (2012), +# reached by way of two ports. Jacob Albrecht (then at Bristol-Myers Squibb) +# ported the JMP add-in to MATLAB in 2015 and released it under the BSD +# 3-Clause licence, Copyright (c) 2015 Jacob Albrecht; the third clause of that +# licence names Bristol-Myers Squibb as the organisation whose name may not be +# used for endorsement. Daniele Ongari translated that MATLAB code to Python +# in the `definitive_screening_design` package (BSD 3-Clause, +# Copyright (c) 2022 Daniele Ongari), retaining Albrecht's header. +# +# Only the numeric table is reused here, and it is verified against +# C.T @ C == 15 * I before use rather than trusted: the order-10 table in the +# same upstream file has a single mistyped entry, which leaves its 9- and +# 10-factor designs non-orthogonal with nothing to catch it. _CONFERENCE_MATRIX_16 = np.array( [ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], From 8437d88e4d2614466ea9260ae38a1ba1b84d7496 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 07:08:06 +0000 Subject: [PATCH 07/16] Move the order-16 table attribution into a docstring The provenance, licensing, and the reason the table is verified rather than trusted now live in an attribute docstring on _CONFERENCE_MATRIX_16 itself, with an Attribution section and a reference to the originating paper. A comment block above an assignment is easy to detach from what it describes; a docstring travels with the object. _conference_matrix points at it. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- .../experiments/designs_response_surface.py | 58 +++++++++++++------ 1 file changed, 40 insertions(+), 18 deletions(-) diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index c4943681..9c035b4f 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -587,23 +587,6 @@ def _paley_conference_matrix(q: int) -> np.ndarray: return matrix -# Conference matrix of order 16. 15 is not a prime power, so Paley's -# construction does not reach this order, and it is the only order below 22 -# that needs a table. -# -# Provenance: the conference-matrix catalogue of Xiao, Lin & Bai (2012), -# reached by way of two ports. Jacob Albrecht (then at Bristol-Myers Squibb) -# ported the JMP add-in to MATLAB in 2015 and released it under the BSD -# 3-Clause licence, Copyright (c) 2015 Jacob Albrecht; the third clause of that -# licence names Bristol-Myers Squibb as the organisation whose name may not be -# used for endorsement. Daniele Ongari translated that MATLAB code to Python -# in the `definitive_screening_design` package (BSD 3-Clause, -# Copyright (c) 2022 Daniele Ongari), retaining Albrecht's header. -# -# Only the numeric table is reused here, and it is verified against -# C.T @ C == 15 * I before use rather than trusted: the order-10 table in the -# same upstream file has a single mistyped entry, which leaves its 9- and -# 10-factor designs non-orthogonal with nothing to catch it. _CONFERENCE_MATRIX_16 = np.array( [ [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], @@ -624,6 +607,43 @@ def _paley_conference_matrix(q: int) -> np.ndarray: [-1, 1, 1, -1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 0], ] ) +"""Skew-symmetric conference matrix of order 16, ``C.T @ C == 15 * I``. + +Needed because 15 is not a prime power, so :func:`_paley_conference_matrix` +cannot reach order 16. It is the only order below 22 that requires a table: +below it every order is covered by Paley's construction, and at 22 no +conference matrix exists at all. + +The table is verified against its defining property on every use, in +:func:`_conference_matrix`, rather than trusted. That is not ceremony. The +order-10 table in the same upstream file carries a single mistyped entry (row +17, column 0, ``+1`` where the foldover requires ``-1``), which leaves that +package's 9- and 10-factor designs correlated at r = 0.11 with nothing in place +to notice. + +Attribution +----------- +The matrix originates in the conference-matrix catalogue of Xiao, Lin & Bai +(2012), and reached this module through two ports: + +- Jacob Albrecht, then at Bristol-Myers Squibb, ported the JMP add-in to MATLAB + in March 2015 and released it under the BSD 3-Clause licence, + Copyright (c) 2015 Jacob Albrecht. The third clause of that licence names + Bristol-Myers Squibb as the organisation whose name may not be used to + endorse derived products. Bristol-Myers Squibb is Albrecht's affiliation + there, not the copyright holder. +- Daniele Ongari translated the MATLAB code to Python in the + ``definitive_screening_design`` package (BSD 3-Clause, + Copyright (c) 2022 Daniele Ongari), retaining Albrecht's header. + +Only the numeric table is reused here; none of the surrounding code was copied. + +References +---------- +.. [1] Xiao, L., Lin, D. K. J. and Bai, F. (2012). "Constructing definitive + screening designs using conference matrices." *Journal of Quality + Technology*, 44(1):2-8. +""" _TABULATED_CONFERENCE_MATRICES: dict[int, np.ndarray] = {16: _CONFERENCE_MATRIX_16} @@ -662,7 +682,9 @@ def _conference_matrix(m: int) -> tuple[np.ndarray, str]: Two constructions are available. Paley's covers every even order *m* whose predecessor ``m - 1`` is an odd prime power, which is most of them. A small - table covers order 16, the one gap below order 22. + table covers order 16, the one gap below order 22; see + :data:`_CONFERENCE_MATRIX_16` for where that table comes from and who holds + the copyright on the ports it travelled through. Parameters ---------- From 9380f06c3c7d8b362fe7634190f13f95019c74f8 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 07:19:15 +0000 Subject: [PATCH 08/16] Make the tool-registry test independent of sibling test execution test_registered_in_registry asserted that "test_dummy_add" and "test_dummy_mul" were present in the global registry, but those names are only registered when the sibling tests that apply the decorator actually run. Under pytest-xdist a sibling can be dispatched to a different worker process, whose registry this test never sees, so the assertion fails depending only on how the run happened to be split. The test now decorates its own function and asserts that name is registered, which is what the test is really about. Pre-existing on main and reproducible there; this branch surfaced it by adding tests, which changed the work split. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- tests/test_tool_spec.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/tests/test_tool_spec.py b/tests/test_tool_spec.py index 81c1f3e9..d4d076ac 100644 --- a/tests/test_tool_spec.py +++ b/tests/test_tool_spec.py @@ -93,9 +93,24 @@ def _plain(spec: _EmptyInput) -> dict: assert _plain._tool_spec["description"] == "Plain description." def test_registered_in_registry(self) -> None: - """Verify decorated tools are added to the global registry.""" - assert "test_dummy_add" in _TOOL_REGISTRY - assert "test_dummy_mul" in _TOOL_REGISTRY + """Verify decorated tools are added to the global registry. + + The tool is decorated here rather than relying on the names registered + by the sibling tests above. Those decorators only run when their own + test body runs, and under pytest-xdist a sibling can be dispatched to a + different worker process, whose registry this one never sees. + """ + + @tool_spec( + name="test_dummy_registered", + description="Add two numbers.", + input_model=_AddInput, + ) + def _registered(spec: _AddInput) -> dict: + return {"result": spec.a + spec.b} + + assert "test_dummy_registered" in _TOOL_REGISTRY + assert _TOOL_REGISTRY["test_dummy_registered"] is _registered def test_rejects_non_basemodel_input_model(self) -> None: """input_model must be a pydantic BaseModel subclass.""" From 94cc489720f67f9ffa5ddf7e9a8ea708d00bf2b2 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 08:14:03 +0000 Subject: [PATCH 09/16] Add two-level categorical factors to definitive screening designs Implements both column-augmentation procedures of Jones and Nachtsheim (2013). generate_design(design_type="dsd") previously raised as soon as a categorical factor was present, because the coded matrix was built as if every factor were continuous and the +/-1/0 values then failed the level check. A categorical factor takes a conference-matrix column like any other, so it costs the same runs as a continuous factor and contributes no quadratic term. The two zeros that column would carry have no meaning for a two-level factor, so each is resolved to a level, and extra centre runs are added because a centre run also needs a categorical setting. The methods differ in exactly that: DSD-augment sets z_j1 = b_j and z_j2 = -b_j, keeping the column balanced and every main effect unbiased by second-order effects; ORTH-augment sets both to +1, giving an orthogonal main-effects plan for up to four categorical factors at the price of aliasing main effects with categorical interactions. DSD-augment is the default. Its sign vectors are chosen by the determinant search the paper prescribes, not a fixed pattern. Verified against every oracle the paper provides: all 36 run sizes in Table 4, the constant-term alias row of Table 3 entry by entry, the information matrix of Equation (2), and Table 4's largest-alias column, which is 2/n throughout. The c = 2 panels of Figure 1 are not used as oracles. Both disagree with the rest of the paper: the DSD-augment panel gives the two categorical columns an inner product of 6, against the 2 of Equation (2) and the 0.1429 of Table 3, and is not the determinant-maximising choice its own Step 3 requires; the ORTH-augment panel prints a minus in run 10 where Step 3 requires a plus, leaving the column unbalanced and the main-effects plan non-orthogonal. Both findings are recorded as tests rather than prose. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 38 +- CITATION.cff | 2 +- docs/doe/coverage.md | 1 + pyproject.toml | 2 +- src/process_improve/experiments/designs.py | 14 +- .../experiments/designs_response_surface.py | 286 ++++++++++++- .../experiments/designs_utils.py | 20 + .../knowledge/data/design_types.yaml | 2 +- tests/test_dsd_categorical.py | 399 ++++++++++++++++++ 9 files changed, 742 insertions(+), 22 deletions(-) create mode 100644 tests/test_dsd_categorical.py diff --git a/CHANGELOG.md b/CHANGELOG.md index dc92c2f5..fa2f8a15 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,41 @@ those changes. ## [Unreleased] +## [1.61.0] - 2026-07-25 + +### Added + +- Definitive screening designs now accept **two-level categorical factors**, via both + column-augmentation procedures of Jones and Nachtsheim (2013). `generate_design(..., + design_type="dsd")` previously raised `ValueError` as soon as a categorical factor was + present. A new `categorical_method` argument selects the procedure: + - `"dsd"` (DSD-augment, the default) keeps the design *definitive*: every main effect + stays unbiased by every second-order effect, and two-factor interactions involving a + categorical factor stay clear of the main effects. The cost is small correlations + among the categorical main-effect columns, so the information matrix is not diagonal. + Adds two centre runs. The sign vectors are chosen to maximise the first-order + information determinant, as the paper prescribes. + - `"orth"` (ORTH-augment) gives an orthogonal linear main-effects plan for up to four + categorical factors, and nearly orthogonal beyond that. The cost is partial aliasing + between main effects and categorical interactions, an unbalanced categorical column, + and two extra runs relative to `"dsd"` when there are two or more categorical factors. + + A categorical factor occupies a conference-matrix column like any other, so it costs + the same runs as a continuous factor while contributing no quadratic term. Factors may + be given in any order; the construction needs the categorical columns trailing and + permutes them back. A categorical factor with three or more levels raises, with the + message pointing at `design_type="d_optimal"`. +- `dsd_centre_runs(n_categorical, categorical_method)` in + `process_improve.experiments.designs_response_surface`, and `dsd_run_count` now takes + `n_categorical` and `categorical_method`, so a design's size can still be predicted + without building it. + +### Fixed + +- `matrix_to_columns` now maps a numerically coded categorical column onto its level + labels. Design families differ in what they return for a categorical factor: the + optimal designs already produce labels, which continue to pass straight through. + ## [1.60.1] - 2026-07-25 ### Fixed @@ -2689,7 +2724,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.1...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.1...v1.61.0 [1.60.1]: https://github.com/kgdunn/process-improve/compare/v1.60.0...v1.60.1 [1.60.0]: https://github.com/kgdunn/process-improve/compare/v1.59.0...v1.60.0 [1.59.0]: https://github.com/kgdunn/process-improve/compare/v1.58.0...v1.59.0 diff --git a/CITATION.cff b/CITATION.cff index 63cba69e..8990021e 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,7 +12,7 @@ authors: repository-code: "https://github.com/kgdunn/process-improve" url: "https://kgdunn.github.io/process-improve/" license: MIT -version: 1.60.1 +version: 1.61.0 date-released: "2026-07-25" keywords: - chemometrics diff --git a/docs/doe/coverage.md b/docs/doe/coverage.md index c8fd0748..a0048831 100644 --- a/docs/doe/coverage.md +++ b/docs/doe/coverage.md @@ -55,6 +55,7 @@ All 14 metrics live in `experiments/evaluate.py` behind the `_METRIC_REGISTRY`: ## Caveats - **DSD conference matrix.** `_conference_matrix` in `designs_response_surface.py` builds order *m* from Paley's construction whenever `m − 1` is an odd prime power (`m ∈ {4, 6, 8, 10, 12, 14, 18, 20, 24, 26, 28, 30, 32, 38, 42, 44, 48, 50, 54, 60, 62, …}`), and from a tabulated matrix at `m = 16`. Every matrix is verified against `Cᵀ C = (m − 1) I` before use. At the remaining orders (`m ∈ {22, 34, 36, 40, 46, 52, 56, 58, …}`) the constructor raises, and `dispatch_dsd` steps up to the next order it can build, dropping the surplus columns. That costs runs but never orthogonality: 21 and 22 factors use order 24, so 49 runs rather than the 45 a minimal design would need. Order 22 is not a gap in this module but a fact about conference matrices: one of order `m ≡ 2 (mod 4)` exists only when `m − 1` is a sum of two squares, and 21 is not. Orders 34 and 58 fail the same test; 36, 40, 46 and 52 are orders this module has no construction for, though a matrix may exist. Use `dsd_run_count(k)` to find the size a given factor count will produce. +- **DSD categorical factors.** Two-level categorical factors are supported through both column-augmentation procedures of Jones & Nachtsheim (2013), selected with `categorical_method`. Three-or-more-level categorical factors are not covered by that family and raise, pointing at `d_optimal`. `ORTH-augment` is exactly orthogonal for up to four categorical factors and only nearly orthogonal beyond that, which is the paper's own limit rather than an implementation shortcut. Note that the c = 2 panels of the paper's Figure 1 disagree with its Table 3, Equation (2) and Section 3 Step 3; we follow the procedures and the tables, and `tests/test_dsd_categorical.py::TestFigure1Errata` records the evidence. - **Optimal designs without `pyoptex`.** D-optimal has a point-exchange fallback; I/A-optimal raise `ImportError` instead. Hard-to-change factors (split-plot structure) are currently ignored with a `logger.warning` when `pyoptex` is not available. - **Taguchi OA auto-selection** picks the smallest standard array that covers the requested factor count and level counts, but the underlying `pyDOE3.taguchi_design` requires the number of `levels_per_factor` entries to match the OA column count exactly, so requesting a Taguchi design with fewer factors than any available OA will fail. Users typically pick *k* to match one of the standard arrays (3, 4, 7, 11, 15, …). diff --git a/pyproject.toml b/pyproject.toml index 65a6a341..cdc0efac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.60.1" +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/designs.py b/src/process_improve/experiments/designs.py index cc79339f..019f7233 100644 --- a/src/process_improve/experiments/designs.py +++ b/src/process_improve/experiments/designs.py @@ -106,7 +106,7 @@ def _dispatch_dsd( ) -> tuple[np.ndarray, dict]: from process_improve.experiments.designs_response_surface import dispatch_dsd # noqa: PLC0415 - return dispatch_dsd(factors) + return dispatch_dsd(factors, categorical_method=kwargs.get("categorical_method", "dsd")) def _dispatch_omars( @@ -295,6 +295,7 @@ def generate_design( # noqa: PLR0913 constraints: list[Constraint] | None = None, hard_to_change: list[str] | None = None, model_type: str = "interactions", + categorical_method: str = "dsd", fixed_runs: pd.DataFrame | None = None, random_seed: int = 42, ) -> DesignResult: @@ -351,6 +352,16 @@ def generate_design( # noqa: PLR0913 the continuous factors only; the categorical enters as a main effect plus its interactions), since a categorical factor has no square. Ignored by the classical (non-optimal) design families. + categorical_method : {"dsd", "orth"} + For ``design_type="dsd"`` with two-level categorical factors, which + column-augmentation procedure of Jones and Nachtsheim (2013) to use. + ``"dsd"`` (default) keeps every main effect unbiased by second-order + effects and leaves categorical interactions clear of the main effects, + at the cost of small correlations among the categorical main-effect + columns. ``"orth"`` gives an orthogonal main-effects plan for up to four + categorical factors, at the cost of partial aliasing between main + effects and categorical interactions and two extra runs. Ignored by + every other design family. fixed_runs : pandas.DataFrame or None Runs to hold fixed while the optimizer fills the rest (design augmentation), for the optimal families only (``"d_optimal"``, ``"i_optimal"``, ``"a_optimal"``, which use @@ -415,6 +426,7 @@ def generate_design( # noqa: PLR0913 "hard_to_change": hard_to_change, "constraints": constraints, "model_type": model_type, + "categorical_method": categorical_method, "fixed_runs": fixed_runs, } diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index 9c035b4f..e2e6727a 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -22,10 +22,21 @@ ccdesign = _MissingExtra("pyDOE3", "expt") # type: ignore[assignment] if TYPE_CHECKING: + from collections.abc import Iterator + from process_improve.experiments.factor import Factor logger = logging.getLogger(__name__) +# Column-augmentation procedures for two-level categorical factors +# (Jones and Nachtsheim, 2013). +_CATEGORICAL_METHODS = frozenset({"dsd", "orth"}) + +# Above this many categorical factors the 2 ** (2c) DSD-augment search stops +# being cheap (2 ** 14 = 16384 determinants at c = 7), so a coordinate-style +# heuristic takes over, as the paper also does for large c. +_MAX_EXHAUSTIVE_CATEGORICAL = 7 + def dispatch_ccd( # noqa: PLR0913 factors: list[Factor], @@ -290,7 +301,7 @@ def dispatch_box_behnken( return coded_matrix, {} -def dispatch_dsd(factors: list[Factor]) -> tuple[np.ndarray, dict]: +def dispatch_dsd(factors: list[Factor], *, categorical_method: str = "dsd") -> tuple[np.ndarray, dict]: """Generate a Definitive Screening Design (DSD). Follows the conference-matrix foldover of Jones & Nachtsheim (2011): for a @@ -312,25 +323,59 @@ def dispatch_dsd(factors: list[Factor]) -> tuple[np.ndarray, dict]: ``C.T @ C == (m - 1) * I`` before it is used, so a degraded design can never reach the caller. + Two-level categorical factors are supported through the column-augmentation + procedures of Jones and Nachtsheim (2013). A categorical factor occupies a + conference-matrix column like any other, so it costs the same runs as a + continuous factor while contributing no quadratic term. The two zeros that + column would carry have no meaning for a two-level factor, so each is + resolved to a level, and extra centre runs are added because a centre run + also needs a categorical setting. The two methods differ in exactly that: + + ``"dsd"`` (DSD-augment, the default) + Sets ``z_{j,1} = b_j`` and ``z_{j,2} = -b_j`` at the foldover pair, so + the column stays balanced. Every main effect remains unbiased by any + second-order effect, which is what makes the design *definitive*, and + all two-factor interactions involving a categorical factor stay clear of + the main effects. The price is small correlations among the categorical + main-effect columns, so the information matrix is not diagonal. Adds two + centre runs. The sign vectors are chosen to maximise the first-order + information determinant, as the paper prescribes. + + ``"orth"`` (ORTH-augment) + Sets ``z_{j,1} = z_{j,2} = 1``, giving an orthogonal linear main-effects + plan for up to four categorical factors (nearly orthogonal beyond that). + The price is partial aliasing between main effects and interactions + involving the categorical factors, an unbalanced categorical column, and + two extra runs relative to ``"dsd"`` when there are two or more + categorical factors. + Parameters ---------- factors : list[Factor] - Continuous factors. At least three are required. + Factors, at least three. Categorical factors must have exactly two + levels; continuous factors are unrestricted. Factor order is preserved + in the returned matrix. + categorical_method : {"dsd", "orth"} + Which column-augmentation procedure to use. Ignored when no categorical + factor is present. Default ``"dsd"``. Returns ------- tuple[np.ndarray, dict] Coded design matrix and metadata. Metadata keys are ``"construction"`` (which conference-matrix construction was used), - ``"conference_order"`` (the order *m* of that matrix) and, when the - minimal order was not constructible, ``"minimal_conference_order"`` - (the order that would have been used had it existed). + ``"conference_order"`` (the order *m* of that matrix), ``"centre_runs"``, + ``"n_categorical"``, ``"categorical_method"`` (only when a categorical + factor is present) and, when the minimal order was not constructible, + ``"minimal_conference_order"`` (the order that would have been used had + it existed). Raises ------ ValueError - If fewer than three factors are supplied, or if no conference matrix - of usable order can be constructed. + If fewer than three factors are supplied, if a categorical factor does + not have exactly two levels, if *categorical_method* is unknown, or if + no conference matrix of usable order can be constructed. References ---------- @@ -340,21 +385,49 @@ def dispatch_dsd(factors: list[Factor]) -> tuple[np.ndarray, dict]: .. [2] Xiao, L., Lin, D. K. J. and Bai, F. (2012). "Constructing definitive screening designs using conference matrices." *Journal of Quality Technology*, 44(1):2-8. + .. [3] Jones, B. and Nachtsheim, C. J. (2013). "Definitive screening + designs with added two-level categorical factors." *Journal of Quality + Technology*, 45(2):121-129. """ k = len(factors) if k < 3: raise ValueError("Definitive Screening Designs require at least 3 factors.") + categorical_positions = [i for i, factor in enumerate(factors) if factor.type.value == "categorical"] + _validate_categorical_factors(factors, categorical_positions) + n_categorical = len(categorical_positions) + minimal_order = k if k % 2 == 0 else k + 1 order = _smallest_constructible_conference_order(minimal_order) conference, construction = _conference_matrix(order) - zero_row = np.zeros((1, order)) - coded_matrix = np.vstack([conference, -conference, zero_row]) + body = np.vstack([conference, -conference]) if order > k: - coded_matrix = coded_matrix[:, :k] + body = body[:, :k] + + if n_categorical == 0: + coded_matrix = np.vstack([body, np.zeros((1, k))]) + centre_runs = 1 + else: + # The construction needs the categorical columns last; permute back afterwards. + continuous_positions = [i for i in range(k) if i not in set(categorical_positions)] + internal_order = continuous_positions + categorical_positions + coded_matrix, centre_runs = _augment_with_categorical_columns( + body[:, internal_order], + n_continuous=k - n_categorical, + method=categorical_method, + ) + inverse = np.argsort(internal_order) + coded_matrix = coded_matrix[:, inverse] - meta = {"construction": construction, "conference_order": order} + meta = { + "construction": construction, + "conference_order": order, + "centre_runs": centre_runs, + "n_categorical": n_categorical, + } + if n_categorical: + meta["categorical_method"] = categorical_method if order > minimal_order: meta["minimal_conference_order"] = minimal_order logger.info( @@ -363,13 +436,146 @@ def dispatch_dsd(factors: list[Factor]) -> tuple[np.ndarray, dict]: minimal_order, k, order, - 2 * order + 1, - 2 * minimal_order + 1, + coded_matrix.shape[0], + 2 * minimal_order + centre_runs, ) return coded_matrix, meta +def _validate_categorical_factors(factors: list[Factor], categorical_positions: list[int]) -> None: + """Reject categorical factors the Jones-Nachtsheim (2013) construction cannot carry. + + The column-augmentation methods are defined for **two-level** categorical + factors only. A factor with three or more levels has no representation as a + single ±1 column, so it needs a different design family. + + Raises + ------ + ValueError + If any categorical factor does not have exactly two levels. + """ + for position in categorical_positions: + factor = factors[position] + n_levels = len(factor.levels or []) + if n_levels != 2: + raise ValueError( + f"Factor {factor.name!r} has {n_levels} levels. Definitive screening designs accept " + "two-level categorical factors only (Jones and Nachtsheim, 2013). For a factor with " + 'three or more levels use design_type="d_optimal" (or "i_optimal"), which places ' + "multi-level categorical factors without this restriction." + ) + + +def _foldover_zero_rows(body: np.ndarray, column: int, half: int) -> tuple[int, int]: + """Locate the foldover pair holding the two zeros of *column*. + + In ``[C; -C]`` a conference-matrix column has exactly one zero in each half, + at mirrored row indices. Those two rows are the foldover pair whose entries + the augmentation replaces with ``z_{j,1}`` and ``z_{j,2}``. + """ + (upper,) = np.flatnonzero(body[:half, column] == 0) + return int(upper), int(upper + half) + + +def _orth_centre_block(n_categorical: int, n_centre: int) -> np.ndarray: + """Build the ORTH-augment centre-run block ``B``. + + Jones and Nachtsheim (2013), Section 3, Step 1: column ``b_j`` has a one in + the ``(5 - j)``th row and negative ones elsewhere, for ``j = 1, ..., c <= 4``. + Beyond four categorical factors the pattern is cycled, which the paper + describes as giving nearly (rather than exactly) orthogonal plans. + """ + block = -np.ones((n_centre, n_categorical)) + if n_centre >= 4: + for j in range(n_categorical): + block[n_centre - 1 - (j % 4), j] = 1.0 + return block + + +def _augment_with_categorical_columns( + body: np.ndarray, + *, + n_continuous: int, + method: str, +) -> tuple[np.ndarray, int]: + """Turn the trailing columns of a folded DSD into two-level categorical columns. + + Implements both procedures of Jones and Nachtsheim (2013). *body* is the + ``[C; -C]`` foldover with the categorical factors already permuted to the + trailing columns; the returned matrix carries the appended centre runs. + + Returns + ------- + tuple[np.ndarray, int] + The augmented design and the number of centre runs appended. + """ + if method not in _CATEGORICAL_METHODS: + raise ValueError(f"categorical_method must be one of {sorted(_CATEGORICAL_METHODS)}; got {method!r}.") + + n_runs, k = body.shape + half = n_runs // 2 + n_categorical = k - n_continuous + zero_rows = [_foldover_zero_rows(body, n_continuous + j, half) for j in range(n_categorical)] + + if method == "orth": + # Step 3: z_{j,1} = z_{j,2} = 1 in every remaining augmented column. + design = body.copy() + for j, (upper, lower) in enumerate(zero_rows): + design[upper, n_continuous + j] = 1.0 + design[lower, n_continuous + j] = 1.0 + n_centre = 2 if n_categorical == 1 else 4 + centre = np.zeros((n_centre, k)) + centre[:, n_continuous:] = _orth_centre_block(n_categorical, n_centre) + return np.vstack([design, centre]), n_centre + + # DSD-augment. Two independent sign choices per categorical factor: z_j at + # the foldover zero pair, and b_j in the two centre runs. Section 2, Step 3 + # searches all 2 ** (2c) combinations and keeps the one maximising the + # first-order information determinant. + best_design = body + best_criterion = -np.inf + for signs in _categorical_sign_candidates(n_categorical): + z_signs, b_signs = signs[:n_categorical], signs[n_categorical:] + design = body.copy() + for j, (upper, lower) in enumerate(zero_rows): + design[upper, n_continuous + j] = z_signs[j] + design[lower, n_continuous + j] = -z_signs[j] + centre = np.zeros((2, k)) + centre[0, n_continuous:] = b_signs + centre[1, n_continuous:] = -b_signs + candidate = np.vstack([design, centre]) + model = np.column_stack([np.ones(candidate.shape[0]), candidate]) + criterion = float(np.linalg.slogdet(model.T @ model)[1]) + if criterion > best_criterion + 1e-12: + best_criterion, best_design = criterion, candidate + + return best_design, 2 + + +def _categorical_sign_candidates(n_categorical: int) -> Iterator[np.ndarray]: + """Yield candidate ``(z, b)`` sign vectors for the DSD-augment search. + + Exhaustive over all ``2 ** (2c)`` combinations while that is cheap. Beyond + :data:`_MAX_EXHAUSTIVE_CATEGORICAL` factors the paper switches to a + coordinate-exchange heuristic; here the search is seeded with the alternating + pattern and each coordinate is offered in both signs, which keeps the cost + linear in *c* at the price of no longer being exhaustive. + """ + length = 2 * n_categorical + if n_categorical <= _MAX_EXHAUSTIVE_CATEGORICAL: + for bits in range(2**length): + yield np.array([1.0 if (bits >> i) & 1 else -1.0 for i in range(length)]) + return + + seed = np.array([1.0 if i % 2 == 0 else -1.0 for i in range(length)]) + yield seed + for i in range(length): + flipped = seed.copy() + flipped[i] = -flipped[i] + yield flipped + + def _is_prime(n: int) -> bool: """Return True iff *n* is a (positive) prime.""" if n < 2: @@ -780,7 +986,7 @@ def dsd_conference_order(n_factors: int) -> int: return _smallest_constructible_conference_order(n_factors if n_factors % 2 == 0 else n_factors + 1) -def dsd_run_count(n_factors: int) -> int: +def dsd_run_count(n_factors: int, n_categorical: int = 0, categorical_method: str = "dsd") -> int: """Return the number of runs in the definitive screening design for *n_factors* factors. This is ``2m + 1`` for the conference-matrix order *m* returned by @@ -791,16 +997,62 @@ def dsd_run_count(n_factors: int) -> int: Parameters ---------- n_factors : int - Number of factors, at least three. + Total number of factors, at least three, counting categorical factors. + n_categorical : int + How many of those are two-level categorical factors. They add centre + runs; see :func:`dsd_centre_runs`. + categorical_method : {"dsd", "orth"} + Which augmentation procedure is in use. Returns ------- int - Run count, including the design's single centre run. + Run count, including the design's centre runs. Examples -------- >>> dsd_run_count(6), dsd_run_count(7), dsd_run_count(22) (13, 17, 49) + + Two-level categorical factors count toward *n_factors* and add centre runs: + + >>> dsd_run_count(6, n_categorical=2) # 4 continuous + 2 categorical + 14 + >>> dsd_run_count(6, n_categorical=2, categorical_method="orth") + 16 + """ + return 2 * dsd_conference_order(n_factors) + dsd_centre_runs(n_categorical, categorical_method) + + +def dsd_centre_runs(n_categorical: int, categorical_method: str = "dsd") -> int: + """Return how many centre runs a definitive screening design carries. + + A design with only continuous factors needs one centre run. A categorical + factor has no centre level, so each centre run must still be assigned a + level, and more than one centre run is needed for the levels to balance. + The counts follow Jones and Nachtsheim (2013) and reproduce every run size + in their Table 4. + + Parameters + ---------- + n_categorical : int + Number of two-level categorical factors. + categorical_method : {"dsd", "orth"} + Which augmentation procedure is in use. Only matters for two or more + categorical factors. + + Returns + ------- + int + 1 with no categorical factor, 2 with exactly one (either method), then + 2 for ``"dsd"`` and 4 for ``"orth"``. """ - return 2 * dsd_conference_order(n_factors) + 1 + if categorical_method not in _CATEGORICAL_METHODS: + raise ValueError( + f"categorical_method must be one of {sorted(_CATEGORICAL_METHODS)}; got {categorical_method!r}." + ) + if n_categorical == 0: + return 1 + if n_categorical == 1: + return 2 + return 2 if categorical_method == "dsd" else 4 diff --git a/src/process_improve/experiments/designs_utils.py b/src/process_improve/experiments/designs_utils.py index 6e54491e..39b6d6bc 100644 --- a/src/process_improve/experiments/designs_utils.py +++ b/src/process_improve/experiments/designs_utils.py @@ -41,6 +41,25 @@ def matrix_to_columns( """ from process_improve.experiments.factor import FactorType # noqa: PLC0415 + def _categorical_labels(values: list, levels: list | None) -> list: + """Translate a numerically coded categorical column into its level labels. + + Design families differ in what they hand back for a categorical factor. + The optimal designs (via pyoptex) already produce labels, which pass + straight through. The definitive screening design produces a two-level + ``-1`` / ``+1`` column, which is mapped onto the factor's two levels in + order. Anything else is left alone for the ``Column`` constructor to + validate and reject if it does not belong. + """ + if not levels or not all(isinstance(v, (int, float)) and not isinstance(v, bool) for v in values): + return values + distinct = sorted({float(v) for v in values}) + if len(levels) == 2 and distinct == [-1.0, 1.0]: + return [levels[0] if float(v) < 0 else levels[1] for v in values] + if distinct and set(distinct) <= set(range(len(levels))): + return [levels[int(v)] for v in values] + return values + columns: list[Column] = [] for i, factor in enumerate(factors): values = matrix[:, i].tolist() @@ -49,6 +68,7 @@ def matrix_to_columns( # from its levels and mark it not-coded so the coded<->actual affine # map (which assumes a numeric low/high range) is skipped and the # labels pass straight through to both the coded and actual designs. + values = _categorical_labels(values, factor.levels) col = c(values, name=factor.name, levels=factor.levels, units=factor.units) col.pi_is_coded = False else: diff --git a/src/process_improve/experiments/knowledge/data/design_types.yaml b/src/process_improve/experiments/knowledge/data/design_types.yaml index 8f36919c..f241dd89 100644 --- a/src/process_improve/experiments/knowledge/data/design_types.yaml +++ b/src/process_improve/experiments/knowledge/data/design_types.yaml @@ -303,7 +303,7 @@ supports_quadratic: "Partial (if <= 3 active factors)" supports_interaction: "Partial (if <= 3 active factors)" min_runs: - formula: "2k + 1 (even k) or 2k + 3 (odd k); more where no conference matrix exists at that order" + formula: "2k + 1 (even k) or 2k + 3 (odd k); more where no conference matrix exists at that order, and one or three extra centre runs when two-level categorical factors are present" examples: {"6": 13, "7": 17, "8": 17, "10": 21, "12": 25, "21": 49, "22": 49} supports_models: - main_effects diff --git a/tests/test_dsd_categorical.py b/tests/test_dsd_categorical.py new file mode 100644 index 00000000..2e09471d --- /dev/null +++ b/tests/test_dsd_categorical.py @@ -0,0 +1,399 @@ +"""Tests for two-level categorical factors in definitive screening designs. + +Every oracle here comes from Jones, B. and Nachtsheim, C. J. (2013), +"Definitive screening designs with added two-level categorical factors", +*Journal of Quality Technology*, 45(2):121-129: + +- **Table 4** lists run sizes for both procedures over 4 to 12 continuous + factors and 1 to 4 categorical factors: 36 published values. +- **Table 4** also lists the largest entry of the constant-term row of the + alias matrix for the DSD-augment designs, which works out to ``2 / n``. +- **Table 3** gives that alias matrix in full for four continuous and two + categorical factors: zero against every continuous-continuous interaction, + 0.1429 against every interaction involving a categorical factor. +- **Section 2** defines the DSD-augment procedure and its defining property, + that main effects stay unbiased by every second-order effect. +- **Section 3** defines the ORTH-augment procedure, which gives an orthogonal + linear main-effects plan for up to four categorical factors. + +Two entries in the paper's Figure 1 disagree with the rest of the paper and are +not used as oracles; see ``TestFigure1Errata`` for the evidence. +""" + +from __future__ import annotations + +import itertools + +import numpy as np +import pytest + +from process_improve.experiments.designs import generate_design +from process_improve.experiments.designs_response_surface import ( + dispatch_dsd, + dsd_centre_runs, + dsd_conference_order, + dsd_run_count, +) +from process_improve.experiments.factor import Factor + +# Jones & Nachtsheim (2013) Table 4: (n_DSD, n_ORTH, m continuous, c categorical). +TABLE_4 = [ + (14, 14, 4, 1), (14, 16, 4, 2), (18, 20, 4, 3), (18, 20, 4, 4), + (14, 14, 5, 1), (18, 20, 5, 2), (18, 20, 5, 3), (22, 24, 5, 4), + (18, 18, 6, 1), (18, 20, 6, 2), (22, 24, 6, 3), (22, 24, 6, 4), + (18, 18, 7, 1), (22, 24, 7, 2), (22, 24, 7, 3), (26, 28, 7, 4), + (22, 22, 8, 1), (22, 24, 8, 2), (26, 28, 8, 3), (26, 28, 8, 4), + (22, 22, 9, 1), (26, 28, 9, 2), (26, 28, 9, 3), (30, 32, 9, 4), + (26, 26, 10, 1), (26, 28, 10, 2), (30, 32, 10, 3), (30, 32, 10, 4), + (26, 26, 11, 1), (30, 32, 11, 2), (30, 32, 11, 3), (34, 36, 11, 4), + (30, 30, 12, 1), (30, 32, 12, 2), (34, 36, 12, 3), (34, 36, 12, 4), +] # fmt: skip + +_TABLE_4_IDS = [f"m{m}c{c}" for _, _, m, c in TABLE_4] + + +def _factors(n_continuous: int, n_categorical: int, *, interleave: bool = False) -> list[Factor]: + """Build a factor list, optionally with the categorical factors not last.""" + continuous = [Factor(name=f"X{i + 1}", low=-1.0, high=1.0) for i in range(n_continuous)] + categorical = [ + Factor(name=f"C{j + 1}", type="categorical", levels=[f"lo{j + 1}", f"hi{j + 1}"]) for j in range(n_categorical) + ] + if not interleave: + return continuous + categorical + merged: list[Factor] = [] + for i in range(max(len(continuous), len(categorical))): + if i < len(categorical): + merged.append(categorical[i]) + if i < len(continuous): + merged.append(continuous[i]) + return merged + + +def _model_matrices(design: np.ndarray, n_continuous: int) -> tuple[np.ndarray, np.ndarray, list[str]]: + """Return the first-order matrix, the interaction matrix, and interaction labels. + + The model is Equation (1) of the paper: intercept, linear main effects for + every factor, all two-factor interactions, and quadratics for the continuous + factors only (a two-level categorical factor has no square). + """ + n, k = design.shape + first_order = np.column_stack([np.ones(n), design]) + labels = [f"{a + 1}{b + 1}" for a, b in itertools.combinations(range(k), 2)] + interactions = np.column_stack([design[:, a] * design[:, b] for a, b in itertools.combinations(range(k), 2)]) + return first_order, interactions, labels + + +def _second_order(design: np.ndarray, n_continuous: int) -> np.ndarray: + """All second-order terms: quadratics on the continuous factors, plus every interaction.""" + _n, k = design.shape + quadratics = [design[:, j] ** 2 for j in range(n_continuous)] + interactions = [design[:, a] * design[:, b] for a, b in itertools.combinations(range(k), 2)] + return np.column_stack(quadratics + interactions) + + +def _alias_matrix(first_order: np.ndarray, terms: np.ndarray) -> np.ndarray: + """Return the alias matrix ``(X1' X1)^-1 X1' X2``.""" + return np.linalg.solve(first_order.T @ first_order, first_order.T @ terms) + + +class TestTable4RunSizes: + """The 36 published run sizes, for both procedures.""" + + @pytest.mark.parametrize(("n_dsd", "_n_orth", "m", "c"), TABLE_4, ids=_TABLE_4_IDS) + def test_dsd_augment_run_size(self, n_dsd: int, _n_orth: int, m: int, c: int) -> None: + design, meta = dispatch_dsd(_factors(m, c), categorical_method="dsd") + assert design.shape == (n_dsd, m + c) + assert meta["categorical_method"] == "dsd" + assert meta["n_categorical"] == c + + @pytest.mark.parametrize(("_n_dsd", "n_orth", "m", "c"), TABLE_4, ids=_TABLE_4_IDS) + def test_orth_augment_run_size(self, _n_dsd: int, n_orth: int, m: int, c: int) -> None: + design, meta = dispatch_dsd(_factors(m, c), categorical_method="orth") + assert design.shape == (n_orth, m + c) + assert meta["categorical_method"] == "orth" + + @pytest.mark.parametrize(("n_dsd", "n_orth", "m", "c"), TABLE_4, ids=_TABLE_4_IDS) + def test_run_count_helper_predicts_the_generated_size(self, n_dsd: int, n_orth: int, m: int, c: int) -> None: + """The planning helper must agree with the constructor, not approximate it.""" + assert dsd_run_count(m + c, n_categorical=c, categorical_method="dsd") == n_dsd + assert dsd_run_count(m + c, n_categorical=c, categorical_method="orth") == n_orth + + def test_orth_costs_two_extra_runs_only_beyond_one_categorical(self) -> None: + """Table 4 shows n_ORTH == n_DSD for every single-categorical row.""" + for n_dsd, n_orth, _m, c in TABLE_4: + if c == 1: + assert n_orth == n_dsd + else: + assert n_orth == n_dsd + 2 + + def test_centre_run_counts(self) -> None: + assert dsd_centre_runs(0) == 1 + assert dsd_centre_runs(1) == 2 + assert dsd_centre_runs(1, "orth") == 2 + assert dsd_centre_runs(3) == 2 + assert dsd_centre_runs(3, "orth") == 4 + + def test_a_categorical_factor_costs_a_conference_column(self) -> None: + """Nine continuous plus three categorical is a twelve-factor design, not a nine-factor one.""" + assert dsd_conference_order(9 + 3) == dsd_conference_order(12) + design, meta = dispatch_dsd(_factors(9, 3)) + assert design.shape[1] == 12 + assert meta["conference_order"] == dsd_conference_order(12) + + +class TestDsdAugmentProperties: + """Section 2: the design stays *definitive*.""" + + @pytest.mark.parametrize(("m", "c"), [(4, 1), (4, 2), (4, 3), (4, 4), (5, 2), (6, 3), (8, 4), (12, 2)]) + def test_main_effects_unbiased_by_every_second_order_effect(self, m: int, c: int) -> None: + """The defining property: main-effect rows of the alias matrix are exactly zero. + + This is what "definitive" means, and it is what separates DSD-augment + from ORTH-augment. + """ + design, _meta = dispatch_dsd(_factors(m, c), categorical_method="dsd") + first_order, _inter, _labels = _model_matrices(design, m) + alias = _alias_matrix(first_order, _second_order(design, m)) + assert np.abs(alias[1:, :]).max() < 1e-12 + + @pytest.mark.parametrize(("n_dsd", "_n_orth", "m", "c"), TABLE_4, ids=_TABLE_4_IDS) + def test_largest_constant_term_alias_is_two_over_n(self, n_dsd: int, _n_orth: int, m: int, c: int) -> None: + """Table 4's rightmost column: the largest entry of the alias matrix's first row. + + Across all 36 published rows this equals ``2 / n`` exactly. It is the + check that is sensitive to the sign vectors the search picks, which the + run-size check is not. + """ + design, _meta = dispatch_dsd(_factors(m, c), categorical_method="dsd") + first_order, interactions, _labels = _model_matrices(design, m) + alias = _alias_matrix(first_order, interactions) + assert np.abs(alias[0]).max() == pytest.approx(2 / n_dsd, abs=5e-5) + + def test_table_3_alias_row_reproduced_entry_by_entry(self) -> None: + """Table 3: the full constant-term alias row for four continuous, two categorical. + + Zero against each of 12, 13, 14, 23, 24 and 34; 0.1429 against each of + the nine interactions that involve a categorical factor, including the + categorical-by-categorical 56. + """ + design, _meta = dispatch_dsd(_factors(4, 2), categorical_method="dsd") + first_order, interactions, labels = _model_matrices(design, 4) + alias = dict(zip(labels, np.abs(_alias_matrix(first_order, interactions)[0]), strict=True)) + + for continuous_pair in ("12", "13", "14", "23", "24", "34"): + assert alias[continuous_pair] == pytest.approx(0.0, abs=1e-12), continuous_pair + for categorical_pair in ("15", "16", "25", "26", "35", "36", "45", "46", "56"): + assert alias[categorical_pair] == pytest.approx(0.1429, abs=5e-5), categorical_pair + + def test_information_matrix_matches_equation_2(self) -> None: + """Equation (2): X'X for four continuous and two categorical factors. + + 14 on the intercept, 10 on each continuous main effect (four of the + fourteen runs are at that factor's centre), 14 on each categorical main + effect, zero between the intercept and every main effect, and entries of + magnitude 2 coupling the categorical columns to the rest. + """ + design, _meta = dispatch_dsd(_factors(4, 2), categorical_method="dsd") + first_order, _inter, _labels = _model_matrices(design, 4) + info = first_order.T @ first_order + + assert info[0, 0] == 14 + assert list(np.diag(info)[1:5]) == [10] * 4 + assert list(np.diag(info)[5:]) == [14] * 2 + assert np.abs(info[0, 1:]).max() == 0 # intercept orthogonal to every main effect + coupling = np.abs(info[1:, 1:] - np.diag(np.diag(info[1:, 1:]))) + assert set(np.unique(coupling)) <= {0.0, 2.0} + + @pytest.mark.parametrize(("m", "c"), [(4, 1), (4, 2), (5, 3), (6, 4), (8, 2)]) + def test_categorical_columns_are_balanced(self, m: int, c: int) -> None: + """z_{j,1} = b_j and z_{j,2} = -b_j keeps the foldover pair balanced.""" + design, _meta = dispatch_dsd(_factors(m, c), categorical_method="dsd") + for column in design[:, m:].T: + assert (column == 1).sum() == (column == -1).sum() + assert set(np.unique(column)) == {-1.0, 1.0} + + @pytest.mark.parametrize(("m", "c"), [(4, 2), (5, 3), (8, 2)]) + def test_information_matrix_is_not_diagonal(self, m: int, c: int) -> None: + """The documented cost of DSD-augment. + + The paper states plainly that "the main effects columns for categorical + factors exhibit small correlations so that the information matrix is not + diagonal". Asserting it keeps the trade-off honest rather than letting a + future change quietly claim orthogonality it does not have. + """ + design, _meta = dispatch_dsd(_factors(m, c), categorical_method="dsd") + first_order, _inter, _labels = _model_matrices(design, m) + info = first_order.T @ first_order + off_diagonal = np.abs(info - np.diag(np.diag(info))) + assert off_diagonal.max() > 0 + assert off_diagonal.max() == 2 # never larger than 2, per Section 2 + + +class TestOrthAugmentProperties: + """Section 3: an orthogonal linear main-effects plan for up to four categorical factors.""" + + @pytest.mark.parametrize(("m", "c"), [(4, 1), (4, 2), (4, 3), (4, 4), (5, 2), (6, 3), (8, 4), (12, 4)]) + def test_main_effects_plan_is_exactly_orthogonal(self, m: int, c: int) -> None: + design, _meta = dispatch_dsd(_factors(m, c), categorical_method="orth") + first_order, _inter, _labels = _model_matrices(design, m) + info = first_order.T @ first_order + assert np.abs(info - np.diag(np.diag(info))).max() < 1e-9 + + @pytest.mark.parametrize(("m", "c"), [(4, 2), (5, 3), (6, 4)]) + def test_main_effects_are_partially_aliased_with_categorical_interactions(self, m: int, c: int) -> None: + """The documented cost of ORTH-augment, and the reason it is not the default. + + The paper: "all of the main effects may be biased by potential two-factor + interactions involving categorical factors", with entries ranging up to + about 0.4. + """ + design, _meta = dispatch_dsd(_factors(m, c), categorical_method="orth") + first_order, interactions, _labels = _model_matrices(design, m) + alias = _alias_matrix(first_order, interactions) + assert np.abs(alias[1:, :]).max() > 0 + assert np.abs(alias[1:, :]).max() <= 0.4 + 1e-9 + + def test_beyond_four_categorical_factors_orthogonality_is_only_approximate(self) -> None: + """Section 3 promises exact orthogonality only up to c = 4, "nearly" beyond. + + The design must still be usable, so the correlation is bounded rather + than absent. + """ + design, _meta = dispatch_dsd(_factors(6, 6), categorical_method="orth") + first_order, _inter, _labels = _model_matrices(design, 6) + info = first_order.T @ first_order + scale = np.sqrt(np.outer(np.diag(info), np.diag(info))) + correlation = np.abs(info / scale) + np.fill_diagonal(correlation, 0.0) + assert correlation.max() < 0.35 + + +class TestFigure1Errata: + """Two entries of the paper's Figure 1 disagree with the rest of the paper. + + Recorded here so a future reader who checks our output against the printed + figure knows why the c = 2 panels differ, and so the reasoning is testable + rather than asserted in prose. + """ + + def test_dsd_c2_panel_is_not_d_optimal(self) -> None: + """Figure 1(a), c = 2 contradicts Table 3, Equation (2), and Step 3. + + As printed, the two categorical columns have inner product 6, giving a + constant-term alias of 6/14 = 0.4286 against the 56 interaction. Table 3 + gives 0.1429 there and Equation (2) gives an inner product of magnitude + 2. Section 2, Step 3 resolves it: the sign vectors are chosen to + maximise the first-order information determinant, and the choice that + does so is the one matching Table 3, not the one printed. + """ + design, _meta = dispatch_dsd(_factors(4, 2), categorical_method="dsd") + ours = float(design[:, 4] @ design[:, 5]) + assert abs(ours) == 2 # Equation (2), not the 6 the figure implies + + # The printed alternative is a valid design, just a worse one: same run + # size, larger aliasing, smaller determinant. + printed = design.copy() + printed[:, 5] = np.where(np.abs(printed[:, 4] - printed[:, 5]) > 0, -printed[:, 5], printed[:, 5]) + model_ours = np.column_stack([np.ones(design.shape[0]), design]) + model_printed = np.column_stack([np.ones(printed.shape[0]), printed]) + assert np.linalg.slogdet(model_ours.T @ model_ours)[1] >= np.linalg.slogdet(model_printed.T @ model_printed)[1] + + def test_orth_c2_panel_breaks_its_own_step_3(self) -> None: + """Figure 1(b), c = 2, run 10 prints a minus in the first categorical column. + + Runs 9 and 10 are the foldover pair carrying z_{1,1} and z_{1,2}, and + Section 3, Step 3 requires both to be +1. With the minus as printed the + column is unbalanced 7 to 9 and the main-effects plan is not orthogonal, + contradicting the procedure's whole purpose. Following Step 3 gives a + balanced column and exact orthogonality. + """ + design, _meta = dispatch_dsd(_factors(4, 2), categorical_method="orth") + for column in design[:, 4:].T: + assert (column == 1).sum() == (column == -1).sum() == 8 + first_order, _inter, _labels = _model_matrices(design, 4) + info = first_order.T @ first_order + assert np.abs(info - np.diag(np.diag(info))).max() < 1e-9 + + +class TestFactorOrderAndValidation: + """API behaviour around the construction.""" + + @pytest.mark.parametrize("method", ["dsd", "orth"]) + def test_categorical_factors_need_not_be_last(self, method: str) -> None: + """The construction needs them trailing; the caller should not have to care.""" + factors = _factors(4, 2, interleave=True) + design, _meta = dispatch_dsd(factors, categorical_method=method) + assert [f.name for f in factors] == ["C1", "X1", "C2", "X2", "X3", "X4"] + for position, factor in enumerate(factors): + column = design[:, position] + if factor.type.value == "categorical": + assert set(np.unique(column)) == {-1.0, 1.0}, f"{factor.name} should be two-level" + else: + assert 0.0 in set(np.unique(column)), f"{factor.name} should carry a centre level" + + @pytest.mark.parametrize("method", ["dsd", "orth"]) + def test_reordering_preserves_the_design(self, method: str) -> None: + """Permuting the factor list permutes the columns and nothing else.""" + ordered, _ = dispatch_dsd(_factors(4, 2), categorical_method=method) + interleaved, _ = dispatch_dsd(_factors(4, 2, interleave=True), categorical_method=method) + # Interleaved order is C1, X1, C2, X2, X3, X4 -> map back to X1..X4, C1, C2. + remapped = interleaved[:, [1, 3, 4, 5, 0, 2]] + assert np.array_equal(np.sort(ordered, axis=0), np.sort(remapped, axis=0)) + + def test_three_level_categorical_is_rejected_with_a_pointer(self) -> None: + factors = [ + *[Factor(name=f"X{i}", low=-1, high=1) for i in range(3)], + Factor(name="Catalyst", type="categorical", levels=["A", "B", "C"]), + ] + with pytest.raises(ValueError, match="two-level categorical factors only"): + dispatch_dsd(factors) + with pytest.raises(ValueError, match="d_optimal"): + dispatch_dsd(factors) + + def test_unknown_method_is_rejected(self) -> None: + with pytest.raises(ValueError, match="categorical_method"): + dispatch_dsd(_factors(4, 2), categorical_method="orthogonal") + with pytest.raises(ValueError, match="categorical_method"): + dsd_centre_runs(2, "nonsense") + + def test_method_is_ignored_without_categorical_factors(self) -> None: + """An all-continuous DSD is unchanged by the argument.""" + a, meta_a = dispatch_dsd(_factors(6, 0), categorical_method="dsd") + b, meta_b = dispatch_dsd(_factors(6, 0), categorical_method="orth") + assert np.array_equal(a, b) + assert "categorical_method" not in meta_a + assert meta_b["centre_runs"] == 1 + + +class TestGenerateDesignIntegration: + """The user-facing path.""" + + def test_categorical_levels_appear_in_the_actual_design(self) -> None: + factors = [ + Factor(name="Temp", low=150, high=200, units="degC"), + Factor(name="Press", low=1, high=5, units="bar"), + Factor(name="Time", low=10, high=60, units="min"), + Factor(name="Solvent", type="categorical", levels=["MeOH", "EtOH"]), + ] + result = generate_design(factors, design_type="dsd", center_points=0) + assert result.n_runs == dsd_run_count(4, n_categorical=1) + solvent = result.design_actual["Solvent"] + assert set(solvent.unique()) == {"MeOH", "EtOH"} + assert set(result.design_actual["Temp"].unique()) == {150.0, 175.0, 200.0} + + def test_method_is_threaded_through_and_recorded(self) -> None: + factors = _factors(4, 2) + default = generate_design(factors, design_type="dsd", center_points=0) + orth = generate_design(factors, design_type="dsd", center_points=0, categorical_method="orth") + assert default.metadata["categorical_method"] == "dsd" + assert orth.metadata["categorical_method"] == "orth" + assert default.n_runs == 14 + assert orth.n_runs == 16 + + def test_three_level_categorical_message_reaches_the_caller(self) -> None: + factors = [ + *[Factor(name=f"X{i}", low=-1, high=1) for i in range(3)], + Factor(name="Catalyst", type="categorical", levels=["A", "B", "C"]), + ] + with pytest.raises(ValueError, match="two-level categorical factors only"): + generate_design(factors, design_type="dsd") From f66d0e5523f7981483396c3e5b62cd52b82074af Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 08:26:23 +0000 Subject: [PATCH 10/16] Cover the new categorical branches, and drop a dead primality helper Codecov flagged eleven uncovered lines in the categorical work. Covering them turned up two things worth having. The large-c heuristic in the DSD-augment sign search had no test, so nothing checked that a design built without the exhaustive search is still definitive. It is, and now that is asserted. The index-coded branch of the categorical label mapping was reachable after all: dispatch_taguchi codes a categorical factor as level indices 0..n-1, and generate_design(design_type="taguchi") with a categorical factor raised on main because those indices reached the Column constructor unmapped. The mapping added for the DSD fixes it, and a regression test pins it. _is_prime became unreachable when _prime_power_factorization replaced it (the smallest divisor above one is always prime, so no separate test is needed). Ruff does not flag unused module-level functions, so it sat there looking load-bearing. Removed. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 4 + .../experiments/designs_response_surface.py | 16 ---- tests/test_dsd_categorical.py | 80 +++++++++++++++++++ 3 files changed, 84 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fa2f8a15..e88fa01d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,10 @@ those changes. - `matrix_to_columns` now maps a numerically coded categorical column onto its level labels. Design families differ in what they return for a categorical factor: the optimal designs already produce labels, which continue to pass straight through. +- `generate_design(..., design_type="taguchi")` with a categorical factor raised + `ValueError: All values must be present in 'levels'`. `dispatch_taguchi` codes a + categorical factor as level indices `0..n-1`, and those indices were handed to the + `Column` constructor without substituting the labels. Fixed by the same mapping. ## [1.60.1] - 2026-07-25 diff --git a/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index e2e6727a..af344aea 100644 --- a/src/process_improve/experiments/designs_response_surface.py +++ b/src/process_improve/experiments/designs_response_surface.py @@ -576,22 +576,6 @@ def _categorical_sign_candidates(n_categorical: int) -> Iterator[np.ndarray]: yield flipped -def _is_prime(n: int) -> bool: - """Return True iff *n* is a (positive) prime.""" - if n < 2: - return False - if n < 4: - return True - if n % 2 == 0: - return False - i = 3 - while i * i <= n: - if n % i == 0: - return False - i += 2 - return True - - def _prime_power_factorization(q: int) -> tuple[int, int] | None: """Factor *q* as ``p ** n`` for a prime *p*, or return None. diff --git a/tests/test_dsd_categorical.py b/tests/test_dsd_categorical.py index 2e09471d..f33ea10e 100644 --- a/tests/test_dsd_categorical.py +++ b/tests/test_dsd_categorical.py @@ -397,3 +397,83 @@ def test_three_level_categorical_message_reaches_the_caller(self) -> None: ] with pytest.raises(ValueError, match="two-level categorical factors only"): generate_design(factors, design_type="dsd") + + +class TestManyCategoricalFactors: + """Beyond seven categorical factors the exhaustive sign search gives way to a heuristic.""" + + def test_heuristic_search_still_produces_a_valid_design(self) -> None: + """2 ** (2c) becomes too many determinants, so the search is seeded and probed. + + The paper switches to a coordinate-exchange algorithm at large c for the + same reason. The design must still be definitive; only optimality of the + sign choice is given up. + """ + design, meta = dispatch_dsd(_factors(4, 8), categorical_method="dsd") + assert meta["n_categorical"] == 8 + assert design.shape[1] == 12 + assert design.shape[0] == dsd_run_count(12, n_categorical=8) + + first_order, _inter, _labels = _model_matrices(design, 4) + alias = _alias_matrix(first_order, _second_order(design, 4)) + assert np.abs(alias[1:, :]).max() < 1e-12 # still definitive + + for column in design[:, 4:].T: + assert set(np.unique(column)) == {-1.0, 1.0} + assert (column == 1).sum() == (column == -1).sum() + + def test_orth_centre_block_cycles_beyond_four(self) -> None: + """The ORTH centre pattern repeats every four columns for c > 4.""" + design, meta = dispatch_dsd(_factors(4, 6), categorical_method="orth") + assert meta["centre_runs"] == 4 + centre = design[-4:, 4:] + assert set(np.unique(centre)) == {-1.0, 1.0} + # Each of the four centre runs raises exactly one or two categorical factors. + assert all((row == 1).sum() >= 1 for row in centre) + + +class TestCategoricalLabelMapping: + """``matrix_to_columns`` translates coded categorical columns into level labels.""" + + def test_labels_pass_through_untouched(self) -> None: + """Families that already emit labels (the optimal designs) must not be re-mapped.""" + from process_improve.experiments.designs_utils import matrix_to_columns + + factors = [ + Factor(name="X1", low=0, high=10), + Factor(name="Cat", type="categorical", levels=["red", "blue"]), + ] + matrix = np.array([[0.0, "red"], [10.0, "blue"]], dtype=object) + columns = matrix_to_columns(matrix, factors) + assert list(columns[1]) == ["red", "blue"] + + def test_taguchi_with_a_categorical_factor(self) -> None: + """Regression: this raised before index-coded categorical columns were mapped. + + ``dispatch_taguchi`` codes a categorical factor as level indices 0..n-1, + which the ``Column`` constructor rejected because the labels were never + substituted. + """ + factors = [ + Factor(name="A", low=0, high=1), + Factor(name="B", low=0, high=1), + Factor(name="Cat", type="categorical", levels=["x", "y"]), + ] + result = generate_design(factors, design_type="taguchi", center_points=0) + assert set(result.design_actual["Cat"].unique()) == {"x", "y"} + + def test_three_level_categorical_indices_are_mapped(self) -> None: + """The index branch is not limited to two levels.""" + from process_improve.experiments.designs_utils import matrix_to_columns + + factors = [Factor(name="Cat", type="categorical", levels=["a", "b", "c"])] + columns = matrix_to_columns(np.array([[0.0], [2.0], [1.0], [0.0]]), factors) + assert list(columns[0]) == ["a", "c", "b", "a"] + + def test_unrecognised_numeric_coding_is_left_for_the_column_to_reject(self) -> None: + """A coding that matches neither convention must not be silently reinterpreted.""" + from process_improve.experiments.designs_utils import matrix_to_columns + + factors = [Factor(name="Cat", type="categorical", levels=["a", "b"])] + with pytest.raises(ValueError, match="All values must be present in `levels`"): + matrix_to_columns(np.array([[7.0], [9.0]]), factors) From f27190327bdb03a2821de979636ce98584ef60d7 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 21:59:51 +0000 Subject: [PATCH 11/16] Reject categorical factors in OMARS designs rather than returning them unverified design_type="omars" shares dispatch_dsd, so adding categorical support to the DSD silently extended it to OMARS as well. The result was a design carrying omars_verified: False, which is the exact silent-degradation failure the rest of this branch removes. On main the same call raised, further downstream, so this was a regression introduced here. OMARS is defined for quantitative factors: the family's properties are stated in terms of quadratic and interaction columns, and is_omars checks them that way. A two-level categorical factor has no quadratic, so a categorical DSD is a different object and belongs to design_type="dsd". The error message says so. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 6 +++++ .../experiments/designs_omars.py | 16 ++++++++++++- tests/test_dsd_categorical.py | 23 +++++++++++++++++++ 3 files changed, 44 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e88fa01d..eb197d5c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,12 @@ those changes. - `matrix_to_columns` now maps a numerically coded categorical column onto its level labels. Design families differ in what they return for a categorical factor: the optimal designs already produce labels, which continue to pass straight through. +- `generate_design(..., design_type="omars")` now rejects categorical factors with a + message pointing at `design_type="dsd"`. OMARS is defined for quantitative factors: + the family's properties are stated in terms of quadratic and interaction columns, and + `is_omars` checks them that way. Because `omars` shares the DSD constructor, adding + categorical support would otherwise have let it return a design carrying + `omars_verified: False` rather than refusing outright. - `generate_design(..., design_type="taguchi")` with a categorical factor raised `ValueError: All values must be present in 'levels'`. `dispatch_taguchi` codes a categorical factor as level indices `0..n-1`, and those indices were handed to the diff --git a/src/process_improve/experiments/designs_omars.py b/src/process_improve/experiments/designs_omars.py index 1d4e2f57..984fd984 100644 --- a/src/process_improve/experiments/designs_omars.py +++ b/src/process_improve/experiments/designs_omars.py @@ -264,13 +264,27 @@ def dispatch_omars(factors: list[Factor], *, verify: bool = True) -> tuple[np.nd Raises ------ ValueError - If fewer than three factors are supplied. + If fewer than three factors are supplied, or if any factor is + categorical. OMARS is defined for quantitative factors: the family's + properties are stated in terms of quadratic and interaction columns, and + :func:`is_omars` checks them that way. A definitive screening design + with a categorical factor is a different object, so it is reached + through ``design_type="dsd"`` rather than silently returned here with + ``omars_verified: False``. """ from process_improve.experiments.designs_response_surface import dispatch_dsd # noqa: PLC0415 if len(factors) < 3: raise ValueError("OMARS designs require at least 3 factors.") + categorical = [factor.name for factor in factors if factor.type.value == "categorical"] + if categorical: + raise ValueError( + f"OMARS designs require quantitative factors; {', '.join(repr(name) for name in categorical)} " + 'is categorical. Use design_type="dsd", which supports two-level categorical factors through ' + "the column-augmentation procedures of Jones and Nachtsheim (2013)." + ) + coded_matrix, dsd_meta = dispatch_dsd(factors) meta = {**dsd_meta, "family": "conference_foldover"} if verify: diff --git a/tests/test_dsd_categorical.py b/tests/test_dsd_categorical.py index f33ea10e..2088c829 100644 --- a/tests/test_dsd_categorical.py +++ b/tests/test_dsd_categorical.py @@ -477,3 +477,26 @@ def test_unrecognised_numeric_coding_is_left_for_the_column_to_reject(self) -> N factors = [Factor(name="Cat", type="categorical", levels=["a", "b"])] with pytest.raises(ValueError, match="All values must be present in `levels`"): matrix_to_columns(np.array([[7.0], [9.0]]), factors) + + +class TestOmarsRejectsCategoricalFactors: + """OMARS shares the DSD constructor but not its categorical support.""" + + def test_omars_rejects_a_categorical_factor(self) -> None: + """A categorical DSD is not an OMARS design, and must not be returned as one. + + ``is_omars`` states the family's properties in terms of quadratic and + interaction columns, which a two-level categorical factor does not have. + Before this guard, ``design_type="omars"`` accepted a categorical factor + and returned a design with ``omars_verified: False`` attached, which is + exactly the silent-degradation failure this module exists to prevent. + """ + factors = _factors(3, 1) + with pytest.raises(ValueError, match="OMARS designs require quantitative factors"): + generate_design(factors, design_type="omars", center_points=0) + with pytest.raises(ValueError, match='design_type="dsd"'): + generate_design(factors, design_type="omars", center_points=0) + + def test_omars_still_works_for_quantitative_factors(self) -> None: + result = generate_design(_factors(6, 0), design_type="omars", center_points=0) + assert result.metadata["omars_verified"] is True From d0c9fe9d67f24480c2ac785c405a323a2a065c3f Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 22:28:55 +0000 Subject: [PATCH 12/16] Teach the planner about categorical runs, and export the DSD sizing helpers estimate_screening_runs had no way to express categorical factors, so its estimate for a definitive screening design under-counted by one to three runs whenever any were present. It now takes n_categorical and categorical_method and asks the constructor. The recommendation engine passes the number of two-level categorical factors specifically, since those are the only ones a DSD can carry; DOEProblemSpec.n_two_level_categorical exposes that count, distinct from n_categorical. The definitive-screening branch moved into its own helper rather than carrying a complexity suppression. dsd_run_count, dsd_conference_order and dsd_centre_runs are now exported from process_improve.experiments. Sizing a study before running it is a normal thing to want, and it should not require importing from a module named for response surfaces. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 8 +++ src/process_improve/experiments/__init__.py | 8 +++ .../experiments/strategy/budget.py | 52 +++++++++++++++--- .../experiments/strategy/engine.py | 9 ++- .../experiments/strategy/models.py | 10 ++++ tests/test_dsd_categorical.py | 55 +++++++++++++++++++ 6 files changed, 131 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 37070871..d757e94b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,6 +35,14 @@ those changes. be given in any order; the construction needs the categorical columns trailing and permutes them back. A categorical factor with three or more levels raises, with the message pointing at `design_type="d_optimal"`. +- `dsd_run_count`, `dsd_conference_order` and `dsd_centre_runs` are exported from + `process_improve.experiments`, so a study can be sized without importing from + `designs_response_surface` directly. +- `estimate_screening_runs` takes `n_categorical` and `categorical_method`, so the + planner's run estimate for a definitive screening design accounts for the centre runs + that categorical factors add rather than under-budgeting by one to three runs. The + recommendation engine passes the count of *two-level* categorical factors, which are + the only ones a DSD can carry; `DOEProblemSpec.n_two_level_categorical` exposes it. - `dsd_centre_runs(n_categorical, categorical_method)` in `process_improve.experiments.designs_response_surface`, and `dsd_run_count` now takes `n_categorical` and `categorical_method`, so a design's size can still be predicted diff --git a/src/process_improve/experiments/__init__.py b/src/process_improve/experiments/__init__.py index 1a1b0e58..8bd42830 100644 --- a/src/process_improve/experiments/__init__.py +++ b/src/process_improve/experiments/__init__.py @@ -6,6 +6,11 @@ from process_improve.experiments.designs_factorial import full_factorial from process_improve.experiments.designs_omars import is_omars, omars_properties from process_improve.experiments.designs_omars_ilp import generate_omars +from process_improve.experiments.designs_response_surface import ( + dsd_centre_runs, + dsd_conference_order, + dsd_run_count, +) from process_improve.experiments.evaluate import evaluate_all, evaluate_design from process_improve.experiments.factor import Constraint, DesignResult, Factor, Response, ResponseGoal from process_improve.experiments.knowledge import doe_knowledge @@ -38,6 +43,9 @@ "augment_design", "c", "doe_knowledge", + "dsd_centre_runs", + "dsd_conference_order", + "dsd_run_count", "evaluate_all", "evaluate_design", "expand_grid", diff --git a/src/process_improve/experiments/strategy/budget.py b/src/process_improve/experiments/strategy/budget.py index 97151ee0..470634eb 100644 --- a/src/process_improve/experiments/strategy/budget.py +++ b/src/process_improve/experiments/strategy/budget.py @@ -47,21 +47,60 @@ # --------------------------------------------------------------------------- -def estimate_screening_runs(n_factors: int, design_type: str) -> int: # noqa: PLR0911 +def _definitive_screening_runs(n_factors: int, n_categorical: int, categorical_method: str) -> int: + """Run size of a definitive screening design, asked of the constructor. + + ``2k + 1`` runs for an even factor count and ``2k + 3`` for an odd one, + except at the counts where the minimal conference matrix does not exist and + a larger one has to be used, plus the centre runs that two-level categorical + factors require. Guessing any of that from a formula gets it wrong. + """ + from process_improve.experiments.designs_response_surface import dsd_run_count # noqa: PLC0415 + + if n_factors < 3: + # Below the constructor's minimum; keep the nominal formula so budget + # allocation still has a number to work with. + return 2 * n_factors + 1 + return dsd_run_count(n_factors, n_categorical=n_categorical, categorical_method=categorical_method) + + +def estimate_screening_runs( # noqa: PLR0911 + n_factors: int, + design_type: str, + n_categorical: int = 0, + categorical_method: str = "dsd", +) -> int: """Estimate the number of runs for a screening design. Parameters ---------- n_factors : int - Number of factors to screen. + Total number of factors to screen, counting categorical factors. design_type : str One of ``"plackett_burman"``, ``"definitive_screening"``, ``"fractional_factorial"``, ``"full_factorial"``. + n_categorical : int + How many of those factors are two-level categorical. Only the + definitive screening design changes size because of them: each one adds + centre runs. Default 0. + categorical_method : {"dsd", "orth"} + Which column-augmentation procedure a definitive screening design would + use, since ORTH-augment costs two extra runs beyond one categorical + factor. Default ``"dsd"``, matching :func:`~process_improve.experiments.designs.generate_design`. Returns ------- int - Estimated run count including center points. + Estimated run count including centre points. + + Examples + -------- + >>> estimate_screening_runs(7, "definitive_screening") + 17 + >>> estimate_screening_runs(6, "definitive_screening", n_categorical=2) + 14 + >>> estimate_screening_runs(6, "definitive_screening", n_categorical=2, categorical_method="orth") + 16 """ if design_type == "plackett_burman": # Next multiple of 4 >= k + 1 @@ -69,12 +108,7 @@ def estimate_screening_runs(n_factors: int, design_type: str) -> int: # noqa: P return int(math.ceil(n / 4) * 4) if design_type == "definitive_screening": - # 2k + 1 runs for an even factor count and 2k + 3 for an odd one, except - # at the counts where the minimal conference matrix does not exist and a - # larger one has to be used. Ask the constructor rather than guessing. - from process_improve.experiments.designs_response_surface import dsd_run_count # noqa: PLC0415 - - return dsd_run_count(n_factors) if n_factors >= 3 else 2 * n_factors + 1 + return _definitive_screening_runs(n_factors, n_categorical, categorical_method) if design_type == "fractional_factorial": # Smallest 2^(k-p) with resolution >= IV diff --git a/src/process_improve/experiments/strategy/engine.py b/src/process_improve/experiments/strategy/engine.py index b7d818b0..4309ed8f 100644 --- a/src/process_improve/experiments/strategy/engine.py +++ b/src/process_improve/experiments/strategy/engine.py @@ -138,6 +138,7 @@ def _classify_problem(spec: DOEProblemSpec) -> dict[str, Any]: "n_factors": n, "n_continuous": spec.n_continuous, "n_categorical": spec.n_categorical, + "n_two_level_categorical": spec.n_two_level_categorical, "n_mixture": spec.n_mixture, "has_mixture": spec.has_mixture, "has_hard_to_change": spec.has_hard_to_change, @@ -211,7 +212,9 @@ def _large_factor_screening_choice( or domain_pref == "definitive_screening" or classification["prior_confidence"] >= 0.6 ): - return "definitive_screening", estimate_screening_runs(n, "definitive_screening") + return "definitive_screening", estimate_screening_runs( + n, "definitive_screening", n_categorical=classification["n_two_level_categorical"] + ) if domain_pref == "plackett_burman" or (n >= 6 and not classification["is_tight_budget"]): return "plackett_burman", estimate_screening_runs(n, "plackett_burman") if domain_pref == "fractional_factorial": @@ -554,7 +557,9 @@ def _build_alternatives(spec: DOEProblemSpec, classification: dict[str, Any]) -> n = classification["n_factors"] if n >= 6: - runs = estimate_screening_runs(n, "definitive_screening") + runs = estimate_screening_runs( + n, "definitive_screening", n_categorical=classification["n_two_level_categorical"] + ) alternatives.append(f"Definitive Screening Design ({runs} runs) to combine screening and curvature detection.") if n <= 5: alternatives.append(f"Full factorial 2^{n} ({2**n} runs) if budget allows complete information.") diff --git a/src/process_improve/experiments/strategy/models.py b/src/process_improve/experiments/strategy/models.py index ff85b6ec..9d374ca0 100644 --- a/src/process_improve/experiments/strategy/models.py +++ b/src/process_improve/experiments/strategy/models.py @@ -223,6 +223,16 @@ def n_categorical(self) -> int: """Number of categorical factors.""" return sum(1 for f in self.factors if f.type.value == "categorical") + @property + def n_two_level_categorical(self) -> int: + """Number of categorical factors with exactly two levels. + + Only these can enter a definitive screening design (Jones and + Nachtsheim, 2013), so run-size estimates for that family key off this + count rather than :attr:`n_categorical`. + """ + return sum(1 for f in self.factors if f.type.value == "categorical" and len(f.levels or []) == 2) + @property def n_mixture(self) -> int: """Number of mixture factors.""" diff --git a/tests/test_dsd_categorical.py b/tests/test_dsd_categorical.py index 2088c829..8cb1c59b 100644 --- a/tests/test_dsd_categorical.py +++ b/tests/test_dsd_categorical.py @@ -500,3 +500,58 @@ def test_omars_rejects_a_categorical_factor(self) -> None: def test_omars_still_works_for_quantitative_factors(self) -> None: result = generate_design(_factors(6, 0), design_type="omars", center_points=0) assert result.metadata["omars_verified"] is True + + +class TestPlannerAndPublicApi: + """The planning helpers, and their exposure on the package.""" + + def test_planner_accounts_for_categorical_centre_runs(self) -> None: + """estimate_screening_runs must match what generate_design produces. + + Categorical factors add centre runs, so a planner that ignores them + under-budgets the study. + """ + from process_improve.experiments.strategy.budget import estimate_screening_runs + + for m, c in [(4, 0), (4, 1), (4, 2), (5, 3), (8, 4)]: + for method in ("dsd", "orth"): + expected = generate_design( + _factors(m, c), design_type="dsd", center_points=0, categorical_method=method + ).n_runs + assert ( + estimate_screening_runs(m + c, "definitive_screening", n_categorical=c, categorical_method=method) + == expected + ), f"m={m} c={c} {method}" + + def test_planner_defaults_are_unchanged_for_continuous_designs(self) -> None: + from process_improve.experiments.strategy.budget import estimate_screening_runs + + assert estimate_screening_runs(7, "definitive_screening") == 17 + assert estimate_screening_runs(8, "definitive_screening") == 17 + assert estimate_screening_runs(21, "definitive_screening") == 49 + + def test_spec_counts_only_two_level_categoricals(self) -> None: + """Only two-level categorical factors can enter a DSD, so only they are counted.""" + from process_improve.experiments.strategy.models import DOEProblemSpec + + spec = DOEProblemSpec( + objective="screen", + factors=[ + Factor(name="X1", low=0, high=1), + Factor(name="C2", type="categorical", levels=["a", "b"]), + Factor(name="C3", type="categorical", levels=["a", "b", "c"]), + ], + responses=[], + ) + assert spec.n_categorical == 2 + assert spec.n_two_level_categorical == 1 + + def test_sizing_helpers_are_importable_from_the_package(self) -> None: + """A user planning a study should not have to reach into a private module.""" + from process_improve import experiments + + assert experiments.dsd_run_count(7) == 17 + assert experiments.dsd_conference_order(7) == 8 + assert experiments.dsd_centre_runs(2, "orth") == 4 + for name in ("dsd_run_count", "dsd_conference_order", "dsd_centre_runs"): + assert name in experiments.__all__ From 46a24f59521c37a9b139668a39af24df267c0b89 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 22:38:22 +0000 Subject: [PATCH 13/16] Test the planner's sub-minimum guard A definitive screening design needs three factors and dsd_run_count raises below that, so the planner keeps the nominal 2k + 1 rather than letting the exception reach budget allocation for a design it would never recommend. That branch had no test; the accompanying assertion that the constructor itself still refuses is the half that matters. budget.py is now fully covered, statements and branches. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- tests/test_dsd_categorical.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/tests/test_dsd_categorical.py b/tests/test_dsd_categorical.py index 8cb1c59b..864a34ba 100644 --- a/tests/test_dsd_categorical.py +++ b/tests/test_dsd_categorical.py @@ -530,6 +530,22 @@ def test_planner_defaults_are_unchanged_for_continuous_designs(self) -> None: assert estimate_screening_runs(8, "definitive_screening") == 17 assert estimate_screening_runs(21, "definitive_screening") == 49 + def test_below_the_constructors_minimum_the_planner_still_returns_a_number(self) -> None: + """A DSD needs three factors, but budget allocation still needs an estimate. + + ``dsd_run_count`` raises below three factors, so the planner keeps the + nominal ``2k + 1`` there rather than propagating the exception into + budget allocation for a design it would never actually recommend. + """ + from process_improve.experiments.strategy.budget import estimate_screening_runs + + assert estimate_screening_runs(1, "definitive_screening") == 3 + assert estimate_screening_runs(2, "definitive_screening") == 5 + + # The constructor itself still refuses, which is what matters. + with pytest.raises(ValueError, match="at least 3 factors"): + dsd_run_count(2) + def test_spec_counts_only_two_level_categoricals(self) -> None: """Only two-level categorical factors can enter a DSD, so only they are counted.""" from process_improve.experiments.strategy.models import DOEProblemSpec From f89b9f35ce159753e90d86dd37b067fa7c0f0e5d Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 22:48:00 +0000 Subject: [PATCH 14/16] Make the dataset tests actually skip when openmv.net is unreachable A 502 from openmv.net failed the 3.11 job on this branch. The test is meant to skip when the host cannot be reached, and its docstring says so, but the guard had never worked: _read_remote_csv deliberately converts OSError and ValueError into a RuntimeError carrying an actionable message, while _load_or_skip caught only URLError, HTTPError and OSError. RuntimeError is not an OSError, so nothing the loader raises was ever caught, and any remote hiccup failed the build. Catching RuntimeError restores the documented behaviour. Nothing is weakened: when the host is up the test fetches and asserts in full, as it does locally. Two tests pin the helper itself, since the mismatch is invisible until an outage happens, which is the worst time to discover it. Neither file is touched by the rest of this branch; the failure was pre-existing and would have hit any PR unlucky enough to build during an outage. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 6 +++++ tests/test_experiments_datasets.py | 40 ++++++++++++++++++++++++++++-- 2 files changed, 44 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d757e94b..2a6dc34f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -53,6 +53,12 @@ those changes. ### Fixed +- The dataset tests' "skip if offline" guard never fired. `_read_remote_csv` converts + network failures into a `RuntimeError` for a readable message, but the test helper + caught only `urllib.error.URLError`, `HTTPError` and `OSError`, none of which can + escape the loader (`RuntimeError` is not an `OSError`). A remote outage at openmv.net + therefore failed the build instead of skipping, which is what the helper and the tests + both claimed would happen. - `matrix_to_columns` now maps a numerically coded categorical column onto its level labels. Design families differ in what they return for a categorical factor: the optimal designs already produce labels, which continue to pass straight through. diff --git a/tests/test_experiments_datasets.py b/tests/test_experiments_datasets.py index eec82f99..dd5c7399 100644 --- a/tests/test_experiments_datasets.py +++ b/tests/test_experiments_datasets.py @@ -9,23 +9,59 @@ from __future__ import annotations +import re import urllib.error from collections.abc import Callable import pandas as pd import pytest +from _pytest.outcomes import Skipped from process_improve.experiments import datasets def _load_or_skip(loader: Callable[[], pd.DataFrame]) -> pd.DataFrame: - """Call the network-backed loader, ``pytest.skip`` on any network error.""" + """Call the network-backed loader, ``pytest.skip`` on any network error. + + ``RuntimeError`` is the one that matters in practice, and it has to be + listed explicitly. ``_read_remote_csv`` deliberately converts the + lower-level ``OSError`` / ``ValueError`` into a ``RuntimeError`` carrying an + actionable message, and ``RuntimeError`` is not an ``OSError``. Catching + only the urllib types therefore skipped nothing: a remote outage (openmv.net + returning 502, say) failed the build rather than skipping, contrary to what + this helper and the tests below claim. The other types are kept in case a + loader ever raises one directly. + """ try: return loader() - except (urllib.error.URLError, urllib.error.HTTPError, OSError) as exc: + except (RuntimeError, urllib.error.URLError, urllib.error.HTTPError, OSError) as exc: pytest.skip(f"could not fetch from openmv.net: {exc}") +def test_load_or_skip_skips_when_the_remote_host_is_unavailable() -> None: + """A remote outage must skip, not fail the build. + + ``_read_remote_csv`` wraps network failures in a ``RuntimeError``, so that + is what the helper actually sees. Pinning it here because the mismatch is + invisible until openmv.net has an outage, at which point every build fails. + """ + + def unavailable() -> pd.DataFrame: + raise RuntimeError( + "Could not download the sample dataset from 'https://openmv.net/file/oil-company-doe.csv': " + "HTTP Error 502: Bad Gateway." + ) + + with pytest.raises(Skipped, match=re.escape("could not fetch from openmv.net")): + _load_or_skip(unavailable) + + +def test_load_or_skip_returns_the_frame_when_the_loader_succeeds() -> None: + """The guard must not swallow a working fetch.""" + frame = pd.DataFrame({"A": [1], "y": [2]}) + assert _load_or_skip(lambda: frame) is frame + + def test_boilingpot_loads() -> None: """``boilingpot()`` returns the documented 11x4 factorial frame.""" df = datasets.boilingpot() From 5b6265f6eb72bca447aa5ad7600724e32105f716 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 26 Jul 2026 05:15:19 +0000 Subject: [PATCH 15/16] Retry a transient failure when fetching the remote sample datasets The sample datasets are fetched live from openmv.net, so a momentary 502 from that host failed an entire CI job on this branch. _read_remote_csv now makes three attempts with 1s and 2s of backoff, which rides out a blip at a cost of at most three seconds on a genuine outage. Only failures that could change on a second attempt are retried: connection errors, and 408, 425, 429 and the 5xx statuses. A 404 or a 403 will say the same thing next time, and a ValueError means the fetch succeeded and the content was wrong, so both are raised immediately. The error raised after the budget is exhausted is unchanged, so the callers that skip on it still skip. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- CHANGELOG.md | 6 ++ src/process_improve/experiments/datasets.py | 80 ++++++++++++++++-- tests/test_experiments_datasets.py | 92 +++++++++++++++++++++ 3 files changed, 170 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a6dc34f..91fd3640 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -53,6 +53,12 @@ those changes. ### Fixed +- `_read_remote_csv` retries a transient failure before giving up: three attempts with + 1s and 2s of backoff. The sample datasets are fetched live from openmv.net, so a + momentary 502 from that host failed a whole CI job. Only failures that can change on a + second attempt are retried (connection errors and 408/425/429/5xx); a 404, a 403 or a + parse error is raised immediately rather than paying the backoff for an answer that + will not change. - The dataset tests' "skip if offline" guard never fired. `_read_remote_csv` converts network failures into a `RuntimeError` for a readable message, but the test helper caught only `urllib.error.URLError`, `HTTPError` and `OSError`, none of which can diff --git a/src/process_improve/experiments/datasets.py b/src/process_improve/experiments/datasets.py index c0a2f8a0..c7ed5116 100644 --- a/src/process_improve/experiments/datasets.py +++ b/src/process_improve/experiments/datasets.py @@ -2,28 +2,92 @@ from __future__ import annotations +import logging +import time +import urllib.error from pathlib import Path import pandas as pd +logger = logging.getLogger(__name__) + _DATASETS_DIR = Path(__file__).resolve().parents[1] / "datasets" / "experiments" +# Retry policy for the live sample-data fetches. Three attempts at 1s and 2s of +# backoff costs at most three seconds on a genuine outage, and rides out the +# momentary 5xx blips that would otherwise fail a whole run. +_FETCH_ATTEMPTS = 3 +_FETCH_BACKOFF_SECONDS = 1.0 + +# Server-side and connection failures worth a second look. A 404 or a 403 will +# say the same thing next time, so those are not retried. +_TRANSIENT_HTTP_STATUS = frozenset({408, 425, 429, 500, 502, 503, 504}) + + +def _is_transient(exc: Exception) -> bool: + """Return True when *exc* looks like a blip rather than a settled answer. + + A parse failure (``ValueError``) means the fetch succeeded and the content + was wrong, so retrying cannot help. + """ + if isinstance(exc, urllib.error.HTTPError): + return exc.code in _TRANSIENT_HTTP_STATUS + return isinstance(exc, OSError) -def _read_remote_csv(url: str) -> pd.DataFrame: + +def _read_remote_csv(url: str, *, attempts: int = _FETCH_ATTEMPTS) -> pd.DataFrame: """Fetch a sample-dataset CSV from a (hard-coded, trusted) remote host. The URL is not user-supplied. Network or parse failures are surfaced as a clear ``RuntimeError`` rather than a lower-level ``URLError`` / parser error, so callers get an actionable message. Note that the content is fetched over the network and is therefore trusted only as far as the remote host is. + + A transient failure is retried with exponential backoff before giving up. + Sample data is fetched live, so a momentary blip at the remote host would + otherwise fail a whole test run or documentation build; a 502 from the host + has done exactly that. Only failures that look transient are retried: a + parse error, or an HTTP status that will not change on a second attempt, is + raised immediately rather than paying the backoff for nothing. + + Parameters + ---------- + url : str + The dataset URL. Hard-coded by the calling loader, never user-supplied. + attempts : int + Maximum number of tries, including the first. Defaults to + :data:`_FETCH_ATTEMPTS`; pass ``1`` to disable retrying. + + Returns + ------- + pandas.DataFrame + The parsed CSV. + + Raises + ------ + RuntimeError + If the dataset could not be fetched or parsed. The message names the URL + and says the data comes from a remote host, so a caller can tell an + outage apart from a genuine bug. """ - try: - return pd.read_csv(url) - except (OSError, ValueError) as exc: - raise RuntimeError( - f"Could not download the sample dataset from {url!r}: {exc}. " - "Check your network connection; this dataset is fetched from a remote host." - ) from exc + last_exc: Exception | None = None + for attempt in range(1, attempts + 1): + try: + return pd.read_csv(url) + except (OSError, ValueError) as exc: # noqa: PERF203 - retry loop; the cost here is the network, not the try + last_exc = exc + if attempt == attempts or not _is_transient(exc): + break + delay = _FETCH_BACKOFF_SECONDS * 2 ** (attempt - 1) + logger.info( + "Fetching %s failed (%s); retrying in %.1fs (attempt %d of %d).", url, exc, delay, attempt + 1, attempts + ) + time.sleep(delay) + + raise RuntimeError( + f"Could not download the sample dataset from {url!r}: {last_exc}. " + "Check your network connection; this dataset is fetched from a remote host." + ) from last_exc def distillateflow() -> pd.DataFrame: diff --git a/tests/test_experiments_datasets.py b/tests/test_experiments_datasets.py index dd5c7399..2a7582db 100644 --- a/tests/test_experiments_datasets.py +++ b/tests/test_experiments_datasets.py @@ -108,3 +108,95 @@ def test_data_dispatch_signature_typed() -> None: # (``from __future__ import annotations``). assert datasets.data.__annotations__["return"] == "pd.DataFrame" assert datasets.data.__annotations__["dataset"] == "str" + + +class TestRemoteFetchRetries: + """A transient blip at the remote host must not fail a whole run. + + Sample data is fetched live, so a momentary 502 from openmv.net has failed + an entire CI job. Retrying rides that out; the tests below pin both that it + retries when it should and that it does not when retrying cannot help. + """ + + def test_transient_error_is_retried_and_succeeds(self, monkeypatch) -> None: + frame = pd.DataFrame({"A": [1]}) + calls = {"n": 0} + + def flaky(_url): + calls["n"] += 1 + if calls["n"] < 3: + raise urllib.error.HTTPError(_url, 502, "Bad Gateway", {}, None) + return frame + + monkeypatch.setattr(datasets.pd, "read_csv", flaky) + monkeypatch.setattr(datasets.time, "sleep", lambda _s: None) + assert datasets._read_remote_csv("https://openmv.net/file/x.csv") is frame + assert calls["n"] == 3 + + def test_gives_up_after_the_attempt_budget(self, monkeypatch) -> None: + calls = {"n": 0} + + def always_down(_url): + calls["n"] += 1 + raise urllib.error.HTTPError(_url, 503, "Service Unavailable", {}, None) + + monkeypatch.setattr(datasets.pd, "read_csv", always_down) + monkeypatch.setattr(datasets.time, "sleep", lambda _s: None) + with pytest.raises(RuntimeError, match="Could not download the sample dataset"): + datasets._read_remote_csv("https://openmv.net/file/x.csv") + assert calls["n"] == datasets._FETCH_ATTEMPTS + + @pytest.mark.parametrize("status", [404, 403, 410]) + def test_settled_http_status_is_not_retried(self, monkeypatch, status: int) -> None: + """A 404 says the same thing next time; do not pay the backoff for it.""" + calls = {"n": 0} + + def missing(_url): + calls["n"] += 1 + raise urllib.error.HTTPError(_url, status, "nope", {}, None) + + monkeypatch.setattr(datasets.pd, "read_csv", missing) + monkeypatch.setattr(datasets.time, "sleep", lambda _s: pytest.fail("should not back off")) + with pytest.raises(RuntimeError, match="Could not download the sample dataset"): + datasets._read_remote_csv("https://openmv.net/file/x.csv") + assert calls["n"] == 1 + + def test_parse_error_is_not_retried(self, monkeypatch) -> None: + """A ValueError means the fetch worked and the content was wrong.""" + calls = {"n": 0} + + def bad_content(_url): + calls["n"] += 1 + raise ValueError("could not parse") + + monkeypatch.setattr(datasets.pd, "read_csv", bad_content) + monkeypatch.setattr(datasets.time, "sleep", lambda _s: pytest.fail("should not back off")) + with pytest.raises(RuntimeError, match="Could not download the sample dataset"): + datasets._read_remote_csv("https://openmv.net/file/x.csv") + assert calls["n"] == 1 + + def test_connection_error_is_retried(self, monkeypatch) -> None: + """A bare OSError (DNS, reset connection) is worth another try.""" + calls = {"n": 0} + + def dns_fail(_url): + calls["n"] += 1 + raise OSError("name resolution failed") + + monkeypatch.setattr(datasets.pd, "read_csv", dns_fail) + monkeypatch.setattr(datasets.time, "sleep", lambda _s: None) + with pytest.raises(RuntimeError): + datasets._read_remote_csv("https://openmv.net/file/x.csv") + assert calls["n"] == datasets._FETCH_ATTEMPTS + + def test_attempts_can_be_disabled(self, monkeypatch) -> None: + calls = {"n": 0} + + def down(_url): + calls["n"] += 1 + raise urllib.error.HTTPError(_url, 502, "Bad Gateway", {}, None) + + monkeypatch.setattr(datasets.pd, "read_csv", down) + with pytest.raises(RuntimeError): + datasets._read_remote_csv("https://openmv.net/file/x.csv", attempts=1) + assert calls["n"] == 1 From 6f17888ee92d594070d5de4e18d3e84a25784cfe Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 26 Jul 2026 05:26:53 +0000 Subject: [PATCH 16/16] Remove the unreachable tail from the fetch retry loop Codecov flagged a partial branch: the for loop's normal exit, which cannot happen because the body either returns or breaks on the last attempt. Rather than leave an uncoverable branch, the loop is now a while-True whose only exits are the return and the raise, so there is no fallthrough and no dead tail to carry. The one case the old shape did reach on that path was attempts < 1, which fell through with last_exc still None and reported "Could not download ...: None" - a failure that never happened. That is now rejected up front as the caller error it is. datasets.py is at 100%, statements and branches. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TeEyMBnkeSKZ8QbSyejvNu --- src/process_improve/experiments/datasets.py | 29 ++++++++++++--------- tests/test_experiments_datasets.py | 5 ++++ 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/process_improve/experiments/datasets.py b/src/process_improve/experiments/datasets.py index c7ed5116..c49d5d84 100644 --- a/src/process_improve/experiments/datasets.py +++ b/src/process_improve/experiments/datasets.py @@ -56,7 +56,8 @@ def _read_remote_csv(url: str, *, attempts: int = _FETCH_ATTEMPTS) -> pd.DataFra The dataset URL. Hard-coded by the calling loader, never user-supplied. attempts : int Maximum number of tries, including the first. Defaults to - :data:`_FETCH_ATTEMPTS`; pass ``1`` to disable retrying. + :data:`_FETCH_ATTEMPTS`; pass ``1`` to disable retrying. Must be at + least 1. Returns ------- @@ -69,26 +70,30 @@ def _read_remote_csv(url: str, *, attempts: int = _FETCH_ATTEMPTS) -> pd.DataFra If the dataset could not be fetched or parsed. The message names the URL and says the data comes from a remote host, so a caller can tell an outage apart from a genuine bug. + ValueError + If *attempts* is less than 1, which would otherwise fetch nothing and + report a failure that never happened. """ - last_exc: Exception | None = None - for attempt in range(1, attempts + 1): + if attempts < 1: + raise ValueError(f"attempts must be at least 1; got {attempts}.") + + attempt = 0 + while True: + attempt += 1 try: return pd.read_csv(url) - except (OSError, ValueError) as exc: # noqa: PERF203 - retry loop; the cost here is the network, not the try - last_exc = exc - if attempt == attempts or not _is_transient(exc): - break + except (OSError, ValueError) as exc: + if attempt >= attempts or not _is_transient(exc): + raise RuntimeError( + f"Could not download the sample dataset from {url!r}: {exc}. " + "Check your network connection; this dataset is fetched from a remote host." + ) from exc delay = _FETCH_BACKOFF_SECONDS * 2 ** (attempt - 1) logger.info( "Fetching %s failed (%s); retrying in %.1fs (attempt %d of %d).", url, exc, delay, attempt + 1, attempts ) time.sleep(delay) - raise RuntimeError( - f"Could not download the sample dataset from {url!r}: {last_exc}. " - "Check your network connection; this dataset is fetched from a remote host." - ) from last_exc - def distillateflow() -> pd.DataFrame: """Return the flow rate of distillate from the top of a distillation column. diff --git a/tests/test_experiments_datasets.py b/tests/test_experiments_datasets.py index 2a7582db..939f5916 100644 --- a/tests/test_experiments_datasets.py +++ b/tests/test_experiments_datasets.py @@ -200,3 +200,8 @@ def down(_url): with pytest.raises(RuntimeError): datasets._read_remote_csv("https://openmv.net/file/x.csv", attempts=1) assert calls["n"] == 1 + + def test_attempts_must_be_at_least_one(self) -> None: + """Zero attempts would report a failure that never happened.""" + with pytest.raises(ValueError, match="attempts must be at least 1"): + datasets._read_remote_csv("https://openmv.net/file/x.csv", attempts=0)