From 9502124c718935ca69437222e64e3e291ecca2be Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 3 Aug 2026 12:19:56 +0000 Subject: [PATCH 1/5] Add sensory screening design: blocked, carry-over balanced serving plans Everything in `process_improve.sensory` so far analyses panel data that already exists. This adds the step before it: turning a long candidate list into a serving sheet a sensory lab can run. A panel screen is not an ordinary designed experiment, and the new module handles the three reasons why: - Session capacity. An assessor can only judge a handful of samples in one sitting, so once the candidate list is longer than a session the design is an incomplete block design with the assessor-session as the block. `cyclic_block_design` builds the blocks by a greedy fill followed by a hill-climb over two move types (substitution to level replication, exchange to level concurrence once replication has settled). Both moves are needed: with replication already level, every substitution makes it worse, so without the exchange the search stalls one step short of an exact balanced incomplete block design. - Carry-over and position. `williams_design` returns serving orders in which every treatment appears once per position and every ordered pair occurs equally often; an odd treatment count needs the mirrored second square. - Sensitivity. `detectable_difference` / `required_panelists` convert a residual SD and a panel size into the smallest resolvable difference and back, with a Bonferroni correction for the comparison family. `sensory_screening_plan` assembles the three into one sheet, optionally anchoring a reference first in every block. The sheet uses the `descriptive_long` column names, so adding scores feeds it straight into `compare_products` with `panelist_id` as the block. A panel too small to cover the candidate list is reported in `warnings` rather than silently dropping candidates. 36 tests, including the Fano plane as a known balanced incomplete block design and the ordered-pair balance of the Williams squares. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TqdSBoXb1DpSobB76pMzFi --- CHANGELOG.md | 25 +- CITATION.cff | 4 +- docs/user_guide/sensory_panel.rst | 71 ++ pyproject.toml | 2 +- src/process_improve/sensory/__init__.py | 18 + src/process_improve/sensory/screening.py | 874 +++++++++++++++++++++++ src/process_improve/sensory/tools.py | 173 ++++- tests/test_sensory_screening.py | 359 ++++++++++ 8 files changed, 1519 insertions(+), 7 deletions(-) create mode 100644 src/process_improve/sensory/screening.py create mode 100644 tests/test_sensory_screening.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 85b64d49..74a8c188 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,28 @@ those changes. ## [Unreleased] +## [1.63.0] - 2026-08-03 + +### Added + +- `process_improve.sensory.screening`: design a descriptive-panel screen before + any tasting happens. `sensory_screening_plan` turns a candidate list plus the + panel's session capacity into a serving sheet: incomplete blocks (one block = + one assessor-session) built by `cyclic_block_design` for near-equal + replication and pairwise concurrence, ordered inside each block by + `williams_design` so first-order carry-over cannot favour any candidate, and + with an optional reference anchored first in every block. The sheet comes back + in the `descriptive_long` shape, so adding scores feeds it straight into + `compare_products`. +- `detectable_difference` / `required_panelists`: the smallest difference a + panel size can resolve, and the panel size a target difference needs, with a + Bonferroni correction for the size of the comparison family. +- `plan_diagnostics` reports replication, pairwise concurrence and whether the + blocks form an exact balanced incomplete block design; a panel too small to + cover the candidate list is reported as a warning rather than silently + dropping candidates. +- Agent tools `sensory_screening_plan` and `sensory_detectable_difference`. + ## [1.62.2] - 2026-07-29 ### Changed @@ -2787,7 +2809,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.62.2...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.63.0...HEAD +[1.63.0]: https://github.com/kgdunn/process-improve/compare/v1.62.2...v1.63.0 [1.62.2]: https://github.com/kgdunn/process-improve/compare/v1.62.1...v1.62.2 [1.62.1]: https://github.com/kgdunn/process-improve/compare/v1.62.0...v1.62.1 [1.62.0]: https://github.com/kgdunn/process-improve/compare/v1.61.0...v1.62.0 diff --git a/CITATION.cff b/CITATION.cff index 17f25a94..ec3889ae 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,8 +12,8 @@ authors: repository-code: "https://github.com/kgdunn/process-improve" url: "https://kgdunn.github.io/process-improve/" license: MIT -version: 1.62.2 -date-released: "2026-07-29" +version: 1.63.0 +date-released: '2026-08-03' keywords: - chemometrics - multivariate analysis diff --git a/docs/user_guide/sensory_panel.rst b/docs/user_guide/sensory_panel.rst index abe59203..63cc1bca 100644 --- a/docs/user_guide/sensory_panel.rst +++ b/docs/user_guide/sensory_panel.rst @@ -406,6 +406,73 @@ with correctly calibrated thresholds. all-pairwise comparisons. *Journal of Computational and Graphical Statistics*, 13(2), 456-466. +Before the panel: designing the screen +------------------------------------- + +Everything above analyses panel data that already exists. When you are at the +other end - a long list of candidate samples and no scores yet - the question is +which samples each assessor should taste, in which order, and whether the panel +is big enough to see the effect at all. :mod:`process_improve.sensory.screening` +covers that step. + +A screen is not an ordinary designed experiment, for three reasons: + +**Session capacity.** An assessor can only judge a handful of samples in one +sitting before fatigue and carry-over dominate. Once the candidate list is +longer than a session, every assessor sees only a *subset*, so the design is an +**incomplete block design** with the assessor-session as the block. +:func:`~process_improve.sensory.screening.cyclic_block_design` builds the blocks +so replication is equal, or as near-equal as the arithmetic allows, and so every +pair of candidates meets inside a block about equally often. That pairwise +**concurrence** is what decides how precisely two candidates can be compared with +each other; when an exact balanced incomplete block design exists for the given +numbers, the construction finds one, and ``diagnostics["balanced"]`` says so. + +**Carry-over and position.** What was tasted just before changes this score, and +the first sample of a session is not scored like the last. +:func:`~process_improve.sensory.screening.williams_design` returns serving orders +in which every treatment appears once in each position and every *ordered* pair +occurs equally often, so no candidate is systematically favoured by its +neighbours. An even number of treatments needs one square; an odd number needs a +mirrored second one. + +**Sensitivity.** Panel scores are noisy, so a screen that is too small returns +"no significant effect" whatever the truth is. +:func:`~process_improve.sensory.screening.detectable_difference` turns a residual +standard deviation and a panel size into the smallest difference the screen can +actually resolve, and +:func:`~process_improve.sensory.screening.required_panelists` inverts it. Correct +for the comparison family with ``n_comparisons``: testing twenty candidates +against one control is twenty comparisons, and pretending otherwise overstates +what the screen can see. + +:func:`~process_improve.sensory.screening.sensory_screening_plan` assembles all +three into one serving sheet, optionally anchoring a reference sample first in +every block so session and assessor drift can be removed later: + +.. code-block:: python + + from process_improve.sensory import sensory_screening_plan + + result = sensory_screening_plan( + [f"Candidate {i:02d}" for i in range(1, 22)], + n_panelists=12, + samples_per_session=6, # 5 candidates + the anchored reference + control="Reference", + replicates=2, + seed=0, + ) + result.plan.head() # panelist_id, session, position, product, role, block + result.diagnostics # replication, concurrence, balanced, control_coverage + result.warnings # e.g. a panel too small to cover the list + +The plan comes back in the same long shape the rest of this page consumes: fill +in the scores next to ``product`` and it feeds straight into +:func:`~process_improve.sensory.compare_products`, with ``panelist_id`` as the +block. A capacity shortfall is reported in ``warnings`` rather than silently +dropping candidates, because quietly screening eighteen of your twenty-one +candidates is the failure mode worth being loud about. + Worked example -------------- @@ -541,6 +608,10 @@ and covariate tables as lists of row-records and returning JSON: option (``"none"`` / ``"align"`` / ``"drop"``), the MAM results, and (unless ``discriminator`` is set false) the cross-validated discriminator in its output. +- ``sensory_screening_plan`` - the blocked, carry-over balanced serving sheet + for a screen that has not been run yet. +- ``sensory_detectable_difference`` - what a given panel size can resolve, or + the panel size a given difference needs. The analyze tool validates first and refuses to run if validation fails, so an agent cannot skip the gate. diff --git a/pyproject.toml b/pyproject.toml index 4aba40f6..7ea2d53b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.62.2" +version = "1.63.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/sensory/__init__.py b/src/process_improve/sensory/__init__.py index f2df1413..671a0378 100644 --- a/src/process_improve/sensory/__init__.py +++ b/src/process_improve/sensory/__init__.py @@ -36,6 +36,16 @@ from process_improve.sensory.mam import MAMResult, align_scores, mixed_assessor_model from process_improve.sensory.panel import PanelScorecard, apply_correction, panel_scorecard from process_improve.sensory.recipes import SENSORY_RECIPES +from process_improve.sensory.screening import ( + SERVING_PLAN_COLUMNS, + ScreeningPlan, + cyclic_block_design, + detectable_difference, + plan_diagnostics, + required_panelists, + sensory_screening_plan, + williams_design, +) from process_improve.sensory.validation import ( DESCRIPTIVE_LONG_COLUMNS, ValidationResult, @@ -45,26 +55,34 @@ __all__ = [ "DESCRIPTIVE_LONG_COLUMNS", "SENSORY_RECIPES", + "SERVING_PLAN_COLUMNS", "AnalysisResult", "ComparisonResult", "MAMResult", "PanelScorecard", + "ScreeningPlan", "ValidationResult", "aggregate_to_product", "align_scores", "analyze_descriptive", "apply_correction", "compare_products", + "cyclic_block_design", + "detectable_difference", "discriminate_observational", "dunnett_vs_control", "factorial_anova", "mixed_assessor_model", "panel_scorecard", "permutation_column_null", + "plan_diagnostics", "product_means", "relate_designed", "relate_observational", + "required_panelists", "reshape_to_long", + "sensory_screening_plan", "tukey_hsd", "validate_descriptive", + "williams_design", ] diff --git a/src/process_improve/sensory/screening.py b/src/process_improve/sensory/screening.py new file mode 100644 index 00000000..9b3e3e45 --- /dev/null +++ b/src/process_improve/sensory/screening.py @@ -0,0 +1,874 @@ +"""(c) Kevin Dunn, 2010-2026. MIT License. + +Designing a descriptive-panel **screening** study: which samples does each +assessor taste, in which order, and is the panel big enough to see the effect? + +The rest of this subpackage analyses panel data that already exists. This module +covers the step before that: turning a long list of candidate samples into a +serving plan a sensory lab can execute. Three constraints make a panel screen +different from an ordinary designed experiment, and each one is handled here: + +* **Session capacity.** An assessor can only judge a handful of samples in one + sitting before fatigue and carry-over dominate. When there are more candidates + than slots, every assessor sees only a *subset*, so the design is an + **incomplete block design** with the assessor-session as the block. + :func:`cyclic_block_design` builds the blocks so that replication is equal (or + as near-equal as the numbers allow) and every pair of candidates is compared + within a block about equally often - the pairwise **concurrence**, which is + what decides how precisely two candidates can be compared against each other. + +* **Carry-over and position.** The sample tasted before this one changes the + score of this one, and the first sample in a session is scored differently + from the last. :func:`williams_design` returns serving orders that are + balanced for first-order carry-over: every ordered pair of samples occurs + equally often, so the carry-over effect cannot bias any one sample's mean. + +* **Sensitivity.** Panel scores are noisy, so a screen that is too small will + return "no significant effect" whatever the truth is. + :func:`detectable_difference` converts a residual standard deviation and a + panel size into the smallest difference the screen can actually resolve, and + :func:`required_panelists` inverts it. + +:func:`sensory_screening_plan` assembles all three into one serving sheet, with +an optional reference (control) sample anchored in every block so drift between +sessions can be removed at the analysis stage. + +The output is deliberately in the same long shape the rest of the subpackage +consumes: once the scores are filled in next to ``product``, the sheet is ready +for :func:`process_improve.sensory.compare_products`. + +References +---------- +Williams, E. J. (1949). Experimental designs balanced for the estimation of +residual effects of treatments. *Australian Journal of Scientific Research*, +2(2), 149-168. + +Cochran, W. G. & Cox, G. M. (1957). *Experimental Designs* (2nd ed.). Wiley. +Chapters 9-11 on incomplete block designs. + +Naes, T., Brockhoff, P. B. & Tomic, O. (2010). *Statistics for Sensory and +Consumer Science*. Wiley. Chapter 2 on designing sensory experiments. + +MacFie, H. J., Bratchell, N., Greenhoff, K. & Vallis, L. V. (1989). Designs to +balance the effect of order of presentation and first-order carry-over effects +in hall tests. *Journal of Sensory Studies*, 4(2), 129-148. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Any + +import numpy as np +import pandas as pd +from scipy.stats import t as _student_t + +#: Columns of the serving sheet returned by :func:`sensory_screening_plan`. +SERVING_PLAN_COLUMNS: tuple[str, ...] = ( + "panelist_id", + "session", + "position", + "product", + "role", + "block", +) + +_MIN_TREATMENTS = 2 + + +@dataclass +class ScreeningPlan: + """Outcome of :func:`sensory_screening_plan`. + + Attributes + ---------- + plan : pandas.DataFrame + The serving sheet: one row per serving, with the columns of + :data:`SERVING_PLAN_COLUMNS`. ``block`` numbers the assessor-session + blocks consecutively; ``role`` is ``"control"`` or ``"test"``. + diagnostics : dict + Balance read-outs for the plan - ``replication`` (min/max/mean servings + per candidate), ``concurrence`` (min/max/mean number of blocks in which + a pair of candidates meet), ``balanced`` (True only for an exact + balanced incomplete block design), ``control_coverage`` (fraction of + blocks containing the reference), ``position_balance`` (largest spread + in how often a candidate lands in any one position) and + ``n_servings`` / ``n_blocks``. + config : dict + The resolved call arguments, for provenance. + warnings : list of str + Practical problems that do not stop the plan being produced - most + importantly a panel too small to cover every candidate. + """ + + plan: pd.DataFrame + diagnostics: dict[str, Any] + config: dict[str, Any] + warnings: list[str] = field(default_factory=list) + + +# --------------------------------------------------------------------------- +# Serving order: carry-over balance +# --------------------------------------------------------------------------- + + +def _williams_first_row(n_treatments: int) -> list[int]: + """Return the generating row ``0, 1, t-1, 2, t-2, ...`` of a Williams square.""" + row = [] + for position in range(n_treatments): + if position % 2 == 0: + row.append(position // 2) + else: + row.append(n_treatments - (position + 1) // 2) + return row + + +def williams_design( + n_treatments: int, + *, + n_subjects: int | None = None, + seed: int | None = None, +) -> pd.DataFrame: + """ + Build serving orders balanced for first-order carry-over (a Williams design). + + Each returned sequence is a permutation of all ``n_treatments`` treatments. + Across the full set of sequences, every treatment appears once in every + position (a Latin square) and every *ordered* pair of treatments occurs + equally often, so the effect of "what was tasted just before" is spread + evenly over the treatments instead of favouring some of them. + + An even ``n_treatments`` achieves this with ``n_treatments`` sequences; an + odd ``n_treatments`` needs a second, mirrored square, so ``2 * n_treatments`` + sequences are returned. + + Parameters + ---------- + n_treatments : int + Number of treatments to order, at least 2. + n_subjects : int, optional + Number of assessors to produce a sequence for. The balanced set of + sequences is recycled in order until this many are produced; balance is + exact only when ``n_subjects`` is a multiple of the number of sequences. + Defaults to the natural size of the design. + seed : int, optional + When given, the *assignment* of sequences to subjects is shuffled. The + balance properties are unaffected. + + Returns + ------- + pandas.DataFrame + Long format with columns ``sequence`` (0-based subject index), + ``position`` (1-based serving order) and ``treatment`` (0-based + treatment index). + + Examples + -------- + >>> design = williams_design(4) + >>> design.pivot(index="sequence", columns="position", values="treatment").to_numpy() + array([[0, 3, 1, 2], + [1, 0, 2, 3], + [2, 1, 3, 0], + [3, 2, 0, 1]]) + """ + if n_treatments < _MIN_TREATMENTS: + msg = f"A carry-over balanced order needs at least 2 treatments; got {n_treatments}." + raise ValueError(msg) + + first_row = _williams_first_row(n_treatments) + square = [[(value + shift) % n_treatments for value in first_row] for shift in range(n_treatments)] + if n_treatments % 2 == 1: + square += [list(reversed(row)) for row in square] + + sequences = square + if n_subjects is not None: + if n_subjects < 1: + msg = f"n_subjects must be at least 1; got {n_subjects}." + raise ValueError(msg) + sequences = [square[i % len(square)] for i in range(n_subjects)] + + if seed is not None: + rng = np.random.default_rng(seed) + order = rng.permutation(len(sequences)) + sequences = [sequences[i] for i in order] + + records = [ + {"sequence": subject, "position": position + 1, "treatment": treatment} + for subject, sequence in enumerate(sequences) + for position, treatment in enumerate(sequence) + ] + return pd.DataFrame.from_records(records) + + +# --------------------------------------------------------------------------- +# Incomplete blocks +# --------------------------------------------------------------------------- + + +def cyclic_block_design( + n_treatments: int, + *, + block_size: int, + n_blocks: int, + seed: int | None = None, +) -> list[list[int]]: + """ + Build ``n_blocks`` incomplete blocks of ``block_size`` treatments each. + + The blocks are filled greedily: each slot takes the treatment that has been + used least so far, breaking ties in favour of the treatment that has met the + treatments already in this block the fewest times. That drives the design + towards equal replication and equal pairwise concurrence, which is what a + balanced incomplete block design achieves exactly. When an exact BIBD exists + for the given ``(n_treatments, block_size, n_blocks)`` this construction + finds one; otherwise it returns the nearest thing the numbers allow, and + :func:`plan_diagnostics` reports how close that is. + + Parameters + ---------- + n_treatments : int + Number of distinct treatments to spread over the blocks. + block_size : int + Treatments per block. Must not exceed ``n_treatments``: a block cannot + contain a treatment twice. + n_blocks : int + Number of blocks to build. + seed : int, optional + Seed for the tie-breaking shuffle, making the result reproducible. + + Returns + ------- + list of list of int + One list of 0-based treatment indices per block. + + Examples + -------- + >>> blocks = cyclic_block_design(7, block_size=3, n_blocks=7, seed=0) + >>> plan_diagnostics(blocks, n_treatments=7)["balanced"] + True + """ + if n_treatments < _MIN_TREATMENTS: + msg = f"Need at least 2 treatments to block; got {n_treatments}." + raise ValueError(msg) + if block_size < 1 or block_size > n_treatments: + msg = f"block_size must be between 1 and n_treatments ({n_treatments}); got {block_size}." + raise ValueError(msg) + if n_blocks < 1: + msg = f"n_blocks must be at least 1; got {n_blocks}." + raise ValueError(msg) + + rng = np.random.default_rng(seed) + replication = np.zeros(n_treatments, dtype=int) + concurrence = np.zeros((n_treatments, n_treatments), dtype=int) + + blocks: list[list[int]] = [] + for _ in range(n_blocks): + block: list[int] = [] + for _slot in range(block_size): + candidates = [t for t in range(n_treatments) if t not in block] + # Least-used first; then least already-met inside this block; then random. + jitter = rng.random(n_treatments) + best = min( + candidates, + key=lambda t: (replication[t], sum(concurrence[t, other] for other in block), jitter[t]), + ) + block.append(best) + for treatment in block: + replication[treatment] += 1 + for i, first in enumerate(block): + for second in block[i + 1 :]: + concurrence[first, second] += 1 + concurrence[second, first] += 1 + blocks.append(block) + + _balance_by_swapping(blocks, n_treatments=n_treatments, replication=replication, concurrence=concurrence) + return [sorted(block) for block in blocks] + + +def _imbalance_penalty(n_treatments: int) -> float: + """Weight making equal replication strictly more important than equal concurrence.""" + return 10.0 * n_treatments + + +def _pair_changes(losers: list[tuple[int, list[int]]], gainers: list[tuple[int, list[int]]]) -> dict: + """Aggregate per-pair concurrence increments for a candidate move. + + ``losers``/``gainers`` are ``(treatment, partners)`` pairs: the treatment + stops (starts) sharing a block with each of its partners. Aggregating first + is what makes overlapping blocks safe - a pair touched by both halves of an + exchange nets out instead of being double-counted. + """ + changes: dict[tuple[int, int], int] = {} + for treatment, partners in losers: + for partner in partners: + key = (min(treatment, partner), max(treatment, partner)) + changes[key] = changes.get(key, 0) - 1 + for treatment, partners in gainers: + for partner in partners: + key = (min(treatment, partner), max(treatment, partner)) + changes[key] = changes.get(key, 0) + 1 + return changes + + +def _apply_move( + pair_changes: dict, + rep_changes: dict, + *, + replication: np.ndarray, + concurrence: np.ndarray, +) -> None: + """Commit the aggregated increments of a move to the running counters.""" + for treatment, change in rep_changes.items(): + replication[treatment] += change + for (first, second), change in pair_changes.items(): + concurrence[first, second] += change + concurrence[second, first] += change + + +def _move_delta( + pair_changes: dict, + rep_changes: dict, + *, + replication: np.ndarray, + concurrence: np.ndarray, + weight: float, +) -> float: + """Change in the imbalance objective if this move were applied.""" + delta = 0.0 + for treatment, change in rep_changes.items(): + current = replication[treatment] + delta += weight * (2 * current * change + change**2) + for (first, second), change in pair_changes.items(): + current = concurrence[first, second] + delta += 2 * current * change + change**2 + return delta + + +def _balance_by_swapping( + blocks: list[list[int]], + *, + n_treatments: int, + replication: np.ndarray, + concurrence: np.ndarray, + max_moves: int = 5_000, +) -> None: + """ + Improve the greedy blocks in place by hill-climbing on the imbalance. + + The objective is ``w * sum (r_i - rbar)^2 + sum_{i tuple[dict, dict, Any] | None: + """Return the first move that lowers the objective, or None when at a local optimum.""" + for index, block in enumerate(blocks): + for out_treatment in block: + others = [t for t in block if t != out_treatment] + for in_treatment in range(n_treatments): + if in_treatment in block: + continue + pair_changes = _pair_changes([(out_treatment, others)], [(in_treatment, others)]) + rep_changes = {out_treatment: -1, in_treatment: +1} + delta = _move_delta( + pair_changes, + rep_changes, + replication=replication, + concurrence=concurrence, + weight=weight, + ) + if delta < 0: + target, leaving, entering = block, out_treatment, in_treatment + + def commit(target: list[int] = target, leaving: int = leaving, entering: int = entering) -> None: + target[target.index(leaving)] = entering + + return pair_changes, rep_changes, commit + + for other_index in range(index + 1, len(blocks)): + other_block = blocks[other_index] + for first in block: + if first in other_block: + continue + rest_here = [t for t in block if t != first] + for second in other_block: + if second in block: + continue + rest_there = [t for t in other_block if t != second] + pair_changes = _pair_changes( + [(first, rest_here), (second, rest_there)], + [(second, rest_here), (first, rest_there)], + ) + delta = _move_delta( + pair_changes, + {}, + replication=replication, + concurrence=concurrence, + weight=weight, + ) + if delta < 0: + here, there, a, b = block, other_block, first, second + + def commit_exchange( + here: list[int] = here, + there: list[int] = there, + a: int = a, + b: int = b, + ) -> None: + here[here.index(a)] = b + there[there.index(b)] = a + + return pair_changes, {}, commit_exchange + return None + + +def plan_diagnostics(blocks: list[list[int]], *, n_treatments: int) -> dict[str, Any]: + """ + Summarise how balanced a set of blocks is. + + Parameters + ---------- + blocks : list of list of int + The blocks, as 0-based treatment indices (the output of + :func:`cyclic_block_design`). + n_treatments : int + Total number of treatments the blocks were drawn from. Treatments that + appear in no block are counted with a replication of zero, so a panel + too small to cover the candidate list shows up here. + + Returns + ------- + dict + ``replication`` and ``concurrence`` (each with ``min``, ``max``, + ``mean``), ``balanced`` (True when replication and off-diagonal + concurrence are both constant, i.e. an exact BIBD), ``n_blocks`` and + ``n_servings``. + """ + replication = np.zeros(n_treatments, dtype=int) + concurrence = np.zeros((n_treatments, n_treatments), dtype=int) + for block in blocks: + for treatment in block: + replication[treatment] += 1 + for i, first in enumerate(block): + for second in block[i + 1 :]: + concurrence[first, second] += 1 + concurrence[second, first] += 1 + + off_diagonal = concurrence[~np.eye(n_treatments, dtype=bool)] + balanced = bool(replication.min() == replication.max() and off_diagonal.min() == off_diagonal.max()) + return { + "replication": { + "min": int(replication.min()), + "max": int(replication.max()), + "mean": float(replication.mean()), + }, + "concurrence": { + "min": int(off_diagonal.min()), + "max": int(off_diagonal.max()), + "mean": float(off_diagonal.mean()), + }, + "balanced": balanced, + "n_blocks": len(blocks), + "n_servings": int(replication.sum()), + } + + +# --------------------------------------------------------------------------- +# The assembled serving plan +# --------------------------------------------------------------------------- + + +def _position_balance(plan: pd.DataFrame) -> float: + """Largest spread (max - min) in how often any candidate lands in one position.""" + tests = plan.loc[plan["role"] == "test"] + if tests.empty: + return 0.0 + counts = tests.pivot_table(index="product", columns="position", aggfunc="size", values="block").fillna(0) + return float((counts.max(axis=1) - counts.min(axis=1)).max()) + + +def _validate_screening_inputs( # noqa: PLR0913 - one guard per screening-plan knob + products: list[str], + *, + n_panelists: int, + samples_per_session: int, + control: str | None, + replicates: int, + n_sessions: int | None, +) -> None: + """Reject inputs that cannot produce a runnable serving plan.""" + if len(products) != len(set(products)): + msg = "The candidate list contains duplicate labels; each sample must appear once." + raise ValueError(msg) + if len(products) < _MIN_TREATMENTS: + msg = f"Screening needs at least 2 candidates; got {len(products)}." + raise ValueError(msg) + if control is not None and control in products: + msg = f"The control {control!r} must not also appear in the candidate list." + raise ValueError(msg) + if n_panelists < 1: + msg = f"n_panelists must be at least 1; got {n_panelists}." + raise ValueError(msg) + if replicates < 1: + msg = f"replicates must be at least 1; got {replicates}." + raise ValueError(msg) + if n_sessions is not None and n_sessions < 1: + msg = f"n_sessions must be at least 1; got {n_sessions}." + raise ValueError(msg) + if samples_per_session - (1 if control is not None else 0) < 1: + msg = ( + f"samples_per_session={samples_per_session} leaves no room for a candidate once the " + "control is anchored; increase it or drop the control." + ) + raise ValueError(msg) + + +def _capacity_warnings(*, n_test: int, n_blocks: int, slots_per_block: int, replicates: int) -> list[str]: + """Flag a panel too small to cover the candidate list at the requested replication.""" + capacity = n_blocks * slots_per_block + if capacity < n_test: + return [ + f"The panel cannot cover the candidate list: {n_blocks} blocks x {slots_per_block} slots = " + f"{capacity} servings for {n_test} candidates. Add sessions or assessors, or shorten the list." + ] + if capacity < n_test * replicates: + return [ + f"Capacity ({capacity} servings) is below the requested {n_test * replicates} " + f"(= {n_test} candidates x {replicates} replicates); realised replication will be lower." + ] + return [] + + +def _serving_records( + blocks: list[list[int]], + *, + products: list[str], + control: str | None, + n_panelists: int, + orders: pd.DataFrame | None, +) -> list[dict[str, Any]]: + """Expand the blocks into one record per serving, in the order they are served.""" + records: list[dict[str, Any]] = [] + for block_index, block in enumerate(blocks): + panelist = f"P{block_index % n_panelists + 1:02d}" + session = block_index // n_panelists + 1 + if orders is None: + ordered = list(block) + else: + sequence = orders.loc[orders["sequence"] == block_index].sort_values("position")["treatment"] + ordered = [block[i] for i in sequence] + + offset = 0 + if control is not None: + records.append( + { + "panelist_id": panelist, + "session": session, + "position": 1, + "product": control, + "role": "control", + "block": block_index + 1, + } + ) + offset = 1 + records.extend( + { + "panelist_id": panelist, + "session": session, + "position": position + 1 + offset, + "product": products[treatment], + "role": "test", + "block": block_index + 1, + } + for position, treatment in enumerate(ordered) + ) + return records + + +def sensory_screening_plan( # noqa: PLR0913 - each argument is a distinct, explicit design knob + products: list[str], + *, + n_panelists: int, + samples_per_session: int, + control: str | None = None, + replicates: int = 1, + n_sessions: int | None = None, + seed: int | None = None, +) -> ScreeningPlan: + """ + Build a blocked, carry-over balanced serving plan for a panel screen. + + Each assessor-session is one **block**. When the candidate list is longer + than a session can hold, the blocks are incomplete: an assessor sees a + subset, chosen by :func:`cyclic_block_design` so that replication and + pairwise concurrence stay as even as the arithmetic allows. Within a block + the serving order comes from :func:`williams_design`, so the sample tasted + beforehand does not systematically favour any candidate. + + Parameters + ---------- + products : list of str + The candidate samples to screen. Must be unique and must not include + ``control``. + n_panelists : int + Number of assessors on the panel. + samples_per_session : int + How many samples one assessor can judge in one session. When a + ``control`` is anchored it occupies one of these slots. + control : str, optional + Label of a reference sample served in every block, at the first + position. Anchoring a reference lets the analysis remove session and + assessor drift; leave it out to use the whole capacity for candidates. + replicates : int + Target number of times each candidate is served across the whole plan. + More sessions are scheduled until this is met, so the realised + replication is at least this and never more than one above the minimum. + n_sessions : int, optional + Fix the number of sessions per assessor instead of deriving it from + ``replicates``. When the fixed number cannot cover every candidate this + is reported in ``warnings`` rather than raising. + seed : int, optional + Seed making the plan reproducible. + + Returns + ------- + ScreeningPlan + The serving sheet plus balance diagnostics; see :class:`ScreeningPlan`. + + Examples + -------- + >>> plan = sensory_screening_plan( + ... [f"C{i}" for i in range(1, 22)], + ... n_panelists=12, + ... samples_per_session=6, + ... control="Base", + ... replicates=2, + ... seed=0, + ... ) + >>> plan.diagnostics["control_coverage"] + 1.0 + """ + _validate_screening_inputs( + products, + n_panelists=n_panelists, + samples_per_session=samples_per_session, + control=control, + replicates=replicates, + n_sessions=n_sessions, + ) + + n_test = len(products) + test_slots = samples_per_session - (1 if control is not None else 0) + slots_per_block = min(test_slots, n_test) + + if n_sessions is None: + blocks_needed = int(np.ceil(n_test * replicates / slots_per_block)) + n_sessions = max(1, int(np.ceil(blocks_needed / n_panelists))) + + n_blocks = n_panelists * n_sessions + warnings = _capacity_warnings( + n_test=n_test, n_blocks=n_blocks, slots_per_block=slots_per_block, replicates=replicates + ) + + blocks = cyclic_block_design(n_test, block_size=slots_per_block, n_blocks=n_blocks, seed=seed) + orders = williams_design(slots_per_block, n_subjects=n_blocks, seed=seed) if slots_per_block > 1 else None + + records = _serving_records(blocks, products=products, control=control, n_panelists=n_panelists, orders=orders) + plan = pd.DataFrame.from_records(records, columns=list(SERVING_PLAN_COLUMNS)) + plan = plan.sort_values(["block", "position"], ignore_index=True) + + diagnostics = plan_diagnostics(blocks, n_treatments=n_test) + diagnostics["control_coverage"] = 1.0 if control is not None else 0.0 + diagnostics["position_balance"] = _position_balance(plan) + diagnostics["n_sessions"] = n_sessions + diagnostics["slots_per_block"] = slots_per_block + + config = { + "n_products": n_test, + "n_panelists": n_panelists, + "samples_per_session": samples_per_session, + "control": control, + "replicates": replicates, + "n_sessions": n_sessions, + "seed": seed, + } + return ScreeningPlan(plan=plan, diagnostics=diagnostics, config=config, warnings=warnings) + + +# --------------------------------------------------------------------------- +# Sensitivity +# --------------------------------------------------------------------------- + + +def detectable_difference( + *, + sd: float, + n_per_product: int, + alpha: float = 0.05, + power: float = 0.80, + n_comparisons: int = 1, +) -> dict[str, Any]: + """ + Smallest difference between two samples the screen can resolve. + + Uses the usual two-sample normal-theory approximation with Student-t + quantiles, ``delta = (t_{1 - alpha', df} + t_{power, df}) * sd * sqrt(2 / n)`` + where ``df = 2 * (n - 1)`` and ``alpha' = alpha / (2 * n_comparisons)`` is + the Bonferroni-adjusted two-sided level. Because it uses the *residual* + standard deviation of the blocked model, the assessor-to-assessor variation + that blocking removes is already excluded, as it should be. + + This is an approximation, and a deliberately honest one: quote it as the + order of magnitude the screen can see, not a guarantee. + + Parameters + ---------- + sd : float + Residual standard deviation of a single score, on the scale the panel + uses. Take it from a previous study on the same attribute and scale. + n_per_product : int + Number of independent scores contributing to each sample's mean + (assessors x replicates), at least 2. + alpha : float + Two-sided significance level before any multiplicity adjustment. + power : float + Probability of detecting a difference of exactly this size. + n_comparisons : int + Size of the comparison family to Bonferroni-correct for; for example the + number of candidates each being compared against a single control. + + Returns + ------- + dict + ``difference`` (the minimum detectable difference), plus the ``sd``, + ``n_per_product``, ``alpha``, ``power`` and ``n_comparisons`` used, so + the number can be quoted with its assumptions attached. + + Examples + -------- + >>> round(detectable_difference(sd=1.5, n_per_product=12)["difference"], 3) + 1.796 + """ + if n_per_product < _MIN_TREATMENTS: + msg = f"n_per_product must be at least 2 to leave error degrees of freedom; got {n_per_product}." + raise ValueError(msg) + if sd <= 0: + msg = f"sd must be positive; got {sd}." + raise ValueError(msg) + if not 0 < alpha < 1: + msg = f"alpha must be strictly between 0 and 1; got {alpha}." + raise ValueError(msg) + if not 0 < power < 1: + msg = f"power must be strictly between 0 and 1; got {power}." + raise ValueError(msg) + if n_comparisons < 1: + msg = f"n_comparisons must be at least 1; got {n_comparisons}." + raise ValueError(msg) + + df = 2 * (n_per_product - 1) + t_alpha = float(_student_t.ppf(1 - alpha / (2 * n_comparisons), df)) + t_beta = float(_student_t.ppf(power, df)) + difference = (t_alpha + t_beta) * sd * np.sqrt(2 / n_per_product) + return { + "difference": float(difference), + "sd": float(sd), + "n_per_product": int(n_per_product), + "alpha": float(alpha), + "power": float(power), + "n_comparisons": int(n_comparisons), + } + + +def required_panelists( # noqa: PLR0913 - explicit sizing knobs, each distinct + *, + sd: float, + difference: float, + alpha: float = 0.05, + power: float = 0.80, + n_comparisons: int = 1, + max_n: int = 10_000, +) -> dict[str, Any]: + """ + Smallest ``n_per_product`` whose detectable difference reaches a target. + + Inverts :func:`detectable_difference` by search, so the two always agree. + + Parameters + ---------- + sd : float + Residual standard deviation of a single score. + difference : float + The difference that must be detectable, on the same scale as ``sd``. + alpha, power, n_comparisons + As in :func:`detectable_difference`. + max_n : int + Search ceiling; exceeding it raises rather than looping forever. + + Returns + ------- + dict + ``n_per_product`` (the answer) and ``achieved_difference`` (what that n + actually resolves, which is at or below ``difference``), plus the + inputs. + + Examples + -------- + >>> required_panelists(sd=1.5, difference=1.0)["n_per_product"] + 37 + """ + if difference <= 0: + msg = f"difference must be positive; got {difference}." + raise ValueError(msg) + + for n in range(2, max_n + 1): + achieved = detectable_difference(sd=sd, n_per_product=n, alpha=alpha, power=power, n_comparisons=n_comparisons)[ + "difference" + ] + if achieved <= difference: + return { + "n_per_product": n, + "achieved_difference": achieved, + "sd": float(sd), + "difference": float(difference), + "alpha": float(alpha), + "power": float(power), + "n_comparisons": int(n_comparisons), + } + + msg = f"No n_per_product up to {max_n} reaches a detectable difference of {difference} with sd={sd}." + raise ValueError(msg) diff --git a/src/process_improve/sensory/tools.py b/src/process_improve/sensory/tools.py index 8f5d2d8f..13e2c314 100644 --- a/src/process_improve/sensory/tools.py +++ b/src/process_improve/sensory/tools.py @@ -21,6 +21,9 @@ from process_improve.sensory.mam import align_scores as _align_scores from process_improve.sensory.mam import mixed_assessor_model as _mixed_assessor_model from process_improve.sensory.panel import panel_scorecard as _panel_scorecard +from process_improve.sensory.screening import detectable_difference as _detectable_difference +from process_improve.sensory.screening import required_panelists as _required_panelists +from process_improve.sensory.screening import sensory_screening_plan as _sensory_screening_plan from process_improve.sensory.validation import DESCRIPTIVE_LONG_COLUMNS from process_improve.sensory.validation import validate_descriptive as _validate_descriptive from process_improve.tool_spec import clean, get_tool_specs, tool_spec @@ -235,9 +238,7 @@ class _AnalyzeInput(BaseModel): "Set false to skip it (faster)." ), ) - n_permutations: int = Field( - 199, ge=1, description="Permutations for the discriminator's selectivity-ratio null." - ) + n_permutations: int = Field(199, ge=1, description="Permutations for the discriminator's selectivity-ratio null.") random_state: int = Field(0, description="Seed for the discriminator's permutations and CV folds.") score_min: float | None = Field(None, description="Optional lower bound for the score scale.") score_max: float | None = Field(None, description="Optional upper bound for the score scale.") @@ -394,10 +395,176 @@ def sensory_panel_check(spec: _PanelCheckInput) -> dict: return clean(out) +class _ScreeningPlanInput(BaseModel): + """Input contract for ``sensory_screening_plan``.""" + + model_config = ConfigDict(extra="forbid") + + products: list[str] = Field( + ..., + min_length=2, + description="Labels of the candidate samples to screen. Must be unique and exclude the control.", + ) + n_panelists: int = Field(..., ge=1, description="Number of assessors on the panel.") + samples_per_session: int = Field( + ..., + ge=1, + description=( + "How many samples one assessor can judge in a single session before fatigue and carry-over " + "dominate. The anchored control, if any, occupies one of these slots." + ), + ) + control: str | None = Field( + None, + description=( + "Label of a reference sample served first in every block, so session and assessor drift can be " + "removed at the analysis stage. Omit to spend the whole capacity on candidates." + ), + ) + replicates: int = Field( + 1, + ge=1, + description="Target number of times each candidate is served across the whole plan.", + ) + n_sessions: int | None = Field( + None, + ge=1, + description="Fix the sessions per assessor instead of deriving them from `replicates`.", + ) + seed: int | None = Field(None, description="Seed making the plan reproducible.") + + +@tool_spec( + name="sensory_screening_plan", + description=( + "Design a blocked, carry-over balanced serving plan for a sensory screen of many candidate " + "samples. Use this BEFORE any tasting, when there are more candidates than one assessor can " + "judge in a session: it splits them into incomplete blocks (one block = one assessor-session) " + "with near-equal replication and near-equal pairwise concurrence, orders each block by a " + "Williams design so the previously tasted sample does not bias any candidate, and optionally " + "anchors a reference sample first in every block. " + "Returns: {ok: true, plan, diagnostics, config, warnings}. 'plan' is the serving sheet, one row " + "per serving (panelist_id, session, position, product, role, block), ready to have scores added " + "and be analysed with sensory_compare_products. 'diagnostics' reports replication and " + "concurrence (min/max/mean), whether the blocks form an exact balanced incomplete block design, " + "control coverage and position balance. 'warnings' flags a panel too small to cover the " + "candidate list rather than silently dropping candidates." + ), + input_model=_ScreeningPlanInput, + category="sensory", +) +def sensory_screening_plan(spec: _ScreeningPlanInput) -> dict: + """Blocked, carry-over balanced serving plan; see tool spec for details.""" + try: + result = _sensory_screening_plan( + list(spec.products), + n_panelists=spec.n_panelists, + samples_per_session=spec.samples_per_session, + control=spec.control, + replicates=spec.replicates, + n_sessions=spec.n_sessions, + seed=spec.seed, + ) + except ValueError as exc: + return clean({"ok": False, "errors": [str(exc)]}) + return clean( + { + "ok": True, + "plan": result.plan.to_dict(orient="records"), + "diagnostics": result.diagnostics, + "config": result.config, + "warnings": result.warnings, + } + ) + + +class _DetectableDifferenceInput(BaseModel): + """Input contract for ``sensory_detectable_difference``.""" + + model_config = ConfigDict(extra="forbid") + + sd: float = Field( + ..., + gt=0, + description=( + "Residual standard deviation of a single score on the panel's scale, from a previous study " + "on the same attribute and scale. This is the within-panel noise after blocking on assessor." + ), + ) + n_per_product: int | None = Field( + None, + ge=2, + description="Scores contributing to each sample's mean (assessors x replicates). Give this OR `difference`.", + ) + difference: float | None = Field( + None, + gt=0, + description="A difference that must be detectable. Give this to size the panel instead of `n_per_product`.", + ) + alpha: float = Field(0.05, gt=0, lt=1, description="Two-sided significance level before multiplicity.") + power: float = Field(0.80, gt=0, lt=1, description="Probability of detecting a difference of that size.") + n_comparisons: int = Field( + 1, + ge=1, + description="Size of the comparison family to Bonferroni-correct for, e.g. every candidate versus one control.", + ) + + +@tool_spec( + name="sensory_detectable_difference", + description=( + "Size a sensory screen, or say what it can see. Given the residual standard deviation of a " + "single score, either (a) pass `n_per_product` to get the smallest difference that panel size " + "can resolve, or (b) pass `difference` to get the number of scores per sample needed to resolve " + "it. Use this before committing to a panel size, and to state honestly what a null result " + "means. Correct for the comparison family with `n_comparisons` (e.g. 20 candidates each tested " + "against one control). " + "Returns: {ok: true, difference, n_per_product, sd, alpha, power, n_comparisons}, plus " + "'achieved_difference' when sizing from a target difference. Approximate: normal-theory " + "two-sample formula with Student-t quantiles, so quote it as an order of magnitude." + ), + input_model=_DetectableDifferenceInput, + category="sensory", +) +def sensory_detectable_difference(spec: _DetectableDifferenceInput) -> dict: + """Minimum detectable difference, or the panel size that reaches one; see tool spec.""" + if (spec.n_per_product is None) == (spec.difference is None): + return clean( + { + "ok": False, + "errors": [ + "Give exactly one of `n_per_product` (what can I see?) or `difference` (how many do I need?)." + ], + } + ) + try: + if spec.n_per_product is not None: + out = _detectable_difference( + sd=spec.sd, + n_per_product=spec.n_per_product, + alpha=spec.alpha, + power=spec.power, + n_comparisons=spec.n_comparisons, + ) + else: + out = _required_panelists( + sd=spec.sd, + difference=spec.difference, + alpha=spec.alpha, + power=spec.power, + n_comparisons=spec.n_comparisons, + ) + except ValueError as exc: + return clean({"ok": False, "errors": [str(exc)]}) + return clean({"ok": True, **out}) + + _register("sensory_reshape_to_long") _register("sensory_validate_descriptive") _register("sensory_analyze_descriptive") _register("sensory_panel_check") +_register("sensory_screening_plan") +_register("sensory_detectable_difference") def get_sensory_tool_specs() -> list[dict]: diff --git a/tests/test_sensory_screening.py b/tests/test_sensory_screening.py new file mode 100644 index 00000000..eedfbbb4 --- /dev/null +++ b/tests/test_sensory_screening.py @@ -0,0 +1,359 @@ +"""(c) Kevin Dunn, 2010-2026. MIT License. + +Tests for :mod:`process_improve.sensory.screening`: carry-over balanced serving +orders, near-balanced incomplete blocks, the assembled panel serving plan, and +the minimum-detectable-difference sizing helpers. +""" + +from __future__ import annotations + +import json +from collections import Counter +from itertools import pairwise + +import pandas as pd +import pytest + +from process_improve.sensory.screening import ( + ScreeningPlan, + cyclic_block_design, + detectable_difference, + plan_diagnostics, + required_panelists, + sensory_screening_plan, + williams_design, +) +from process_improve.sensory.tools import ( + _DetectableDifferenceInput, + _ScreeningPlanInput, + get_sensory_tool_specs, +) +from process_improve.sensory.tools import ( + sensory_detectable_difference as detectable_difference_tool, +) +from process_improve.sensory.tools import ( + sensory_screening_plan as screening_plan_tool, +) + +# --------------------------------------------------------------------------- +# williams_design +# --------------------------------------------------------------------------- + + +def _ordered_pairs(sequences: list[list[int]]) -> Counter: + """Count every ordered (predecessor, successor) pair across the sequences.""" + pairs: Counter = Counter() + for seq in sequences: + for before, after in pairwise(seq): + pairs[before, after] += 1 + return pairs + + +def _sequences(design: pd.DataFrame) -> list[list[int]]: + return [grp.sort_values("position")["treatment"].tolist() for _, grp in design.groupby("sequence", sort=True)] + + +@pytest.mark.parametrize("n_treatments", [2, 4, 6, 8]) +def test_williams_design_even_is_carryover_balanced(n_treatments: int) -> None: + """For an even number of treatments each ordered pair follows exactly once.""" + design = williams_design(n_treatments) + sequences = _sequences(design) + + assert len(sequences) == n_treatments + assert all(sorted(seq) == list(range(n_treatments)) for seq in sequences) + + pairs = _ordered_pairs(sequences) + expected = {(i, j) for i in range(n_treatments) for j in range(n_treatments) if i != j} + assert set(pairs) == expected + assert set(pairs.values()) == {1} + + +@pytest.mark.parametrize("n_treatments", [3, 5, 7]) +def test_williams_design_odd_uses_two_squares_and_balances(n_treatments: int) -> None: + """An odd treatment count needs 2t sequences; each ordered pair then follows twice.""" + design = williams_design(n_treatments) + sequences = _sequences(design) + + assert len(sequences) == 2 * n_treatments + + pairs = _ordered_pairs(sequences) + expected = {(i, j) for i in range(n_treatments) for j in range(n_treatments) if i != j} + assert set(pairs) == expected + assert set(pairs.values()) == {2} + + +def test_williams_design_every_treatment_appears_once_per_position() -> None: + """The design is a Latin square: each treatment occupies each position once.""" + design = williams_design(6) + counts = design.pivot_table(index="treatment", columns="position", aggfunc="size", values="sequence") + assert (counts.to_numpy() == 1).all() + + +def test_williams_design_n_subjects_cycles_the_sequences() -> None: + """Asking for more subjects than sequences recycles them in order.""" + design = williams_design(4, n_subjects=10) + assert design["sequence"].nunique() == 10 + assert len(design) == 40 + + +def test_williams_design_rejects_too_few_treatments() -> None: + """Fewer than two treatments has no ordering to balance.""" + with pytest.raises(ValueError, match="at least 2"): + williams_design(1) + + +# --------------------------------------------------------------------------- +# cyclic_block_design +# --------------------------------------------------------------------------- + + +def test_cyclic_block_design_recovers_a_known_bibd() -> None: + """t=7, k=3, b=7 is the Fano plane: r=3 and every pair concurs exactly once.""" + blocks = cyclic_block_design(7, block_size=3, n_blocks=7, seed=0) + assert len(blocks) == 7 + assert all(len(b) == 3 for b in blocks) + assert all(len(set(b)) == 3 for b in blocks) + + diag = plan_diagnostics(blocks, n_treatments=7) + assert diag["replication"]["min"] == diag["replication"]["max"] == 3 + assert diag["concurrence"]["min"] == diag["concurrence"]["max"] == 1 + assert diag["balanced"] is True + + +def test_cyclic_block_design_keeps_replication_near_equal_when_not_a_bibd() -> None: + """An arbitrary (t, k, b) has no exact BIBD; replication must still be near-equal.""" + blocks = cyclic_block_design(20, block_size=6, n_blocks=10, seed=1) + diag = plan_diagnostics(blocks, n_treatments=20) + + assert sum(len(b) for b in blocks) == 60 + assert diag["replication"]["max"] - diag["replication"]["min"] <= 1 + assert diag["balanced"] is False + + +def test_cyclic_block_design_is_deterministic_for_a_seed() -> None: + """The same seed reproduces the same blocks; a different seed may differ.""" + assert cyclic_block_design(12, block_size=4, n_blocks=9, seed=7) == cyclic_block_design( + 12, block_size=4, n_blocks=9, seed=7 + ) + + +def test_cyclic_block_design_rejects_block_larger_than_treatments() -> None: + """A block cannot hold more distinct treatments than exist.""" + with pytest.raises(ValueError, match="block_size"): + cyclic_block_design(4, block_size=5, n_blocks=3) + + +# --------------------------------------------------------------------------- +# sensory_screening_plan +# --------------------------------------------------------------------------- + +CANDIDATES = [f"C{i:02d}" for i in range(1, 22)] + + +def test_screening_plan_covers_every_candidate_at_least_the_requested_replicates() -> None: + """Every candidate is served at least `replicates` times, and near-equally.""" + result = sensory_screening_plan( + CANDIDATES, + n_panelists=12, + samples_per_session=6, + control="Base", + replicates=2, + seed=0, + ) + assert isinstance(result, ScreeningPlan) + served = result.plan.loc[result.plan["role"] == "test", "product"].value_counts() + assert set(served.index) == set(CANDIDATES) + assert served.min() >= 2 + assert served.max() - served.min() <= 1 + assert result.diagnostics["replication"]["min"] == served.min() + + +def test_screening_plan_puts_the_control_in_every_block() -> None: + """The reference anchors every panelist-session block, at the first position.""" + result = sensory_screening_plan( + CANDIDATES, + n_panelists=10, + samples_per_session=6, + control="Base", + seed=0, + ) + blocks = result.plan.groupby(["panelist_id", "session"]) + for _, block in blocks: + assert (block["product"] == "Base").sum() == 1 + assert block.loc[block["product"] == "Base", "position"].iloc[0] == 1 + assert result.diagnostics["control_coverage"] == pytest.approx(1.0) + + +def test_screening_plan_respects_the_session_capacity() -> None: + """No panelist ever sees more than `samples_per_session` samples in a session.""" + result = sensory_screening_plan(CANDIDATES, n_panelists=8, samples_per_session=5, control="Base", seed=2) + sizes = result.plan.groupby(["panelist_id", "session"]).size() + assert sizes.max() <= 5 + assert result.plan["position"].max() <= 5 + + +def test_screening_plan_never_repeats_a_product_within_a_block() -> None: + """A panelist does not taste the same sample twice in one session.""" + result = sensory_screening_plan(CANDIDATES, n_panelists=9, samples_per_session=7, control="Base", seed=3) + for _, block in result.plan.groupby(["panelist_id", "session"]): + assert block["product"].is_unique + + +def test_screening_plan_without_a_control_uses_the_whole_capacity() -> None: + """With no control every slot in the block is a test sample.""" + result = sensory_screening_plan(CANDIDATES, n_panelists=7, samples_per_session=6, seed=4) + assert (result.plan["role"] == "test").all() + assert result.diagnostics["control_coverage"] == 0.0 + + +def test_screening_plan_is_reproducible_for_a_seed() -> None: + """The same seed gives an identical plan.""" + kwargs = {"n_panelists": 6, "samples_per_session": 5, "control": "Base", "seed": 11} + first = sensory_screening_plan(CANDIDATES, **kwargs).plan + second = sensory_screening_plan(CANDIDATES, **kwargs).plan + pd.testing.assert_frame_equal(first, second) + + +def test_screening_plan_warns_when_the_panel_is_too_small_for_the_candidates() -> None: + """A capacity shortfall is reported rather than silently dropping candidates.""" + result = sensory_screening_plan( + CANDIDATES, + n_panelists=2, + samples_per_session=4, + control="Base", + n_sessions=1, + seed=5, + ) + assert result.warnings + assert any("cannot cover" in w or "not every" in w.lower() for w in result.warnings) + + +def test_screening_plan_rejects_a_capacity_of_one_when_a_control_is_anchored() -> None: + """A control that fills the only slot leaves no room for a test sample.""" + with pytest.raises(ValueError, match="samples_per_session"): + sensory_screening_plan(CANDIDATES, n_panelists=5, samples_per_session=1, control="Base") + + +def test_screening_plan_rejects_duplicate_candidates() -> None: + """Duplicated labels would silently double a candidate's replication.""" + with pytest.raises(ValueError, match="duplicate"): + sensory_screening_plan(["A", "B", "A"], n_panelists=4, samples_per_session=3) + + +def test_screening_plan_carries_the_long_format_columns() -> None: + """The plan is a serving sheet, ready to be joined to scores.""" + result = sensory_screening_plan(CANDIDATES, n_panelists=6, samples_per_session=6, control="Base", seed=6) + assert list(result.plan.columns) == [ + "panelist_id", + "session", + "position", + "product", + "role", + "block", + ] + assert result.config["n_panelists"] == 6 + + +# --------------------------------------------------------------------------- +# detectable_difference / required_panelists +# --------------------------------------------------------------------------- + + +def test_detectable_difference_shrinks_as_the_panel_grows() -> None: + """More assessors per product means a smaller difference is detectable.""" + small = detectable_difference(sd=1.5, n_per_product=8)["difference"] + large = detectable_difference(sd=1.5, n_per_product=32)["difference"] + assert large < small + + +def test_detectable_difference_grows_with_more_comparisons() -> None: + """Correcting for a family of comparisons costs sensitivity.""" + one = detectable_difference(sd=1.5, n_per_product=16, n_comparisons=1)["difference"] + many = detectable_difference(sd=1.5, n_per_product=16, n_comparisons=21)["difference"] + assert many > one + + +def test_detectable_difference_scales_linearly_with_the_noise() -> None: + """Doubling the residual SD doubles the minimum detectable difference.""" + base = detectable_difference(sd=1.0, n_per_product=20)["difference"] + doubled = detectable_difference(sd=2.0, n_per_product=20)["difference"] + assert doubled == pytest.approx(2 * base, rel=1e-9) + + +def test_detectable_difference_reports_its_inputs() -> None: + """The result is self-describing so it can be quoted in a report.""" + out = detectable_difference(sd=1.2, n_per_product=12, alpha=0.05, power=0.8) + assert out["sd"] == 1.2 + assert out["n_per_product"] == 12 + assert out["power"] == 0.8 + assert out["difference"] > 0 + + +def test_detectable_difference_rejects_a_single_assessor() -> None: + """With one assessor per product there are no error degrees of freedom.""" + with pytest.raises(ValueError, match="n_per_product"): + detectable_difference(sd=1.0, n_per_product=1) + + +def test_required_panelists_delivers_the_requested_difference() -> None: + """The returned n actually reaches the target difference, and n-1 does not.""" + target = 0.8 + n = required_panelists(sd=1.5, difference=target)["n_per_product"] + assert detectable_difference(sd=1.5, n_per_product=n)["difference"] <= target + assert detectable_difference(sd=1.5, n_per_product=n - 1)["difference"] > target + + +def test_required_panelists_rejects_a_non_positive_difference() -> None: + """A zero or negative target difference is never reachable.""" + with pytest.raises(ValueError, match="difference"): + required_panelists(sd=1.0, difference=0.0) + + +# --------------------------------------------------------------------------- +# Agent-callable tool wrappers +# --------------------------------------------------------------------------- + + +def test_screening_plan_tool_returns_json_serialisable_rows() -> None: + """The tool hands back plain records the agent can pass straight on.""" + result = screening_plan_tool( + _ScreeningPlanInput(products=CANDIDATES, n_panelists=12, samples_per_session=6, control="Base", seed=0) + ) + assert result["ok"] is True + assert isinstance(result["plan"], list) + assert set(result["plan"][0]) == {"panelist_id", "session", "position", "product", "role", "block"} + assert result["diagnostics"]["control_coverage"] == 1.0 + json.dumps(result) + + +def test_screening_plan_tool_reports_bad_input_instead_of_raising() -> None: + """A design that cannot be built comes back as ok=False so the agent can relay it.""" + result = screening_plan_tool( + _ScreeningPlanInput(products=CANDIDATES, n_panelists=5, samples_per_session=1, control="Base") + ) + assert result["ok"] is False + assert any("samples_per_session" in error for error in result["errors"]) + + +def test_detectable_difference_tool_answers_both_questions() -> None: + """One argument asks what the panel can see; the other asks how big it must be.""" + seeing = detectable_difference_tool(_DetectableDifferenceInput(sd=1.5, n_per_product=24)) + sizing = detectable_difference_tool(_DetectableDifferenceInput(sd=1.5, difference=1.0)) + assert seeing["ok"] is True + assert seeing["difference"] > 0 + assert sizing["ok"] is True + assert sizing["achieved_difference"] <= 1.0 + + +def test_detectable_difference_tool_refuses_both_or_neither() -> None: + """Giving both (or neither) is ambiguous and is rejected, not guessed at.""" + both = detectable_difference_tool(_DetectableDifferenceInput(sd=1.5, n_per_product=24, difference=1.0)) + neither = detectable_difference_tool(_DetectableDifferenceInput(sd=1.5)) + assert both["ok"] is False + assert neither["ok"] is False + + +def test_screening_tools_are_registered_for_the_agent() -> None: + """Both tools appear in the sensory tool specs so an LLM can call them.""" + names = {spec["name"] for spec in get_sensory_tool_specs()} + assert {"sensory_screening_plan", "sensory_detectable_difference"} <= names From e7d5b0e0de1be6381b01077e9fd9bcdb45c61120 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 3 Aug 2026 12:22:38 +0000 Subject: [PATCH 2/5] Narrow the optional sizing arguments so mypy can check the tool The "exactly one of n_per_product / difference" guard was written as a single equality on two `is None` tests, which mypy cannot use to narrow either field, so `difference: float | None` reached `required_panelists`. Written as two positive tests instead, each branch narrows at its own call site and the ambiguous case falls through to the same error response. Also restores a Field() call that the formatter had collapsed; it is not part of this change. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TqdSBoXb1DpSobB76pMzFi --- src/process_improve/sensory/tools.py | 29 ++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/process_improve/sensory/tools.py b/src/process_improve/sensory/tools.py index 13e2c314..68ae03d1 100644 --- a/src/process_improve/sensory/tools.py +++ b/src/process_improve/sensory/tools.py @@ -238,7 +238,9 @@ class _AnalyzeInput(BaseModel): "Set false to skip it (faster)." ), ) - n_permutations: int = Field(199, ge=1, description="Permutations for the discriminator's selectivity-ratio null.") + n_permutations: int = Field( + 199, ge=1, description="Permutations for the discriminator's selectivity-ratio null." + ) random_state: int = Field(0, description="Seed for the discriminator's permutations and CV folds.") score_min: float | None = Field(None, description="Optional lower bound for the score scale.") score_max: float | None = Field(None, description="Optional upper bound for the score scale.") @@ -528,17 +530,10 @@ class _DetectableDifferenceInput(BaseModel): ) def sensory_detectable_difference(spec: _DetectableDifferenceInput) -> dict: """Minimum detectable difference, or the panel size that reaches one; see tool spec.""" - if (spec.n_per_product is None) == (spec.difference is None): - return clean( - { - "ok": False, - "errors": [ - "Give exactly one of `n_per_product` (what can I see?) or `difference` (how many do I need?)." - ], - } - ) try: - if spec.n_per_product is not None: + # Written as two positive tests rather than one "exactly one is set" guard so the + # optional fields are narrowed for the type checker at each call site. + if spec.n_per_product is not None and spec.difference is None: out = _detectable_difference( sd=spec.sd, n_per_product=spec.n_per_product, @@ -546,7 +541,7 @@ def sensory_detectable_difference(spec: _DetectableDifferenceInput) -> dict: power=spec.power, n_comparisons=spec.n_comparisons, ) - else: + elif spec.difference is not None and spec.n_per_product is None: out = _required_panelists( sd=spec.sd, difference=spec.difference, @@ -554,6 +549,16 @@ def sensory_detectable_difference(spec: _DetectableDifferenceInput) -> dict: power=spec.power, n_comparisons=spec.n_comparisons, ) + else: + return clean( + { + "ok": False, + "errors": [ + "Give exactly one of `n_per_product` (what can I see?) or `difference` " + "(how many do I need?)." + ], + } + ) except ValueError as exc: return clean({"ok": False, "errors": [str(exc)]}) return clean({"ok": True, **out}) From 8a2f3c5e99c89836b16e3389167385d84fb1d3fb Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 3 Aug 2026 12:26:18 +0000 Subject: [PATCH 3/5] Fix the title underline length in the new docs section Sphinx builds with -W, so the one-character-short underline under "Before the panel: designing the screen" failed the docs job. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TqdSBoXb1DpSobB76pMzFi --- docs/user_guide/sensory_panel.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user_guide/sensory_panel.rst b/docs/user_guide/sensory_panel.rst index 63cc1bca..2a93e6b5 100644 --- a/docs/user_guide/sensory_panel.rst +++ b/docs/user_guide/sensory_panel.rst @@ -407,7 +407,7 @@ with correctly calibrated thresholds. 13(2), 456-466. Before the panel: designing the screen -------------------------------------- +-------------------------------------- Everything above analyses panel data that already exists. When you are at the other end - a long list of candidate samples and no scores yet - the question is From b6a39a12a822def1f1d7a66f6f499ebc70a8b61d Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 3 Aug 2026 12:36:53 +0000 Subject: [PATCH 4/5] Cover the screening guard clauses, and de-couple a tool_spec test Two fixes for CI on this branch. 1. `test_registered_in_registry` asserted that dummies registered by two OTHER tests in the class were in the global registry. That only holds when all three run in the same process, and pytest-xdist distributes them across workers, so it was one scheduling change away from failing. Adding tools to the registry was that change. It now registers its own tool and asserts on that, so it is independently runnable. 2. Screening coverage 88% -> 99%. The uncovered lines were the guard clauses, and those are worth testing individually: the message is the deliverable, because a scientist who asks for an impossible design needs to know which number to change. Also covers the tool's ValueError path (an unreachable target difference comes back as ok=False rather than raising). Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TqdSBoXb1DpSobB76pMzFi --- tests/test_sensory_screening.py | 105 ++++++++++++++++++++++++++++++++ tests/test_tool_spec.py | 21 ++++++- 2 files changed, 123 insertions(+), 3 deletions(-) diff --git a/tests/test_sensory_screening.py b/tests/test_sensory_screening.py index eedfbbb4..c383bcb3 100644 --- a/tests/test_sensory_screening.py +++ b/tests/test_sensory_screening.py @@ -357,3 +357,108 @@ def test_screening_tools_are_registered_for_the_agent() -> None: """Both tools appear in the sensory tool specs so an LLM can call them.""" names = {spec["name"] for spec in get_sensory_tool_specs()} assert {"sensory_screening_plan", "sensory_detectable_difference"} <= names + + +# --------------------------------------------------------------------------- +# Guard clauses +# +# Every one of these is a way to ask for a design that cannot exist. They are +# tested individually because the message is the deliverable: a scientist who +# gets these back needs to know which number to change. +# --------------------------------------------------------------------------- + + +def test_williams_design_rejects_a_non_positive_subject_count() -> None: + """Zero subjects is not a design.""" + with pytest.raises(ValueError, match="n_subjects"): + williams_design(4, n_subjects=0) + + +def test_cyclic_block_design_rejects_fewer_than_two_treatments() -> None: + """There is nothing to block with a single treatment.""" + with pytest.raises(ValueError, match="at least 2 treatments"): + cyclic_block_design(1, block_size=1, n_blocks=3) + + +def test_cyclic_block_design_rejects_an_empty_block() -> None: + """A block of size zero would serve nothing.""" + with pytest.raises(ValueError, match="block_size"): + cyclic_block_design(6, block_size=0, n_blocks=3) + + +def test_cyclic_block_design_rejects_a_non_positive_block_count() -> None: + """Asking for no blocks is a caller error, not an empty design.""" + with pytest.raises(ValueError, match="n_blocks"): + cyclic_block_design(6, block_size=3, n_blocks=0) + + +def test_screening_plan_rejects_an_empty_panel() -> None: + """No assessors means no plan.""" + with pytest.raises(ValueError, match="n_panelists"): + sensory_screening_plan(CANDIDATES, n_panelists=0, samples_per_session=5) + + +def test_screening_plan_rejects_non_positive_replicates() -> None: + """Every candidate must be served at least once.""" + with pytest.raises(ValueError, match="replicates"): + sensory_screening_plan(CANDIDATES, n_panelists=6, samples_per_session=5, replicates=0) + + +def test_screening_plan_rejects_a_non_positive_session_count() -> None: + """Fixing the sessions at zero contradicts asking for a plan.""" + with pytest.raises(ValueError, match="n_sessions"): + sensory_screening_plan(CANDIDATES, n_panelists=6, samples_per_session=5, n_sessions=0) + + +def test_screening_plan_rejects_a_single_candidate() -> None: + """One candidate is not a screen.""" + with pytest.raises(ValueError, match="at least 2 candidates"): + sensory_screening_plan(["Only one"], n_panelists=6, samples_per_session=5) + + +def test_screening_plan_rejects_a_control_that_is_also_a_candidate() -> None: + """The reference cannot double as a treatment; its replication would be wrong.""" + with pytest.raises(ValueError, match="must not also appear"): + sensory_screening_plan([*CANDIDATES, "Base"], n_panelists=6, samples_per_session=5, control="Base") + + +def test_screening_plan_handles_a_single_slot_per_block() -> None: + """With one test slot per block there is no order to balance, and that is fine.""" + result = sensory_screening_plan(CANDIDATES, n_panelists=6, samples_per_session=2, control="Base", seed=12) + + for _, block in result.plan.groupby("block"): + assert (block["role"] == "test").sum() == 1 + assert result.diagnostics["position_balance"] == 0.0 + + +def test_detectable_difference_rejects_a_non_positive_sd() -> None: + """A zero or negative noise estimate is not a noise estimate.""" + with pytest.raises(ValueError, match="sd"): + detectable_difference(sd=0.0, n_per_product=10) + + +@pytest.mark.parametrize(("alpha", "power"), [(0.0, 0.8), (1.0, 0.8), (0.05, 0.0), (0.05, 1.0)]) +def test_detectable_difference_rejects_out_of_range_rates(alpha: float, power: float) -> None: + """Both alpha and power are probabilities strictly inside (0, 1).""" + with pytest.raises(ValueError, match=r"alpha|power"): + detectable_difference(sd=1.0, n_per_product=10, alpha=alpha, power=power) + + +def test_detectable_difference_rejects_an_empty_comparison_family() -> None: + """There is always at least one comparison being made.""" + with pytest.raises(ValueError, match="n_comparisons"): + detectable_difference(sd=1.0, n_per_product=10, n_comparisons=0) + + +def test_required_panelists_refuses_an_unreachable_target() -> None: + """A target the search ceiling cannot reach raises rather than looping or lying.""" + with pytest.raises(ValueError, match="No n_per_product up to"): + required_panelists(sd=10.0, difference=0.001, max_n=20) + + +def test_detectable_difference_tool_relays_an_unreachable_target() -> None: + """A target no panel size can reach comes back as ok=False, not as an exception.""" + result = detectable_difference_tool(_DetectableDifferenceInput(sd=1e6, difference=1e-6)) + + assert result["ok"] is False + assert any("No n_per_product up to" in error for error in result["errors"]) diff --git a/tests/test_tool_spec.py b/tests/test_tool_spec.py index 81c1f3e9..59abb389 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. + + Registers its own tool rather than relying on the two tests above having + already run: pytest-xdist distributes tests across worker processes, so + the registry this assertion reads is not necessarily the one those tests + wrote to. + """ + + @tool_spec( + name="test_dummy_registered", + description="Registered.", + input_model=_AddInput, + ) + def _registered(spec: _AddInput) -> dict: + return {"result": spec.a + spec.b} + + assert "test_dummy_registered" in _TOOL_REGISTRY + assert _TOOL_REGISTRY["test_dummy_registered"] is _registered def test_rejects_non_basemodel_input_model(self) -> None: """input_model must be a pydantic BaseModel subclass.""" From 20e9b2dc130f86f9fbbaaaf4944af73260c080ff Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 3 Aug 2026 18:45:31 +0000 Subject: [PATCH 5/5] Drop unverified chapter pointers from the screening references The Cochran & Cox and Naes et al. citations carried specific chapter numbers I could not verify. A precise-looking wrong pointer is worse than no pointer, so they are replaced with what each source actually contributes - and the Naes entry now says plainly that the block construction is NOT in it, which is the thing a reader would otherwise go looking for there and not find. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TqdSBoXb1DpSobB76pMzFi --- src/process_improve/sensory/screening.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/process_improve/sensory/screening.py b/src/process_improve/sensory/screening.py index 9b3e3e45..fab5144d 100644 --- a/src/process_improve/sensory/screening.py +++ b/src/process_improve/sensory/screening.py @@ -44,10 +44,12 @@ 2(2), 149-168. Cochran, W. G. & Cox, G. M. (1957). *Experimental Designs* (2nd ed.). Wiley. -Chapters 9-11 on incomplete block designs. +The classical catalogue of balanced incomplete block designs. Naes, T., Brockhoff, P. B. & Tomic, O. (2010). *Statistics for Sensory and -Consumer Science*. Wiley. Chapter 2 on designing sensory experiments. +Consumer Science*. Wiley. Covers the assessor as a blocking factor and the +practical design of a panel session; the block *construction* below is not in +it, and comes from Cochran & Cox and from Williams. MacFie, H. J., Bratchell, N., Greenhoff, K. & Vallis, L. V. (1989). Designs to balance the effect of order of presentation and first-order carry-over effects