diff --git a/CHANGELOG.md b/CHANGELOG.md index e78e309a..91fd3640 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,101 @@ those changes. ## [Unreleased] +## [1.62.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_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 + without building it. +- `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. + +### 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 + 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. +- `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 + `Column` constructor without substituting the labels. Fixed by the same mapping. +- 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. + + ## [1.61.0] - 2026-07-25 ### Fixed @@ -2707,7 +2802,8 @@ this entry records them together. - Reworked the README with a sharper value proposition and a "Why not scikit-learn?" comparison table. -[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.61.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.62.0...HEAD +[1.62.0]: https://github.com/kgdunn/process-improve/compare/v1.61.0...v1.62.0 [1.61.0]: https://github.com/kgdunn/process-improve/compare/v1.60.0...v1.61.0 [1.60.0]: https://github.com/kgdunn/process-improve/compare/v1.59.0...v1.60.0 [1.59.0]: https://github.com/kgdunn/process-improve/compare/v1.58.0...v1.59.0 diff --git a/CITATION.cff b/CITATION.cff index 8990021e..ca5b16c7 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.61.0 +version: 1.62.0 date-released: "2026-07-25" keywords: - chemometrics diff --git a/docs/doe/coverage.md b/docs/doe/coverage.md index fec9d443..a0048831 100644 --- a/docs/doe/coverage.md +++ b/docs/doe/coverage.md @@ -54,7 +54,8 @@ 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. +- **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 cdc0efac..b8b154c9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.61.0" +version = "1.62.0" description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.' readme = "README.md" license = "MIT" diff --git a/src/process_improve/experiments/__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/datasets.py b/src/process_improve/experiments/datasets.py index c0a2f8a0..c49d5d84 100644 --- a/src/process_improve/experiments/datasets.py +++ b/src/process_improve/experiments/datasets.py @@ -2,28 +2,97 @@ 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. Must be at + least 1. + + 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. + ValueError + If *attempts* is less than 1, which would otherwise fetch nothing and + report a failure that never happened. """ - 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 + 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: + 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) def distillateflow() -> pd.DataFrame: 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_omars.py b/src/process_improve/experiments/designs_omars.py index e6d06ebb..984fd984 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 ------- @@ -260,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/src/process_improve/experiments/designs_response_surface.py b/src/process_improve/experiments/designs_response_surface.py index cc0e58bc..af344aea 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,34 +301,81 @@ 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-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. + + 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. + 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 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), ``"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, 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 ---------- @@ -327,54 +385,367 @@ 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.") - # 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) + 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) - 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: - coded_matrix = coded_matrix[:, :k] + body = np.vstack([conference, -conference]) + if order > k: + body = body[:, :k] - return coded_matrix, {"construction": construction} + 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, + "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( + "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, + coded_matrix.shape[0], + 2 * minimal_order + centre_runs, + ) -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 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 _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 + # 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]: + """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 +753,290 @@ 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) - n = q + 1 - c_matrix = np.zeros((n, n), dtype=int) - c_matrix[0, 1:] = 1 + # Jacobsthal matrix Q[a, b] = chi(b - a), where the subtraction is in the field. + jacobsthal = character[subtraction.T] + + 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 + matrix[1:, 0] = -1 + matrix[1:, 1:] = jacobsthal + return matrix + + +_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], + ] +) +"""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} + +# 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 _cyclic_conference_matrix(k: int) -> np.ndarray: - """Legacy cyclic approximation of a conference matrix. +def _validate_conference_matrix(matrix: np.ndarray, order: int) -> None: + """Check the defining property ``C.T @ C == (m - 1) * I`` 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. + 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; see + :data:`_CONFERENCE_MATRIX_16` for where that table comes from and who holds + the copyright on the ports it travelled through. 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, 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 + :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 + 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 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"``. + """ + 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 b7de854e..f241dd89 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, 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 - "quadratic (limited active factors)" diff --git a/src/process_improve/experiments/strategy/budget.py b/src/process_improve/experiments/strategy/budget.py index 80ffa333..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,7 +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": - return 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 c3cba685..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": @@ -226,7 +229,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 +557,10 @@ 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", 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.") if n >= 4: 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_conference_matrices.py b/tests/test_conference_matrices.py new file mode 100644 index 00000000..16ca4ba0 --- /dev/null +++ b/tests/test_conference_matrices.py @@ -0,0 +1,349 @@ +"""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 ( + _MAX_CONFERENCE_ORDER, + _conference_matrix, + _gf_multiplication_table, + _irreducible_polynomial, + _is_constructible_conference_order, + _paley_conference_matrix, + _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)) + + 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.""" + 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 + + 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.""" + + 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_dsd_categorical.py b/tests/test_dsd_categorical.py new file mode 100644 index 00000000..864a34ba --- /dev/null +++ b/tests/test_dsd_categorical.py @@ -0,0 +1,573 @@ +"""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") + + +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) + + +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 + + +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_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 + + 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__ diff --git a/tests/test_experiments_datasets.py b/tests/test_experiments_datasets.py index eec82f99..939f5916 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() @@ -72,3 +108,100 @@ 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 + + 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) 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") 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."""