DSD: fix the conference-matrix construction, and add two-level categorical factors

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

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

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, …).
 

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"

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",

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:

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,
     }