# DOE Pro+ — Algorithm Sources & License-Safe Implementation Spec

**Created:** 2026-06-27 · For the **DOE Pro+ (Tier 2)** features: no-code Bayesian/sequential optimization + design augmentation, in `doe-pro/app/index.html`.

## Why this document exists
DOE Pro is a **commercial** product. These features are non-trivial math, so they are **implemented clean-room from the canonical published algorithms** (mathematics is not copyrightable) and only **validated against permissively-licensed open-source** (BSD/MIT/Apache). No GPL/AGPL or closed-source (e.g. Atinary SDLabs) code or algorithms are used. Every source below was verified to exist (CrossRef / arXiv / GitHub) and its license checked.

---

## Feature 1 — Bayesian / sequential optimization

**Approach:** an exact small-n Gaussian-Process surrogate + analytic Expected Improvement, optimized over a finite candidate set (no gradients — ideal for a few coded factors). Default kernel = anisotropic **Matérn-5/2**; default acquisition = **Expected Improvement**. GP-UCB offered as an alternative (no erf, never underflows). All math in coded ([-1,1]) / standardized space; convert to real units only at the I/O boundary.

**Pipeline:** code factors to [-1,1] and z-score the response → build K → `A = K + (σ_n² + jitter)·I`, symmetrize → Cholesky `A = L Lᵀ` (escalate jitter ×10, ≤5 retries) → `α = Lᵀ \ (L \ y)` → per candidate `x*`: `mean = k*·α`, `v = L \ k*`, `var = k(x*,x*) − v·v` (clamp ≥0) → de-standardize → score by EI (rank via log-EI for stability) → argmax = next run. Batch suggestions via **Constant-Liar** (lie = min(y)).

**Exact formulas**
- Matérn-5/2 (R&W eq. 4.17): `k = σ_f²·(1 + √5·r + 5r²/3)·exp(−√5·r)`, `r = √(Σ_d (x_d−x'_d)²/l_d²)`
- RBF alternative (R&W eq. 4.9): `k = σ_f²·exp(−½ Σ_d (x_d−x'_d)²/l_d²)`
- GP solve (R&W **Algorithm 2.1**): never invert K; use Cholesky + triangular solves; `var = k(x*,x*) − vᵀv`
- Log marginal likelihood (R&W eq. 2.30): `−½ yᵀα − Σ log L_ii − (n/2) log 2π`
- EI, maximization (Jones 1998 / Brochu 2010): `d = μ − f_best − ξ; z = d/σ; EI = d·Φ(z) + σ·φ(z)` (ξ≈0.01; if σ≤1e-9 → EI=0)
- Stable log-EI ranking (Ament 2023): `EI = σ·h(z)`, `h(z)=z·Φ(z)+φ(z)`, `log h(z) ≈ −z²/2 − log|z| − ½log2π` for z≪0
- GP-UCB (Srinivas 2010): `μ + κ·σ` (κ≈2)
- erf (Abramowitz & Stegun **7.1.26**, max abs err 1.5e-7): `t=1/(1+0.3275911x)`, `erf = 1−(((((1.061405429t−1.453152027)t+1.421413741)t−0.284496736)t+0.254829592)t)·e^{−x²}`; `Φ(z)=½(1+erf(z/√2))`, `φ(z)=e^{−z²/2}/√(2π)`
- Constant-Liar batch (Ginsbourger 2010): pick argmax EI → append fake (x, L=min(y)) → refit → repeat to q → discard lies

