Skills › Scientific Agent Skills › experimental-design › doe_designs.py
doe_designs.py
227 lines · skills/experimental-design/scripts/doe_designs.py at 5f2bc03
1"""Design-of-experiments (DOE) matrices as labeled, decoded pandas DataFrames.
3pydoe returns designs in *coded* units (-1/+1, or 0..k-1). Researchers want the
4design in *real* factor units (temperature in C, concentration in mM) with named
5columns, randomized run order, and a clear sense of what each design is for. This
6module wraps pydoe 1.5.0 to do exactly that.
8A `factors` spec maps factor names to their real-world levels:
9 - two-level / continuous: {"temp": (20, 60), "conc": (1, 10)} # (low, high)
10 - multi-level categorical: {"catalyst": ["A", "B", "C"]}
12Functions:
13 full_factorial every combination of given levels (cost grows fast)
14 two_level_factorial 2^k full factorial (screening + interactions)
15 fractional_factorial 2^(k-p) fraction (screening many factors cheaply)
16 plackett_burman very economical main-effects-only screening
17 central_composite response-surface design (curvature / optimization)
18 box_behnken response-surface design, no extreme corners
19 latin_hypercube space-filling sample for simulation / computer experiments
21Each returns a DataFrame in real units; pass randomize=True to also get a
22randomized 'run_order' (default except for Latin hypercube).
23Requires: pydoe, numpy, pandas.
24"""
26from __future__ import annotations
28import operator
30import numpy as np
31import pandas as pd
34def _factor_names(factors):
35 if not factors or "run_order" in factors:
36 raise ValueError("provide at least one factor; 'run_order' is reserved")
37 return list(factors)
40def _ranges(factors):
41 """Reject malformed ranges rather than silently discarding a third level."""
42 _factor_names(factors)
43 out = {}
44 for name, levels in factors.items():
45 try:
46 values = np.asarray(levels, dtype=float)
47 valid = (values.shape == (2,) and np.isfinite(values).all()
48 and values[0] < values[1])
49 except (TypeError, ValueError, OverflowError):
50 valid = False
51 if not valid:
52 raise ValueError(f"{name!r} needs exactly two finite numeric levels, low < high")
53 out[name] = tuple(values)
54 return out
57def _decode_two_level(coded, factors):
58 """Map a -1/+1 coded matrix to real (low/high) units per factor."""
59 names = list(factors)
60 out = {}
61 for j, name in enumerate(names):
62 lvl = factors[name]
63 low, high = lvl[0], lvl[1]
64 mid, half = (high + low) / 2.0, (high - low) / 2.0
65 out[name] = mid + coded[:, j] * half
66 return pd.DataFrame(out)
69def _randomize(df, randomize, seed):
70 if not randomize:
71 return df.reset_index(drop=True)
72 rng = np.random.default_rng(seed)
73 order = rng.permutation(len(df)) + 1
74 df = df.copy()
75 df.insert(0, "run_order", order)
76 return df.sort_values("run_order").reset_index(drop=True)
79def full_factorial(factors, randomize=True, seed=0):
80 """Every combination of the listed levels.
82 factors values are explicit level lists, e.g.
83 {"temp": [20, 40, 60], "catalyst": ["A", "B"]} -> 3*2 = 6 runs.
84 Runs = product of level counts, so this explodes quickly with many factors.
85 """
86 from pydoe import fullfact
87 names = _factor_names(factors)
88 levels = [list(factors[n]) for n in names]
89 if any(not level or pd.Series(level, dtype=object).isna().any()
90 or not pd.Series(level, dtype=object).is_unique for level in levels):
91 raise ValueError("each factor needs nonempty, non-missing, distinct levels")
92 counts = [len(l) for l in levels]
93 coded = fullfact(counts).astype(int)
94 data = {n: [levels[j][coded[i, j]] for i in range(len(coded))]
95 for j, n in enumerate(names)}
96 return _randomize(pd.DataFrame(data), randomize, seed)
99def two_level_factorial(factors, randomize=True, seed=0):
100 """Full 2^k factorial: all main effects and all interactions, estimable.
102 Each factor needs a (low, high) pair. Use for k up to ~5; beyond that the
103 run count (2^k) gets expensive — switch to fractional_factorial or
104 plackett_burman for screening.
105 """
106 from pydoe import ff2n
107 factors = _ranges(factors)
108 coded = ff2n(len(factors))
109 return _randomize(_decode_two_level(coded, factors), randomize, seed)
112def fractional_factorial(factors, generator, randomize=True, seed=0):
113 """2^(k-p) fractional factorial from a generator string.
115 `generator` is pydoe's Yates notation, e.g. for 4 factors in 8 runs:
116 "a b c abc" (I=ABCD, resolution IV). Each token defines a column; multi-letter tokens
117 alias a factor with an interaction (this is the tradeoff — fewer runs, some
118 effects confounded). Choose a higher-resolution generator if you need to
119 separate main effects from two-factor interactions.
120 """
121 from pydoe import fracfact
122 factors = _ranges(factors)
123 coded = fracfact(generator)
124 if coded.shape[1] != len(factors):
125 raise ValueError(f"generator defines {coded.shape[1]} factors but "
126 f"{len(factors)} were named")
127 return _randomize(_decode_two_level(coded, factors), randomize, seed)
130def plackett_burman(factors, randomize=True, seed=0):
131 """Plackett-Burman screening design: main effects only, very few runs.
133 The target run count is the next multiple of 4 above k. pydoe 1.5.0 supports
134 counts of form 2^a, 12*2^a, or 20*2^a; other targets (e.g. 28 runs for k=25)
135 raise ValueError. Interactions can confound main effects; inspect aliasing.
136 """
137 from pydoe import pbdesign
138 factors = _ranges(factors)
139 coded = pbdesign(len(factors))
140 return _randomize(_decode_two_level(coded, factors), randomize, seed)
143def central_composite(factors, center=(0, 1), alpha="orthogonal",
144 face="circumscribed", randomize=True, seed=0):
145 """Central composite design (CCD) for response-surface / optimization work.
147 Adds axial ("star") points and center points to a 2^k factorial so you can fit
148 a quadratic model and locate an optimum. With face='circumscribed' the axial
149 points sit OUTSIDE the (low, high) box (so real levels exceed your stated
150 range); use face='inscribed' or 'faced' to keep everything within range.
151 `center` = (n center pts in factorial block, n in axial block). The default
152 supplies just one center run, so it does not estimate pure error. Add true
153 replicate runs when required. Output is globally shuffled: it does not
154 retain block membership for an experiment run in separate batches.
155 """
156 from pydoe import ccdesign
157 factors = _ranges(factors)
158 coded = ccdesign(len(factors), center=center, alpha=alpha, face=face)
159 return _randomize(_decode_two_level(coded, factors), randomize, seed)
162def box_behnken(factors, center=1, randomize=True, seed=0):
163 """Box-Behnken response-surface design (needs >= 3 factors).
165 Like a CCD it fits a quadratic, but it never uses the extreme corner
166 combinations (all-low or all-high), which is useful when those corners are
167 unsafe or infeasible. Compare actual run counts; it is not always smaller.
168 The default one center run does not supply pure-error replication.
169 """
170 from pydoe import bbdesign
171 factors = _ranges(factors)
172 if len(factors) < 3:
173 raise ValueError("box_behnken requires at least 3 factors")
174 coded = bbdesign(len(factors), center=center)
175 return _randomize(_decode_two_level(coded, factors), randomize, seed)
178def latin_hypercube(factors, n_samples, criterion="maximin", seed=0,
179 randomize=False):
180 """Space-filling Latin hypercube sample over continuous factor ranges.
182 For computer experiments / simulations where you want even coverage of a
183 high-dimensional space with relatively few points. Each factor needs a
184 (low, high) range. `criterion`: 'maximin' spreads points apart;
185 'center'/'centermaximin'/'correlation' are alternatives.
186 """
187 from pydoe import lhs
188 factors = _ranges(factors)
189 try:
190 n_samples = operator.index(n_samples)
191 except TypeError:
192 raise ValueError("n_samples must be a positive integer") from None
193 if n_samples < 1:
194 raise ValueError("n_samples must be a positive integer")
195 # pydoe draws from its own default_rng, so seeding numpy's global RNG has
196 # no effect -- the seed has to be handed to lhs itself.
197 names = list(factors)
198 unit = lhs( # in [0,1]
199 len(names), samples=n_samples, criterion=criterion, seed=seed
200 )
201 out = {}
202 for j, n in enumerate(names):
203 low, high = factors[n][0], factors[n][1]
204 out[n] = low + unit[:, j] * (high - low)
205 df = pd.DataFrame(out)
206 return _randomize(df, randomize, seed)
209if __name__ == "__main__":
210 f2 = {"temp": (20, 60), "conc": (1, 10), "ph": (6, 8)}
212 print("== 2^3 full factorial ==")
213 print(two_level_factorial(f2, seed=1).to_string(index=False))
215 print("\n== fractional 2^(4-1), generator 'a b c abc' ==")
216 f4 = {"A": (-1, 1), "B": (-1, 1), "C": (-1, 1), "D": (-1, 1)}
217 print(fractional_factorial(f4, "a b c abc", seed=1).to_string(index=False))
219 print("\n== Plackett-Burman screening, 5 factors ==")
220 f5 = {f"x{i}": (0, 1) for i in range(1, 6)}
221 print(f"runs = {len(plackett_burman(f5, randomize=False))}")
223 print("\n== central composite (2 factors) ==")
224 print(central_composite({"temp": (20, 60), "conc": (1, 10)}, seed=1).round(2).to_string(index=False))
226 print("\n== Latin hypercube, 8 samples over 3 factors ==")
227 print(latin_hypercube(f2, 8, seed=1).round(2).to_string(index=False))