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..2a93e6b5 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..fab5144d --- /dev/null +++ b/src/process_improve/sensory/screening.py @@ -0,0 +1,876 @@ +"""(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. +The classical catalogue of balanced incomplete block designs. + +Naes, T., Brockhoff, P. B. & Tomic, O. (2010). *Statistics for Sensory and +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 +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..68ae03d1 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 @@ -394,10 +397,179 @@ 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.""" + try: + # 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, + alpha=spec.alpha, + power=spec.power, + n_comparisons=spec.n_comparisons, + ) + elif spec.difference is not None and spec.n_per_product is None: + out = _required_panelists( + sd=spec.sd, + difference=spec.difference, + alpha=spec.alpha, + 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}) + + _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..c383bcb3 --- /dev/null +++ b/tests/test_sensory_screening.py @@ -0,0 +1,464 @@ +"""(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 + + +# --------------------------------------------------------------------------- +# 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."""