**Critical pitfalls (from the research):** never invert K (Cholesky only); clamp `var≥0` before sqrt; branch `σ≤1e-9 → 0` before dividing (the #1 NaN bug); rank by log-EI once the GP is confident (raw EI underflows to 0 everywhere); pick ONE sign convention (internal maximization; minimize by optimizing −f); standardize y / code factors so ξ, κ, length-scales transfer and K stays conditioned; with <~5-8 points bound/fix length-scales (l_i ∈ [0.05, 2] coded) or the GP overfits; symmetrize A before Cholesky; log-det from `Σ log L_ii`, never a determinant call.

**Cite in code + sources:**
- Rasmussen & Williams (2006), *Gaussian Processes for Machine Learning*, MIT Press. ISBN 978-0-262-18253-9. Free full text: gaussianprocess.org/gpml (Alg. 2.1; eqs 2.23–2.26, 2.30, 4.9, 4.17). *Freely readable, not open-licensed — cite, don't redistribute.*
- Jones, Schonlau & Welch (1998), *Efficient Global Optimization of Expensive Black-Box Functions*, J. Global Optim. 13(4):455-492. doi:10.1023/A:1008306431147
- Snoek, Larochelle & Adams (2012), *Practical Bayesian Optimization of Machine Learning Algorithms*, NeurIPS 25. arXiv:1206.2944 (Matérn-5/2 default)
- Brochu, Cora & de Freitas (2010), *A Tutorial on Bayesian Optimization…*, arXiv:1012.2599 (EI/PI/UCB, ξ/κ, maximization)
- Srinivas, Krause, Kakade & Seeger (2010), *Gaussian Process Optimization in the Bandit Setting (GP-UCB)*, ICML 2010. arXiv:0912.3995
- Ginsbourger, Le Riche & Carraro (2010), *Kriging Is Well-Suited to Parallelize Optimization*, Springer. doi:10.1007/978-3-642-10701-6_6 (Constant Liar / Kriging Believer)
- Ament et al. (2023), *Unexpected Improvements to Expected Improvement for Bayesian Optimization*, NeurIPS 2023. arXiv:2310.20708 (LogEI)
- Abramowitz & Stegun (1964/1972), *Handbook of Mathematical Functions*, NBS AMS 55, eq. 7.1.26 p.299 (public domain)

**Validate behavior against (permissive OSS, reference only):** scikit-learn `GaussianProcessRegressor` (BSD-3) for GP mean/std; scikit-optimize `skopt/acquisition.py` (BSD-3) and BoTorch `analytic.py` (MIT, incl. LogEI) for EI/PI/UCB; scikit-optimize `cl_min` for Constant-Liar.

---

## Feature 2 — Design augmentation

Two paths, both reusing the app's existing engine.

**Path 1 (default, pure JS) — screening → RSM growth.** Take the existing two-level factorial/screen ("cube" block), append the **axial/star block** (2k points) + center points to form a **Central Composite Design**, reusing all original runs. `α = (n_f)^{1/4}` for rotatability, or `α = 1` (face-centered CCF) when levels beyond ±1 are infeasible (common in bioprocess — flag if rotatable α pushes runs outside the validated window). Optional **foldover** (reverse all coded signs, or one column) to de-alias a low-resolution screen first.

**Path 2 (general) — D-optimal augmentation, existing runs FIXED.** Adapt the app's existing `generateDOptimal()` coordinate-exchange (it already does random-restart cyclic exchange over levels, incremental X'X via buildXtX/addOuter, det via matDet, model via getModelTerms/termValue). The only changes: (a) seed X with the fixed existing runs' model rows first, then the new random rows; (b) restrict the coordinate-exchange row loop to indices ≥ n_fixed (never touch existing runs); (c) maximize `det(X_combᵀX_comb)`; (d) report combined D-efficiency. Raise restarts from 4 to 10–50 for augmentation.

**Pitfalls:** optimize **log-determinant** (2·Σ log diag of Cholesky), not raw det (over/underflow); `cond(X'X)=cond(X)²` so watch conditioning; guard singularity (`n_fixed+n_new ≥ p`, no duplicate support); the augmentation-specific bug is iterating over ALL rows (re-optimizes existing runs — defeats augmentation); Fedorov delta must include the cross term `d_ij`; many restarts (local optima); relative improvement tolerance (avoid float cycling); scale to coded [-1,1] before building X; flag rotatable α that extrapolates beyond the screened range.

**Cite in code + sources:**
- Meyer & Nachtsheim (1995), *The Coordinate-Exchange Algorithm for Constructing Exact Optimal Experimental Designs*, Technometrics 37(1):60-69. doi:10.1080/00401706.1995.10485889
- Montgomery (2013), *Design and Analysis of Experiments*, 8th ed., Wiley, Ch. 11 (CCD; α=(n_f)^{1/4}; orthogonal blocking). ISBN 978-1-118-14692-7
- Box, Hunter & Hunter (2005), *Statistics for Experimenters*, 2nd ed., Wiley (foldover; sequential RSM assembly). ISBN 978-0-471-71813-0
- Matrix determinant lemma / Sherman-Morrison (standard linear algebra)

**Validate against (permissive OSS, reference only):** pyDOE2 `ccdesign`/`fracfact` (BSD-3) for Path 1; dexpy `optimal.py` `delta()`/`update()` (Apache-2.0) for Path 2 — note stock dexpy bootstraps fresh and does NOT hold runs fixed, so the augmentation change is ours. The app's existing `generateDOptimal()` is first-party.

---

## License summary
**Safe to build.** Implement clean-room from the papers/textbooks above (math not copyrightable). Reference only permissive OSS (BSD-3: scikit-learn, scikit-optimize, pyDOE2; MIT: BoTorch, bayesian-optimization, ml-matrix, jStat; Apache-2.0: BayBE, dexpy). If any small snippet is ported, retain the upstream copyright/LICENSE (and NOTICE for Apache-2.0). Do **not** reference Atinary's closed product or proprietary algorithms. Augmentation extends the project's own existing `generateDOptimal()`.