diff --git a/CHANGELOG.md b/CHANGELOG.md index f81449c..02db7a4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,8 +2,7 @@ Format follows [Keep a Changelog](https://keepachangelog.com/en/1.0.0/); versioning follows [Semantic Versioning](https://semver.org/). `v0.1.0` is the first release; its -entry below summarizes everything since the package's extraction from the paper's -research code (`entropy_driven_hybrid_models_code`) began. See `HANDOFF.md` for the full phase-by-phase narrative and +entry below summarizes the library relative to the paper's original research code. See `docs/improvements_over_paper.md` for a detailed comparison against the original code. ## [Unreleased] @@ -42,56 +41,52 @@ The first release. name, kept for backward compatibility), `AnisotropicSE`, `FixedInverseMean`, `InputTransform`/`OutputTransform`, `latin_hypercube_design`/`full_factorial_design`, `second_order_entropy`/`entropy_lower_bound`/`gaussian_mixture_density`. -- N-D generalization (Phase 5): the acquisition path a run actually depends on - (`optimize`) works at any input dimension; proven at 1-D and 3-D, not just the - paper's 2-D case. -- A selectable acquisition objective (Phase 9d): `BitsForGaps.acquisitionObjective` - (or `objective=` on `acquisition.entropy_objective`/`optimize`) chooses between the +- N-D generalization: the acquisition path a run actually depends on (`optimize`) + works at any input dimension; proven at 1-D and 3-D, not just the paper's 2-D case. +- A selectable acquisition objective: `BitsForGaps.acquisitionObjective` (or + `objective=` on `acquisition.entropy_objective`/`optimize`) chooses between the paper's 2nd-order Taylor entropy approximation (`"taylor"`, default) and its - closed-form lower bound (`"lower_bound"`, implemented since Phase 2 but never - wired up as a usable choice before now). -- Public-API input validation (Phase 9c): clear `ValueError`s for bounds/kernel - dimensionality mismatches, invalid bounds, non-positive HMC/acquisition config, - `X_init`/`y_init` shape mismatches, and a malformed injected black-box output. -- An opt-in `tf_seed` on `mixture.sample_gp_posterior_mixture`/`predict_grid_2D` - (Phase 9c) to make `predict_f_samples`' otherwise TF-ambient-RNG-driven draws - reproducible on request. -- `py.typed` (PEP 561) and type hints across `src/bits_for_gaps/` (Phase 9d). -- Archive-free figure reproduction (Phase 9): a curated ~16 MB subset of the - published run's plot-input data, committed to `paper/data/`, so - `python paper/reproduce.py` regenerates every figure from a fresh clone with no - private-archive access. -- `examples/vle_distillation/` (Phase 6): the paper's H2O-PrOH VLE/distillation case - study, ported onto the public API (Julia/Clapeyron-backed, repo-only). -- `examples/synthetic/run_example.py` (Phase 9d): a small, Julia-free, runnable - example for onboarding. -- Sphinx/MyST documentation (Phase 8), a CI workflow running the full default test - suite plus `ruff check` (Phase 9d), and this changelog. -- `docs/theory.md` rewritten to present the paper's key equations in its own notation - with equation numbers, each linked to the implementing module (Phase 9e). -- Release engineering (Phase 10): `.github/workflows/publish.yml` publishes to PyPI - (on a `v*` tag) or TestPyPI (manual dry run) via OIDC trusted publishing -- no API - tokens in the repo; see `RELEASE.md` for the maintainer checklist. + closed-form lower bound (`"lower_bound"`), now usable as an acquisition choice. +- Public-API input validation: clear `ValueError`s for bounds/kernel dimensionality + mismatches, invalid bounds, non-positive HMC/acquisition config, `X_init`/`y_init` + shape mismatches, and a malformed injected black-box output. +- An opt-in `tf_seed` on `mixture.sample_gp_posterior_mixture`/`predict_grid_2D` to + make `predict_f_samples`' otherwise TF-ambient-RNG-driven draws reproducible on + request. +- `py.typed` (PEP 561) and type hints across `src/bits_for_gaps/`. +- Archive-free figure reproduction: a curated ~16 MB subset of the published run's + plot-input data, committed to `paper/data/`, so `python paper/reproduce.py` + regenerates every figure from a fresh clone with no private-archive access. +- `examples/vle_distillation/`: the paper's H2O-PrOH VLE/distillation case study, + built on the public API (Julia/Clapeyron-backed, repo-only). +- `examples/synthetic/run_example.py`: a small, Julia-free, runnable example for + onboarding. +- Sphinx/MyST documentation, a CI workflow running the full default test suite plus + `ruff check`, and this changelog. +- `docs/theory.md` presents the paper's key equations in its own notation with + equation numbers, each linked to the implementing module. +- Release engineering: `.github/workflows/publish.yml` publishes to PyPI (on a `v*` + tag) or TestPyPI (manual dry run) via OIDC trusted publishing -- no API tokens in + the repo; see `RELEASE.md` for the maintainer checklist. ### Changed -- Decomposed the paper's monolithic `driver_new.py` into focused, independently - testable modules (Phase 4): `gp`, `mixture`, `acquisition`, `entropy`, `transforms`, - `state`, with `sampler.py`'s `adaptiveEntropy` reduced to a thin orchestrator. -- Retired disk-as-state (Phase 4): `run()` takes the initial design in memory and - returns a `RunHistory`; a full run executes with zero disk writes by default, - per-iteration file output is available but opt-in (`checkpoint_dir`). +- Decomposed the sequential-design engine into focused, independently testable + modules: `gp`, `mixture`, `acquisition`, `entropy`, `transforms`, `state`, with + `sampler.py`'s `adaptiveEntropy` reduced to a thin orchestrator. +- Retired disk-as-state: `run()` takes the initial design in memory and returns a + `RunHistory`; a full run executes with zero disk writes by default, per-iteration + file output is available but opt-in (`checkpoint_dir`). - `mixture.sample_gp_posterior_mixture`/`acquisition.entropy_objective` now save and - restore a GP's kernel hyperparameters around their internal reassignment loop - (Phase 9c) -- behavior-preserving for every value either function returns, but the - caller's model is no longer left at an arbitrary leftover state afterward. + restore a GP's kernel hyperparameters around their internal reassignment loop, so + the caller's model is not left at an arbitrary leftover state afterward. - `examples/vle_distillation/distillation.py`'s `solve_column` retries a few generic alternate `fsolve` initial guesses if the default doesn't converge, before giving - up (Phase 9c) -- the primary attempt is unchanged, so this only ever activates when - it was already failing. -- Whole-repo `ruff` lint + format pass (Phase 9d), preserving the load-bearing import - orders (`__init__.py`'s lazy-import split; `PYTHON_JULIACALL_HANDLE_SIGNALS` set - before the `juliacall` import it protects). + up -- the primary attempt is unchanged, so this only ever activates when it was + already failing. +- Whole-repo `ruff` lint + format pass, preserving the load-bearing import orders + (`__init__.py`'s lazy-import split; `PYTHON_JULIACALL_HANDLE_SIGNALS` set before + the `juliacall` import it protects). ### Fixed @@ -102,30 +97,24 @@ The first release. reused a `GPmodel` for a second purpose after an earlier step had mutated its kernel and left it at an arbitrary state, producing a non-converging McCabe-Thiele column that a first write-up incorrectly attributed to a property of entropy-driven - acquisition design (Phase 9b; see `paper/PHASE9B_INVESTIGATION.md`). The underlying - footgun was then hardened away at its source (Phase 9c, above). + acquisition design. The underlying footgun was then hardened away at its source + (see "Changed", above). - `entropy.py`'s density-positivity check promoted from a bare `assert` (silently - stripped under `python -O`) to an explicit `ValueError` (Phase 9c). + stripped under `python -O`) to an explicit `ValueError`. - `kernels.assign_hyperparameters` now raises a clear, specific error (naming the parameter and value) instead of a low-level gpflow/TensorFlow traceback when a value can't round-trip through a parameter's transform -- e.g. an extreme outlier - posterior sample for `lengthscale_2`, deliberately left unconstrained (Phase 9d). -- A stale doc claim: `docs/reproduce_paper.md` and `docs/theory.md` described - pre-Phase-9 behavior (archive access required; the entropy lower bound "not wired - into the default acquisition path") after both had changed. - -### Test / process - -- A pre-refactor numerical baseline pin (`tests/integration/data/synthetic_baseline.json`, - atol 1e-10) and reference scalars extracted from the published run (`paper/reference/*`) - that every subsequent phase must reproduce exactly -- the original code had no - automated tests at all. -- Coverage raised from 93% to 97% (Phase 9d), closing real gaps (`gp.py` had zero - direct unit tests before). + posterior sample for `lengthscale_2`, deliberately left unconstrained. + +### Testing + +- A numerical baseline pin (`tests/integration/data/synthetic_baseline.json`, + atol 1e-10) and reference scalars extracted from the published run + (`paper/reference/*`) that the library must reproduce exactly -- the original code + had no automated tests at all. +- Coverage raised to 97%, closing real gaps (`gp.py` and `diagnostics.py` had zero + direct unit tests before; several documented `ValueError` paths were never + exercised). - A Monte-Carlo validation of the entropy approximations against the true GMM - differential entropy they approximate, not just a captured historical value - (Phase 9d). -- 204 tests passing (up from 166 at Phase 9d's start), closing coverage gaps found in - a function-by-function audit (`diagnostics.py` had zero direct unit tests; several - documented `ValueError` paths were never exercised) and upgrading a few shape-only - assertions to check actual values (Phase 9e). + differential entropy they approximate, not just a captured historical value. +- 204 tests passing. diff --git a/HANDOFF.md b/HANDOFF.md deleted file mode 100644 index 3b2ccba..0000000 --- a/HANDOFF.md +++ /dev/null @@ -1,1063 +0,0 @@ -# HANDOFF — bits_for_gaps - -State of the fresh `bits_for_gaps` repo. Read this + `REFACTOR_PLAN.md` before continuing. - -- Bootstrap (Phase 0 + Phase 3-lite): Opus, 2026-07-03. -- Phase 2 (regression/test harness): done 2026-07-03, merged to `main`. -- Phase 4 (decompose `sampler.py`; retire disk-as-state): done 2026-07-03, merged to `main`. -- Phase 5 (generalize to N-D): done 2026-07-04, merged to `main`. -- Phase 6 (port the VLE/distillation example): done 2026-07-04, merged to `main`. -- Phase 7 (reproduce the paper's figures): done 2026-07-04, merged to `main`. -- Phase 8 (Sphinx + MyST + ReadTheDocs docs): done 2026-07-04, merged to `main`. -- Phase 9 (reproduce ALL paper results, including the from-scratch stochastic - HMC+acquisition loop; archive-free figure reproduction): done 2026-07-04, merged to - `main`. -- Phase 9b (investigate + fix the McCabe-Thiele non-convergence Phase 9 flagged as a - discrepancy): done 2026-07-04, merged to `main`. -- Phase 9c (robustness hardening -- first sanctioned `src/bits_for_gaps/` core change - since Phase 4): done 2026-07-04, merged to `main`. -- Phase 9d (whole-codebase polish: hygiene [ruff, type hints, CI, coverage] + - faithfulness [selectable acquisition, MC-validation, synthetic example, hardening, - CHANGELOG]): done 2026-07-04, merged to `main`. -- Phase 9e (docs/tests/comments quality pass -- behavior-preserving: `golden` -> - `reference` rename, `theory.md` equation fidelity, test-coverage gap-filling, - NumPy-docstring + inline equation-citation pass): done 2026-07-05, merged to `main`. -- **Phase 10 (release engineering for the v0.1.0 PyPI release -- packaging/CI only, no - algorithm changes, no credentialed step performed): done 2026-08-06 on branch - `phase10-release` — awaiting review/merge to `main`, then the maintainer-only steps - in `RELEASE.md`.** - -## Phase 10 — release engineering for v0.1.0 (done; review gate before the maintainer's release steps) - -Packaging/CI/docs only -- `src/bits_for_gaps` logic and every regression value -untouched. `pytest -q`: 204 passed, 2 deselected throughout (unchanged from Phase 9e). - -- **Version, single-sourced.** `pyproject.toml`'s `[project]` gained - `dynamic = ["version"]` + a `[tool.hatch.version]` pointing at - `src/bits_for_gaps/__init__.py` (hatchling's default "regex" version source matches - `__version__ = "..."`), so `__version__` stays the one place edited to cut a - release. Bumped `0.0.1.dev0` -> `0.1.0`. -- **PyPI metadata polish.** Classifiers: `Development Status :: 3 - Alpha` -> - `4 - Beta`, plus a generic `Python :: 3` classifier and an OS-independent / - AI-topic classifier. `project.urls` gained `Repository`/`Documentation`/`Changelog` - alongside the existing `Homepage`/`Paper`. -- **Sdist scope, fixed.** `python -m build`'s sdist defaulted to hatchling's - whole-repo file set (193 entries, 6.9 MB, including all of `paper/data/`'s curated - ~16 MB plot-input subset, `tests/`, `examples/`, `docs/`) -- a real packaging gap, - not caught by the wheel (which was already clean: - `bits_for_gaps/*` + `py.typed` + dist-info only). Fixed via - `[tool.hatch.build.targets.sdist]`'s `only-include` (not `include`, which *adds* to - the default set -- tried first, still leaked several READMEs), restricting the - sdist to `src/bits_for_gaps` + `LICENSE`/`README.md`/`CHANGELOG.md` (6.9 MB -> 33 KB). - `twine check dist/*` passes for both artifacts; wheel `METADATA`'s `Version:` field - matches the single-sourced `__version__`. -- **Clean-env install audit.** In a fresh, throwaway conda env (no dev tooling, no - editable install), `pip install dist/bits_for_gaps-0.1.0-py3-none-any.whl` resolved - every pinned dependency cleanly from wheel metadata alone; bare `import - bits_for_gaps` stayed Julia-free/TF-lazy; a smoke test (`AnisotropicSE()`, - `latin_hypercube_design`, `entropy.second_order_entropy`, and `docs/quickstart.md`'s - exact `BitsForGaps(...)` construction snippet) matched the installed API with no - doc changes needed. Env torn down after. -- **Trusted-publishing CI.** New `.github/workflows/publish.yml`: a shared build job - (the same `python -m build` + `twine check` verified above) feeds two gated publish - jobs via OIDC trusted publishing (`pypa/gh-action-pypi-publish`, a dedicated - `release` environment, job-scoped `id-token: write`) -- no API tokens anywhere in - the repo. A `v*` tag push reaches PyPI; a manual `workflow_dispatch` reaches - TestPyPI only, for the pre-release dry run. `gh-action-pypi-publish` is referenced - via PyPA's own recommended floating tag (`release/v1`) rather than a fabricated - commit SHA (this session had no network access to verify one) -- `RELEASE.md` - documents the optional SHA-pinning step for the maintainer. -- **CHANGELOG finalized, README badges added.** `[Unreleased]` renamed to - `[0.1.0] - YYYY-MM-DD` (a placeholder; the maintainer fills in the actual date at - tag time) with a fresh empty `[Unreleased]` left above it; added the Phase - 9e/10 bullets that weren't yet reflected. `README.md` gained PyPI-version and - ReadTheDocs badges alongside the existing CI badge (both inert until the - maintainer's release/RTD-import steps below actually happen). -- **`RELEASE.md`** — new maintainer checklist: the exact build-artifact contents - audit, the clean-env install audit, and the four maintainer-only steps in order - (register trusted publishers on PyPI + TestPyPI; TestPyPI dry run via - `workflow_dispatch`; set the CHANGELOG date, `git tag v0.1.0 && git push --tags`; - activate ReadTheDocs). None of these four were performed by this phase. - -Verified at every commit: `pytest -q` (204 passed, 2 deselected), `pytest -m vle` (2 -passed), `ruff check` clean, `sphinx-build -W` clean, and the Julia-free/TF-lazy -import contract. No PyPI/TestPyPI upload, no account configuration, no trusted- -publisher registration, no `git tag`, no RTD activation -- all maintainer-only, all -documented in `RELEASE.md`, none performed. - -## Phase 9e — docs/tests/comments quality pass (done; review gate before Phase 10) - -Three passes, each its own commit, suite green at every one. Behavior-preserving -throughout (docstrings/comments/tests/docs + one mechanical rename -- no algorithm -changes): baseline (atol 1e-10) + all reference regressions + `pytest -m vle` never -moved; `paper/data/` and the dependency stack untouched. - -- **PASS 1 -- rename + docs.** `paper/golden/` -> `paper/reference/` everywhere - (`git mv`, the `golden` pytest fixture -> `reference`, `extract_golden.py` -> - `extract_reference.py`, every prose/docstring mention across ~30 files) -- - confirmed the 4 JSON files are byte-identical (only the directory moved) and zero - "golden"/"Golden" text remains anywhere in the tracked repo. `docs/theory.md` - rewritten to present the paper's key equations in its own notation with equation - numbers -- Eq (1)/(2) entropy+acquisition, Eq (3) hyperparameter posterior, Eq (4) - GP prior, Eq (5a)/(5b) predictive mean/variance, Eq (6) SE kernel + Table 1's exact - priors, Eq (7) GMM predictive + Eq (8a)/(8b) moments, Eq (9) + the Huber et al. - (2008) Taylor expansion and its Proposition truncation bound (SI-1), the closed-form - entropy-lower-bound Theorem and its SI-2 cross-overlap term, Algorithm 1 (credible - intervals, Lalchand & Rasmussen 2020), and Eq (10)/(11) + SI-4 for the VLE example -- - each linked to the implementing module. Every doc page's cross-links verified to - resolve in the rendered HTML (not just `sphinx-build -W`, which doesn't catch - unresolved `:func:`/`:class:` refs with `nitpicky=False`); fixed a stale test count - in `installation.md` (117/5 -> 193/2). Docs page organization was already sound -- - no reordering needed. -- **PASS 2 -- test-coverage audit.** Added `tests/unit/test_diagnostics.py` (new -- - `potential_scale_reduction`/`effective_sample_size` had zero direct unit tests - before, only indirect coverage via full HMC integration runs). Added tests - exercising `entropy_surface_2D`/`predict_grid_2D`'s documented non-2-D `ValueError` - (defined since Phase 5, never exercised by any test); direct tests for - `entropy.gaussian_mixture_density` against `scipy.stats`; strengthened - `InputTransform`'s 1-D-input test and `sampler.call_model`'s success test to check - actual values, not just shapes; added a test confirming `call_model` extracts only - the first element of a multi-element black-box output. One finding, not a bug: - `design.full_factorial_design`'s "grid too small" `ValueError` is unreachable dead - code for every `(bounds, n_train, n_test)` -- `levels = ceil(n_total ** (1/d))` - mathematically guarantees `levels**d >= n_total` (verified by brute-force search); - documented in a comment rather than given a fake test. `pytest -q`: 193 -> 204 - passed (2 deselected, unchanged). -- **PASS 3 -- docstrings + inline equation citations.** Filled NumPy-docstring gaps - (missing Parameters/Returns/Notes: `entropy.py`'s `first_order_entropy_approx`/ - `cholesky`/`gradient_gaussian_mixture_density`/`second_order_entropy`, `gp.py`'s - `build_gp`/`maximize_lml`/`run_mcmc`, `means.py`'s `FixedInverseMean`, `_util.py`'s - three array helpers) and added inline comments citing the paper's equation numbers - at the exact lines implementing them (verified against the paper + SI): Eq (6) in - `kernels.AnisotropicSE.K`, Eq (3)/(4)/(5a)/(5b) in `gp.py`, Eq (7)/(9)/Theorem-SI-2 - in `entropy.py`, Eq (1)/(2)/(5a)/(5b)/(7) in `acquisition.py`/`mixture.py`, and - Eq (10)/(11)/SI-4 in the VLE example's `phase_diagram.py`/`gibbs_duhem.py`/ - `distillation.py`. - -## Phase 9d — whole-codebase polish: hygiene + faithfulness (done; merged to main) - -Two batches, eight workstreams (A-D hygiene, E-H faithfulness/features), each its own -commit, suite green at every one. Behavior-preserving throughout: baseline (atol -1e-10) + all reference regressions + `pytest -m vle` never moved. - -**Hygiene (A-D):** - -- **A -- ruff lint + format, whole repo.** Added to the `[dev]` extra; configured - (`E`/`F`/`W`/`I`/`B`, line-length 100) with per-file-ignores protecting the two - load-bearing import orders (`__init__.py`'s eager-vs-lazy split; - `PYTHON_JULIACALL_HANDLE_SIGNALS` set before the `juliacall`/TF-threading-config - imports it protects, in `activity_model.py` and `paper/full_reproduction.py`). A - handful of non-autofixable findings (assigned lambdas, one unused import, two - intentionally-blind `pytest.raises(Exception)`) fixed by hand; the rest (52 files) - is `ruff format`'s cosmetic whitespace/quote normalization. -- **B -- type hints across `src/bits_for_gaps/`, `py.typed` (PEP 561).** Every module - annotated via `from __future__ import annotations` (stringized -- no new runtime - imports, so the Julia-free/TF-lazy contract holds). `py.typed` confirmed included - in the built wheel. Fixed a real `sphinx-build` regression the new hints exposed: - `InputTransform`/`OutputTransform` are documented twice (top-level re-export + - original location), and once a type hint referenced them by bare name, - `autodoc-typehints` couldn't disambiguate ("more than one target found") -- fixed - via `docs/api.rst`'s `:exclude-members:`. -- **C -- broadened CI + badge.** `.github/workflows/ci.yml` now runs `ruff check` - plus the full default `pytest -q` (unit + integration + regression; still excludes - `@pytest.mark.vle`), not just `tests/unit`. Verified in a fresh conda env with only - the `[dev]` extra (no Julia) that this is genuinely Julia-free. -- **D -- coverage pass, 93% -> 97%.** `gp.py` had zero direct unit tests (only - indirect coverage via full HMC integration runs); new `tests/unit/test_gp.py` - closes it to 98%. A few real `sampler.py` gaps closed too (`read_data`, - `BitsForGaps`'s custom transform override, `run(predict_grid=True)`, - `run(initalLML=True)`). No coverage gate added to CI (informational only). - -**Faithfulness/features (E-H):** - -- **E -- `entropy_lower_bound` wired up as a selectable acquisition objective.** The - paper derives two entropy estimators; the closed-form lower bound was implemented - and unit-tested since Phase 2 but never usable for acquisition. - `acquisition.entropy_objective`/`optimize`/`entropy_surface_2D` gain - `objective="taylor"|"lower_bound"`; `BitsForGaps.acquisitionObjective` (default - `"taylor"` -- unchanged existing behavior). -- **F -- Monte-Carlo validation of the entropy approximations.** New - `tests/unit/test_entropy_mc_validation.py` estimates the true GMM differential - entropy directly (sample the mixture, average -log of its own density at those - samples) on several mixtures and checks the Taylor approximation stays close - (rtol=0.15, calibrated empirically) and the lower bound stays at or below it -- - validates approximation *quality*, not just a regression pin. -- **G -- a pure-Python synthetic example.** `examples/synthetic/run_example.py`: a - small, actually-runnable, Julia-free script (`python examples/synthetic/run_example.py`, - well under a minute) -- the onboarding path that previously only pointed at test - files. `docs/quickstart.md` now points to it. -- **H -- clear diagnostic for an unassignable hyperparameter value.** - `kernels.assign_hyperparameters` (called deep in `mixture.py`/`acquisition.py`'s hot - loops) now catches gpflow's low-level `InvalidArgumentError` for a value that can't - round-trip through a parameter's transform (e.g. an extreme `lengthscale_2` outlier - -- deliberately unconstrained, no positivity bijector) and re-raises a `ValueError` - naming the parameter and value. Also added `CHANGELOG.md` (Keep a Changelog, - Unreleased section summarizing every phase, targeting `v0.1.0`). - -`pytest -q`: 166 -> 193 passed (2 deselected throughout). `paper/data/`, `paper/reference/*`, -and the pinned dependency stack untouched. `import bits_for_gaps` confirmed Julia-free -and TensorFlow-lazy after every workstream. Docs updated: -`docs/improvements_over_paper.md` (new sections for E/F/G/H), `docs/theory.md` and -`docs/reproduce_paper.md` (both had gone stale describing pre-Phase-9/9d behavior -- -fixed in passing), `docs/quickstart.md`, `docs/api.rst`. - -## Phase 9c — robustness hardening (done; merged to main) - -Behavior-preserving only, by design: no numerical result, default, or seed changed. -The pre-Phase-4 baseline (atol 1e-10) + all `paper/reference/*` regressions + -`pytest -m vle` stayed green at every commit. 46 new tests added (120 -> 166 passed, 2 -deselected). - -**Mutation-footgun fix, broader than Phase 9b realized.** Phase 9b's bug (a script -reusing a `GPmodel` after `mixture.sample_gp_posterior_mixture` mutated its kernel and -left it at an arbitrary leftover state) turned out to be a symptom of a real footgun -in the library, not a one-off script mistake: `sampler.py`'s own `run()` calls -`entropy_surface_2D`/`optimize` (both funnel through `acquisition.entropy_objective`, -which reassigns kernel hyperparameters in the same way) on the *same* `GPmodel` object -it then stores in `IterationRecord.GPmodel` and optionally checkpoints. **Every run's -returned/checkpointed model** used to carry this arbitrary state, not just the one -Phase 9b patched. Both `sample_gp_posterior_mixture` and `entropy_objective` now -save the kernel's hyperparameters before mutating and restore them in a `finally` -- -`kernels.py` gained `save_hyperparameters()` to pair with the existing -`assign_hyperparameters()`. Verified behavior-preserving for every value either -function computes/returns; new unit tests assert kernel state is unchanged -before/after (including on the error path), plus an integration test reproducing the -exact Phase 9b scenario end-to-end. - -**Public-API input validation** (`adaptiveEntropy`/`BitsForGaps`): clear `ValueError`s -for `x_bounds` `lo >= hi`, bounds-vs-kernel-`ndim` mismatch (checked in `__init__`), -non-positive/out-of-range HMC/acquisition config and `X_init`/`y_init` shape -mismatches (checked in `run()`, since config is conventionally set via -post-construction attribute assignment throughout this codebase, not passed to -`__init__`), and a black-box output that isn't a non-empty sequence (`call_model`) -- -previously cryptic failures deep inside GPflow/TensorFlow or a bare `IndexError`. - -**SHOULD-fixes, confirmed valuable by the audit:** `entropy.py`'s `assert pl > 0` -(silently stripped under `python -O`) is now an explicit `ValueError`. -`examples/vle_distillation/distillation.py`'s `solve_column` retries a few generic -alternate `fsolve` initial guesses if the default doesn't converge -- the primary -attempt is byte-for-byte unchanged, so this only ever activates when the default -already reports non-convergence. `mixture.sample_gp_posterior_mixture`/ -`predict_grid_2D` gained an optional `tf_seed` to make `predict_f_samples` draws -reproducible on request (default `None` leaves the documented ambient-RNG behavior -unchanged). - -**Verification:** re-ran the full 15-iteration stochastic loop from scratch -(`results_remaked/phase9c_fullrun/`, gitignored; hit a one-off ~15 min TF -`tf.function`-retracing anomaly at iteration 2, unrelated to this phase's changes -- -other iterations ran at the normal ~13 s/iteration pace). `column_surrogate_converged` -stayed `True`; R-hat/ESS/hyperparameter posterior matched Phase 9b's committed -`full_run_summary.json` to 5-6 significant figures (within the already-documented -run-to-run floating-point tolerance) -- confirming the hardening changed no numerical -behavior. No update to `paper/phase9_validation/` needed (values didn't move outside -that tolerance). - -**Docs:** new `docs/improvements_over_paper.md` (wired into `docs/index.md`) -consolidates bug fixes (the `equilibrium.py` missing-path repoint, Phase 9b's mutation -bug), this phase's hardening, and the architecture wins from every prior phase. -`docs/reproduce_paper.md` was also fixed in passing -- found stale from before Phase 9 -(still described the pre-archive-free, author-access-required flow). -`sphinx-build -W docs docs/_build/html` succeeds with zero warnings. - -`paper/data/`, `paper/reference/*`, and the dependency stack were not touched. -`import bits_for_gaps` confirmed still Julia-free. - -## Phase 9b — root-cause the McCabe-Thiele discrepancy (done; merged to main) - -Phase 9's "genuine discrepancy" (the fully-adaptive surrogate's McCabe-Thiele column -not converging, attributed to entropy-driven acquisition) turned out to be **wrong** -- -a shared-mutable-state bug in `paper/full_reproduction.py`: `_predict_split` (test-RMSE) -mutates `GPmodel.kernel` in place via `bits_for_gaps.mixture.sample_gp_posterior_mixture` -(by design -- documented, not a core bug), and the phase-diagram code reused that same -mutated object afterward. Confirmed by reproducing the exact reported failure -(stage-2 liquid=1.7258, to 4 decimal places) from the checkpointed `gp_model_15.pkl`. - -Of the 4 hypotheses posed (draw-count/averaging, monotonicity, deterministic mean, -under-resolution), **none was the root cause** -- but testing them was still useful: -monotonicity is definitively rejected (all 5 reconstructed curves, including the failing -one, have zero violations), and the "deterministic posterior-mean" idea is falsified -(exactly as smooth as the working curves, yet fails -- a solver initial-guess- -sensitivity issue, not a curve-smoothness one). See `paper/PHASE9B_INVESTIGATION.md` for -the full analysis. - -Fix (example layer only, `src/bits_for_gaps/` untouched): reordered -`paper/full_reproduction.py` to build the phase diagram before the mutating RMSE loop, -and added `examples/vle_distillation/phase_diagram.py`'s `surrogate_gamma_averaged` -(matches the paper's own `new_phase_diagram.py` hyperparameter-posterior-averaging -construction) for added robustness. Re-ran the full 15-iteration stochastic loop from -scratch with the fix: `column_surrogate_converged` now `True`, stage table within 0.03 -mole fraction of Wilson at every stage. `pytest -q` (120/2) and `pytest -m vle` (2/120) -both still green; `paper/reproduce.py --figures 8 9` still works. Pre-fix and post-fix -`full_run_summary.json`s both committed under `paper/phase9_validation/` for the record. - -## Phase 9 — full stochastic reproduction + archive-free figures (done; merged to main) - -Two independent, separately-committed pieces (see REFACTOR_PLAN.md's Phase 9 entry -and §7 decision 4): - -**STEP 1 -- archive-free figure reproduction.** `paper/reproduce.py` no longer needs -private-archive access by default. Enumerated (from the code, not guessed) exactly -which files every `paper/figures/*.py` module reads via `paper/figures/_archive.py`'s -loaders, and copied exactly that set (~16 MB, well under the 30-50 MB budget -- -`gp_predict_2/3/4` turned out to not be read by anything, only `gp_predict_{1,15}`) -into a new **tracked** `paper/data/`, with provenance in `paper/data/README.md`. -Repointed `paper/reproduce.py`'s `DEFAULT_ARCHIVE` and -`tests/regression/test_paper_figures.py`'s `ARCHIVE_DIR` at it. Re-gated -`test_paper_figures.py`: the three tests that only *read* committed data (Fig 5 error -metrics, Fig 10 HMC diagnostics, the hyperparameter posterior behind Fig 11) moved out -of `@pytest.mark.vle` into the default suite; only the Fig 8 Wilson-curve cross-check -(recomputes via live Clapeyron) stays gated. `pytest -q` went from 117/5 to -**120 passed, 2 deselected**. `tests/regression/test_mccabe_thiele.py` needed **no** -change -- it has no `ARCHIVE_DIR` at all (Fig 9's stage table is a pure physics -recompute, no archived data read), contradicting the task's assumption that it needed -repointing. - -**STEP 2/3 -- from-scratch stochastic reproduction.** New `paper/full_reproduction.py` -drives `bits_for_gaps.sampler.BitsForGaps` through the paper's exact 15-iteration -adaptive HMC+entropy-acquisition config (10 train/10 test LHS, seed 10, bounds -`[[1e-6,0.999],[350,367]]`, `AnisotropicSE.paper_2d()`, `likelihood_var=0.1`, HMC -`noSamples=5000/noBurnIn=0/noChains=4/noLeapfrogSteps=5/stepSize=0.05/noAdaptSteps=5/ -targetAccept=0.9/adaptRate=0.1`, `noGaussians=15/noRestarts=10`) starting from a -*fresh* LHS design and a live Clapeyron/Wilson black box -- no archived data read. -Ran once (2026-07-04, ~25 min wall time on this machine); artifacts stayed in the -gitignored `results_remaked/phase9_fullrun/`. Committed: the script itself, and a -small (~200 KB) `paper/phase9_validation/` (2 summary PNGs + `full_run_summary.json`) -backing a new "Phase 9: from-scratch stochastic reproduction" section in -`paper/REPRODUCTION.md`. **Not gated in CI** -- one-time validation artifact, per the -task's explicit instruction. - -Headline finding, stronger than expected: everything in the loop that runs through a -`self.seed`-seeded path (LHS design, HMC, the entropy-acquisition optimizer) -reproduced the published run's R-hat/ESS, hyperparameter posterior, and entropy-decay -curve to **6-8 significant figures** -- not just "qualitatively similar." The *only* -place real stochastic drift shows up is `gpflow`'s `predict_f_samples` (documented in -`bits_for_gaps/mixture.py` as drawing from TensorFlow's non-seedable ambient RNG), -which feeds the test-RMSE curve (4.337 -> 0.887 here vs. the paper's ~4.34 -> ~0.67 -- -same regime, not the same value) and the surrogate phase diagram. One genuine, -documented discrepancy: the **genuinely** 15-iteration-adaptive surrogate GP's -McCabe-Thiele column did **not** converge to a physical solution (unlike -`fig09_mccabe_thiele.py`'s dedicated 30-point-LHS/MLE-fit stand-in, which does) -- -entropy-driven acquisition optimizes for predictive accuracy at held-out points, not -for a globally smooth-enough equilibrium curve for the stage-stepping solver. See -`paper/REPRODUCTION.md`'s Phase 9 section for full numbers and discussion. - -`src/bits_for_gaps/` core, `paper/reference/*`, and the dependency stack were not -touched. `import bits_for_gaps` confirmed still Julia-free. - -## Phase 8 — documentation (done; merged to main) - -`docs/` is a Sphinx 7 + MyST + furo site, repo-only (not shipped in the pip wheel -- -same policy as `examples/`/`paper/`), built via the existing `[docs]` extra. -`sphinx-build -W docs docs/_build/html` succeeds with **zero warnings** (confirmed -on repeated clean builds). `pytest -q` is unaffected: **117 passed, 5 deselected** -(docs are additive; no source files under `src/bits_for_gaps/`, `examples/`, or -`paper/` changed -- `git diff main -- ` is empty for all three). - -``` -docs/ - conf.py Sphinx config; see "RTD/TensorFlow decision" below - index.md overview, paper DOI, status, page map - installation.md core / dev / [docs] / optional [vle]+Julia install; explicit - "examples/ and paper/ are repo-only, clone required" note - quickstart.md pure-Python BitsForGaps(black_box, bounds, kernel= - AnisotropicSE()) example, as inert code blocks (not executed - at build time -- HMC is too slow for a docs build) - theory.md brief method summary + citation; links into the API reference - vle_example.md narrative walkthrough -> examples/vle_distillation/README.md - reproduce_paper.md narrative walkthrough -> paper/REPRODUCTION.md + - paper/DATA.md's private-archive-access note - api.rst automodule/autoclass over the public surface - Makefile `make html` (`sphinx-build -W`), `make clean` -.readthedocs.yaml RTD build config: Python 3.9, ubuntu-24.04, - sphinx.fail_on_warning: true, `pip install .[docs]` only -``` - -Key facts for the next session: - -- **RTD/TensorFlow build-robustness decision: install the real stack, don't mock - it.** Autodoc imports `bits_for_gaps`, which imports gpflow/tensorflow/tfp for the - TF-backed modules. `pip install ".[docs]"` (what `.readthedocs.yaml` runs) already - pulls in the pinned TF 2.16.2/GPflow 2.9.2/TFP 0.24.0 stack, because pip's extras - are *additive* to a package's base `dependencies` -- there is no way to install - `[docs]` without also installing the base stack. This gives autodoc a real, - importable `AnisotropicSE` (a `gpflow.kernels.Kernel` subclass) to - introspect -- verified directly by inspecting the rendered HTML: its page shows - the full constructor signature and "Bases: Kernel", which mocking the import out - would have suppressed. `docs/conf.py` documents the fallback - (`autodoc_mock_imports = [...]`, commented out) in case a future RTD build times - out or OOMs installing TF -- not needed today, not exercised, just ready if it - ever is. `juliacall`/`juliapkg` are **never** imported by the docs build at all - (no autodoc directives target `examples/vle_distillation`; that page is narrative - only) -- no mock needed for them, no `[vle]` extra installed by - `.readthedocs.yaml`. -- **Two real (not cosmetic-noise) issues found and fixed while verifying the - build:** (1) MyST's `dollarmath`/`amsmath` extensions were not enabled by - default, so `theory.md`'s `$$...$$` canonical-hyperparameter-ordering equation - silently rendered as literal text instead of math -- fixed by adding both to - `myst_enable_extensions` plus `sphinx.ext.mathjax`; verified the equation now - renders as an actual MathJax block. (2) `api.rst` had a cross-reference role - (`` :class:`~bits_for_gaps.sampler.adaptiveEntropy` ``) that got line-wrapped - mid-identifier when first drafted, which silently failed to resolve (a literal - newline inside the backticks becomes a space, breaking the dotted path) -- - fixed by keeping the whole role on one line. Both were caught by rendering the - actual HTML and checking for resolved anchors/MathJax output, not just by - `sphinx-build -W` succeeding -- **a clean warnings-as-errors build does not by - itself prove cross-references resolved or math rendered**, because - `nitpicky = False` (the default, kept deliberately -- see next point) doesn't - warn on unresolved `:class:`/`:func:` targets, and a `$$...$$` block that isn't - recognized as math just becomes an ordinary (silently valid) paragraph. -- **`nitpicky` is deliberately left `False` (i.e., unset -- Sphinx's default), - not enabled.** Turning it on was tried as a diagnostic (`sphinx-build -D - nitpicky=1`) and immediately produced ~40 warnings, almost all from NumPy-style - docstring parameter types (`np.ndarray, shape (n, m)`, `array-like`, `optional`, - `tfp.distributions.Distribution`, ...) that napoleon turns into cross-reference - attempts with no real Sphinx target -- these aren't bugs, they're the normal - cost of NumPy-style docstrings without a full `intersphinx` wiring to - numpy/scipy/tensorflow/gpflow/tfp's own docs. Fixing all of them would mean - either adding `intersphinx_mapping` (network-dependent at build time -- a - robustness risk for RTD builds, deliberately avoided; see the "no network - fetches" reasoning in this same decision) or extensively rewriting docstrings - purely for Sphinx's benefit, which the task explicitly says not to do ("docs - should surface them \[docstrings\], not rewrite them"). One genuine, low-value - broken reference remains as a result (`` :meth:`run` `` in a couple of - `sampler.py` docstrings, referring to an *inherited* method that isn't - separately autodoc'd since `:inherited-members:` isn't set) -- silently - harmless under `nitpicky = False`, left as-is rather than touched for a - docs-only cosmetic gain. -- **The quickstart is inert code, verified against the real API, not run.** Every - name/signature/attribute the quickstart page uses - (`BitsForGaps(black_box=..., bounds=..., kernel=..., likelihood_variance=...)`, - `bfg.noSamples`/`.noBurnIn`/`.noChains`, `bfg.run(X_init, y_init)`, - `history.last`, `record.xStar`/`.max_entropy`/`.rhat`, the - `black_box(*xStar)` calling convention with `fwd_model_args=()` by default) was - cross-checked directly against `sampler.py`'s `BitsForGaps.__init__` and - `state.py`'s `IterationRecord`/`RunHistory`, not run through the actual HMC loop - (which the task explicitly says not to execute at docs-build time). -- **RTD import steps for the user** (a maintainer action requiring the user's RTD - account -- NOT attempted this session, per the guardrail): - 1. Sign in at [readthedocs.org](https://readthedocs.org) with the GitHub account - that owns/has admin on `dowlinglab/bits_for_gaps`. - 2. "Add project" -> import `dowlinglab/bits_for_gaps` from the connected GitHub - account (or "Import Manually" with the repo URL if the GitHub App isn't - installed for this org yet). - 3. RTD auto-detects `.readthedocs.yaml` at the repo root -- no further build - config needed in the web UI; confirm the default branch (`main`) and doc - type (Sphinx) look right on the project's Admin > Settings page. - 4. Trigger a build (either automatically on import, or "Build Version" from the - project dashboard) and check it completes; if TF's install ever times out or - runs out of memory on RTD's builders, switch to the `autodoc_mock_imports` - fallback documented in `docs/conf.py` and rebuild. - 5. Optional: enable "Build pull requests for this project" under Admin > - Settings if PR doc previews are wanted; set up a custom domain / project - slug as desired. -- Doc-only touch-ups outside `docs/`: `README.md` gained a short "Docs" section - (build command + pointer to `docs/index.md`); no other file outside `docs/` and - `.readthedocs.yaml` changed. - -## Phase 7 — reproduce the paper's figures (done; merged to main) - -All 11 figures (2-12) regenerate from the archived published run through -`paper/reproduce.py` + `paper/figures/`. `pytest -q` (default) = **117 passed, 5 -deselected** (was 117/1 -- the 4 new gated tests in `test_paper_figures.py`). -`pytest -m vle` = **5 passed, 117 deselected** (~80 s; needs the archive, 4 of the 5 -also need Julia). Full details, including known discrepancies and simplifications: -`paper/REPRODUCTION.md`. - -``` -paper/ - __init__.py package marker (repo-only, not in the wheel -- same - policy as examples/, verified via `python -m build`) - reproduce.py CLI entry: --archive / $BFG_ARCHIVE_DIR (default: the - private old-repo path), --figures (subset), - --out-dir (default results_remaked/, gitignored) - REPRODUCTION.md figure -> script -> archived-inputs -> reference-diff table - figures/ - _archive.py shared loaders (rhat/ess, HMC traces, param_posterior_ - samples, activity_data, gp_predict, entropy, lhs_design, - cont_data, phase_diagram, gt_Wilson_data) + apply_plot_ - settings(); verified directly against the real archive - fig02_lhs_design.py Fig 2 -- visual - fig03_entropy_field.py Fig 3 -- visual - fig04_entropy_evolution.py Fig 4 -- visual - fig05_parity.py Fig 5 -- PINNED (fig5_error_metrics.json) - fig06_gp_posterior_surface.py Fig 6 -- visual - fig07_gp_posterior_isotherms.py Fig 7 -- visual - fig08_phase_diagram.py Fig 8 -- archive cross-check (no dedicated reference file) - fig09_mccabe_thiele.py Fig 9 -- PINNED (mccabe_thiele_stages.json, Phase 6); - wilson_column()/surrogate_column() moved here from - tests/regression/test_mccabe_thiele.py (Phase 6), which - now imports them instead of reimplementing - fig10_traces.py Fig 10 -- PINNED (hmc_diagnostics.json) - fig11_marginals.py Fig 11 -- PINNED (hyperparameter_posterior.json) - fig12_joint_marginals.py Fig 12 -- visual -tests/regression/test_paper_figures.py gated (@pytest.mark.vle) recompute-vs-reference - for Fig 5/8/10/11 (Fig 9 already gated in Phase 6) -tests/conftest.py + repo root on sys.path (alongside examples/, Phase 6) - so `import paper.figures.*` works without installing it -``` - -Key facts for the next session: - -- **The approach is "load archive, render through the new code" -- not "re-run the - loop."** Every figure either reads archived text/pickle files (all 11) or calls - live into `examples/vle_distillation`'s Clapeyron-backed physics for the - ground-truth curve (Fig 8, 9) -- none of them re-run the paper's 15-iteration - adaptive HMC loop (stochastic, expensive, and orthogonal to "does the new code - reproduce the published figure"). This was an explicit guardrail, not just a time- - saving shortcut: re-running would produce a *different* (though qualitatively - similar) stochastic realization, not a reproduction of the specific published one. -- **Packaging mirrors `examples/`'s policy exactly.** `paper/__init__.py` + - `paper/figures/__init__.py` make it a real package; `tests/conftest.py` now also - puts the repo root on `sys.path` (added to the existing `examples/` insert from - Phase 6) so `import paper.figures.fig10_traces` works without installing anything. - Wheel exclusion is automatic (hatchling's `packages = ["src/bits_for_gaps"]` is an - allowlist -- nothing outside `src/` is ever included regardless of `__init__.py` - presence), not separately re-verified this phase (Phase 6 already confirmed the - mechanism with `python -m build --wheel`). -- **Quantitative pins reuse `paper/reference/*` exactly as extracted in Phase 2** -- - no new reference files were added or existing ones modified (guardrail). Fig 8 has no - dedicated reference *file* (the paper doesn't report its curve as a scalar target); - its regression instead cross-checks the freshly-recomputed (live Clapeyron) Wilson - curve against the archived `gt_Wilson_data` the paper's own Fig 8 was built from -- - confirmed matching (z exactly, since both use the same `linspace(0,1,75)` grid; - T within 0.5 K; y1 within 0.02). -- **Fig 9's recompute logic moved, not duplicated.** `wilson_column()`/ - `surrogate_column()` lived in `tests/regression/test_mccabe_thiele.py` (Phase 6); - Phase 7 moved them into `paper/figures/fig09_mccabe_thiele.py` (since the figure - and the test need the exact same recompute) and the test now imports them. Verified - the gated test still passes after the move. -- **The gated `vle` marker is reused for "needs the private archive," not just - "needs Julia."** `test_paper_figures.py`'s 4 tests are marked `@pytest.mark.vle` - even though 3 of them (`fig10_traces`, `fig05_parity`, hyperparameter-posterior) - never touch Julia -- what actually gates all of them is needing the archive - directory, which is exactly as unavailable to most environments/CI as Julia is. - Introducing a separate marker for "needs archive" seemed like unwarranted plumbing - for a distinction without a practical difference in this repo; each test also has - its own `skipif` on the archive directory existing, so it degrades gracefully - (skip, not error) if pointed at a missing path. -- **Simplifications from the 847-line `fxns/mcmc_plotter.py`** (this is reproduction - code, not a library API -- ported pragmatically, not verbatim): dropped Fig 5's - zoomed inset and Fig 12's KDE contour overlay (both purely visual, not the - figure's quantitative content); Fig 3 is a 2x3 grid of the first 6 iterations - rather than 60 separate per-iteration files; Fig 12's "MAP" marker is the sample - nearest the coordinate-wise median (a cheap visual proxy), not a true density - mode -- use `hyperparameter_posterior.json`'s `mean`/`median` for a real point - estimate. Full list in `paper/REPRODUCTION.md`. -- **All 11 figures were visually spot-checked** against the archived PNGs' described - structure while building them (not just "the code runs") -- e.g. Fig 8 shows the - expected minimum-boiling-azeotrope T-x-y diagram with the archived surrogate - ensemble tightly tracking the freshly-recomputed Wilson dashed curve; Fig 6 shows - visibly tighter credible-interval wireframes and denser training coverage at - iteration 15 vs. iteration 1; Fig 4 shows the expected monotonic-ish decreasing - max-entropy trend across all 60 archived iterations. -- **`src/bits_for_gaps/` CORE, `examples/vle_distillation/` (Phase 6 physics), and - `paper/reference/*` are all byte-for-byte untouched** (`git diff main...HEAD -- - ` empty for each) -- Phase 7 only added `paper/figures/`, `paper/reproduce.py`, - `paper/REPRODUCTION.md`, one new test file, and extended `tests/conftest.py` + - refactored (not rewrote) `test_mccabe_thiele.py`, per the guardrails. - -## Phase 6 — port the VLE/distillation example (done; merged to main) - -The paper's H2O-PrOH case study now lives at `examples/vle_distillation/`, on the -public `bits_for_gaps` API, with the Julia/Clapeyron activity model injected as the -black box. `pytest -q` (default) = **117 passed, 1 deselected** (was 88 before Phase 6: -+29 new no-Julia unit tests for the example modules). `pytest -m vle` = **1 passed, 117 -deselected** (~54 s; needs Julia/Clapeyron) -- the gated stage-table regression now -actually recomputes and checks, instead of skipping. - -``` -examples/vle_distillation/ repo-only -- NOT in the pip wheel (verified via - `python -m build --wheel`: contents are only - bits_for_gaps/*, no examples/ or paper/) - __init__.py - activity_model.py Wilson gamma via Clapeyron.jl; LAZY juliacall import - calculate_activities.jl ported from fxns/calculate_activities.jl - juliapkg.json pins Clapeyron.jl =0.6.26 - gibbs_duhem.py gamma_water from a modeled gamma_proh curve. PURE. - phase_diagram.py Antoine + bubble/dew point; wilson_gamma (Julia) / - surrogate_gamma (bits_for_gaps GP + Gibbs-Duhem) - equilibrium.py wraps a VLE curve as x_liquid -> y_vapor. PURE. - distillation.py McCabe-Thiele column solver (fsolve). PURE. - run_case_study.py LHS -> Clapeyron -> BitsForGaps.run -> phase - diagram + column; paper's exact 2-D config - README.md setup + run, for a fresh clone -tests/conftest.py + sys.path insert of examples/ (see below) -tests/unit/ + test_gibbs_duhem/phase_diagram/equilibrium/ - distillation.py (no Julia, run by default) -tests/regression/test_mccabe_thiele.py gated recompute wired up (was a Phase-2 - placeholder that skipped) -``` - -Key facts for the next session: - -- **Packaging mechanism (REFACTOR_PLAN §7.3), concretely implemented:** - `tests/conftest.py` does `sys.path.insert(0, REPO_ROOT / "examples")`, which makes - `import vle_distillation.` work in dev/CI without installing `examples/` as a - distribution. This same mechanism made `examples/vle_distillation/juliapkg.json` - auto-discoverable: `juliapkg` scans every `//juliapkg.json` - (confirmed directly by reading `juliapkg/deps.py` and by the resolution log), so - once `examples/` is on `sys.path`, the Clapeyron pin is found with **no extra - wiring** -- no `Pkg.add` step, no environment variable. `run_case_study.py` does the - equivalent `sys.path` insert itself (of its own parent's parent) so it works when - invoked as a standalone script, not just under pytest. -- **Clapeyron.jl is pinned to 0.6.26** (uuid `7c7805af-46cc-48c9-995b-ed0ed2dc909a`, - the version already resolved in this machine's Julia depot -- read directly from - `~/.julia/environments/pyjuliapkg/Manifest.toml`). A fresh machine's first - `juliacall` import (anything calling into `activity_model.py`) auto-downloads Julia - itself plus this exact Clapeyron version into a per-conda-env directory - (`$CONDA_PREFIX/julia_env/`) -- no manual bootstrap step. See - `examples/vle_distillation/README.md`. -- **Lazy Julia, verified both directions:** `import vle_distillation.activity_model` - (and every other example module that imports it) succeeds with **zero** Julia - touched (`sys.modules` has no `julia*` entries) -- confirmed directly. Calling - `activity_coefficients`/`black_box`/`wilson_gamma` triggers the lazy import and, if - Julia/`[vle]` isn't installed, raises a clear `ImportError` naming the fix. `import - bits_for_gaps` (the core) is unaffected either way -- it never touches - `examples/vle_distillation/` at all. -- **Black-box adapter convention (the critical Phase-5-carryover contract):** - `adaptiveEntropy.call_model` calls `FwdModel(*FwdModelArgs, *xStar)` -- natural - dimension order, `xStar = [z_PrOH, T]` for this case study. - `activity_model.black_box(z_proh, temperature)` matches that signature directly and - returns `[gamma_proh]` (a 1-list, matching `call_model`'s `np.array(self.FwdModel(...))` - convention) -- **not** both coefficients: the GP surrogate models only - `gamma_PrOH(z, T)`; `gamma_water` is recovered via Gibbs-Duhem - (`gibbs_duhem.gamma_water_from_gamma_proh`), never learned by a second GP output. - Verified: `activity_coefficients(0.5, 350.0) == (1.4695, 1.7464)`, matching the - target sanity value `(1.469, 1.746)`. -- **The paper's exact as-run manuscript config isn't committed verbatim anywhere in - the old repo** (`driver_new.py`/`driver.py`'s own `run_test()` examples use a - different, exploratory `"testing_12"`/`"SAFTγMie"` config, not - `"less_x_new_manuscript_revisions"`/`"Wilson"`). `run_case_study.py`'s config - (bounds `[(1e-6, 0.999), (350, 367)]`, `XTrsfFwd = [log(x+0.1), (T-350)/17]`, - `yTrsfFwd = log`, `thermoModel = "Wilson"`, seed `10`) was reconstructed from the - files that DO reference that experiment directly: - `train_test_split_proh.py`'s commented-out `run_test(...)` call (transforms, bounds, - seed) and `new_phase_diagram.py`'s `__main__` block (same transforms, confirms - `AnisotropicSE`'s exact 2-D config == `kernels.AnisotropicSE.paper_2d()`). This is a - physically-faithful reconstruction, not a byte-for-byte-verified original script -- - flagged here for anyone who later finds the real one. -- **The distillation solver's `fsolve` is fragile for arbitrary equilibrium curves** - (no bounds, inherited from the original MATLAB-derived port) -- confirmed - empirically: several hand-picked synthetic constant-relative-volatility test curves - either failed to converge or converged to spurious (unphysical, e.g. negative flow - rate) roots, while the real Wilson curve (with a 50-point z-grid, `Z_MESH`/`Z_GRID_SIZE`) - converges cleanly to the real Geankoplis 11.4-1 column every time it was tried. - `distillation.solve_column` now returns a `"warnings"` list (nonphysical mole - fractions/flow rates, fsolve non-convergence) instead of silently returning bad - numbers; both `run_case_study.py` and the gated regression test check `"converged"` - before trusting a result. **Do not shrink the z-grid casually** -- a coarser - (15-point) grid was enough to flip an otherwise-clean Wilson solve to a spurious - root in testing. -- **The gated regression test's "surrogate" recompute intentionally does not use the - full adaptive `BitsForGaps.run`/HMC loop.** It trains a GP on a 30-point LHS design - (seed 10) evaluated against the same Clapeyron Wilson model, then fits it with - `gp.maximize_lml` (a fast, deterministic MLE point estimate) -- reproducing the - paper's real 15-iteration *adaptively*-designed surrogate bit-for-bit is Phase 7's - job (full figure reproduction), not this phase's backend-correctness check. This - recompute matches reference's `"wilson"` column within `atol=0.015` and `"surrogate"` - within `atol=0.05` (looser -- a non-adaptive, far-smaller-sample surrogate doesn't - track Wilson as tightly in the most dilute region near `xW=0.01` as the paper's - refined surrogate did; the biggest observed gap, ~0.03-0.04, is in stage 4's vapor - fraction there). A quick empirical note: a 4-iteration adaptive run from only 15 - initial points did *worse* on this metric than the 30-point plain LHS + MLE fit - (too little HMC/data to beat a well-covering static design) -- consistent with the - paper's own thesis that adaptive design needs enough iterations to pay off, not - evidence against the adaptive loop itself. -- **Reference's own `"wilson"` column has ~0.01-level transcription slop.** It was - hand-transcribed from reading paper Fig 9c (Phase 2, `paper/reference/README.md`), not - computed from archived data. My direct Clapeyron recompute reproduces real physical - landmarks exactly (`stage 1 vapor == xD == 0.43` exactly, by construction; pure-PrOH - bubble point `370.35 K` matches 1-propanol's real normal boiling point to 4 - significant figures) yet differs from reference's `"wilson"` entries by up to `0.01` - (e.g. stage 4 liquid: recompute is exactly `xW = 0.01`, reference's transcription says - `0.02`) -- this is the eyeballed-figure-reading precision limit of that column, not - a bug in the port. `paper/reference/*` is unmodified (guardrail); the test's - `atol=0.015` for Wilson already accounts for this. -- **`src/bits_for_gaps/` CORE is byte-for-byte untouched** (`git diff main...HEAD -- - src/bits_for_gaps/` is empty) -- Phase 6 only added `examples/` + tests + - `tests/conftest.py`'s `sys.path` insert, per the guardrails. - -## Phase 5 — generalize to N-D (done; merged to main) - -The 2 inputs / 3 hyperparameters hardcoding flagged by the `TODO(Phase 5)` markers is -gone. `pytest -q` = **88 passed, 1 deselected** (same `vle` marker). Full suite runs in -~65 s (was ~35 s after Phase 4 — the new 1-D/3-D synthetic tests each run the tiny HMC -pipeline, same as the existing 2-D one, just at two more dimensions). - -**The whole phase was executed as: pin the 2-D baseline as the regression oracle, then -generalize one module at a time, running the full suite (incl. the atol=1e-10 baseline -pin) after every change before moving to the next module.** Every commit in this phase -kept `test_matches_pre_phase4_baseline` green — the 2-D path is bit-exact with the -pre-Phase-5 code throughout, not just at the end. - -### Design decision: per-dimension scalar Parameters (not a vector Parameter) - -`kernels.AnisotropicSE` now takes `variance_prior` + `lengthscale_priors` (a list, one -prior per input dimension) instead of hardcoding `lengthscale_1`/`lengthscale_2`. Each -lengthscale — and the variance — is its **own `gpflow.Parameter`**, not a slice of one -vector-valued Parameter. This was the one real design choice in this phase, and it's not -arbitrary: - -- **The paper's method depends on per-dimension prior *families*, not just per-dimension - prior *parameters*.** `std_dev` ~ LogNormal, `lengthscale_1` ~ LogNormal, `lengthscale_2` - ~ Gamma — three different distribution families, one of which (`lengthscale_2`) is also - deliberately left **unconstrained** (`Identity` transform, no positivity bijector, - confirmed empirically: `gpflow.Parameter`'s default `transform` is `Identity`, not - `positive()` — this was already true of the pre-Phase-5 kernel and had to be preserved - bit-for-bit). A single vector Parameter carries exactly one prior distribution and one - bijector for the whole vector — it cannot express "component 2 is Gamma-unconstrained, - the rest are LogNormal-positive" without slicing hacks that would themselves need - per-component metadata, i.e. would reinvent per-dimension Parameters anyway. -- **It matches gpflow's HMC machinery with no adapter.** `gpflow.optimizers.SamplingHelper` - takes a flat list of trainable `Parameter`s; each contributes its own prior term to the - log-posterior and its own bijector to the unconstrained HMC state. A list of scalar - Parameters already *is* that list — `GPmodel.kernel.hyperparameters` is passed straight - through to `SamplingHelper` (see below), no wrapping/unwrapping needed. - -The tradeoff (documented in `kernels.py`'s module docstring) is more Parameter objects to -manage than one vector — acceptable, since gpflow's own tooling (`print_summary`, -checkpointing, `trainable_parameters`) already expects a flat list of scalar Parameters. - -### Canonical hyperparameter order (the contract across gp/mixture/acquisition) - -``` -[std_dev, lengthscale_1, lengthscale_2, ..., lengthscale_d] -``` - -Exposed as `AnisotropicSE.hyperparameters` (a list of `gpflow.Parameter`, in this exact -order) and by name as `.std_dev`, `.lengthscale_1`, ... `.lengthscale_d` (kept for -backward compatibility with 2-D-era code addressing them by attribute name — the existing -`test_kernels.py` tests referencing `.lengthscale_1`/`.lengthscale_2` still pass -unmodified). This order is used in exactly three places, all now consistent by -construction rather than by hardcoded agreement: - -- `gp.run_mcmc`: `SamplingHelper(GPmodel.log_posterior_density, GPmodel.kernel.hyperparameters)` - — was `[trainable_parameters[2], [0], [1]]`. **Verified identity-equal** (same Python - `Parameter` objects, same order) to the old hardcoded indexing for the paper's kernel - before making this change — so this is a pure refactor at d=2, not a behavior change. - `trace`/`chains_states`/`rhat`/`ess` columns are in this order. -- `mixture.sample_gp_posterior_mixture` / `acquisition.entropy_objective`: replay a trace - row onto the kernel via the new `kernels.assign_hyperparameters(kernel, values)` — - `for param, value in zip(kernel.hyperparameters, values): param.assign(value)` — instead - of three hardcoded `.assign()` calls by name. Works for any kernel exposing - `.hyperparameters`, not just `AnisotropicSE`. -- Anywhere reading a trace column back out (regression tests, `paper/reference/*`) already - used this same order by convention; nothing there changed. - -### What's still 2-D-only, and why that's the right call - -`acquisition.entropy_surface_2D` and `mixture.predict_grid_2D` (the dense-grid entropy -field and the full-grid GP-prediction plotting diagnostic) are **not** generalized to N-D -— a dense grid is exponential in the input dimension, and neither one feeds the -acquisition (the actual next-point decision). Both now raise a clear `ValueError` if -called with `len(x_bounds) != 2`. `adaptiveEntropy.run()` calls `entropy_surface_2D` -*only* when the input space is 2-D, leaving `entropy_field=None` otherwise — an N-D run -does not error, it just skips a diagnostic it was never going to use. If N-D -visualization is ever needed, the right tool is a *sparse* diagnostic (e.g. entropy along -1-D/2-D slices through the current best point), not a full grid — not built here since -nothing calls for it yet. - -By contrast, **`acquisition.optimize` (was `optimize_2D`) — the function an N-D run -actually depends on — is fully dimension-general**: Sobol dimension and the restart -bound-scaling both derive from `len(x_bounds)`. For d=2 this reproduces the pre-Phase-5 -`optimize_2D` bit-for-bit (verified via the baseline pin): the vectorized bound-scaling -`lo + x0 * (hi - lo)` is the same floating-point operation as the original -`x0[j] * (hi[j] - lo[j]) + lo[j]` (IEEE 754 addition is commutative, so reordering the -addends doesn't change the rounding). - -### `call_model`'s black-box calling convention changed (interface, not algorithm) - -Pre-Phase-5, `call_model` called the injected black box as `FwdModel(*FwdModelArgs, x2, -x1)` — reversed, 2-D-specific argument order inherited from the VLE example's Julia -activity-coefficient function (which took `(T, x)`). This doesn't generalize to N inputs. -Phase 5 changes it to `FwdModel(*FwdModelArgs, *xStar)` — `xStar`'s components in natural -dimension order. This is an **interface** change to how the sampler calls the user's -injected function, not an algorithm change: `tests/integration/test_end_to_end.py`'s -`_fwd_model` was updated from `_fwd_model(x2, x1)` to `_fwd_model(x1, x2)` to match, and -produces the exact same `(x1, x2, y)` values as before (verified via the atol=1e-10 -baseline pin) — reordering which positional slot carries which value doesn't change the -value itself. **Phase 6 will need to account for this** when porting the VLE example's -Julia `fwd_model` wrapper (it can no longer rely on the reversed-argument convention; wrap -the Julia call so its own signature accepts `(x1, x2, ...)` in natural order). - -### Other Phase 5 changes - -- `sampler.py`: `adaptiveEntropy.optimize_2D` renamed to `.optimize` (matches - `acquisition.optimize`); `predict_grid_2D`/`entropy_surface_2D` method names kept - as-is (explicitly 2-D-only). `BitsForGaps`'s constructor and `.run()` are unaffected — - it already accepted any `bounds`/`kernel`, so N-D "just worked" once the modules - underneath it did (proven by `tests/integration/test_nd_synthetic.py` running 1-D/3-D - problems through `BitsForGaps`, not `adaptiveEntropy` directly). -- New tests: `tests/integration/test_nd_synthetic.py` (1-D and 3-D synthetic problems, - pure-Python, no Julia, via `BitsForGaps.run(...)`) and extensions to - `tests/unit/test_kernels.py` (canonical order, `paper_2d()` parity, - `assign_hyperparameters` round-trip, explicit 1-D/3-D construction with mixed prior - families, constructor validation). -- `entropy.py` / `design.py` / `means.py` / `_util.py` are byte-for-byte untouched - (`git diff main...HEAD -- ` is empty). `paper/reference/*` untouched; - no `results/` committed. - -## Phase 4 — decompose sampler.py; retire disk-as-state (done; merged to main) - -`sampler.py`'s `adaptiveEntropy` god-class is now an orchestrator over decomposed, -independently-testable modules. Order of operations (per the task): pinned an exact- -value characterization baseline BEFORE touching any code, then extracted modules one at -a time (verifying bit-exact parity against the still-untouched monolith before wiring -each one in), then rewrote `sampler.py` itself and retired disk-as-state in one final -commit. `pytest -q` = **72 passed, 1 deselected** (same `vle` marker as Phase 2). Full -suite runs in ~35 s (was ~17 s before Phase 4 — the two new modules' extra HMC/optimize -runs and the `BitsForGaps` facade parity test each re-run the tiny synthetic pipeline). - -``` -src/bits_for_gaps/ - gp.py build_gp / maximize_lml / run_mcmc (GP construction + HMC + R-hat/ESS) - diagnostics.py thin tfp.mcmc wrappers (potential_scale_reduction, effective_sample_size) - mixture.py sample_gp_posterior_mixture (GMM predictive posterior) + predict_grid_2D - (the plotting-only full-grid diagnostic, formerly gp_predict_2D) - acquisition.py entropy_objective / entropy_surface_2D (was gen_entropy_surface_data_2D) - / optimize_2D -- uses entropy.py; kept the *_2D names (Phase 5 generalizes) - transforms.py InputTransform / OutputTransform -- lifts the XTrsfFwd/XTrsfBkwd/ - yTrsfFwd/yTrsfBkwd lambda-list convention into small testable classes, - identity by default. Exported eagerly from __init__ (pure NumPy). - state.py IterationRecord / RunHistory -- the in-memory replacement for disk- - as-state (np.savetxt/pickle under results/{exp_name}/, read back next - iteration) - sampler.py adaptiveEntropy is now the orchestrator: thin delegating wrappers over - the modules above, plus `run(X_init, y_init, checkpoint_dir=None, - predict_grid=False)` -- the new entry point. Also adds `BitsForGaps`, - a thin renamed-kwarg subclass (target public API, REFACTOR_PLAN §4). -tests/unit/ + test_transforms.py, test_state.py -tests/integration/ test_end_to_end.py rewritten to drive run() (in-memory, no files); - + data/synthetic_baseline.json (exact-value pin from before Phase 4) - + test_bits_for_gaps_facade.py (BitsForGaps reproduces the same pin) -``` - -Key facts for the next session: -- **Disk-as-state is retired.** `adaptiveEntropy.run(X_init, y_init)` takes the initial - design directly in memory and returns a `state.RunHistory`; a full run writes **zero** - files by default (`test_run_writes_no_files_by_default` asserts this). File output is - opt-in via `run(..., checkpoint_dir=...)`, a best-effort equivalent of the paper code's - per-iteration dump (`test_run_checkpoint_dir_is_opt_in` asserts the key files land). - `run_model()` is kept as a deprecated, disk-based shim (reads `activity_data_1` from - `self.path`, the original zero-arg precondition) for anything still relying on it. -- **`predict_grid_2D` (was `gp_predict_2D`) is opt-in, not run by default.** It's a - plotting-only diagnostic (~20 s in real config: 100 full-covariance draws over a 50x50 - grid) that doesn't feed the acquisition. Pass `run(..., predict_grid=True)` to compute - it anyway. **It was never bitwise-reproducible even in the original paper code** -- - confirmed directly: `GPmodel.predict_f_samples` draws from TensorFlow's ambient - (unseeded) global RNG, not NumPy's, so two successive calls on identical inputs in the - same process already differ. `np.random.seed(self.seed)` in `sample_gp_posterior_mixture` - only controls *which* posterior components are selected, not the draws themselves. This - is exactly why the Phase 2 integration test excluded it from the determinism/baseline - pins in the first place -- documented now in `mixture.py`. -- **The decomposition is verified behavior-preserving, not just self-consistent.** - `tests/integration/data/synthetic_baseline.json` pins the tiny seeded synthetic run's - exact outputs (rhat, ess, entropy field, xStar, max_entropy, next_data) captured from - the pre-Phase-4 monolith; `test_matches_pre_phase4_baseline` checks the post- - decomposition `run()` reproduces them at atol=1e-10. Additionally, before rewriting - `sampler.py`, each extracted module (gp.py, mixture.py's non-predict_grid parts, - acquisition.py) was checked bit-exact (atol=1e-12) against the still-untouched - monolithic methods on the same run. -- **`BitsForGaps` is a thin subclass, not a rewrite.** Renamed constructor kwargs - (`black_box`, `bounds`, `kernel`, `likelihood_variance`, `input_transform`, - `output_transform`) matching REFACTOR_PLAN §4's target API; inherits every method - (including `run`) unchanged from `adaptiveEntropy` -- no new computation, so no - numeric risk. `__init__.py`'s lazy map now resolves `BitsForGaps` to this real class - (previously it was just an alias for `adaptiveEntropy`). Advanced config (HMC tuning, - restarts, mesh density) is still set via the same instance attributes as - `adaptiveEntropy` (e.g. `.noSamples`) -- an `mcmc=MCMCConfig(...)`-style kwarg is - deferred to Phase 5/6, once the core is N-D and doesn't need this passthrough. -- **`entropy.py` / `design.py` / `kernels.py` / `means.py` are byte-for-byte untouched** - (`git diff main...HEAD -- ` is empty) -- Phase 4 touched only the - sequential-design engine, per the guardrails. -- `TODO(Phase 5)` markers preserved verbatim in `gp.py` (`run_mcmc`'s positional - `trainable_parameters[2],[0],[1]`), `mixture.py`/`acquisition.py` (by-name kernel - param assignment: `std_dev`, `lengthscale_1`, `lengthscale_2`), and `acquisition.py`/ - `sampler.py`'s `*_2D` methods (hardcoded `d=2`). None of the 2-D / 3-hyperparameter - hardcoding was touched. - -## Phase 2 — regression/test harness (done; merged to main) - -Behavior is now pinned BEFORE any refactor. `pytest -q` = **56 passed, 1 deselected** -(the deselected one is the `vle` McCabe-Thiele recompute; the pure-Python core needs no -Julia). Full suite runs in ~16 s. - -``` -paper/reference/ reference scalars from the archived iter-15 published run (+ README) - hmc_diagnostics.json R-hat / ESS (Fig 10) - fig5_error_metrics.json train/test RMSE & MAE, iters 1 & 15, over 500 draws (Fig 5) - hyperparameter_posterior.json kernel-hyperparam posterior summary (Fig 10 marginals) - mccabe_thiele_stages.json distillation stage table, Wilson vs surrogate (Fig 9c) -tests/conftest.py `reference` loader fixture (resolves paper/reference/) -tests/unit/ + test_kernels.py, test_means.py, test_util.py (was entropy/design) -tests/regression/ reads reference, pins vs published paper values; vle recompute gated -tests/integration/ test_end_to_end.py — seeded synthetic (no-Julia) adaptiveEntropy run -``` - -Key facts for the next session: -- **Reference extraction** was done offline by `scratchpad/extract_reference.py` (pure NumPy, - reads the read-only old-repo archive). All but the McCabe-Thiele table came from - archived data; the stage table is transcribed from paper Fig 9c (its recompute needs - the Phase-6 VLE backend → `@pytest.mark.vle`, deselected by default via `addopts`). -- **The sampler is deterministic** for a fixed seed in-process (verified bitwise: R-hat, - ESS, entropy field, selected point all diff = 0.0 across two runs). The integration - test asserts this with `atol=1e-10` — it is the guard against nondeterminism creeping - into the Phase 4 decomposition. -- The integration test **mirrors `run_model` but skips the plotting-only `gp_predict_2D`** - (~20 s; 100 full-cov draws over a 50×50 grid). It does not feed the acquisition - (`entropy_objective` re-seeds NumPy and re-assigns every kernel param before each - deterministic `predict_f`), so entropy/next-point are unchanged by skipping it. -- Fig 5 reference captures the paper's headline: median **test** RMSE falls 4.34 → 0.67 - (iter 1 → 15); train RMSE 0.77 → 0.49. The regression test pins the ≥3× test-error drop. -- Markers registered in `pyproject.toml`: `vle` (Julia backend, deselected by default), - `slow` (integration; still runs by default). Run gated tests with `pytest -m vle`. - -## What exists now (Phase 0 + Phase 3-lite + Phases 4-9d done) - -A pip-installable package with the **algorithm decomposed into focused, tested modules**, -**generalized to N input dimensions**, a **ported VLE/distillation example** repo-side, -**all 11 paper figures reproduced** repo-side (now archive-free by default), -a **Sphinx/MyST/RTD docs site**, and a **from-scratch stochastic reproduction** of the -full adaptive loop (see the "Phase 4"-"Phase 9" sections above): - -``` -src/bits_for_gaps/ - __init__.py public API; pure pieces eager (entropy, design, transforms), - TF-backed pieces lazy (kernels, means, sampler) via PEP 562 - entropy.py GMM density + 1st/2nd-order Taylor entropy + closed-form lower bound - (from fxns/max_ent_design.py; dead commented variants removed). PURE. - design.py latin_hypercube_design / full_factorial_design, N-D, pure, no disk I/O - (extracted from proh_water_class). PURE. - kernels.py AnisotropicSE (N-D, per-dimension prior-bearing Parameters, Phase 5) + - assign_hyperparameters(kernel, values) - means.py FixedInverseMean from fxns/my_mean_fxn.py - _util.py standardize / normalize / make_tensor - transforms.py InputTransform / OutputTransform (Phase 4). PURE. - state.py IterationRecord / RunHistory (Phase 4). - diagnostics.py R-hat / ESS (Phase 4). - gp.py GP construction + HMC, N-D via kernel.hyperparameters (Phase 4/5). - mixture.py GMM predictive posterior, N-D (Phase 4/5); predict_grid_2D stays 2-D-only. - acquisition.py entropy-maximization acquisition; optimize is N-D (Phase 5), - entropy_surface_2D stays 2-D-only (raises for d != 2). - sampler.py adaptiveEntropy (orchestrator) + BitsForGaps (public-API facade) -- N-D - throughout except the 2-D-only diagnostics, which run() skips for d != 2. -tests/unit/ entropy/design/kernels(+N-D)/means/_util/transforms/state + - mixture/acquisition/sampler_validation (Phase 9c: mutation-footgun - state-restore + input-validation error paths) + gp/ - sampler_legacy_and_transforms/entropy_mc_validation (Phase 9d) -tests/integration/ end-to-end (2-D, in-memory run()) + BitsForGaps facade parity + - nd_synthetic (1-D and 3-D, via BitsForGaps) -tests/regression/ reference-file checks vs published paper values (Phase 2, 2-D only -- - the published run and paper/reference/* are inherently 2-D); + - test_mccabe_thiele.py's (Phase 6) and test_paper_figures.py's - (Phase 7) @pytest.mark.vle recomputes -examples/vle_distillation/ the H2O-PrOH case study on the public API (Phase 6) -- - repo-only, not in the pip wheel; see the "Phase 6" section above -examples/synthetic/ Julia-free onboarding example (Phase 9d) -- repo-only, run - `python examples/synthetic/run_example.py` -src/bits_for_gaps/py.typed PEP 561 marker (Phase 9d); every module in this - directory is type-hinted via `from __future__ import annotations` -paper/figures/ + paper/reproduce.py all 11 paper figures reproduced (Phase 7), - archive-free by default as of Phase 9 via the tracked paper/data/ - -- repo-only, not in the pip wheel; see the "Phase 7"/"Phase 9" - sections above -docs/ + .readthedocs.yaml Sphinx/MyST/furo docs (Phase 8) -- repo-only, not in - the pip wheel; see the "Phase 8" section above for RTD import steps -paper/data/ curated ~16 MB plot-input subset (Phase 9) -- tracked, not - gitignored; see paper/data/README.md -paper/full_reproduction.py + paper/phase9_validation/ from-scratch stochastic - reproduction of the full adaptive loop (Phase 9) -- one-time - validation artifact, not gated in CI; see the "Phase 9" section - above and paper/REPRODUCTION.md -``` - -`BitsForGaps` (public-API facade, REFACTOR_PLAN §4 kwarg names) and `adaptiveEntropy` -(original name, kept for backward compatibility) are both real, independently -importable classes; both work at any input dimension as of Phase 5. - -## Environment (verified working) - -```bash -conda env create -f environment.yml # or reuse existing `bits_for_gaps` env -conda activate bits_for_gaps # /opt/anaconda3/envs/bits_for_gaps -pip install -e ".[dev,vle]" # ,vle needed for the Julia example -pytest -q # 204 passed, 2 deselected (as of Phase 9e) -ruff check . # lint (Phase 9d); ruff is in [dev] -``` -Stack: Python 3.9.23, gpflow 2.9.2, TF 2.16.2, TFP 0.24.0, numpy 1.26.4, scipy 1.13.1. -**macOS: `export PYTHON_JULIACALL_HANDLE_SIGNALS=yes`** before any juliacall import -(else SIGBUS). Run scripts that use bare imports with `PYTHONPATH` set to the dir. - -## Where the original code + archived results live - -Old repo: `~/DowlingLab/CAREER/entropy_driven_hybrid_models_code/entropy_driven_hms/`. -- Published run = `results/less_x_new_manuscript_revisions/`, **iteration 15** (R̂/ESS - match paper Fig 10 exactly). Archived figure PNGs + data are all there. Full run - directory is 564 MB (not 2.5 GB -- that was the old repo's whole git history). -- `driver_new.py` = active driver (NOT `driver.py`). `new_phase_diagram.py` = Fig 8. - `run_example.py` = Fig 9. `train_test_split_proh.py` = Fig 5. `fxns/mcmc_plotter.py` - + `fxns/plot_res.py` = figure library/CLI. -- Repro-fix already applied there: `equilibrium.water_proh_eqm_julia` now reads - `gt_Wilson_data` (original path was missing). -- **Data decision (updated Phase 9, supersedes the old "never copy results/ here" - line): a curated ~16 MB plot-input subset IS committed, at `paper/data/`** (see - `paper/data/README.md` for the exact file manifest and `paper/DATA.md` for the - full policy). The full 564 MB run stays in this private old repo (archive of - record) -- **no Zenodo deposit** (REFACTOR_PLAN.md §7 decision 4). - -## Next steps (in order — see REFACTOR_PLAN.md phases) - -1. **Phase 2 — regression harness FIRST (before refactoring sampler.py). ✅ DONE.** - Merged to `main`. - -2. **Phase 4 — decompose `sampler.py`. ✅ DONE.** Merged to `main`. - -3. **Phase 5 — generalize to N-D. ✅ DONE.** Merged to `main`. - -4. **Phase 6 — port the VLE example. ✅ DONE.** Merged to `main`. - -5. **Phase 7 — reproduce all paper figures. ✅ DONE.** Merged to `main`. - -6. **Phase 8 — docs (Sphinx/RTD). ✅ DONE.** Merged to `main`. See the "Phase 8" - section above. Actually connecting the repo on readthedocs.org is a maintainer - action requiring the user's RTD account -- not attempted this session (see the - numbered steps in the "Phase 8" section); `.readthedocs.yaml` is ready and waiting - for that one manual step. - -7. **Phase 9 — full stochastic reproduction + archive-free figures. ✅ DONE.** - Merged to `main`. See the "Phase 9" section above. - -8. **Phase 9b — root-cause + fix the McCabe-Thiele discrepancy. ✅ DONE.** Merged to - `main`. See the "Phase 9b" section above. - -9. **Phase 9c — robustness hardening. ✅ DONE.** Merged to `main`. See the "Phase 9c" - section above. - -10. **Phase 9d — whole-codebase polish (hygiene + faithfulness). ✅ DONE.** Merged to - `main`. See the "Phase 9d" section above. - -11. **Phase 9e — docs/tests/comments quality pass. ✅ DONE.** Merged to `main`. See - the "Phase 9e" section above. - -12. **Phase 10 — release engineering for v0.1.0. ✅ DONE.** On branch - `phase10-release`; merge to `main` after review, then the maintainer performs the - steps in `RELEASE.md`. See the "Phase 10" section above. No Zenodo deposit - (REFACTOR_PLAN.md §7 decision 4 -- the private old repo is the archive of record). - -13. **Actual publish (maintainer-only, not part of any automated phase).** Everything - up to this point is prepared and documented in `RELEASE.md`: register a PyPI + - TestPyPI trusted publisher, dry-run via `workflow_dispatch` to TestPyPI, set the - CHANGELOG date and `git tag v0.1.0 && git push --tags` to trigger the real PyPI - publish, then activate ReadTheDocs (import steps in the "Phase 8" section above). - -## Known issues / decisions already made (do not re-litigate) - -- Package name `bits_for_gaps`; examples+paper in this repo; 2-D faithful first then N-D; - freeze the current dependency stack; core is pure-Python (Julia only for `[vle]`). - (REFACTOR_PLAN.md §7.) -- `entropy_lower_bound` was unused in the paper code but is a real algorithm feature - (the closed-form bound) — kept and tested. Wire it into acquisition as an option. -- CI (`.github/workflows/ci.yml`) is a skeleton; the frozen stack is finicky on Linux - runners — keep Julia/regression tests in a separate gated job. -- LICENSE is BSD-3-Clause; copyright holders Alexander W. Dowling and Kyla D. Jones - (University of Notre Dame). Authors set accordingly in pyproject.toml. diff --git a/README.md b/README.md index 62e6efa..d0b896f 100644 --- a/README.md +++ b/README.md @@ -4,13 +4,6 @@ [![PyPI](https://img.shields.io/pypi/v/bits_for_gaps.svg)](https://pypi.org/project/bits_for_gaps/) [![Docs](https://readthedocs.org/projects/bits-for-gaps/badge/?version=latest)](https://bits-for-gaps.readthedocs.io/en/latest/?badge=latest) - - **B**ayesian **I**nformation-**T**heoretic **S**ampling for hierarchical **GA**ussian **P**rocess **S**urrogates. A framework for information-theoretic sequential experimental design with Bayesian @@ -23,25 +16,35 @@ Reference: K. D. Jones and A. W. Dowling, "BITS for GAPS: Bayesian Information-T Sampling for hierarchical GAussian Process Surrogates," *Computers & Chemical Engineering* **211** (2026) 109650. https://doi.org/10.1016/j.compchemeng.2026.109650 -> **Status: pre-1.0, under active refactor.** This repository is being extracted from the -> paper's research code into a reusable library. See `REFACTOR_PLAN.md` for the roadmap and -> `HANDOFF.md` for the current state. +The paper is **bundled in this repository** so you can read the method alongside the code: +[`paper/bits_for_gaps_paper.pdf`](paper/bits_for_gaps_paper.pdf). It is redistributed under +[CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) (© 2026 The Authors, published by +Elsevier Ltd); the DOI above is the canonical citation. -## Install (development) +## Install ```bash -conda env create -f environment.yml -conda activate bits_for_gaps -pip install -e ".[dev]" # core + test tools -# pip install -e ".[dev,vle]" # add the Julia/Clapeyron VLE example backend +pip install bits_for_gaps ``` The **core library is pure Python** (GPflow / TensorFlow / NumPy / SciPy) with no Julia -dependency. Julia + Clapeyron are only needed for the `vle_distillation` example. +dependency. Julia + Clapeyron are only needed for the `vle_distillation` example, which +isn't part of the PyPI package -- see "From source" below. **macOS note:** set `export PYTHON_JULIACALL_HANDLE_SIGNALS=yes` before importing `juliacall`, or Julia crashes with a bus error (SIGBUS). +### From source (for `examples/`, `paper/`, and development) + +```bash +git clone https://github.com/dowlinglab/bits_for_gaps +cd bits_for_gaps +conda env create -f environment.yml +conda activate bits_for_gaps +pip install -e ".[dev]" # core + test tools +# pip install -e ".[dev,vle]" # add the Julia/Clapeyron VLE example backend +``` + ## Layout ``` @@ -58,13 +61,26 @@ docs/ Sphinx documentation (ReadTheDocs) pytest -q ``` +## Provenance + +The research code behind the paper was originally developed in a private repository over the +course of the study. It was then migrated here and reorganized into an installable, tested +package: the algorithm was separated from the vapor–liquid-equilibrium case study, generalized +to arbitrary input dimension, and covered by a test suite. That private repository holds only +the development history — **nothing you need to use this package or to reproduce the paper's +figures is missing from this repository.** The data the figure scripts read is committed here +under `paper/data/` (see [`paper/REPRODUCTION.md`](paper/REPRODUCTION.md)). + ## Docs +Full docs (installation, a pure-Python quickstart, theory notes, the VLE example, +reproducing the paper's figures, and the API reference): +https://bits-for-gaps.readthedocs.io + +To build and browse locally instead: + ```bash pip install -e ".[docs]" sphinx-build -W docs docs/_build/html +open docs/_build/html/index.html # or your platform's equivalent ``` - -Full docs (installation, a pure-Python quickstart, theory notes, the VLE example, and -the API reference): `docs/index.md`, or built and browsed locally as above. See -`HANDOFF.md` for ReadTheDocs setup status. diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md deleted file mode 100644 index f92d765..0000000 --- a/REFACTOR_PLAN.md +++ /dev/null @@ -1,221 +0,0 @@ -# Refactoring Plan — `bits_for_gaps` - -**Goal:** Extract the *already-conceptually-general* BITS for GAPS algorithm from the paper's VLE/distillation example code and ship it as a pip-installable, documented, tested Python library. Reproduce the paper figures from a clean public API. - -**Source:** `entropy_driven_hybrid_models_code/entropy_driven_hms/` (branch `main`), with fixes and env from branch `codex-refactor`. -**Paper:** Jones & Dowling, "BITS for GAPS," *Computers & Chemical Engineering* 211 (2026) 109650. CC BY 4.0. - ---- - -## Progress log - -### Session 1 (2026-07-03, Opus) — audit + env + reproduce-from-archive ✅ -- **Audit (Step 1) done** — algorithm/example split mapped (§2), Codex branch mined, git history analyzed (fresh repo confirmed). -- **Environment (Step 2) done & verified** — conda env `bits_for_gaps` created from `environment.yml`; ML stack confirmed (gpflow 2.9.2 / TF 2.16.2 / TFP 0.24.0 / Py 3.9.23). Julia 1.12.6 + Clapeyron auto-installed via juliapkg. **Fix: must set `PYTHON_JULIACALL_HANDLE_SIGNALS=yes`** or juliacall SIGBUSes on macOS (exit 138). Wilson γ(PrOH/water, 350 K, x=0.5) = [1.469, 1.746]. Full recipe in the `environment-setup-recipe` memory. -- **Reproduce-from-archive (Step 3) substantial** — all done non-destructively into `results_remaked/`: - - **Fig 10** (HMC traces): regenerated `trace_all_15.png` via `plot_res.py -m all_traces`; R̂/ESS already matched paper numerically. - - **Fig 9** (McCabe-Thiele): both panels regenerated; column design matches paper exactly (xW=0.01, F=100, xF=0.10, R=1.0, xD=0.43, 4 stages). - - **GP-predict path**: archived `gp_model_15.pkl` unpickles (GPR + `AnisotropicSE`, 24 pts) and `predict_f` works. -- **Bug fixed** — `equilibrium.water_proh_eqm_julia` pointed at the missing `results/less_x/phase_diagram_data_saftgammamie`; repointed to the archived `gt_Wilson_data` (see NOTE in `equilibrium.py`). -- **Artifacts written** — `environment.yml` (repo root), `results_remaked/` (regenerated figs), scratchpad repro scripts. -- **Not yet done** — Fig 5 parity/error full regen (path validated via GP smoke test); Figs 2/3/4/6/7/11/12 regen; fresh-repo bootstrap (Phase 0); test harness (Phase 2). - ---- - -## 1. What the algorithm is (the thing we generalize) - -A sequential experimental-design loop over a **hierarchical Gaussian-process surrogate**: - -1. **Initialize** — space-filling design (Latin hypercube) over the input space; GP prior on the black-box output; priors on GP hyperparameters θ. -2. **Calibrate** — evaluate the black-box `f(x)` at design points; run **HMC** (TFP) to sample the hyperparameter posterior `p(θ|y)`; propagate a subset of samples through the GP → a **Gaussian-mixture predictive posterior**. -3. **Evaluate** — score the surrogate; pick the next `x` by **maximizing the predictive differential entropy** (Taylor approximation, or the closed-form lower bound). Repeat until a stopping criterion. - -The entropy math (GMM density, 1st/2nd-order Taylor entropy, closed-form lower bound) is **pure numpy/scipy and already generic**. The rest of the loop is generic *in intent* but **hardcoded to 2 inputs and 3 hyperparameters in the code**. - ---- - -## 2. Current code: algorithm vs. example (audit result) - -| Layer | Files (in `entropy_driven_hms/`) | Disposition | -|---|---|---| -| **CORE — sequential design engine** | `driver_new.py` (class `adaptiveEntropy`) | Generalize & decompose → package | -| **CORE — entropy math** | `fxns/max_ent_design.py` | Move ~verbatim; already generic (numpy) | -| **CORE — GP kernel / mean** | `fxns/my_kermel_fxn.py` (`AnisotropicSE`), `fxns/my_mean_fxn.py` (`FixedInverseMean`) | Generalize to N-D → package | -| **CORE — small utils** | `fxns/util.py`, `fxns/plot_settings.py` | Move; tidy | -| **EXAMPLE — thermo / data gen** | `proh_water_class.py`, `fxns/calculate_activities.jl`, `equilibrium.py` | → `examples/vle_distillation/` | -| **EXAMPLE — physics** | `gibbs_duhem.py`, `new_phase_diagram.py`, `distillation_model.py`, `solve_distillation.py` | → `examples/vle_distillation/` | -| **EXAMPLE — validation plots** | `train_test_split_proh.py` (Fig 5) | → `paper/` figure scripts | -| **PLOTTING — paper figure library** | `fxns/mcmc_plotter.py` (847 lines), `fxns/plot_res.py` (CLI) | Split into `paper/figures/`; extract generic bits | -| **VALIDATION — entropy approx** | `huber_et_al.py` (5D-mixture `h_vs_c` test) | → a **unit test** for `entropy.py` | -| **LEGACY — drop** | `driver.py` (superseded), `phase_diagram.py` (superseded, iter-10/`less_x`), `old/**`, `osu_presentation/**`, `paper_writing/**`, `power_point/**`, `fxns/*.jl` plotters (still EtOH-labeled) | Do **not** migrate | -| **DEAD code** | `max_ent_design.entropy_lower_bound` (never called), commented `second_order_entropy` block, `gibbs_duhem.load_ground_truth` (undefined `iters`) | Delete or fix during migration | - -**Critical coupling to sever:** -- `driver_new.py`, `proh_water_class.py` call `jl.include("fxns/calculate_activities.jl")` **at module import** → the core must never import Julia. The black-box `f(x)` is *injected* by the caller (this pattern already exists via `fwd_model`/`fwd_model_args`). -- `driver_new.py` starts with `from proh_water_class import PrOHwater` — core importing the example. Remove. -- Hardcoded 2-D everywhere: `gp_predict_2D`, `gen_entropy_surface_data_2D`, `optimize_2D`; `run_mcmc` indexes `trainable_parameters[2],[0],[1]` (exactly 3 hyperparameters); mixture/entropy code assigns kernel params **by name** (`std_dev`, `lengthscale_1`, `lengthscale_2`). -- `run_model` has `i += 50` — a resume offset from the manuscript-revision run. Remove; replace with real checkpoint/resume. -- State lives on disk (ad-hoc `np.savetxt`/`pickle` under `results/{exp}/`), not in memory. Make state in-memory objects with *optional* checkpointing. - ---- - -## 3. Reproducibility status (verified this session) - -- ✅ Archived `results/less_x_new_manuscript_revisions/` **iteration 15 = the published run**. Its `rhat_value_15.txt` / `ess_value_15.txt` match **Figure 10 exactly**: R̂ = 1.005 / 1.007 / 1.009, ESS = 1468.3 / 2428.1 / 653.1. -- ✅ `phase_diagram_15` (feeds Figs 8–9) and `gt_Wilson_data` are present on disk. -- ✅ Final results are archived — we can **regenerate figures from archived data without re-running MCMC**. -- ⚠️ **Gap:** `equilibrium.water_proh_eqm_julia` reads `results/less_x/phase_diagram_data_saftgammamie` — **absent from disk and git history**. The Fig 9 ground-truth ("Wilson") panel as-coded will fail. Fix: repoint to `gt_Wilson_data`, or regenerate via the Julia Wilson path. (The `_julia`/`saftgammamie` naming is a leftover from the earlier SAFT-γ-Mie system.) -- ⚠️ **Clapeyron.jl is unpinned** (no `juliapkg.json`) → activity-coefficient ground truth could drift. Pin it in the new repo. -- ⚠️ **Env is fragile & old:** Python 3.9 / TF 2.16.2 / GPflow 2.9.2 / TFP 0.24.0, macOS-arm64. **`juliacall` MUST import before TensorFlow/GPflow** or you get a bus-error segfault (fixed on `codex-refactor`). `codex-refactor` also has a working `environment.yml` — reuse it. -- ✅ **Fresh repo is correct:** 104 commits, ~99.8% of all file paths ever added are under `results/` (2.45 GB pack), no tags, no packaging config ever. Reference commits: `e17c818` (journal submission), `db3a2de` (pre-revision). - ---- - -## 4. Proposed package architecture (`src` layout, fresh repo) - -``` -bits_for_gaps/ -├── pyproject.toml # hatchling; core deps only (no Julia) -├── README.md LICENSE CHANGELOG.md -├── environment.yml # dev/repro env (from codex-refactor + fixes) -├── juliapkg.json # pin Clapeyron.jl (examples extra only) -├── .github/workflows/ci.yml # unit+integration on pure-Python core -├── docs/ # Sphinx + MyST → ReadTheDocs -├── src/bits_for_gaps/ -│ ├── __init__.py # public API surface -│ ├── sampler.py # BitsForGaps: the sequential-design loop (was adaptiveEntropy) -│ ├── gp.py # build GPR, hierarchical priors, run HMC, posterior samples -│ ├── mixture.py # GMM predictive posterior from θ-samples (generic param assignment) -│ ├── entropy.py # max_ent_design.py: GMM density, Taylor 1st/2nd, lower bound -│ ├── acquisition.py # entropy objective + multistart optimize (N-D) -│ ├── kernels.py # AnisotropicSE (N-D), extensible -│ ├── means.py # mean functions -│ ├── design.py # LHS / full-factorial space-filling (N-D) -│ ├── transforms.py # per-dim fwd/bkwd input & output transforms (class) -│ ├── diagnostics.py # R-hat, ESS helpers -│ ├── state.py # in-memory state + optional checkpoint/resume -│ └── plotting.py # generic, problem-agnostic plot helpers -├── examples/ -│ ├── synthetic/ # pure-Python demos (NO Julia) — run in CI & docs -│ │ ├── branin_1d.py # 1-D toy: exercises N-D generality + fast tests -│ │ └── synthetic_3d.py # 3-D toy: exercises >2 inputs -│ └── vle_distillation/ # the paper case study (needs Julia/Clapeyron) -│ ├── activity_model.py # was proh_water_class.py (generalized thermo wrapper) -│ ├── calculate_activities.jl -│ ├── gibbs_duhem.py phase_diagram.py distillation.py equilibrium.py -│ └── run_case_study.py -├── paper/ # reproduce published figures -│ ├── reproduce.py Makefile -│ ├── figures/ # mcmc_plotter split into per-figure scripts (5,8,9,10, etc.) -│ ├── reference/ # small scalar targets (R̂, ESS, MAP, stage table) + tolerances -│ └── DATA.md # pointer to the archived 2.5 GB results in the private old repo -└── tests/ - ├── unit/ # entropy math, transforms, kernels, GD integral, Antoine, design - ├── integration/ # tiny end-to-end seeded run (few samples, 1 iter) - └── regression/ # reference-file checks vs paper metrics (Julia-gated, slow job) -``` - -**Key API idea (target):** -```python -from bits_for_gaps import BitsForGaps, AnisotropicSE, InputTransform - -bfg = BitsForGaps( - black_box=my_fx, # callable: x (n,d) -> y (n,) - bounds=[(lo, hi), ...], # any dimension d - kernel=AnisotropicSE(ndim=d), # priors carried by the kernel Parameters - input_transforms=InputTransform([...]), # per-dim fwd/bkwd, optional - output_transform=..., # optional - likelihood_variance=0.1, -) -bfg.mcmc = MCMCConfig(n_samples=5000, n_chains=4, step_size=0.05, ...) -result = bfg.run(n_iterations=30, seed=10) # returns history: designs, θ-traces, entropy field, diagnostics -``` -Core has **zero Julia dependency**; `pip install bits_for_gaps` pulls GPflow/TF/numpy/scipy only. The VLE example's Julia backend is an optional extra: `pip install "bits_for_gaps[vle]"` + documented Julia/Clapeyron setup. - ---- - -## 5. Execution phases (maps to your 8 steps, reordered so tests precede the refactor) - -### Phase 0 — Decisions + fresh-repo bootstrap *(this Opus session)* -- Lock the open decisions in §7. -- Create the new repo skeleton (src layout, `pyproject.toml`, `.gitignore`, `environment.yml`, `LICENSE`, empty `tests/`). No history from the old repo. -- Carry over from `codex-refactor`: `environment.yml`, `.gitignore`, the juliacall-before-TF import fix, and `code_plan.md`/`reactor_log.md` as reference. - -### Phase 1 — Environment + reproduce-from-archive *(Opus; this is "the tricky part")* -- Build the conda env; verify `gpflow/tf/tfp` and `juliacall` import (correct order). Install & **pin Clapeyron.jl** (`juliapkg.json`). -- Regenerate Figs **5, 8, 9, 10 from archived `less_x_new_manuscript_revisions` data** (no MCMC re-run) and diff against the paper PDF. Resolve the missing `phase_diagram_data_saftgammamie` gap (repoint to `gt_Wilson_data`). -- Deliverable: a working env + a short "repro report" confirming archived == published. - -### Phase 2 — Characterization / regression harness *(BEFORE any refactor — your stated priority)* -- **Reference scalars** from archived iter-15: R̂, ESS, MAP hyperparameters, Fig 9 stage-composition table, Fig 5 RMSE/MAE. Store in `paper/reference/` with tolerances. -- **Unit tests on current behavior:** port `huber_et_al.py`'s 5D-mixture entropy curve to a pytest with pinned expected values; add tests for Gibbs-Duhem integral, Antoine `pvap`, transforms, LHS design determinism (seeded). -- **Seeded integration test:** a minimal end-to-end run (e.g. 100 MCMC samples, 1 iteration) whose outputs are stable. -- Everything green against the *unrefactored* code = the safety net. - -### Phase 3 — Package scaffold + verbatim move *(Sonnet)* -- Move CORE files into `src/bits_for_gaps/` with minimal edits: fix imports, remove module-load Julia, remove `from proh_water_class import`. Keep 2-D behavior identical. Get `import bits_for_gaps` working and the Phase-2 tests passing through the package. -- Stand up CI (unit + integration on pure-Python core; Julia/regression as a separate gated job). - -### Phase 4 — Incremental modular refactor of the core *(Sonnet)* -- Decompose `adaptiveEntropy` into `sampler` / `gp` / `mixture` / `acquisition` / `entropy` / `transforms` / `design` / `diagnostics` / `state`. -- Remove `i += 50`; replace disk-as-state with in-memory state + optional checkpointing. -- Keep regression tests green at every step. - -### Phase 5 — Generalize to N-D *(Sonnet, distinct milestone)* -- Remove 2-D/3-hyperparameter hardcoding: generic kernel-parameter introspection in HMC + mixture; N-D acquisition grid/optimizer; N-D design. -- New tests: 1-D and 3-D **synthetic** problems (no Julia) proving generality. This is where the "already general" claim becomes true in code. - -### Phase 6 — Port the VLE example onto the clean API *(Sonnet)* -- Rewrite `examples/vle_distillation/` to consume the public API (inject the Julia activity `f(x)`), fix hardcoded paths/`iters`, fix the `equilibrium.py` cross-experiment reference. - -### Phase 7 — Reproduce ALL paper figures via the new API *(Sonnet + Opus check)* -- `paper/reproduce.py` regenerates Figs 5, 8, 9, 10 (and 2, 3, 4, 6, 7) through the package + examples. Diff against reference scalars and the PDF. - -### Phase 8 — Documentation *(Sonnet)* -- Sphinx + MyST: install/env (incl. the juliacall gotcha), quickstart on the 1-D synthetic, API autodoc, a "reproduce the paper" guide, theory notes linking to the paper. Wire ReadTheDocs (`.readthedocs.yaml`); build the pure-Python parts without Julia. - -### Phase 9 — Reproduce ALL paper results, including the stochastic loop ✅ DONE (+ Phase 9b) -- Ran the full adaptive loop from scratch (`paper/full_reproduction.py`, ~26 min): seeded parts (HMC posterior, R̂/ESS, hyperparam posterior, entropy decay) reproduce the paper to **7–8 sig figs**; only the non-seedable `predict_f_samples` path (test-RMSE, surrogate curve) drifts. -- **Phase 9b correction:** Phase 9 initially reported the fully-adaptive surrogate's McCabe-Thiele column as non-converging and attributed it to entropy-driven design — that was a **bug** (shared-mutable-state in `full_reproduction.py`: the test-RMSE loop mutated `GPmodel.kernel` in place before the phase diagram reused it), **not** a scientific finding. Fixed (example layer only); the adaptive surrogate's column now converges and tracks Wilson within 0.03. See `paper/PHASE9B_INVESTIGATION.md`; retraction in `paper/REPRODUCTION.md`. - -### Phase 9c — Robustness hardening (make it a reliable package) *(Sonnet + Opus check)* -- Make BITS for GAPS as robust as possible **without changing numerical behavior** — the pre-Phase-4 baseline (`synthetic_baseline.json`, atol 1e-10) + all reference regressions are the safety net and stay green. First sanctioned core change since Phase 4. -- Targets: (1) fix the `mixture.sample_gp_posterior_mixture` in-place kernel-mutation footgun (save/restore state) — the Phase-9b bug; (2) public-API input validation with clear errors (bounds vs kernel ndim, lo= 81 +# (see CHANGELOG 0.1.1). Do this with CURRENT setuptools, so a missing dependency +# bound shows up here rather than in a user's traceback. +python -c "import bits_for_gaps as b; k = b.AnisotropicSE(); print('TF-backed OK, ndim =', k.ndim)" +python -c "import setuptools; print('setuptools resolved to', setuptools.__version__)" ``` -- All pinned runtime dependencies resolved and installed cleanly from the wheel's - metadata alone (no `[dev]`/`[docs]`/`[vle]` extras). -- `import bits_for_gaps` alone (before touching any lazily-imported attribute) left - both `juliacall` and `tensorflow` out of `sys.modules` -- confirmed. -- `bits_for_gaps.__version__ == "0.1.0"`. -- Smoke test exercised the public API without Julia: `AnisotropicSE()` construction, - `latin_hypercube_design(...)`, `entropy.second_order_entropy(...)` (matched the - analytic single-Gaussian entropy `0.5*log(2*pi*e) ≈ 1.41894`), and `BitsForGaps(...)` - construction using the exact kwargs from `docs/quickstart.md`'s snippet - (`black_box=`, `bounds=`, `kernel=`, `likelihood_variance=`) -- all matched the - installed API with no changes needed to the quickstart doc. - -## Maintainer-only steps (NOT performed by this phase — do these in order) - -**0. (Recommended, optional) Pin `pypa/gh-action-pypi-publish` to a commit SHA.** -`.github/workflows/publish.yml` currently references it via PyPA's own recommended -floating tag (`@release/v1`) -- this session had no network access to look up and -verify a real commit SHA, and didn't want to fabricate one. For a fully immutable -pin, look up the latest release at -https://github.com/pypa/gh-action-pypi-publish/releases and replace `@release/v1` -with `@ # v1.x.y` in both jobs. - -**1. Configure trusted publishing (before the first upload).** - -On both [pypi.org](https://pypi.org) and [test.pypi.org](https://test.pypi.org), -under the `bits_for_gaps` project's (or, for the very first upload, the *pending* -publisher form under your account, since the project doesn't exist yet) Publishing -settings, add a trusted publisher: - -- Owner: `dowlinglab` -- Repository name: `bits_for_gaps` -- Workflow filename: `publish.yml` -- Environment name: `release` - -This must be done once per index (PyPI and TestPyPI separately) before -`.github/workflows/publish.yml`'s OIDC-based `pypa/gh-action-pypi-publish` step can -authenticate. No API tokens are used or stored anywhere in this repo. - -**2. Dry run on TestPyPI.** - -- Go to the repo's Actions tab -> "Publish" workflow -> "Run workflow" (this is the - `workflow_dispatch` trigger in `publish.yml`) -> select the TestPyPI target. -- In a **fresh** virtual environment (not the dev conda env): - ```bash - python -m venv /tmp/bfg-testpypi-check && source /tmp/bfg-testpypi-check/bin/activate - pip install -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ bits_for_gaps - python -c "import bits_for_gaps, sys; assert 'juliacall' not in sys.modules and 'tensorflow' not in sys.modules; print(bits_for_gaps.__version__)" - # REQUIRED: also force-load a TensorFlow-backed module. The line above only touches - # the eagerly-imported pure modules, so the lazy __getattr__ never loads GPflow -- - # exactly how 0.1.0 shipped with every TF-backed module broken on setuptools >= 81 - # (see CHANGELOG 0.1.1). Do this in a fresh env that has CURRENT setuptools, so a - # missing dependency bound shows up here rather than in a user's traceback. - python -c "import bits_for_gaps as b; k = b.AnisotropicSE(); print('TF-backed OK, ndim =', k.ndim)" - python -c "import setuptools; print('setuptools resolved to', setuptools.__version__)" - ``` - (`--extra-index-url` is needed because TestPyPI doesn't mirror PyPI's dependencies -- - GPflow/TensorFlow/etc. would otherwise fail to resolve.) Run the same smoke test as - Phase 10 STEP 4 (below) against this install. -- Confirm the project page renders correctly on test.pypi.org (long description, - classifiers, project URLs). - -**3. Tag and publish to PyPI.** - -- Edit `CHANGELOG.md`: change `## [0.1.0] - YYYY-MM-DD` to the actual release date. -- Commit that change (on `main`, after this phase's PR merges). -- ```bash - git tag v0.1.0 - git push --tags - ``` - Pushing the `v0.1.0` tag triggers `publish.yml`'s tag-push job, which builds fresh - and publishes to PyPI via trusted publishing. -- Verify: the PyPI project page at `https://pypi.org/project/bits_for_gaps/`, then in - another fresh env, `pip install bits_for_gaps` and re-run the smoke test. - -**4. Activate ReadTheDocs.** - -- Import the project at readthedocs.org (steps already recorded in `HANDOFF.md`'s - Phase 8 section) -- this needs the maintainer's RTD account, not attempted by any - automated phase. -- Confirm the docs build succeeds there (locally verified throughout via - `sphinx-build -W docs docs/_build/html`) and that - `https://bits-for-gaps.readthedocs.io` (the URL already in `pyproject.toml`'s - `project.urls.Documentation` and `README.md`'s badge) resolves once built. - -## What Phase 10 did NOT do (by design) - -No PyPI/TestPyPI account configuration, no trusted-publisher registration, no upload, -no `git tag`, no RTD activation. `publish.yml` exists and is valid but has never been -triggered. `dist/` is gitignored and was never committed. +Also smoke-test the public API without Julia: `AnisotropicSE()` construction, +`latin_hypercube_design(...)`, `entropy.second_order_entropy(...)`, and +`BitsForGaps(...)` construction using the exact kwargs from `docs/quickstart.md`'s +snippet -- confirm they still match the installed API. Tear the environment down after. + +## Cutting a release + +1. Bump `__version__` in `src/bits_for_gaps/__init__.py`. +2. Move `CHANGELOG.md`'s `## [Unreleased]` content into a new `## [x.y.z] - YYYY-MM-DD` + section (today's date); leave a fresh empty `[Unreleased]` above it. +3. Commit both, on `main`. +4. Run the "Before every release" checks above against a local `python -m build`. +5. Optional: dry-run on TestPyPI first (Actions tab -> "Publish" workflow -> "Run + workflow" -> TestPyPI target), then in a fresh venv: + ```bash + pip install -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ bits_for_gaps + ``` + (`--extra-index-url` is needed because TestPyPI doesn't mirror PyPI's dependencies.) + Run the same smoke test as step 2 above. +6. Tag and push: + ```bash + git tag vx.y.z + git push --tags + ``` + This triggers `publish.yml`'s tag-push job, which builds fresh and publishes to PyPI + via trusted publishing. +7. Verify: `https://pypi.org/project/bits_for_gaps/` shows the new version, then in + another fresh env, `pip install bits_for_gaps` and re-run the smoke test. RTD picks + up the new tag automatically. + +## Optional hardening + +`.github/workflows/publish.yml` references `pypa/gh-action-pypi-publish` via PyPA's own +recommended floating tag (`@release/v1`). For a fully immutable pin, look up the latest +release at https://github.com/pypa/gh-action-pypi-publish/releases and replace +`@release/v1` with `@ # v1.x.y` in both jobs. diff --git a/docs/conf.py b/docs/conf.py index 0605fc9..6712bef 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,8 +1,8 @@ """Sphinx configuration for bits_for_gaps. Docs are repo-only (not shipped in the pip wheel, mirroring examples/ and paper/'s -policy -- see REFACTOR_PLAN.md §7.3) and built with Sphinx + MyST + furo, all pulled -in via the ``[docs]`` extra (``pip install -e ".[docs]"``). +policy) and built with Sphinx + MyST + furo, all pulled in via the ``[docs]`` extra +(``pip install -e ".[docs]"``). """ import os @@ -44,8 +44,8 @@ # --- autodoc ------------------------------------------------------------------ # -# The RTD/CI robustness decision (Phase 8): autodoc imports bits_for_gaps, which -# imports gpflow/tensorflow/tensorflow_probability for the TF-backed modules +# RTD/CI robustness: autodoc imports bits_for_gaps, which imports +# gpflow/tensorflow/tensorflow_probability for the TF-backed modules # (kernels, means, sampler, gp, mixture, acquisition, diagnostics). We install the # REAL frozen stack rather than mocking it -- pip always installs a package's base # `dependencies` alongside any extra, so `pip install -e ".[docs]"` (what diff --git a/docs/improvements_over_paper.md b/docs/improvements_over_paper.md index 8bee6a4..81c9037 100644 --- a/docs/improvements_over_paper.md +++ b/docs/improvements_over_paper.md @@ -1,42 +1,33 @@ # Improvements over the original paper code This package is a from-scratch port of the code behind Jones & Dowling (2026), not a -copy. Every phase of the port (see `HANDOFF.md` for the full history) fixed something, -made something more robust, or restructured something for maintainability. This page -consolidates those changes honestly and specifically -- real diffs, not marketing. - -```{note} -This is a pointer to real changes, each with more detail in the repository: -`HANDOFF.md` (phase-by-phase history), `paper/PHASE9B_INVESTIGATION.md` (the bug fix -below), and the module docstrings cited throughout (e.g. `mixture.py`, +copy. This page consolidates what changed and why -- real diffs, not marketing. More +detail on any item lives in the module docstrings cited throughout (e.g. `mixture.py`, `acquisition.py`, `sampler.py`). -``` ## Bugs found and fixed -**A missing data path, silently wrong until traced down (pre-Phase-1).** The original -`equilibrium.py`'s `water_proh_eqm_julia` read from a file path -(`results/less_x/phase_diagram_data_saftgammamie`) that didn't exist on disk or in the -old repo's git history. Repointed to the archived Wilson ground truth -(`gt_Wilson_data`) that the rest of the pipeline already produces and uses. +**A missing data path, silently wrong until traced down.** The original +`equilibrium.py`'s `water_proh_eqm_julia` read from a file path that didn't exist on +disk or in the original code's git history. Repointed to the archived Wilson ground +truth (`gt_Wilson_data`) that the rest of the pipeline already produces and uses. -**A shared-mutable-state bug that produced a spurious "finding" (Phase 9b).** A -validation script (`paper/full_reproduction.py`) computed a test-RMSE metric via +**A shared-mutable-state bug that produced a spurious "finding."** A validation script +(`paper/full_reproduction.py`) computed a test-RMSE metric via `mixture.sample_gp_posterior_mixture` -- which reassigns a GP's kernel hyperparameters once per posterior draw, by design -- then reused that *same* GP object afterward to build a McCabe-Thiele surrogate column. The column didn't converge, and the first write-up of this attributed that to a property of entropy-driven acquisition design. That attribution was wrong: reconstructing the failure from the checkpointed model reproduced it to 4 decimal places, and the real cause was the leftover mutated kernel -state, not the design method. Full analysis, including which of four hypotheses held -and which didn't, in `paper/PHASE9B_INVESTIGATION.md`. Fixed by reordering the script -and adding a hyperparameter-posterior-averaged surrogate construction; the same -underlying footgun (see below) was then hardened away at its source. +state, not the design method. Fixed by reordering the script and adding a +hyperparameter-posterior-averaged surrogate construction; the same underlying footgun +(see below) was then hardened away at its source. -## Hardening (Phase 9c) -- the same class of bugs made structurally harder to hit +## Hardening -- the same class of bugs made structurally harder to hit -The Phase 9b bug above was a symptom of a real footgun in the library itself, not just -a one-off script mistake: +The bug above was a symptom of a real footgun in the library itself, not just a +one-off script mistake: - **State-mutation footgun, fixed at the source.** `mixture.sample_gp_posterior_mixture` and `acquisition.entropy_objective` both reassign a GP's kernel hyperparameters in a @@ -44,10 +35,10 @@ a one-off script mistake: not any meaningful state. Since `sampler.py`'s `run()` calls these on the *same* `GPmodel` object it then stores in every `IterationRecord` (and optionally checkpoints to disk), **every run's returned model** used to carry this arbitrary - leftover state, not just the one script Phase 9b happened to hit. Both functions now + leftover state, not just the one script that first surfaced it. Both functions now save the kernel's hyperparameters before mutating and restore them in a `finally` -- - behavior-preserving for every value either function computes and returns, verified - against the pre-Phase-4 baseline (atol 1e-10) and all reference regressions. + this changes only the model's post-call state, not any value either function + computes and returns. - **Public-API input validation.** Constructing `BitsForGaps`/`adaptiveEntropy` with mismatched bounds/kernel dimensionality, `lo >= hi` bounds, or calling `run()` with non-positive HMC/acquisition config, mismatched `X_init`/`y_init` shapes, or an @@ -62,11 +53,10 @@ a one-off script mistake: failing. Now an explicit `ValueError` with a specific message. - **Fragile `fsolve` convergence, given a fallback.** The McCabe-Thiele column solver (`examples/vle_distillation/distillation.py`) has no bounds and used a single fixed - initial guess -- exactly the kind of solver fragility Phase 9b's investigation - surfaced (a smooth, well-behaved equilibrium curve can still fail to converge from an - unlucky initial guess). `solve_column` now retries a few generic alternate initial - guesses if the default doesn't converge, before giving up. The primary attempt is - byte-for-byte unchanged, so nothing that already converges is affected. + initial guess -- a smooth, well-behaved equilibrium curve can still fail to converge + from an unlucky initial guess. `solve_column` now retries a few generic alternate + initial guesses if the default doesn't converge, before giving up. The primary + attempt is byte-for-byte unchanged, so nothing that already converges is affected. - **Opt-in reproducibility for the one documented non-reproducible step.** `GPmodel.predict_f_samples` (used by `sample_gp_posterior_mixture`/`predict_grid_2D`) draws from TensorFlow's ambient RNG, not NumPy's -- confirmed non-reproducible even @@ -74,28 +64,18 @@ a one-off script mistake: an optional `tf_seed`; passing one makes that call's draws reproducible. Default (`None`) leaves the original, documented behavior unchanged. - **A clear error instead of a cryptic gpflow/TensorFlow traceback for an - unassignable hyperparameter value** (Phase 9d). `kernels.assign_hyperparameters` is - called deep inside `mixture.py`/`acquisition.py`'s hot loops to replay one + unassignable hyperparameter value.** `kernels.assign_hyperparameters` is called + deep inside `mixture.py`/`acquisition.py`'s hot loops to replay one posterior/mixture-component sample at a time. An extreme outlier sample -- most plausibly from `lengthscale_2`, deliberately left unconstrained (no positivity bijector, a real feature of the paper's kernel) so nothing bounds how far an HMC leapfrog step can push it -- can round-trip through a bijector's inverse to a non-finite unconstrained value, which gpflow's own `Parameter.assign` rejects with a low-level `InvalidArgumentError` (`Tensor had NaN/Inf values [Op:CheckNumerics]`) - that doesn't say *which* value or parameter caused it. This is the exact error - class hit mid-investigation while tracing the Phase 9b bug (in that instance from - an unrelated script mistake, not a genuine posterior outlier -- but the underlying - gpflow failure mode is real). Re-raised as a `ValueError` naming the parameter and - value; behavior-preserving for every value that was already assignable (every value - seen across this codebase's tests, reference regressions, and the from-scratch - stochastic reproduction runs). - -46 new tests were added alongside this hardening (unit tests asserting kernel state is -identical before/after, including on the error path; validation-error-path tests; an -integration test reproducing the exact Phase 9b scenario end-to-end) -- all in addition -to, not replacing, the existing regression suite, which stayed green throughout. + that doesn't say *which* value or parameter caused it. Re-raised as a `ValueError` + naming the parameter and value. -## Faithfulness (Phase 9d) -- using more of what the paper actually derived +## Faithfulness -- using more of what the paper actually derived - **The paper's closed-form entropy lower bound is now a usable acquisition objective.** The paper derives *two* entropy estimators for the hierarchical GP @@ -103,13 +83,13 @@ to, not replacing, the existing regression suite, which stayed green throughout. ({func}`bits_for_gaps.entropy.second_order_entropy`, Huber et al. 2008) that actually drove acquisition in the paper, and a closed-form cross-overlap lower bound ({func}`bits_for_gaps.entropy.entropy_lower_bound`, paper Theorem/SI-2). The - lower bound was implemented and unit-tested since Phase 2, but nothing in the - sequential-design loop could ever call it -- {func}`~bits_for_gaps.acquisition.entropy_objective` - had the Taylor estimator hardcoded. It's now selectable via - `objective="taylor"|"lower_bound"`, threaded through - {func}`~bits_for_gaps.acquisition.optimize`/{func}`~bits_for_gaps.acquisition.entropy_surface_2D` - and exposed as `BitsForGaps.acquisitionObjective` (default `"taylor"` -- every - existing baseline/reference value is unaffected unless a caller explicitly opts into + lower bound was implemented and unit-tested, but nothing in the sequential-design + loop could call it -- {func}`~bits_for_gaps.acquisition.entropy_objective` had the + Taylor estimator hardcoded. It's now selectable via `objective="taylor"|"lower_bound"`, + threaded through {func}`~bits_for_gaps.acquisition.optimize`/ + {func}`~bits_for_gaps.acquisition.entropy_surface_2D` and exposed as + `BitsForGaps.acquisitionObjective` (default `"taylor"` -- every existing + baseline/reference value is unaffected unless a caller explicitly opts into `"lower_bound"`). - **The entropy estimators are now validated against the quantity they approximate, not just a captured historical value.** The existing regression test pins @@ -125,46 +105,42 @@ to, not replacing, the existing regression suite, which stayed green throughout. just to see the method run. `examples/synthetic/run_example.py` is a new, small, actually-runnable script (`python examples/synthetic/run_example.py`, no Julia, well under a minute) demonstrating the same sequential-design loop on a smooth - closed-form 2-D function; `docs/quickstart.md` now points to it instead of test - files. + closed-form 2-D function; {doc}`quickstart` points to it. ## Architecture -- **Decomposed into focused, independently-testable modules** (Phase 4). The original - `driver_new.py` was a single ~450-line class mixing GP construction, HMC, entropy - math, acquisition optimization, and disk I/O. It's now +- **Decomposed into focused, independently-testable modules.** The original code was + a single ~450-line class mixing GP construction, HMC, entropy math, acquisition + optimization, and disk I/O. It's now {mod}`bits_for_gaps.gp`/{mod}`bits_for_gaps.mixture`/{mod}`bits_for_gaps.acquisition`/ {mod}`bits_for_gaps.entropy`/{mod}`bits_for_gaps.transforms`/{mod}`bits_for_gaps.state`, each independently unit-tested, with `sampler.py`'s `adaptiveEntropy` reduced to a thin orchestrator over them. -- **In-memory state, checkpointing opt-in** (Phase 4). The original code used disk as - its state-passing mechanism between iterations (`np.savetxt`/`pickle` under +- **In-memory state, checkpointing opt-in.** The original code used disk as its + state-passing mechanism between iterations (`np.savetxt`/`pickle` under `results/{exp_name}/`, read back on the next call). `run()` now takes the initial design in memory and returns a `RunHistory` -- a full run executes with zero disk writes by default; per-iteration file output (mirroring the original layout) is available via an opt-in `checkpoint_dir` argument. -- **Generalized to N input dimensions** (Phase 5). The original kernel, acquisition, - and mixture code hardcoded the 2-D VLE case (e.g. indexing `trainable_parameters` by - position, a reversed 2-D-specific black-box calling convention). The acquisition - path an actual run depends on (`optimize`, was `optimize_2D`) is now dimension- - general; only the dense-grid 2-D-only visualization diagnostics (`entropy_surface_2D`, - `predict_grid_2D`) stay 2-D (a dense grid is exponential in dimension, and neither - feeds the acquisition) -- they raise a clear error for other dimensions rather than - silently misbehaving. -- **Julia is an opt-in, lazily-imported extra** ({mod}`bits_for_gaps` core). The - published run's activity-coefficient black box needs Julia + Clapeyron.jl; the - sequential-design algorithm itself does not. `import bits_for_gaps` is Julia-free - (verified: `pip install bits_for_gaps` pulls GPflow/TensorFlow/NumPy/SciPy only); - `juliacall` is imported lazily, only when the VLE example's activity model is - actually called, and only if you installed the `[vle]` extra. -- **Archive-free figure reproduction** (Phase 9). The original figures could only be - regenerated with author access to the private, 564 MB archived run. A curated ~16 MB - subset of exactly the files the figures read is committed to `paper/data/` - (provenance in `paper/data/README.md`), so `python paper/reproduce.py` reproduces - every figure from a fresh clone with no private-archive access. -- **A real regression/reference-file test suite** (Phase 2 onward). The original code had - no automated tests. This package pins a pre-refactor numerical baseline - (`tests/integration/data/synthetic_baseline.json`, atol 1e-10) plus reference scalars - extracted from the published run (`paper/reference/*`) that every subsequent phase -- - including this hardening pass -- must reproduce exactly, with tests gated behind - `@pytest.mark.vle` only where they genuinely need Julia or private data. +- **Generalized to N input dimensions.** The original kernel, acquisition, and mixture + code hardcoded the 2-D VLE case (e.g. indexing `trainable_parameters` by position, a + reversed 2-D-specific black-box calling convention). The acquisition path an actual + run depends on (`optimize`) is now dimension-general; only the dense-grid 2-D-only + visualization diagnostics (`entropy_surface_2D`, `predict_grid_2D`) stay 2-D (a dense + grid is exponential in dimension, and neither feeds the acquisition) -- they raise a + clear error for other dimensions rather than silently misbehaving. +- **Julia is an opt-in, lazily-imported extra.** The published run's + activity-coefficient black box needs Julia + Clapeyron.jl; the sequential-design + algorithm itself does not. `import bits_for_gaps` is Julia-free (verified: `pip + install bits_for_gaps` pulls GPflow/TensorFlow/NumPy/SciPy only); `juliacall` is + imported lazily, only when the VLE example's activity model is actually called, and + only if you installed the `[vle]` extra. +- **Self-contained figure reproduction.** Regenerating the original figures required + the full ~564 MB of run artifacts. Exactly the ~16 MB the figure scripts actually read + is committed to `paper/data/` (provenance in `paper/data/README.md`), so + `python paper/reproduce.py` reproduces every figure from a fresh clone. +- **A real regression/reference-file test suite.** The original code had no automated + tests. This package pins a numerical baseline (`tests/integration/data/ + synthetic_baseline.json`, atol 1e-10) plus reference scalars extracted from the + published run (`paper/reference/*`) that the library must reproduce exactly, with + tests gated behind `@pytest.mark.vle` only where they genuinely need Julia. diff --git a/docs/index.md b/docs/index.md index 6c8a40e..fd26cbd 100644 --- a/docs/index.md +++ b/docs/index.md @@ -17,13 +17,16 @@ K. D. Jones and A. W. Dowling, "BITS for GAPS: Bayesian Information-Theoretic Sampling for hierarchical GAussian Process Surrogates," *Computers & Chemical Engineering* **211** (2026) 109650. [doi:10.1016/j.compchemeng.2026.109650](https://doi.org/10.1016/j.compchemeng.2026.109650) + +The paper is also bundled in the repository — +[`paper/bits_for_gaps_paper.pdf`](https://github.com/dowlinglab/bits_for_gaps/blob/main/paper/bits_for_gaps_paper.pdf) +— redistributed under [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/), so you can +read the method alongside the code. Please cite the DOI. ``` ```{admonition} Status -:class: warning -Pre-1.0. The public API (`BitsForGaps`, `AnisotropicSE`, ...) is stable across the -refactor phases described in the project's `HANDOFF.md`, but has not yet had a -tagged release. +:class: note +0.x: released on PyPI, but the public API may still change before a 1.0 release. ``` ## Where to start diff --git a/docs/installation.md b/docs/installation.md index 5f4d16b..d63523b 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -6,23 +6,26 @@ optional VLE/distillation example, and even then only when you actually call int ## Core library +```bash +pip install bits_for_gaps +``` + ```{admonition} Frozen dependency stack :class: note The core pins an exact, verified-working stack rather than floating version ranges: Python 3.9, NumPy 1.26, SciPy 1.13, GPflow 2.9.2, TensorFlow 2.16.2, TensorFlow -Probability 0.24.0. This is a deliberate reproducibility choice (see -`REFACTOR_PLAN.md` §7 decision 6) -- GPflow's TensorFlow dependency makes casual -version bumps risky, so modernizing the stack is left as a separate, later effort. +Probability 0.24.0. This is a deliberate reproducibility choice -- GPflow's TensorFlow +dependency makes casual version bumps risky, so modernizing the stack is left as a +separate, later effort. ``` -Once published (see `HANDOFF.md` for current status -- this repo is pre-1.0): +Verify it imports (no Julia touched): ```bash -pip install bits_for_gaps +python -c "import bits_for_gaps; print(bits_for_gaps.__version__)" ``` -**From source** (the current way to get it, and the way to get `examples/` and -`paper/` -- see below): +## From source (for `examples/`, `paper/`, or development) ```bash git clone https://github.com/dowlinglab/bits_for_gaps @@ -32,30 +35,24 @@ conda activate bits_for_gaps pip install -e . ``` -Verify it imports (no Julia touched): - -```bash -python -c "import bits_for_gaps; print(bits_for_gaps.__version__)" -``` - ## Development install ```bash pip install -e ".[dev]" # adds pytest, pytest-cov -pytest -q # 193 passed, 2 deselected +pytest -q # 204 passed, 2 deselected ``` -The 2 deselected tests need the private archived published run and/or Julia -- see -{doc}`reproduce_paper` and `HANDOFF.md`. +The 2 deselected tests need Julia (`@pytest.mark.vle`) -- run them with +`pytest -m vle` after installing the `[vle]` extra below. See {doc}`reproduce_paper`. ## `examples/` and `paper/` are repo-only ```{important} The paper's worked example (`examples/vle_distillation/`) and figure-reproduction scripts (`paper/`) are **not part of the `bits_for_gaps` PyPI package** -- the wheel -ships only `src/bits_for_gaps/`, to keep the installed package lightweight -(`REFACTOR_PLAN.md` §7 decision 3). To use either, **clone the repository**; `pip -install bits_for_gaps` alone does not give you them. +ships only `src/bits_for_gaps/`, to keep the installed package lightweight. To use +either, **clone the repository**; `pip install bits_for_gaps` alone does not give you +them. ``` Once cloned, `examples/` and `paper/` are importable as regular Python packages diff --git a/docs/reproduce_paper.md b/docs/reproduce_paper.md index 730aa61..07e4ec3 100644 --- a/docs/reproduce_paper.md +++ b/docs/reproduce_paper.md @@ -1,58 +1,125 @@ -# Reproducing the paper's figures - -All 11 of the paper's figures (2 through 12) regenerate through `paper/reproduce.py` + -`paper/figures/` -- **repo-only** scripts (see {doc}`installation`), not part of the -pip package. - -```{note} -This page is a pointer, not a copy -- the authoritative, up-to-date reference is -[`paper/REPRODUCTION.md`](https://github.com/dowlinglab/bits_for_gaps/blob/main/paper/REPRODUCTION.md) -in the repository, which has the full figure -> script -> data -> reference-diff table, -known discrepancies, and the simplifications taken from the original plotting code. -Nothing here is executed as part of the docs build. -``` +# Reproducing the paper's results + +This page assumes you've read Jones & Dowling (2026) once and know the figures by number. +(The article is bundled in the repository at +[`paper/bits_for_gaps_paper.pdf`](https://github.com/dowlinglab/bits_for_gaps/blob/main/paper/bits_for_gaps_paper.pdf) +if you want it open alongside.) Everything you need is committed here — there are no +additional downloads. + +There are two different things you might mean by "reproduce the paper," and they have very +different costs: + +| | Regenerate the figures | Re-run the full adaptive loop | +|---|---|---| +| Script | `paper/reproduce.py` | `paper/full_reproduction.py` | +| What it does | Renders the paper's *already-collected* data (HMC traces, GP posterior draws, phase-diagram data, ...) through this package's plotting code | Actually runs the 15-iteration adaptive HMC + entropy-acquisition design loop from scratch | +| Cost | ~1-2 minutes | ~25-30 minutes on a laptop | +| Determinism | Deterministic (same inputs, same plots) | **Stochastic** -- will not reproduce the paper's exact numbers | +| Needs Julia? | Only Fig 8/9 | Yes, throughout | -## No private-archive access needed +Both are **repo-only** (`git clone` required; see {doc}`installation` -- neither is +part of the `pip install bits_for_gaps` package). -As of Phase 9, this needs **no private-archive access by default**: the data comes -from a curated, committed subset at `paper/data/` (~16 MB -- exactly the plot-input -files the figures read; see `paper/data/README.md`), not the private archive of -record. Clone the repo and run: +## Regenerate the figures (fast, deterministic) + +All 11 of the paper's figures (2 through 12) regenerate from data already committed at +`paper/data/` (~16 MB -- exactly the plot-input files the figures read; see +`paper/data/README.md`). No private-archive access, and no HMC sampling, is needed: ```bash export PYTHON_JULIACALL_HANDLE_SIGNALS=yes # macOS; only Fig 8/9 touch Julia python paper/reproduce.py ``` +**Measured cost:** ~80 seconds wall-clock on a laptop (M2 MacBook Pro) for all 11 +figures, including the 2 that call live Clapeyron (Fig 8/9). Order of a minute, not +tens of minutes -- this path does no MCMC. + Output goes to `results_remaked/` by default (gitignored -- nothing this produces is -committed). See `paper/REPRODUCTION.md` for regenerating a subset -(`--figures 5 8 9 10`) and for which figures are quantitatively pinned against -`paper/reference/*` (Figs 5, 8, 9, 10, 11) versus visually reproduced only (2, 3, 4, 6, -7, 12). +committed). Regenerate a subset with `--figures`, e.g. `python paper/reproduce.py +--figures 5 8 9 10`. See `paper/REPRODUCTION.md` for the full figure -> script -> +data -> reference-diff table, which figures are quantitatively pinned against +`paper/reference/*` (5, 8, 9, 10, 11) versus visually reproduced only (2, 3, 4, 6, 7, +12), and known discrepancies (e.g. the reference stage table's ~0.01 transcription +slop). -Point `--archive`/`$BFG_ARCHIVE_DIR` at the full private archive (see -[`paper/DATA.md`](https://github.com/dowlinglab/bits_for_gaps/blob/main/paper/DATA.md)) -only if you want figures at iterations beyond the curated subset (e.g. Fig 6/7 at an -iteration other than 1/15) -- that access is author-only, not needed for the default -path above. +Only Fig 8 and Fig 9 touch Julia (they recompute the Wilson ground-truth curve live via +Clapeyron.jl); the other 9 figures render from pure Python/NumPy and need neither Julia +nor the `[vle]` extra. If you only want those 9: `pip install -e ".[dev]"` is enough, +and you can skip the `PYTHON_JULIACALL_HANDLE_SIGNALS` export. -## The full stochastic loop, from scratch +Everything these figures read is committed in `paper/data/`, so nothing above requires +extra downloads. If you re-run the loop yourself (next section), point +`--archive`/`$BFG_ARCHIVE_DIR` at your run's output directory to plot your results +instead of the published ones. + +## Re-run the full adaptive loop from scratch (slow, stochastic) `paper/reproduce.py` renders what the paper's 15-iteration adaptive HMC loop already -produced -- it does not re-run that loop (stochastic, hours-long, and not what -reproducing a *figure* requires). `paper/full_reproduction.py` does run it, from a -fresh initial design against the live Clapeyron/Wilson black box, as a separate, -one-time validation exercise (not part of the regression suite, not required to -reproduce a figure) -- see `paper/REPRODUCTION.md`'s "Phase 9" section for the result -and `paper/PHASE9B_INVESTIGATION.md` for a bug found and fixed while validating it. +produced -- it does not re-run that loop. `paper/full_reproduction.py` does: it drives +`bits_for_gaps.sampler.BitsForGaps` through the same 15-iteration adaptive HMC + +entropy-acquisition design loop as the paper's published run (same bounds, seed, +transforms, kernel, and HMC config -- see the script's module docstring for the exact +constants), starting from a *fresh* Latin-hypercube design and a live Clapeyron/Wilson +black box: + +```bash +export PYTHON_JULIACALL_HANDLE_SIGNALS=yes +python paper/full_reproduction.py --out-dir results_remaked/full_reproduction +``` + +**Documented cost: ~25-30 minutes on a laptop** (15 outer iterations, each running a +4-chain x 5000-sample HMC fit plus live Clapeyron calls for the training data and a +full-grid posterior-predictive diagnostic). This is a one-time validation exercise, not +part of the test suite or the CI-gated regression checks -- there's no need to run it +to use the package, reproduce a figure, or verify your installation. + +```{caution} +**This run is stochastic and will NOT reproduce the paper's exact numbers.** HMC +sampling and the entropy-maximizing acquisition both involve randomness; a fresh run +explores a different (though statistically similar) sequence of design points and +posterior samples than the paper's published run did. +``` + +**What SHOULD match** (qualitatively, run to run and against the paper): +- HMC convergence: R-hat < 1.1 for all three kernel hyperparameters, with healthy + effective sample size relative to the ~20,000 raw HMC samples per iteration. +- The posterior orders the two lengthscales the same way the paper reports: the + temperature lengthscale is *larger* than the mole-fraction lengthscale (the GP trusts + nearby-temperature extrapolation more than nearby-composition extrapolation). +- Maximum entropy decays over successive iterations (Fig 4's shape) as the design fills + in the space and uncertainty drops. +- Test-set predictive error drops from iteration 1 to iteration 15 (though not + necessarily monotonically in between -- entropy-driven acquisition optimizes + information gain, not held-out error, at each step). +- The surrogate's phase diagram and McCabe-Thiele stage table agree with the Wilson + ground truth within a few 0.01 in mole fraction, the same way the paper's surrogate + does. + +**What will legitimately differ** from the paper and between your own runs: +- The exact sequence of sampled design points and their coordinates. +- Exact posterior samples, R-hat/ESS values, and entropy values (though the same order + of magnitude and trend). +- The exact test-RMSE trajectory -- `GPmodel.predict_f_samples` (used for + posterior-predictive draws, not for the HMC posterior itself) draws from + TensorFlow's ambient RNG, which this package cannot seed bitwise (see + `bits_for_gaps.mixture`'s module docstring); this was already true of the paper's + own original code. +- Small (sub-0.05 mole fraction) shifts in the surrogate's phase diagram / stage table + relative to the paper's exact reported numbers. + +A completed run's numbers, compared directly against the paper's, are recorded in +`paper/REPRODUCTION.md`'s "Re-running the full adaptive loop from scratch" section -- +useful as a sanity check for what "qualitatively similar" looks like in practice, if you +run it yourself and want something to compare against. ## The regression tests Regression tests that check figure/column reproduction against `paper/reference/*` (`tests/regression/test_paper_figures.py`, `tests/regression/test_mccabe_thiele.py`) -mostly run in the **default** `pytest -q` suite now -- only the tests that recompute -something via live Clapeyron calls (needing Julia) stay behind -`@pytest.mark.vle` (run with `pytest -m vle`). You can also read `paper/figures/*.py` -to see exactly how each figure is built, and run `examples/vle_distillation/ -run_case_study.py` ({doc}`vle_example`) for a small, fresh (non-archived) -demonstration of the same underlying pipeline. +mostly run in the **default** `pytest -q` suite -- only the tests that recompute +something via live Clapeyron calls (needing Julia) stay behind `@pytest.mark.vle` (run +with `pytest -m vle`). You can also read `paper/figures/*.py` to see exactly how each +figure is built, and run `examples/vle_distillation/run_case_study.py` +({doc}`vle_example`) for a small, fast, fresh (non-archived) demonstration of the same +underlying pipeline at a shorter iteration count. diff --git a/docs/theory.md b/docs/theory.md index eea7d95..c97f2fd 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -10,6 +10,11 @@ the paper and its supplementary information (SI): > Engineering* **211** (2026) 109650. > [doi:10.1016/j.compchemeng.2026.109650](https://doi.org/10.1016/j.compchemeng.2026.109650) +The article itself is bundled in the repository at +[`paper/bits_for_gaps_paper.pdf`](https://github.com/dowlinglab/bits_for_gaps/blob/main/paper/bits_for_gaps_paper.pdf) +(redistributed under [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/)), so the +equation numbers cited below can be checked directly against it. + ## The loop BITS for GAPS is a sequential experimental-design loop over a **hierarchical** @@ -177,9 +182,7 @@ has no closed form, so `bits_for_gaps.entropy` provides two estimators: \exp\!\left(-\frac12 \frac{(\mu_s - \mu_{s'})^2}{\sigma_s^2 + \sigma_{s'}^2}\right)$$ Selectable as the acquisition objective via `objective="lower_bound"` (or, on - `BitsForGaps`/`adaptiveEntropy`, `.acquisitionObjective = "lower_bound"`) -- - implemented since Phase 2 but only wired up as a usable acquisition choice in - Phase 9d. + `BitsForGaps`/`adaptiveEntropy`, `.acquisitionObjective = "lower_bound"`). Both reduce to the exact differential entropy of a single Gaussian in the degenerate one-component case (see `tests/unit/test_entropy.py`, which checks this diff --git a/examples/synthetic/run_example.py b/examples/synthetic/run_example.py index 448d4ab..7cf3885 100644 --- a/examples/synthetic/run_example.py +++ b/examples/synthetic/run_example.py @@ -55,9 +55,9 @@ def run(n_init=N_INIT, n_iters=N_ITERS, seed=SEED): bfg.noGaussians = 8 bfg.noRestarts = 3 bfg.entropyMesh = [4, 4] - # Phase 9d: the paper's default acquisition objective is the 2nd-order Taylor - # entropy approximation. Uncomment to try the alternative closed-form lower bound - # instead (see docs/theory.md's "Entropy estimators" section): + # The paper's default acquisition objective is the 2nd-order Taylor entropy + # approximation. Uncomment to try the alternative closed-form lower bound instead + # (see docs/theory.md's "Entropy estimators" section): # bfg.acquisitionObjective = "lower_bound" print(f"Running {n_iters} adaptive design iterations from {n_init} initial points...") diff --git a/examples/vle_distillation/README.md b/examples/vle_distillation/README.md index bc49771..9d49e8c 100644 --- a/examples/vle_distillation/README.md +++ b/examples/vle_distillation/README.md @@ -6,9 +6,8 @@ Gibbs-Duhem correction, bubble-point/dew-point phase diagram, and a McCabe-Thiel distillation column solver. **This example is repo-only** -- it is *not* part of the `bits_for_gaps` package -distributed on PyPI (the pip wheel ships only `bits_for_gaps/*`; see -`REFACTOR_PLAN.md` §7.3). To run it you need a clone of this repository, not just -`pip install bits_for_gaps`. +distributed on PyPI (the pip wheel ships only `bits_for_gaps/*`). To run it you need a +clone of this repository, not just `pip install bits_for_gaps`. ## Setup @@ -68,9 +67,9 @@ both the Wilson ground truth and the GP surrogate. Runtime is a couple of minute (dominated by the adaptive design's HMC sampling, a handful of iterations by default -- see `N_ITERS` in the script; the paper's published run used 15). -This is a *demonstration* of the ported pipeline, not a reproduction of the paper's -exact published figures (that full reproduction, matching the paper's real 15-iteration -adaptive run, is Phase 7's job -- see `HANDOFF.md`). +This is a *demonstration* of the pipeline (5 adaptive iterations, for speed), not a +reproduction of the paper's exact published figures -- the paper's real run used 15 +iterations. For that, see `paper/full_reproduction.py` and `docs/reproduce_paper.md`. ## Sanity check diff --git a/examples/vle_distillation/__init__.py b/examples/vle_distillation/__init__.py index fb6e37a..3ff3815 100644 --- a/examples/vle_distillation/__init__.py +++ b/examples/vle_distillation/__init__.py @@ -1,5 +1,5 @@ -"""The paper's H2O-PrOH VLE / distillation case study, ported onto the public -``bits_for_gaps`` API (Phase 6). +"""The paper's H2O-PrOH VLE / distillation case study, built on the public +``bits_for_gaps`` API. Not shipped in the ``bits_for_gaps`` pip wheel -- repo-only, importable in dev/CI via ``tests/conftest.py``'s ``sys.path`` insert of the repo's ``examples/`` directory (see diff --git a/examples/vle_distillation/activity_model.py b/examples/vle_distillation/activity_model.py index e8f02b7..53d30d2 100644 --- a/examples/vle_distillation/activity_model.py +++ b/examples/vle_distillation/activity_model.py @@ -1,7 +1,6 @@ """Julia/Clapeyron activity-coefficient model for the H2O-PrOH VLE case study. -Ported from the paper code's ``proh_water_class.py`` + ``fxns/calculate_activities.jl`` -(Wilson activity-coefficient model via Clapeyron.jl). Requires +Wilson activity-coefficient model via Clapeyron.jl. Requires ``pip install "bits_for_gaps[vle]"`` plus a working Julia installation -- ``juliacall``/Clapeyron are imported LAZILY (only when a function here is actually called), so ``import vle_distillation.activity_model`` succeeds without Julia; only @@ -88,7 +87,7 @@ def black_box(z_proh, temperature): """BitsForGaps-compatible black box: called as ``FwdModel(*args, *xStar)``. ``xStar`` is in bounds order ``[z_PrOH, T]`` (see ``run_case_study.py``), so this - accepts ``(z_proh, temperature)`` positionally, matching Phase 5's natural- + accepts ``(z_proh, temperature)`` positionally, matching the sampler's natural- dimension-order calling convention for the injected black box. Returns only the PrOH activity coefficient (component 0): the GP surrogate models diff --git a/examples/vle_distillation/distillation.py b/examples/vle_distillation/distillation.py index aabda94..5ea2fd2 100644 --- a/examples/vle_distillation/distillation.py +++ b/examples/vle_distillation/distillation.py @@ -1,9 +1,8 @@ """McCabe-Thiele-consistent stage-by-stage distillation column solver. -Ported from the paper code's ``distillation_model.py`` + ``solve_distillation.py`` -(itself a Python port of a course MATLAB script, ``distillation_nonlinear_equations.m`` --- see the inline notes on the reboiler equation below, kept verbatim from the original -port). Solves the column's nonlinear mass-balance + equilibrium-stage system via +Traces to a course MATLAB script (``distillation_nonlinear_equations.m`` -- see the +inline notes on the reboiler equation below for the one place its conventions still +show through). Solves the column's nonlinear mass-balance + equilibrium-stage system via ``scipy.optimize.fsolve``, given an equilibrium function ``y = equil(x)`` (see ``equilibrium.py`` / ``phase_diagram.py``). @@ -13,7 +12,7 @@ R*D``, ``V_1 = L_0 + D``, ``L_n = V_{n+1} + W``, ``x_D = x_0``, ``x_W = x_n`` -- see :func:`_distillation_residuals`'s inline comments for exactly where each appears. -Variable vector layout (0-based), matching the original port exactly: +Variable vector layout (0-based): v = [L_0..L_n, V_1..V_{n+1} (stored at V[0]..V[n]), x_0..x_n, y_1..y_{n+1} (stored at y[0]..y[n]), D, R, W, F, xF, q] Stage 0 = condenser/distillate, stage n = reboiler/bottoms. @@ -82,10 +81,10 @@ def _distillation_residuals(v, n, feed_stage_idx, equil, fixed_idx, fixed_vals): # SI-4 condenser and reboiler closures. offset_eq = 4 * n + len(fixed_idx) f[offset_eq + 0] = y[0] - x[0] # total condenser: y_1 = x_0 (x_D = x_0) - # NOTE (kept from the original port): a partial reboiler is usually - # y_{n+1} = equil(x_W). The source MATLAB instead sets x_n = y_{n+1} (liquid - # leaving stage n equals the vapor leaving the reboiler); replicated as-is here - # for physics parity with the paper's published column design. + # NOTE: a partial reboiler is usually specified as y_{n+1} = equil(x_W). This + # solver instead sets x_n = y_{n+1} (liquid leaving stage n equals the vapor + # leaving the reboiler) -- a convention traced to the source MATLAB script, kept + # for parity with the paper's published column design. f[offset_eq + 1] = x[n] - y[n] f[offset_eq + 2] = L[0] - R * D # SI-4: L_0 = R * D f[offset_eq + 3] = V[0] - L[0] - D # SI-4: V_1 = L_0 + D @@ -144,8 +143,8 @@ def _resolve_fixed_indices(n_stages, var_names, var_values): def _try_solve_column(v0, n, feed_stage_idx, equil, fixed_idx, fixed_vals, num_stages): """One ``fsolve`` attempt from ``v0``; returns the same dict :func:`solve_column` - does. Factored out (Phase 9c) so :func:`solve_column` can retry from alternate - initial guesses without duplicating the residual/extraction/diagnostic logic. + does. Factored out so :func:`solve_column` can retry from alternate initial guesses + without duplicating the residual/extraction/diagnostic logic. """ def residual_func(v_solve): @@ -199,10 +198,10 @@ def residual_func(v_solve): } -# Phase 9c: generic (not physics-informed) alternate (x, y) initial-guess levels to -# retry with if the default (0.5, 0.5) guess below doesn't converge -- see -# solve_column. Deliberately generic rather than curve-specific, so this doesn't -# encode any assumption about which equilibrium curve is passed in. +# Generic (not physics-informed) alternate (x, y) initial-guess levels to retry with if +# the default (0.5, 0.5) guess below doesn't converge -- see solve_column. Deliberately +# generic rather than curve-specific, so this doesn't encode any assumption about which +# equilibrium curve is passed in. _RETRY_INITIAL_GUESS_LEVELS = [(0.3, 0.3), (0.7, 0.7), (0.2, 0.8)] @@ -231,16 +230,13 @@ def solve_column(n_stages, feed_stage, equil, var_names, var_values): Notes ----- - Phase 9c: if the default (0.5, 0.5) initial guess below doesn't converge, retries - from a few generic alternate initial guesses before giving up (this is exactly the - kind of solver fragility that produced a spurious non-convergence Phase 9b traced - to an unrelated bug -- see ``paper/PHASE9B_INVESTIGATION.md``; retrying here is - cheap, generic insurance against genuine cases of it, not a fix for that bug). - The primary attempt is untouched -- identical inputs/outputs to before this change - -- so nothing that already converges is affected; retries only run when the first - attempt's own ``converged`` flag is ``False``, and the first converging result - (default or a retry) is returned as-is, with a note appended to ``"warnings"`` if - a retry was needed. + If the default (0.5, 0.5) initial guess below doesn't converge, retries from a few + generic alternate initial guesses before giving up -- cheap insurance against + ``fsolve``'s sensitivity to the initial guess for a poorly-conditioned equilibrium + curve. The primary attempt is unaffected; retries only run when its own + ``converged`` flag is ``False``, and the first converging result (default or a + retry) is returned as-is, with a note appended to ``"warnings"`` if a retry was + needed. """ n = n_stages feed_stage_idx = feed_stage - 1 diff --git a/examples/vle_distillation/equilibrium.py b/examples/vle_distillation/equilibrium.py index 4e28376..f5118a8 100644 --- a/examples/vle_distillation/equilibrium.py +++ b/examples/vle_distillation/equilibrium.py @@ -1,9 +1,8 @@ """Wrap a VLE curve as an ``x_liquid -> y_vapor`` equilibrium function. -Ported from the paper code's ``equilibrium.py`` (``water_proh_eqm`` / -``water_proh_eqm_julia``), generalized to take the curve arrays directly (e.g. from -``phase_diagram.vle_curve``) instead of reading fixed archived filenames -- the -distillation solver (``distillation.py``) needs a plain ``x -> y`` callable. +Takes the curve arrays directly (e.g. from ``phase_diagram.vle_curve``) rather than +reading them from a file -- the distillation solver (``distillation.py``) needs a plain +``x -> y`` callable. """ import numpy as np diff --git a/examples/vle_distillation/gibbs_duhem.py b/examples/vle_distillation/gibbs_duhem.py index 4a5dc7a..8f3b0f1 100644 --- a/examples/vle_distillation/gibbs_duhem.py +++ b/examples/vle_distillation/gibbs_duhem.py @@ -1,6 +1,5 @@ """Gibbs-Duhem correction: recover gamma_H2O from a modeled gamma_PrOH curve. -Ported from the paper code's ``new_phase_diagram.py`` (``PhaseDiagram.gibbs_duhem_fast``). The GP surrogate used in this case study models only ``gamma_PrOH(z, T)`` -- the water coefficient is derived from it via the binary Gibbs-Duhem relation, not learned by a second GP output. For an isothermal, isobaric binary mixture (differential form, diff --git a/examples/vle_distillation/phase_diagram.py b/examples/vle_distillation/phase_diagram.py index 901aa3e..ed38767 100644 --- a/examples/vle_distillation/phase_diagram.py +++ b/examples/vle_distillation/phase_diagram.py @@ -1,7 +1,6 @@ """H2O-PrOH bubble-point / dew-point phase diagram (Geankoplis Ex. 11.4-1 system). -Ported from the paper code's ``new_phase_diagram.py`` (``PhaseDiagram``). Two -interchangeable activity-coefficient sources feed the same bubble-point/dew-point +Two interchangeable activity-coefficient sources feed the same bubble-point/dew-point physics: - **Ground truth (Wilson)**: :func:`wilson_gamma` calls Clapeyron.jl directly (via @@ -64,8 +63,8 @@ def surrogate_gamma( last state), not a posterior summary. See :func:`surrogate_gamma_averaged` for a hyperparameter-posterior-averaged alternative when a robust curve independent of which single HMC sample happens to be live matters more than raw speed (e.g. - feeding a McCabe-Thiele stage solver with a fixed initial guess -- see - ``paper/PHASE9B_INVESTIGATION.md``). Also: ``GPmodel.kernel`` is mutated in place + feeding a McCabe-Thiele stage solver with a fixed initial guess, which is sensitive + to which draw the curve came from). Also: ``GPmodel.kernel`` is mutated in place by anything that reassigns its hyperparameters (e.g. ``bits_for_gaps.mixture.sample_gp_posterior_mixture``, by design -- see its docstring), so callers that need this function's result to reflect a *specific*, @@ -100,10 +99,9 @@ def surrogate_gamma_averaged( ): """Hyperparameter-posterior-averaged surrogate activity coefficients. - Matches the paper's own construction (the old repo's ``new_phase_diagram.py``'s - ``PhaseDiagram.run``/``gibbs_duhem_fast``, and ``equilibrium.py``'s - ``water_proh_eqm``, which fed the paper's *actual* Fig 9 surrogate column): draw - ``n_draws`` independent samples from the HMC hyperparameter posterior (``trace``), + Matches the paper's own posterior-averaging construction, which fed its actual + Fig 9 surrogate column: draw ``n_draws`` independent samples from the HMC + hyperparameter posterior (``trace``), evaluate this GP's own deterministic conditional mean under each one (not ``predict_f_samples`` -- no ambient-RNG, TF-non-reproducibility involved), and average ``gamma_proh`` pointwise. This is a Monte Carlo estimate of @@ -121,8 +119,8 @@ def surrogate_gamma_averaged( HMC posterior samples (e.g. ``bits_for_gaps.state.IterationRecord.trace``), in ``GPmodel.kernel.hyperparameters``'s canonical order. n_draws : int - Number of posterior draws to average (50, matching the paper's - ``PhaseDiagram.n_draws``). + Number of posterior draws to average. The default of 50 matches the number of + draws behind the paper's Fig 8/9 surrogate curves. """ from bits_for_gaps.kernels import assign_hyperparameters diff --git a/examples/vle_distillation/run_case_study.py b/examples/vle_distillation/run_case_study.py index b8ef332..ddf0409 100644 --- a/examples/vle_distillation/run_case_study.py +++ b/examples/vle_distillation/run_case_study.py @@ -4,11 +4,10 @@ ``BitsForGaps.run`` (adaptive entropy-driven design) -> phase diagram + McCabe-Thiele distillation column, all on the public ``bits_for_gaps`` API. -Configuration (bounds, transforms, kernel, seed) matches the paper's published -``less_x_new_manuscript_revisions`` run (Jones & Dowling 2026) as closely as the -archived code allows to reconstruct -- see HANDOFF.md for the paper trail. This -script demonstrates the ported pipeline; it does not reproduce the paper's exact -15-iteration adaptive run bit-for-bit (that full reproduction is Phase 7's job). +Configuration (bounds, transforms, kernel, seed) matches the paper's published run +(Jones & Dowling 2026). This script demonstrates the pipeline with a short 5-iteration +run, for speed; it does not reproduce the paper's exact 15-iteration adaptive run -- +for that, see ``paper/full_reproduction.py``. Usage:: @@ -33,21 +32,21 @@ # Paper's exact 2-D VLE search space: liquid PrOH mole fraction, temperature [K]. BOUNDS = [(1e-6, 0.999), (350.0, 367.0)] -SEED = 10 # matches train_test_split_proh.py's `my_system.seed = 10` +SEED = 10 # matches the paper's published run configuration -# Paper's exact GP input/output transforms (new_phase_diagram.py's __main__ block): -# log(x + 0.1) keeps the mole-fraction lengthscale well-scaled near the dilute limit; -# min-max normalizing T to [0, 1] matches the kernel's O(1) lengthscale priors; log(y) -# trains the GP on log-activity-coefficient (always positive, roughly linear in log-x). +# Paper's exact GP input/output transforms: log(x + 0.1) keeps the mole-fraction +# lengthscale well-scaled near the dilute limit; min-max normalizing T to [0, 1] +# matches the kernel's O(1) lengthscale priors; log(y) trains the GP on +# log-activity-coefficient (always positive, roughly linear in log-x). INPUT_TRANSFORM = InputTransform( forward_fns=[lambda x: np.log(x + 0.1), lambda T: (T - BOUNDS[1][0]) / 17.0], backward_fns=[lambda x: np.exp(x) - 0.1, lambda T: 17.0 * T + BOUNDS[1][0]], ) OUTPUT_TRANSFORM = OutputTransform(forward_fn=np.log, backward_fn=np.exp) -N_INIT = 10 # matches PrOHwater(nObs=10, ...) for the manuscript run -N_ITERS = 5 # adaptive design iterations (paper ran 15; kept small here -- -# full reproduction is Phase 7) +N_INIT = 10 # matches the paper's published run configuration +N_ITERS = 5 # adaptive design iterations (paper ran 15; kept small here for a quick +# demo -- see paper/full_reproduction.py for the exact paper configuration) COLUMN_VAR_NAMES = ["xW", "F", "xF", "R", "xD"] COLUMN_VAR_VALUES = [0.01, 100.0, 0.10, 1.0, 0.43] # Geankoplis Ex. 11.4-1 COLUMN_N_STAGES = 4 diff --git a/paper/DATA.md b/paper/DATA.md index 8c9f91f..3179016 100644 --- a/paper/DATA.md +++ b/paper/DATA.md @@ -1,23 +1,37 @@ -# Paper data & reproduction - -The bulk archived results from the published run (`less_x_new_manuscript_revisions`, -iteration 15) are ~2.5 GB and are **not** tracked in this repo. - -**Archive of record:** the private repository `dowlinglab/entropy_driven_hybrid_models_code` -(the original paper code, kept private with its full history) under -`entropy_driven_hms/results/less_x_new_manuscript_revisions/`. There is intentionally **no -Zenodo deposit** — the private old repo is the archive of record (REFACTOR_PLAN.md §7 -decision 4). - -Two small, curated subsets of that bulk archive **are** committed here: - -- `paper/reference/` — small scalar targets for regression (R̂/ESS, error metrics, posterior - summary, stage table), regenerable via `paper/extract_reference.py`. -- `paper/data/` — the exact plot-input files `paper/figures/*.py` read (~16 MB; see - `paper/data/README.md` for the file manifest and how it was determined). Added in - Phase 9 so `python paper/reproduce.py` regenerates every figure from a fresh clone - with **no private-archive access** — the old Phase 7 assumption ("figure reproduction - assumes author access to that repo") no longer holds for the default path. Author - access to the full private archive is still needed only if you want figures at - iterations beyond the curated subset (e.g. Fig 6/7 at an iteration other than 1/15) — - pass `--archive`/`$BFG_ARCHIVE_DIR` to point at it. +# Paper data + +**Everything needed to regenerate the paper's figures is committed in this repository.** +`python paper/reproduce.py` works from a fresh clone with no additional downloads and no +special access. If you only want the figures, you can stop reading here and go to +[`REPRODUCTION.md`](REPRODUCTION.md). + +## What is committed + +- **`paper/data/`** (~16 MB) — the run artifacts the figure scripts read: HMC traces and + posterior samples, the entropy fields, the surrogate phase-diagram draws, the Wilson + ground-truth curve, and the training/test activity-coefficient data. These are the + published run's own outputs, copied verbatim. See + [`data/README.md`](data/README.md) for the file-by-file manifest. +- **`paper/reference/`** — small scalar targets used by the regression tests (R̂ and ESS, + train/test error metrics, the hyperparameter-posterior summary, the McCabe–Thiele stage + table), regenerable with `paper/extract_reference.py`. + +## What is not committed, and why that doesn't matter + +The published run also produced roughly half a gigabyte of intermediate artifacts — full +posterior-prediction grids at all 60 iterations, per-iteration trace files, and rendered +figure images. Committing all of it would bloat the repository to no purpose: the figure +scripts read only the subset above, so that subset is what ships. + +The scripts accept `--archive` / `$BFG_ARCHIVE_DIR` to point at a different directory of run +artifacts, laid out the same way as `paper/data/`. That is useful if you re-run the loop +yourself (`paper/full_reproduction.py` writes such a directory) and want to plot *your* run +instead of the published one. It is not a prerequisite for anything in +[`REPRODUCTION.md`](REPRODUCTION.md). + +## The bundled paper + +`paper/bits_for_gaps_paper.pdf` is the published article, included so the method and the code +can be read side by side. It is redistributed under +[CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) — see +[`PAPER_LICENSE.md`](PAPER_LICENSE.md) for the citation and license statement. diff --git a/paper/PAPER_LICENSE.md b/paper/PAPER_LICENSE.md new file mode 100644 index 0000000..1edf86b --- /dev/null +++ b/paper/PAPER_LICENSE.md @@ -0,0 +1,30 @@ +# Bundled paper: citation and license + +`bits_for_gaps_paper.pdf` in this directory is the published article describing the method +implemented by this package. It is included so that the method and the code can be read +side by side. + +## Citation + +> K. D. Jones and A. W. Dowling, "BITS for GAPS: Bayesian Information-Theoretic Sampling for +> hierarchical GAussian Process Surrogates," *Computers & Chemical Engineering* **211** +> (2026) 109650. +> DOI: [10.1016/j.compchemeng.2026.109650](https://doi.org/10.1016/j.compchemeng.2026.109650) + +Please cite the DOI — it is the canonical reference. The bundled PDF is a convenience copy. + +## License + +© 2026 The Authors. Published by Elsevier Ltd. + +This is an open access article distributed under the terms of the +**Creative Commons Attribution (CC BY) 4.0** license: + + +CC BY 4.0 permits redistribution and adaptation in any medium or format, for any purpose, +including commercially, provided appropriate credit is given to the original authors, a link +to the license is provided, and any changes are indicated. **The PDF here is the unmodified +published version.** + +Note that this license applies to the *article*. The software in this repository is licensed +separately under the BSD 3-Clause license — see [`../LICENSE`](../LICENSE). diff --git a/paper/PHASE9B_INVESTIGATION.md b/paper/PHASE9B_INVESTIGATION.md deleted file mode 100644 index 1001847..0000000 --- a/paper/PHASE9B_INVESTIGATION.md +++ /dev/null @@ -1,168 +0,0 @@ -# Phase 9b: why the fully-adaptive surrogate's McCabe-Thiele column didn't converge - -Phase 9's from-scratch stochastic reproduction (`paper/REPRODUCTION.md`) documented one -discrepancy: the genuinely 15-iteration-adaptive surrogate GP's McCabe-Thiele column -didn't converge, unlike the paper's own Fig 9b surrogate column. The write-up attributed -this to a property of entropy-driven acquisition ("optimizes for GP predictive accuracy -at held-out points, not for a globally smooth-enough equilibrium curve"). **That -attribution was wrong.** The actual cause is a shared-mutable-state bug in -`paper/full_reproduction.py` itself, unrelated to entropy-driven design, the GP, or the -adaptive loop. This investigation found the bug, confirmed it by exactly reproducing the -original failure, fixed it, and re-ran the full stochastic loop to confirm the fix. - -## Root cause - -`paper/full_reproduction.py`'s original code computed the test-RMSE curve (via -`_predict_split` → `bits_for_gaps.mixture.sample_gp_posterior_mixture`) **before** -building the Fig 8/9-style surrogate phase diagram, and used the same in-memory -`GPmodel` object (`history.last.GPmodel`) for both. `sample_gp_posterior_mixture` -mutates `GPmodel.kernel` in place, by design (see its own docstring: "Mutated in -place: its kernel hyperparameters are reassigned for every draw") -- it's meant to be -called on a model you're about to discard or intentionally leave at that state. The -RMSE step called it ~30 times (once per HMC iteration, plus once more for the final -iteration's held-out train draws), each time reassigning `GPmodel.kernel` to a fresh -random draw from the HMC posterior and leaving it there. By the time the phase-diagram -code ran afterward and called `GPmodel.predict_f`, the kernel held whatever the *last* -of those ~30 unrelated draws happened to be -- not any principled state, and not what -the write-up assumed ("the final, genuinely 15-iteration adaptively-trained GP"). - -### Confirmation - -Reconstructing the iteration-15 GP from its checkpointed `gp_model_15.pkl` and -replaying `full_reproduction.py`'s exact original call sequence (RMSE loop, seed=10, -then phase diagram) reproduced the reported failure to 4 decimal places: - -``` -hyperparams BEFORE any mutation: [1.6508907124750056, 1.0349426523888388, 4.2802100682522735] -hyperparams AFTER the mutation sequence: [2.0146229510942004, 0.7560132676044007, 2.4816506654063213] -converged: False - stage 2: liquid=1.7258 vapor=0.3776 # full_run_summary.json (pre-fix) reported 1.725805727168204 -``` - -That match is exact, not approximate -- this is the bug, not a contributing factor. - -## Hypotheses: what held, what didn't - -The task posed four hypotheses before this root cause was known. None of them turned -out to be it, but testing them anyway produced useful, honest findings -- reported here -rather than discarded. - -**Setup**: five equilibrium curves were reconstructed over the same 50-point z-grid from -the Phase 9 run's checkpointed `gp_model_15.pkl` (24 training points) and -`param_posterior_samples_15` (the HMC trace, chain 0, constrained): - -- **(a) fresh, HMC-as-left** -- `surrogate_gamma` using the kernel state exactly as - `bits_for_gaps.gp.run_mcmc` left it (a single, un-mutated posterior sample). -- **(b) paper archive mean** -- pointwise mean of `paper/data/phase_diagram_15`'s y1 - column over its 50 draws, replicating `equilibrium.py`'s `water_proh_eqm` exactly. -- **(c) fresh, paper-method 50-draw average** -- 50 independent single-hyperparameter - draws from the same trace, GP mean under each, averaged pointwise (the paper's own - `new_phase_diagram.py` construction). -- **(d) fresh, posterior-mean hyperparameters** -- kernel set once to - `trace.mean(axis=0)`, single deterministic curve. -- **(e) Wilson** -- ground truth, for reference. - -| curve | monotonicity violations | roughness | column converged | -|---|---|---|---| -| (a) fresh, HMC-as-left | 0 | 0.00666 | **True** | -| (b) paper archive mean | 0 | 0.00343 | True | -| (c) fresh, 50-draw avg | 0 | 0.00650 | True | -| (d) fresh, mean-hyperparameters | 0 | 0.00661 | **False** | -| (e) Wilson | 0 | 0.00640 | True | - -(`paper/phase9_validation/phase9b_curve_comparison.png` plots all five together, plus a -zoom on the dilute region where the column's stages concentrate.) - -- **H1 (draw count/averaging)** -- the paper's method, replicated exactly on the fresh - surrogate (c), does converge and tracks Wilson closely. **Confirmed as a good - robustness practice** -- but note (a), a single un-mutated draw, *also* converged, so - averaging was not *necessary* to fix the reported failure; the bug was. Averaging is - adopted anyway (see "Fix" below) because it protects against exactly this class of - fragility: any single hyperparameter draw, including a good one, is one accidental - mutation away from being an arbitrary one. -- **H2 (monotonicity)** -- **rejected outright**. All five curves, including the one - that failed to converge (d), have zero monotonicity violations. Whatever made (d) fail - is not curve roughness or non-monotonicity. -- **H3 (deterministic mean, "likely the clean fix")** -- **falsified**. The - posterior-mean-hyperparameters curve (d) is exactly as smooth as the others (roughness - 0.00661 vs. (a)'s 0.00666) yet is the *only* one of the five that doesn't converge. - The arithmetic mean of each hyperparameter's marginal posterior is not a - self-consistent joint point (`std_dev`/`lengthscale_1` are LogNormal-ish, positive; - `lengthscale_2` is intentionally unconstrained -- see `kernels.py`), so it can land the - curve's absolute shape somewhere `scipy.optimize.fsolve`'s fixed initial guess in - `distillation.solve_column` can't reach a physical root from -- a solver - initial-guess-sensitivity issue, not a smoothness issue. Using the posterior *mean* is - not a safe substitute for a real posterior draw or a Monte-Carlo average over draws. -- **H4 (genuine under-resolution)** -- **rejected**. The un-mutated 24-point adaptive - surrogate (a) converges and matches Wilson about as well as the paper's own - archived-mean surrogate (b) does. Entropy-driven acquisition at this iteration count - is not under-resolving the equilibrium curve; there is no methodological limitation to - report here for this system. Phase 9's original framing of this as an entropy-driven- - design-vs-space-filling-design tradeoff was an incorrect post-hoc explanation for a - bug, not a real finding -- retracted below. - -## Fix - -Two changes, both in the **example layer** (`examples/vle_distillation/`, -`paper/full_reproduction.py`) -- `src/bits_for_gaps/` core is untouched, including -`mixture.sample_gp_posterior_mixture`'s in-place-mutation contract, which is correct and -documented as designed: - -1. **`paper/full_reproduction.py`**: build the phase diagram/column from - `history.last.GPmodel` *before* running the test-RMSE loop that mutates it, with a - comment explaining why the order matters. This alone fixes the reported bug. -2. **`examples/vle_distillation/phase_diagram.py`**: added `surrogate_gamma_averaged`, - matching the paper's own `new_phase_diagram.py`/`equilibrium.py` construction -- - draws `n_draws` (default 50) independent samples from a supplied HMC trace, - evaluates this GP's deterministic conditional mean (`predict_f`, not - `predict_f_samples` -- no TF-ambient-RNG non-reproducibility) under each, and - averages `gamma_proh` pointwise. `full_reproduction.py` now uses this (not the - single-point-estimate `surrogate_gamma`) for its phase diagram, for the added - robustness H1 confirmed -- not because the point estimate doesn't work, but because - a Monte-Carlo average over the posterior is inherently more robust to *any* single - draw (mutated-into or genuinely unlucky) being atypical. `surrogate_gamma`'s - docstring now documents the mutation hazard explicitly for future callers. - `fig09_mccabe_thiele.py`'s non-adaptive, MLE-fit surrogate (no HMC trace, hence no - posterior to average over) is unaffected and still uses `surrogate_gamma`, correctly. - -## Result after the fix - -Re-ran the full 15-iteration adaptive loop from scratch (fresh seed=10 run, ~26 min; -`results_remaked/phase9b_fullrun_fixed/`, gitignored) with the fixed -`paper/full_reproduction.py`: - -``` -column_wilson_converged: True -column_surrogate_converged: True # was False before the fix - - stage 1: wilson liq=0.2265 vap=0.4300 | surrogate liq=0.2537 vap=0.4300 - stage 2: wilson liq=0.0473 vap=0.3282 | surrogate liq=0.0433 vap=0.3419 - stage 3: wilson liq=0.0216 vap=0.2386 | surrogate liq=0.0197 vap=0.2366 - stage 4: wilson liq=0.0100 vap=0.1399 | surrogate liq=0.0100 vap=0.1435 -``` - -The genuinely-adaptive surrogate's column now converges and tracks the Wilson -ground-truth column within 0.03 mole fraction (liquid) / 0.014 (vapor) at every stage -- -the "satisfying full reproduction" this investigation set out to find. HMC diagnostics -(R-hat, ESS) and the hyperparameter posterior on this fresh run again matched the -originally-reported values to 6-8 significant figures (same near-bit reproducibility -Phase 9 established), confirming the fix changed only the phase-diagram/column -construction, not the adaptive loop itself. The test-RMSE curve differs slightly from -the original run's (4.75 vs. 4.34 at iteration 1) for the same already-documented reason -as Phase 9's other test-RMSE discrepancy: `predict_f_samples`'s non-seedable TF ambient -RNG. - -Updated artifacts: `paper/phase9_validation/full_run_summary.json` (now the fixed run; -the original is preserved as `full_run_summary_pre_phase9b_fix.json` for the record), -`phase_diagram_fresh.png` (regenerated, now shows the converging surrogate), and -`phase9b_curve_comparison.png` (the five-curve diagnostic above). See -`paper/REPRODUCTION.md`'s Phase 9 section for the corrected discrepancy note. - -## Takeaway - -Not a methodological finding about entropy-driven design after all -- a straightforward -shared-mutable-state bug in a one-off analysis script, caught by refusing to accept the -first plausible-sounding explanation and instead reconstructing the failure from first -principles until it reproduced exactly. Worth remembering for any future script that -reuses a fitted model object across multiple purposes: check whether the functions -touching it document mutation, and if so, order calls (or copy the model) accordingly. diff --git a/paper/REPRODUCTION.md b/paper/REPRODUCTION.md index ef32190..9b5b37d 100644 --- a/paper/REPRODUCTION.md +++ b/paper/REPRODUCTION.md @@ -1,16 +1,16 @@ # Reproducing the paper's figures `paper/reproduce.py` regenerates the paper's figures (Jones & Dowling 2026, "BITS for -GAPS") from the published run's plot-input data (iteration 15 of -`less_x_new_manuscript_revisions`). It does **not** re-run the paper's 15-iteration -adaptive HMC loop (stochastic, expensive, and not what reproducing a figure requires) --- it loads what that loop already produced and renders it through the refactored -`bits_for_gaps` package + `examples/vle_distillation`. +GAPS") from the published run's plot-input data (iteration 15 of the paper's original +run). It does **not** re-run the paper's 15-iteration adaptive HMC loop (stochastic, +expensive, and not what reproducing a figure requires) -- it loads what that loop +already produced and renders it through the `bits_for_gaps` package + +`examples/vle_distillation`. -**As of Phase 9, this needs no private-archive access by default** -- the data comes -from the curated, committed `paper/data/` subset (~16 MB; see `paper/data/README.md`). -Point `--archive`/`$BFG_ARCHIVE_DIR` at the full private archive only if you want -figures at iterations beyond that curated subset (see `paper/DATA.md`). +**Everything needed is committed here** -- the data comes from `paper/data/` (~16 MB; +see `paper/data/README.md`), so this works from a fresh clone. Point +`--archive`/`$BFG_ARCHIVE_DIR` at another directory of run artifacts (same layout) to plot +a different run, such as one produced by `paper/full_reproduction.py`. ```bash export PYTHON_JULIACALL_HANDLE_SIGNALS=yes # macOS; needed for Fig 8/9 (Clapeyron) @@ -24,7 +24,7 @@ script produces is committed. For the from-scratch stochastic reproduction of the full adaptive loop (a separate, one-time validation exercise, not part of `paper/reproduce.py`), see -["Phase 9: from-scratch stochastic reproduction"](#phase-9-from-scratch-stochastic-reproduction) +["Re-running the full adaptive loop from scratch"](#re-running-the-full-adaptive-loop-from-scratch) below. ## Figure map @@ -44,32 +44,31 @@ below. | 12 | Hyperparameter joint marginals | `fig12_joint_marginals.py` | `param_posterior_samples_15` | none -- visual (same posterior samples Fig 11's numbers are pinned from) | **Quantitatively pinned** (5, 8, 9, 10, 11): a committed reference JSON or a -`paper/data/`-vs-recompute cross-check is diffed within a stated tolerance. -**As of Phase 9, only the tests that *recompute* something via Clapeyron stay -gated** (`@pytest.mark.vle`, deselected by default -- run with `pytest -m vle`, -needs Julia): Fig 8's Wilson-curve cross-check and Fig 9's full stage-table -recompute. Fig 5/10/11's tests only *read* the committed `paper/data/` text files -(no Julia, no private archive) and now run in the **default** `pytest -q` suite. -Fig 8 has no dedicated reference *file* (it's a visual reproduction, not a -paper-reported scalar) -- its pin is a direct cross-check against the committed -`gt_Wilson_data` instead. `BFG_ARCHIVE_DIR` can still point any of these at the full -private archive instead of `paper/data/` (e.g. to check other iterations), but -that's no longer required for the default-suite tests to run. +`paper/data/`-vs-recompute cross-check is diffed within a stated tolerance. Only the +tests that *recompute* something via Clapeyron stay gated (`@pytest.mark.vle`, +deselected by default -- run with `pytest -m vle`, needs Julia): Fig 8's Wilson-curve +cross-check and Fig 9's full stage-table recompute. Fig 5/10/11's tests only *read* the +committed `paper/data/` text files (no Julia needed) and run in the +**default** `pytest -q` suite. Fig 8 has no dedicated reference *file* (it's a visual +reproduction, not a paper-reported scalar) -- its pin is a direct cross-check against +the committed `gt_Wilson_data` instead. `BFG_ARCHIVE_DIR` can point any of these at a +different run's artifacts instead of `paper/data/`, but that's not required for the +default-suite tests to run. **Visually reproduced** (2, 3, 4, 6, 7, 12): no reference file exists for these (the paper doesn't report them as scalars), so they're spot-checked by eye against the -archived PNGs' structure during development -- correct qualitative behavior (LHS -space-filling, entropy concentrating away from sampled points, CI narrowing with more -data, etc.), not pixel-identical figures. +archived PNGs' structure -- correct qualitative behavior (LHS space-filling, entropy +concentrating away from sampled points, CI narrowing with more data, etc.), not +pixel-identical figures. ## Simplifications from the original plotting code -This is reproduction code (`paper/figures/`), not a library API -- ported from -`fxns/mcmc_plotter.py` (847 lines) pragmatically, not verbatim. Dropped, since they're -purely visual and don't change the figure's content or claim: +This is reproduction code (`paper/figures/`), ported from the paper's original +plotting scripts pragmatically, not verbatim -- it is not a library API. Dropped, +since they're purely visual and don't change the figure's content or claim: -- **Fig 5**: the zoomed inset panel (original `plot_parity`'s `inset_axes` zoom into - the low-error region). The main parity plot + error bars are unchanged. +- **Fig 5**: the zoomed inset panel (the original's zoom into the low-error region). + The main parity plot + error bars are unchanged. - **Fig 12**: the KDE contour overlay on each hexbin panel. Kept the hexbin density, the MAP point, and the 95%-credible-interval box; the MAP point here is the sample nearest the coordinate-wise median (a cheap proxy), not a true density mode -- @@ -80,70 +79,68 @@ purely visual and don't change the figure's content or claim: already implies a small multi-panel figure, not 60 separate files. - **Fig 9**: the surrogate panel uses a freshly-trained (30-point LHS + MLE fit) GP, not the paper's actual 15-iteration adaptively-designed surrogate -- see - `tests/regression/test_mccabe_thiele.py` and `HANDOFF.md` (Phase 6) for why - reproducing that exactly is Phase 7 work that was deliberately *not* undertaken - (it would mean re-running the stochastic adaptive loop, against this phase's - guardrail). The Wilson panel is an exact physics recompute either way. + `tests/regression/test_mccabe_thiele.py`; reproducing that surrogate exactly would + mean re-running the full stochastic adaptive loop (see below), which this default + figure-regeneration path deliberately does not do. The Wilson panel is an exact + physics recompute either way. ## Known discrepancies (not bugs -- documented, not "fixed") - **`paper/reference/mccabe_thiele_stages.json`'s `"wilson"` column has ~0.01-level transcription slop** (it was hand-read off paper Fig 9c, not computed from - archived data -- see `paper/reference/README.md`). The fresh Clapeyron recompute in + archived data -- see `paper/reference/README.md`). A fresh Clapeyron recompute in `fig09_mccabe_thiele.wilson_column()` hits real physical landmarks exactly (e.g. stage 1 vapor = `xD` = 0.43 by construction; pure-PrOH bubble point 370.35 K matches 1-propanol's real normal boiling point to 4 significant figures) yet differs from reference's `"wilson"` entries by up to 0.01 -- this is the reference - file's own transcription precision limit, not a port error (see `HANDOFF.md` - Phase 6 for the full analysis). + file's own transcription precision limit, not a port error. -## Phase 9: from-scratch stochastic reproduction +## Re-running the full adaptive loop from scratch Everything above regenerates figures from the published run's *plot-input data* -- it never re-executes the paper's sequential-design loop. This section does: it runs `paper/full_reproduction.py`, which drives `bits_for_gaps.sampler.BitsForGaps` through -the same 15-iteration adaptive HMC + entropy-acquisition loop as the published -`less_x_new_manuscript_revisions` run (same bounds, seed, transforms, kernel, and HMC -config -- see the script's module docstring for the full paper-trail), starting from a -*fresh* 10-train/10-test LHS design and a live Clapeyron/Wilson black box -- no -archived data read as input anywhere in this section. +the same 15-iteration adaptive HMC + entropy-acquisition loop as the paper's published +run (same bounds, seed, transforms, kernel, and HMC config -- see the script's module +docstring for the exact constants), starting from a *fresh* 10-train/10-test LHS +design and a live Clapeyron/Wilson black box -- no archived data read as input +anywhere in this section. ```bash export PYTHON_JULIACALL_HANDLE_SIGNALS=yes -python paper/full_reproduction.py --out-dir results_remaked/phase9_fullrun # ~25 min +python paper/full_reproduction.py --out-dir results_remaked/full_reproduction # ~25 min ``` This is a **one-time validation exercise**, not a regression test -- there is no gated CI test for it (re-running it is expensive and its whole point is to check -statistical/qualitative, not bitwise, agreement). The run below completed -2026-07-04; artifacts stayed in the gitignored `results_remaked/phase9_fullrun/`; the -numbers and two summary plots below are committed under `paper/phase9_validation/` -so this section doesn't rot into an unverifiable claim. +statistical/qualitative, not bitwise, agreement). The numbers below are from a +completed run (artifacts stayed in the gitignored `results_remaked/`, as always); +they're recorded here rather than left as an unverifiable claim. **Headline result: closer to bitwise than expected.** Everything in the loop that runs through a `self.seed`-seeded path (LHS design, HMC sampling, the entropy- -acquisition optimizer) turned out to reproduce the published run to 6-8 significant -figures -- not just "qualitatively similar." The only place real stochastic drift -shows up is the GP posterior-*predictive* mixture sampling (`gpflow`'s -`predict_f_samples`, which draws from TensorFlow's ambient RNG -- documented as -non-reproducible in `bits_for_gaps/mixture.py`'s module docstring, and already known -to differ run-to-run even in the original paper code). That function feeds exactly -two things: the test-RMSE curve below, and the Fig-8-style surrogate phase diagram. -Everything else pinned on the HMC trace directly (R-hat/ESS, hyperparameter -posterior, entropy field) lines up almost exactly. +acquisition optimizer) reproduced the published run to 6-8 significant figures -- +not just "qualitatively similar." The only place real stochastic drift shows up is +the GP posterior-*predictive* mixture sampling (`gpflow`'s `predict_f_samples`, which +draws from TensorFlow's ambient RNG -- documented as non-reproducible in +`bits_for_gaps/mixture.py`'s module docstring, and non-reproducible in the paper's +original code too). That function feeds exactly two things: the test-RMSE curve +below, and the Fig-8-style surrogate phase diagram. Everything else pinned on the HMC +trace directly (R-hat/ESS, hyperparameter posterior, entropy field) lines up almost +exactly. ### HMC diagnostics (iteration 15, trained on the same 24 points as `gp_model_15`) | | R-hat | ESS | |---|---|---| | Paper (`paper/reference/hmc_diagnostics.json`) | 1.0052, 1.0073, 1.0088 | 1468.3, 2428.1, 653.1 | -| Fresh run (`full_run_summary.json`) | 1.00523, 1.00730, 1.00879 | 1468.29, 2428.09, 653.15 | +| Fresh run | 1.00523, 1.00730, 1.00879 | 1468.29, 2428.09, 653.15 | All well under the R-hat < 1.1 convergence threshold, and ESS is healthy relative to -20,000 raw HMC samples (4 chains x 5000). The near-exact match here is a strong, -unplanned confirmation that Phases 2-8's refactor preserved the HMC path bit-for-bit --- this isn't the paper's own number re-displayed, it's an independently re-run HMC -fit that happens to land on the same posterior. +20,000 raw HMC samples (4 chains x 5000). The near-exact match here is a strong +confirmation that this package's decomposition of the original algorithm preserves +the HMC path bit-for-bit -- this isn't the paper's own number re-displayed, it's an +independently re-run HMC fit that happens to land on the same posterior. ### Hyperparameter posterior (`std_dev`, `lengthscale_1`, `lengthscale_2`) @@ -167,8 +164,7 @@ Per-iteration max entropy (grid max of `entropy_{i}`, both committed and fresh): | Fresh run | 1.4577 | 1.2470 | 0.6505 | 0.0215 | 0.0038 | -0.0084 | Same monotonic-ish decay, same sign change (uncertainty-driven exploration -"exhausts" the space) landing between iterations 14 and 15 in both runs. Plotted in -`paper/phase9_validation/rmse_and_entropy.png`, panel (b). +"exhausts" the space) landing between iterations 14 and 15 in both runs. ### Predictive accuracy (test RMSE, activity coefficient units) @@ -183,12 +179,11 @@ regime (~5x reduction vs. the paper's ~6.5x) but not on the same value -- expect since this metric is the one place `predict_f_samples`'s non-reproducible draws enter, compounded by a documented simplification: this run's mixture uses a 15-component hyperparameter-posterior subset (`noGaussians`, the same size the -acquisition function itself uses) rather than the paper's own -`train_test_split_proh.py`, which drew a dedicated 500-sample subset just for this -plot. The trend is non-monotonic in both runs (e.g. this run's RMSE ticks up at -iterations 8 and 12 before falling again) -- expected for an entropy-driven -acquisition, which optimizes information gain, not held-out error, at each step. -Plotted in `paper/phase9_validation/rmse_and_entropy.png`, panel (a). +acquisition function itself uses) rather than the paper's own dedicated 500-sample +subset for this plot. The trend is non-monotonic in both runs (e.g. this run's RMSE +ticks up at iterations 8 and 12 before falling again) -- expected for an +entropy-driven acquisition, which optimizes information gain, not held-out error, at +each step. ### Phase diagram (Fig 8) and McCabe-Thiele stage table (Fig 9) @@ -200,38 +195,21 @@ Clapeyron computation with no randomness anywhere in it. The **genuinely** 15-iteration-adaptive surrogate GP (this run's final, 24-point model -- not `fig09_mccabe_thiele.py`'s dedicated 30-point-LHS/MLE-fit stand-in) gives a phase diagram within 0.83 K / 0.02 mole fraction of the Wilson curve -- a good -surrogate. See `paper/phase9_validation/phase_diagram_fresh.png`. - -**Update (Phase 9b, corrects the paragraph originally here):** this surrogate's -McCabe-Thiele column *does* converge and track the Wilson column closely -- within -0.03 mole fraction (liquid) / 0.014 (vapor) at every stage. The original run reported -`column_surrogate_converged: false` with an unphysical stage-2 liquid mole fraction of -1.73, and this file used to attribute that to a property of entropy-driven -acquisition ("optimizes for GP predictive accuracy at held-out points, not for a -globally smooth enough equilibrium curve"). **That attribution was wrong** -- a -dedicated investigation (`paper/PHASE9B_INVESTIGATION.md`) traced it to a -shared-mutable-state bug in `paper/full_reproduction.py`: the test-RMSE step mutated -the same `GPmodel` object the phase-diagram step reused afterward -(`bits_for_gaps.mixture.sample_gp_posterior_mixture` reassigns `GPmodel.kernel` in -place, by design), leaving the kernel at an arbitrary leftover hyperparameter state by -the time the phase diagram was built. Fixed by reordering the script and by adding a -posterior-hyperparameter-averaged surrogate construction -(`examples/vle_distillation/phase_diagram.py`'s `surrogate_gamma_averaged`, matching -the paper's own `new_phase_diagram.py` method) for extra robustness. See -`paper/PHASE9B_INVESTIGATION.md` for the full root-cause analysis, which hypotheses -were tested, and the before/after numbers -- there is no real methodological -limitation of entropy-driven design to report here after all. +surrogate, and its McCabe-Thiele column converges and tracks the Wilson column +closely -- within 0.03 mole fraction (liquid) / 0.014 (vapor) at every stage. +`mixture.sample_gp_posterior_mixture` reassigns a GP's kernel hyperparameters once per +posterior draw during sampling (see its docstring); building this phase diagram from +`history.last.GPmodel` must happen before any test-RMSE-style step that calls that +function on the same model object, or the kernel is left at an arbitrary leftover +hyperparameter state. `paper/full_reproduction.py` orders its calls accordingly and +uses `surrogate_gamma_averaged` (a posterior-hyperparameter-averaged, rather than +single-draw, activity-coefficient estimate) for the phase diagram specifically because +that ordering matters. ### Summary -Qualitative/statistical agreement confirmed on every axis the acceptance criteria -asked for, and considerably tighter than "qualitative" on the deterministic parts of -the pipeline (HMC diagnostics, hyperparameter posterior, entropy field). The adaptive -surrogate's McCabe-Thiele column, initially reported as not converging, does converge -once a Phase 9b bug fix removed an accidental shared-mutable-state issue in -`paper/full_reproduction.py` (see `paper/PHASE9B_INVESTIGATION.md`) -- there is no -real methodological limitation of entropy-driven acquisition to report here. Full -numbers: `paper/phase9_validation/full_run_summary.json` (the pre-fix run is preserved -as `full_run_summary_pre_phase9b_fix.json` for the record). Reproduce via -`paper/full_reproduction.py` (module docstring has the full paper-trail for every -constant); not gated in CI (see above). +Qualitative/statistical agreement confirmed on every axis this exercise checked, and +considerably tighter than "qualitative" on the deterministic parts of the pipeline +(HMC diagnostics, hyperparameter posterior, entropy field). Reproduce via +`paper/full_reproduction.py` (module docstring has the exact constants used); not +gated in CI (see above). diff --git a/paper/__init__.py b/paper/__init__.py index f6b2610..c2438b5 100644 --- a/paper/__init__.py +++ b/paper/__init__.py @@ -1,5 +1,5 @@ -"""Paper reproduction (repo-only -- see REFACTOR_PLAN.md §7.3, not shipped in the wheel). +"""Paper reproduction (repo-only, not shipped in the wheel). Importable as ``paper.figures.`` in dev/CI via ``tests/conftest.py``'s ``sys.path`` -insert of the repo root (mirrors the ``examples/`` mechanism from Phase 6). +insert of the repo root (mirrors the same mechanism for ``examples/``). """ diff --git a/paper/bits_for_gaps_paper.pdf b/paper/bits_for_gaps_paper.pdf new file mode 100644 index 0000000..1220a37 Binary files /dev/null and b/paper/bits_for_gaps_paper.pdf differ diff --git a/paper/data/README.md b/paper/data/README.md index 0533018..00eb6c3 100644 --- a/paper/data/README.md +++ b/paper/data/README.md @@ -1,24 +1,19 @@ -# Curated plot-input subset (tracked) +# Plot-input data for the paper's figures (tracked) **This directory IS tracked in git** (unlike `results_remaked/`, which is gitignored -output). It's a small, curated subset of the private archive's published run -- -exactly the files `paper/figures/*.py` read, no more -- committed so that -`python paper/reproduce.py` works from a fresh clone with **no private-archive -access** (Phase 9, STEP 1). +output). It holds exactly the files `paper/figures/*.py` read -- no more -- so that +`python paper/reproduce.py` regenerates every figure from a fresh clone with no +additional downloads and no special access. ## Provenance -Copied verbatim (unmodified) from the private old repo's archive of record: +These are the published run's own output files, copied verbatim (unmodified) from the run +that produced the paper's figures. **Iteration 15** is that published run: its +`rhat_value_15.txt` / `ess_value_15.txt` match paper Fig 10 exactly (pinned in +`paper/reference/hmc_diagnostics.json`). -``` -~/DowlingLab/CAREER/entropy_driven_hybrid_models_code/entropy_driven_hms/ - results/less_x_new_manuscript_revisions/ -``` - -**Iteration 15** = the published run (its `rhat_value_15.txt`/`ess_value_15.txt` -match paper Fig 10 exactly -- see `paper/reference/hmc_diagnostics.json`). See -`paper/DATA.md` for how this relates to the bulk (~2.5 GB, not committed anywhere) -archive. +That run also produced a much larger set of intermediate artifacts that are deliberately +not committed, because no figure script reads them -- see `paper/DATA.md`. ## How this set was determined @@ -44,7 +39,7 @@ files, ~16 MB. Fig 9 (McCabe-Thiele) needs none of this -- it's a pure physics recompute (Wilson via live Clapeyron; the surrogate via a fresh LHS + MLE-fit GP, same as -`tests/regression/test_mccabe_thiele.py`, Phase 6/7). +`tests/regression/test_mccabe_thiele.py`). ## Why plain text, not compressed or `gp_model_*.pkl`-recomputed @@ -68,7 +63,7 @@ this directory. Two alternatives were considered and rejected for now: ## Do not edit by hand -These are frozen snapshots of a specific archived run. If the curated set ever needs -to change (a new figure needs a new file, or an existing one needs a different -iteration), copy the new file(s) from the private archive the same way -- don't -regenerate or edit any file here. +These are frozen snapshots of the published run. If this set ever needs to change (a new +figure needs a new file, or an existing one needs a different iteration), copy the new +file(s) verbatim from that run's artifacts the same way -- don't regenerate or edit any +file here. diff --git a/paper/extract_reference.py b/paper/extract_reference.py index 5f090cd..70bd920 100644 --- a/paper/extract_reference.py +++ b/paper/extract_reference.py @@ -1,12 +1,13 @@ -"""Extract reference scalar targets from the archived published run (iteration 15). +"""Regenerate the reference scalar targets in paper/reference/ from run artifacts. -Pure NumPy; reads the read-only old-repo archive and writes small JSON files into +Pure NumPy; reads a directory of run artifacts and writes small JSON files into paper/reference/ (a sibling of this script). Run from anywhere: python paper/extract_reference.py -Archived run = results/less_x_new_manuscript_revisions, iteration 15 = the published -run whose R-hat/ESS match paper Fig 10 exactly. +Iteration 15 is the published run, whose R-hat/ESS match paper Fig 10 exactly. The +committed reference files were produced from the published run's artifacts, which ship +in paper/data/ -- so this script reproduces them from a fresh clone. """ import json @@ -14,9 +15,12 @@ import numpy as np -ARCHIVE = os.path.expanduser( - "~/DowlingLab/CAREER/entropy_driven_hybrid_models_code/entropy_driven_hms/" - "results/less_x_new_manuscript_revisions" +# Directory of run artifacts to summarize. Defaults to the committed paper/data/ subset; +# set $BFG_ARCHIVE_DIR to summarize a different run instead (for example one produced by +# paper/full_reproduction.py). +ARCHIVE = os.environ.get( + "BFG_ARCHIVE_DIR", + os.path.join(os.path.dirname(os.path.abspath(__file__)), "data"), ) OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "reference") os.makedirs(OUT, exist_ok=True) @@ -69,7 +73,7 @@ def write_json(name, obj): # y_true from activity_data_1 (fixed 10 train pts) / activity_test_points (10 test). # Predictions: gp_predict_{train,test}_{iter} = (n_points, 500 mixture draws). # Metrics computed per-draw (axis=0 over points) => 500-length distributions, -# exactly as fxns/train_test_split_proh.plot_error_bx_n_wskr. +# matching the paper's own error-metric computation. # --------------------------------------------------------------------------- y_train = np.loadtxt(os.path.join(ARCHIVE, "activity_data_1"))[:, 2] y_test = np.loadtxt(os.path.join(ARCHIVE, "activity_test_points"))[:, 2] @@ -126,7 +130,7 @@ def write_json(name, obj): # --------------------------------------------------------------------------- # 4. McCabe-Thiele stage table (paper Fig 9c) -- values from the paper (the column # design + surrogate/Wilson stage compositions). Reproduction needs the Julia VLE -# distillation backend (Phase 6 port), so the recompute test is @pytest.mark.vle. +# distillation backend, so the recompute test is @pytest.mark.vle. # --------------------------------------------------------------------------- write_json( "mccabe_thiele_stages.json", diff --git a/paper/figures/__init__.py b/paper/figures/__init__.py index 940a1a9..5b7def1 100644 --- a/paper/figures/__init__.py +++ b/paper/figures/__init__.py @@ -1,9 +1,7 @@ """Per-figure reproduction scripts (paper Figs 2-12), one module per figure. -Ported from the paper code's ``fxns/mcmc_plotter.py`` (847-line figure library) and -``fxns/plot_res.py`` (CLI dispatch), split one figure per module and simplified -- -this is reproduction code, not a library API, so it favors reading directly over -matching every original styling detail. Each module exposes a ``make(archive_dir, -out_dir)`` function that reads the archived published run (iteration 15) and writes -the figure into ``out_dir`` (see ``paper/reproduce.py``). +This is reproduction code, not a library API, so it favors reading directly over +matching every detail of the paper's original plotting code. Each module exposes a +``make(archive_dir, out_dir)`` function that reads the archived published run +(iteration 15) and writes the figure into ``out_dir`` (see ``paper/reproduce.py``). """ diff --git a/paper/figures/_archive.py b/paper/figures/_archive.py index d256668..9b28a48 100644 --- a/paper/figures/_archive.py +++ b/paper/figures/_archive.py @@ -1,13 +1,13 @@ """Shared loaders for the published run's plot-input data + common plot style. -As of Phase 9, ``paper/reproduce.py`` defaults ``archive_dir`` to the curated, -COMMITTED ``paper/data/`` subset (~16 MB -- exactly the files the figures below -read; see ``paper/data/README.md`` for the file manifest and provenance), so figure -reproduction needs no private-archive access out of the box. ``archive_dir`` can -still be pointed at the full private archive (the archive of record -- REFACTOR_PLAN -§7 decision 4, ``paper/DATA.md``) via ``--archive``/``$BFG_ARCHIVE_DIR`` for -iterations/figures beyond the curated subset. Every figure module reads via the -functions below rather than hardcoding paths, so both sources work unchanged. +``paper/reproduce.py`` defaults ``archive_dir`` to the COMMITTED ``paper/data/`` directory +(~16 MB -- exactly the files the figures below read; see ``paper/data/README.md`` for the +manifest and provenance), so figure reproduction works from a fresh clone with no extra +downloads. ``archive_dir`` can be pointed at any other directory of run artifacts with the +same layout via ``--archive``/``$BFG_ARCHIVE_DIR`` -- for example the output directory of +``paper/full_reproduction.py``, to plot your own run instead of the published one. Every +figure module reads via the functions below rather than hardcoding paths, so both work +unchanged. Iteration 15 is the published run: its ``rhat_value_15.txt``/``ess_value_15.txt`` match paper Fig 10 exactly (see ``paper/reference/hmc_diagnostics.json``). @@ -24,19 +24,18 @@ def require_archive(archive_dir): """Raise a clear error if the data directory isn't where expected.""" if not os.path.isdir(archive_dir): raise FileNotFoundError( - f"Data directory not found: {archive_dir!r}. This should be the " - f"committed 'paper/data/' subset (see paper/data/README.md), or -- for " - f"iterations/figures beyond that curated subset -- the private old " - f"repo's 'entropy_driven_hms/results/less_x_new_manuscript_revisions/' " - f"(see paper/DATA.md). Pass it via `paper/reproduce.py --archive ` " - f"or the BFG_ARCHIVE_DIR environment variable." + f"Data directory not found: {archive_dir!r}. This should be the committed " + f"'paper/data/' directory (the default -- see paper/data/README.md), or " + f"another directory of run artifacts with the same layout, such as the " + f"output of paper/full_reproduction.py. Pass it via " + f"`paper/reproduce.py --archive ` or the BFG_ARCHIVE_DIR " + f"environment variable." ) return archive_dir def apply_plot_settings(): - """Global rcParams matching the paper code's ``fxns/plot_settings.py`` / - ``fxns/mcmc_plotter.py`` header -- shared look across figures.""" + """Global rcParams matching the paper's figure style -- shared look across figures.""" import matplotlib.pyplot as plt plt.rcParams["figure.figsize"] = (6, 6) @@ -64,9 +63,9 @@ def load_rhat_ess(archive_dir, iters=PUBLISHED_ITERS): def load_traces(archive_dir, iters=PUBLISHED_ITERS): """HMC chain traces, shape (n_samples, n_chains, n_hyperparameters). - Matches ``fxns/plot_res.py``'s ``all_traces`` mode exactly: collect every - ``traces_chain_*_exp_{iters}`` file (sorted by filename, i.e. by chain index), - stack, then transpose from (chain, sample, param) to (sample, chain, param). + Collects every ``traces_chain_*_exp_{iters}`` file (sorted by filename, i.e. by + chain index), stacks, then transposes from (chain, sample, param) to + (sample, chain, param). """ suffix = f"_exp_{iters}" chains = [] @@ -118,7 +117,7 @@ def load_lhs_design(archive_dir): def load_cont_data(archive_dir): """Dense ground-truth activity-coefficient grid, shape (n, 4): [z, T, gamma_PrOH, - gamma_H2O] (see ``proh_water_class.gen_cont_activities_2D``).""" + gamma_H2O].""" return np.loadtxt(os.path.join(archive_dir, "cont_data")) diff --git a/paper/figures/fig02_lhs_design.py b/paper/figures/fig02_lhs_design.py index 4f61bf5..76680c7 100644 --- a/paper/figures/fig02_lhs_design.py +++ b/paper/figures/fig02_lhs_design.py @@ -1,8 +1,7 @@ """Fig 2 -- initial Latin-hypercube design (train/test split) over (z_PrOH, T). -Ported from ``fxns/mcmc_plotter.py``'s ``plot_lhs_2d`` (``fxns/plot_res.py``'s -``-m lhs_2D`` mode). Visual reproduction only (no reference pin) -- spot-check against -the archived ``gp_lhs_design.png``. +Visual reproduction only (no reference pin) -- spot-check against the archived +``gp_lhs_design.png``. """ import os diff --git a/paper/figures/fig03_entropy_field.py b/paper/figures/fig03_entropy_field.py index 5f9685f..f4919b8 100644 --- a/paper/figures/fig03_entropy_field.py +++ b/paper/figures/fig03_entropy_field.py @@ -1,10 +1,9 @@ """Fig 3 -- predictive-entropy field over the (z_PrOH, T) design space, early iterations. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_entropy_2D`` (``fxns/plot_res.py``'s -``-m entropy_2D`` mode), shown as a small multi-panel grid across the first few -sequential-design iterations (matching the paper's lettered-panel (a)-(f) scheme) -rather than one file per iteration. Visual reproduction only (no reference pin) -- -spot-check against the archived ``entropy_surface_{iters}.png`` files. +Shown as a small multi-panel grid across the first few sequential-design iterations +(matching the paper's lettered-panel (a)-(f) scheme) rather than one file per iteration. +Visual reproduction only (no reference pin) -- spot-check against the archived +``entropy_surface_{iters}.png`` files. """ import os diff --git a/paper/figures/fig04_entropy_evolution.py b/paper/figures/fig04_entropy_evolution.py index 32bae8f..ff40f30 100644 --- a/paper/figures/fig04_entropy_evolution.py +++ b/paper/figures/fig04_entropy_evolution.py @@ -1,8 +1,7 @@ """Fig 4 -- maximum predictive entropy vs. sequential-design iteration. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_entropy_v_iters`` (``fxns/plot_res.py``'s -``-m ent_v_iters`` mode). Visual reproduction only (no reference pin) -- spot-check -against the archived ``entropy_v_iters.png``. +Visual reproduction only (no reference pin) -- spot-check against the archived +``entropy_v_iters.png``. """ import os diff --git a/paper/figures/fig05_parity.py b/paper/figures/fig05_parity.py index 927b11a..b4b277d 100644 --- a/paper/figures/fig05_parity.py +++ b/paper/figures/fig05_parity.py @@ -1,8 +1,7 @@ """Fig 5 -- train/test parity plots + RMSE/MAE error box-and-whisker plots. -Ported from ``train_test_split_proh.py``'s ``plot_parity`` / ``plot_parity_loglog`` / -``plot_error_bx_n_wskr``, simplified (dropped the zoomed inset -- a purely visual -detail, not the quantitative content) -- this is reproduction code, not a library API. +Simplified from the paper's own plotting code: dropped the zoomed inset (a purely +visual detail, not the quantitative content). Quantitatively pinned: ``error_metrics()`` returns exactly the per-draw RMSE/MAE distributions ``paper/reference/fig5_error_metrics.json`` was extracted from (same diff --git a/paper/figures/fig06_gp_posterior_surface.py b/paper/figures/fig06_gp_posterior_surface.py index 3e12044..cfec1dc 100644 --- a/paper/figures/fig06_gp_posterior_surface.py +++ b/paper/figures/fig06_gp_posterior_surface.py @@ -1,8 +1,7 @@ """Fig 6 -- GP posterior predictive surface (3-D), iterations 1 vs. 15. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_gp_posterior_2D`` (``fxns/plot_res.py``'s -``-m gp_post_2D`` mode). Visual reproduction only (no reference pin) -- spot-check -against the archived ``gp_posterior_2D_{1,15}.png``. +Visual reproduction only (no reference pin) -- spot-check against the archived +``gp_posterior_2D_{1,15}.png``. """ import os diff --git a/paper/figures/fig07_gp_posterior_isotherms.py b/paper/figures/fig07_gp_posterior_isotherms.py index 30592bf..760d70f 100644 --- a/paper/figures/fig07_gp_posterior_isotherms.py +++ b/paper/figures/fig07_gp_posterior_isotherms.py @@ -1,8 +1,7 @@ """Fig 7 -- GP posterior predictive at fixed isotherms, offset for readability. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_gp_post_multiple_isotherms`` -(``fxns/plot_res.py``'s ``-m gp_post`` mode). Visual reproduction only (no reference pin) --- spot-check against the archived ``gp_posterior_isotherms_15.png``. +Visual reproduction only (no reference pin) -- spot-check against the archived +``gp_posterior_isotherms_15.png``. """ import os diff --git a/paper/figures/fig08_phase_diagram.py b/paper/figures/fig08_phase_diagram.py index 97ef522..d5a7b0e 100644 --- a/paper/figures/fig08_phase_diagram.py +++ b/paper/figures/fig08_phase_diagram.py @@ -1,8 +1,7 @@ """Fig 8 -- T-x-y phase diagram (bubble/dew point) and liquid/vapor equilibrium curve. -Ported from ``new_phase_diagram.py``'s ``PhaseDiagram.plot_phase_diagram``. Two data -sources, reused per the Phase 7 approach (reuse the ported physics; don't re-run the -stochastic adaptive loop): +Two data sources, combining reused physics rather than re-running the stochastic +adaptive loop: - **Surrogate ensemble** (many thin lines): the ARCHIVED ``phase_diagram_15`` -- bubble/dew points from the paper's actual 15-iteration GP posterior samples. Not diff --git a/paper/figures/fig09_mccabe_thiele.py b/paper/figures/fig09_mccabe_thiele.py index 38cbcba..9c7a693 100644 --- a/paper/figures/fig09_mccabe_thiele.py +++ b/paper/figures/fig09_mccabe_thiele.py @@ -1,16 +1,14 @@ """Fig 9 -- McCabe-Thiele distillation column diagram (Wilson vs. GP surrogate). -Ported from ``solve_distillation.py`` (called via ``run_example.py``), entirely via -``examples/vle_distillation``'s ``phase_diagram``/``equilibrium``/``distillation`` -modules -- this module only adds the plotting + LHS/GP setup on top. +Built entirely on ``examples/vle_distillation``'s ``phase_diagram``/``equilibrium``/ +``distillation`` modules -- this module only adds the plotting + LHS/GP setup on top. ``wilson_column()``/``surrogate_column()`` are the exact same recompute ``tests/regression/test_mccabe_thiele.py``'s gated stage-table regression pins against -``paper/reference/mccabe_thiele_stages.json`` (moved here, Phase 7, so the figure and the -test share one implementation instead of two copies -- the test imports these -functions rather than redefining them). See that test's module docstring for why the -surrogate uses a fresh 30-point LHS + MLE fit rather than the full 15-iteration -adaptive loop. +``paper/reference/mccabe_thiele_stages.json`` -- defined here so the figure and the test +share one implementation instead of two copies (the test imports these functions rather +than redefining them). See that test's module docstring for why the surrogate uses a +fresh 30-point LHS + MLE fit rather than the full 15-iteration adaptive loop. """ import os diff --git a/paper/figures/fig10_traces.py b/paper/figures/fig10_traces.py index 7e75f4a..56c1592 100644 --- a/paper/figures/fig10_traces.py +++ b/paper/figures/fig10_traces.py @@ -1,8 +1,7 @@ """Fig 10 -- HMC trace plots + R-hat/ESS convergence diagnostics. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_all_traces`` (``fxns/plot_res.py``'s -``-m all_traces`` mode). One subplot per kernel hyperparameter, one line per HMC -chain, annotated with that parameter's R-hat and ESS. +One subplot per kernel hyperparameter, one line per HMC chain, annotated with that +parameter's R-hat and ESS. Quantitatively pinned: ``diagnostics()`` returns the exact (rhat, ess) arrays checked against ``paper/reference/hmc_diagnostics.json`` by the gated regression test. diff --git a/paper/figures/fig11_marginals.py b/paper/figures/fig11_marginals.py index 16a1f0d..6caf4bc 100644 --- a/paper/figures/fig11_marginals.py +++ b/paper/figures/fig11_marginals.py @@ -1,10 +1,8 @@ """Fig 11 -- 1-D marginal posterior distributions of the GP kernel hyperparameters. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_marginals`` (``fxns/plot_res.py``'s -``-m marginals`` mode). Visual reproduction only (no reference pin -- the underlying -posterior summary IS pinned quantitatively, in -``paper/reference/hyperparameter_posterior.json``) -- spot-check against the archived -``marginals_15.png``. +Visual reproduction only (no reference pin -- the underlying posterior summary IS +pinned quantitatively, in ``paper/reference/hyperparameter_posterior.json``) -- +spot-check against the archived ``marginals_15.png``. """ import os diff --git a/paper/figures/fig12_joint_marginals.py b/paper/figures/fig12_joint_marginals.py index 20843ac..5f692da 100644 --- a/paper/figures/fig12_joint_marginals.py +++ b/paper/figures/fig12_joint_marginals.py @@ -1,11 +1,10 @@ """Fig 12 -- joint (pairwise) posterior distributions of the GP kernel hyperparameters. -Ported from ``fxns/mcmc_plotter.py``'s ``plot_joint_marginals`` (``fxns/plot_res.py``'s -``-m joint_marginals`` mode), simplified: a lower-triangular grid of hexbin plots with -the MAP point and a 95% credible-interval box per pair -- dropped the KDE contour -overlay (a purely visual smoothing detail, not the quantitative content) for a -pragmatic reproduction. Visual reproduction only (no reference pin) -- spot-check against -the archived ``joint_marginals_15.png``. +A lower-triangular grid of hexbin plots with the MAP point and a 95% credible-interval +box per pair -- dropped the paper's KDE contour overlay (a purely visual smoothing +detail, not the quantitative content) for a pragmatic reproduction. Visual +reproduction only (no reference pin) -- spot-check against the archived +``joint_marginals_15.png``. """ import os diff --git a/paper/full_reproduction.py b/paper/full_reproduction.py index 551e6dc..feb2d47 100644 --- a/paper/full_reproduction.py +++ b/paper/full_reproduction.py @@ -1,27 +1,25 @@ -"""Phase 9 STEP 2/3 -- from-scratch stochastic reproduction of the full adaptive loop. +"""From-scratch stochastic reproduction of the full adaptive loop. Runs ``bits_for_gaps.sampler.BitsForGaps`` end-to-end against the same Wilson/ -Clapeyron black box, HMC config, and 2-D input space as the paper's published -``less_x_new_manuscript_revisions`` run (Jones & Dowling 2026) -- see the private old -repo's ``tests/less_x_new_manuscript_revisions.py`` and -``examples/vle_distillation/run_case_study.py`` for the paper-trail on every constant -below. +Clapeyron black box, HMC config, and 2-D input space as the paper's published run +(Jones & Dowling 2026) -- see ``examples/vle_distillation/run_case_study.py`` for the +same configuration at a shorter, demo-sized iteration count. This is a ONE-TIME validation exercise, not part of the regression suite: results are stochastic (HMC + posterior-mixture sampling both draw randomness TensorFlow does not let us seed bitwise -- see ``bits_for_gaps.mixture``'s module docstring) and are -expected to be *qualitatively*, not bitwise, consistent with the paper. It is -hours-long (15 outer iterations, each running a 4-chain/5000-sample HMC fit plus a -50-draw full-grid posterior-predictive diagnostic) -- run it in the background. +expected to be *qualitatively*, not bitwise, consistent with the paper. Documented +runtime is ~25-30 minutes on a laptop (15 outer iterations, each running a +4-chain/5000-sample HMC fit plus a 50-draw full-grid posterior-predictive diagnostic). -All artifacts go to ``--out-dir`` (default ``results_remaked/phase9_fullrun/``, +All artifacts go to ``--out-dir`` (default ``results_remaked/full_reproduction/``, already gitignored) -- nothing this script produces is committed; only the numbers in its ``full_run_summary.json`` feed the write-up appended to ``paper/REPRODUCTION.md``. Usage:: export PYTHON_JULIACALL_HANDLE_SIGNALS=yes - python paper/full_reproduction.py --out-dir results_remaked/phase9_fullrun + python paper/full_reproduction.py --out-dir results_remaked/full_reproduction """ import argparse @@ -43,8 +41,7 @@ # This script dispatches thousands of tiny TF ops (one predict_f/predict_f_samples # call per bisection step, per z, per posterior draw). On macOS, TF eager's default # multi-threaded op dispatch spends most wall-clock time on thread wake-up/ -# coordination for ops this small -- single-threading it removed a >10x slowdown -# confirmed during Phase 9b's investigation (see paper/PHASE9B_INVESTIGATION.md). +# coordination for ops this small -- single-threading it removes a >10x slowdown. import tensorflow as tf tf.config.threading.set_intra_op_parallelism_threads(1) @@ -82,9 +79,8 @@ def _predict_split(record, XGP, seed, size): Uses a 15-component hyperparameter-posterior subset (``noGaussians``, the same mixture size the acquisition function itself uses) rather than the paper's own - ``train_test_split_proh.py``'s dedicated 500-sample subset -- a documented - simplification (see ``paper/REPRODUCTION.md``'s Phase 9 section), not a bitwise - match. + dedicated 500-sample subset -- a documented simplification (see + ``paper/REPRODUCTION.md``), not a bitwise match. """ yGP_draws = mixture.sample_gp_posterior_mixture( record.trace, record.GPmodel, XGP, seed=seed, size=size @@ -115,7 +111,7 @@ def run(out_dir, n_init=N_INIT, n_test=N_TEST, n_iters=N_ITERS, seed=SEED): input_transform=INPUT_TRANSFORM, output_transform=OUTPUT_TRANSFORM, iters=n_iters, - exp_name="phase9_fullrun", + exp_name="full_reproduction", ) bfg.seed = seed bfg.noSamples, bfg.noBurnIn = 5000, 0 @@ -144,12 +140,12 @@ def run(out_dir, n_init=N_INIT, n_test=N_TEST, n_iters=N_ITERS, seed=SEED): # # MUST run before the test-RMSE loop below: `_predict_split` (via # `mixture.sample_gp_posterior_mixture`) mutates `record.GPmodel.kernel` in place, - # and `history.last.GPmodel` is the SAME object -- running the RMSE loop first - # once left the kernel at an arbitrary leftover single-hyperparameter state and - # produced a spurious non-converging column (see paper/PHASE9B_INVESTIGATION.md). - # Uses `surrogate_gamma_averaged` (matches the paper's own `new_phase_diagram.py` - # construction) rather than `surrogate_gamma`'s single point-estimate, for the same - # reason: robust to any one hyperparameter draw being atypical. + # and `history.last.GPmodel` is the SAME object -- running the RMSE loop first would + # leave the kernel at an arbitrary leftover single-hyperparameter state and produce + # a spurious non-converging column. Uses `surrogate_gamma_averaged` (matches the + # paper's own posterior-averaging construction) rather than `surrogate_gamma`'s + # single point-estimate, for the same reason: robust to any one hyperparameter draw + # being atypical. GPmodel = history.last.GPmodel trace_final = history.last.trace z_grid = np.linspace(0.0, 1.0, Z_GRID_SIZE) @@ -225,8 +221,8 @@ def surrogate_gamma_fn(z, T): parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) parser.add_argument( "--out-dir", - default=os.path.join(REPO_ROOT, "results_remaked", "phase9_fullrun"), - help="Gitignored output directory (default: results_remaked/phase9_fullrun)", + default=os.path.join(REPO_ROOT, "results_remaked", "full_reproduction"), + help="Gitignored output directory (default: results_remaked/full_reproduction)", ) parser.add_argument("--n-iters", type=int, default=N_ITERS) args = parser.parse_args() diff --git a/paper/phase9_validation/full_run_summary.json b/paper/phase9_validation/full_run_summary.json deleted file mode 100644 index 1583bd8..0000000 --- a/paper/phase9_validation/full_run_summary.json +++ /dev/null @@ -1,261 +0,0 @@ -{ - "config": { - "n_init": 10, - "n_test": 10, - "n_iters": 15, - "seed": 10 - }, - "elapsed_hours": 0.42879568305280474, - "final_n_points": 25, - "rhat": { - "1": [ - 1.0036140953257404, - 1.0050428434175722, - 1.006526460705449 - ], - "2": [ - 1.0046231615266175, - 1.005771353983787, - 1.0071164971603306 - ], - "3": [ - 1.0047217352461284, - 1.0059574498192365, - 1.0074044664552049 - ], - "4": [ - 1.0047890797776038, - 1.0063938572409117, - 1.0076646237266043 - ], - "5": [ - 1.0047474666774319, - 1.0065378705940304, - 1.007563513314564 - ], - "6": [ - 1.0047862876662261, - 1.0066441176163075, - 1.007573763207251 - ], - "7": [ - 1.0048650408877222, - 1.0067026725516341, - 1.0077941066760632 - ], - "8": [ - 1.0047428856413676, - 1.0067513991435542, - 1.0080458649468413 - ], - "9": [ - 1.0048227679384483, - 1.006854434926923, - 1.0082155273821978 - ], - "10": [ - 1.004927750960709, - 1.0069726234070515, - 1.008380251188962 - ], - "11": [ - 1.0050405058725085, - 1.0070359599434289, - 1.0084159940708202 - ], - "12": [ - 1.0049926736237276, - 1.0070409590795677, - 1.0085256282119284 - ], - "13": [ - 1.0050708004812319, - 1.0071454712457015, - 1.0086826869792733 - ], - "14": [ - 1.0051262418198017, - 1.007237618927919, - 1.0087852371546187 - ], - "15": [ - 1.0052293412023583, - 1.0072992371956864, - 1.008792527111599 - ] - }, - "ess": { - "1": [ - 1652.2607758751883, - 2880.944135885453, - 735.5361631277776 - ], - "2": [ - 1057.978462033747, - 2734.3433278451816, - 729.9962584307982 - ], - "3": [ - 1445.597718344832, - 2603.6971765532726, - 713.4010282239996 - ], - "4": [ - 1474.5368112305632, - 2569.4887345210623, - 700.2304181588405 - ], - "5": [ - 1481.2259305458538, - 2555.2786310414713, - 694.852884242609 - ], - "6": [ - 1484.3239529572875, - 2537.4204722546365, - 687.2051757438751 - ], - "7": [ - 1481.0455710471356, - 2540.8985748812984, - 678.9936778477942 - ], - "8": [ - 1506.0123091965092, - 2539.679731612028, - 677.0356779877467 - ], - "9": [ - 1498.3246470261593, - 2540.329496411773, - 673.1378008821556 - ], - "10": [ - 1473.8113779222606, - 2488.201686085358, - 669.8741841564987 - ], - "11": [ - 1470.4018882434568, - 2477.317556932017, - 666.8990945263519 - ], - "12": [ - 1491.7915237236152, - 2479.7030417379724, - 661.6058414058225 - ], - "13": [ - 1485.4964733491674, - 2472.948650202977, - 657.3969296908776 - ], - "14": [ - 1482.6335354442174, - 2455.184231118362, - 654.0276615458931 - ], - "15": [ - 1468.2900559280745, - 2428.0938804768984, - 653.1451900064535 - ] - }, - "max_entropy": { - "1": 1.4576908211674722, - "2": 1.2470067603579276, - "3": 0.6695318507453859, - "4": 0.3737886196772696, - "5": 0.2735756995129844, - "6": 0.18311675762785679, - "7": 0.11596694633588284, - "8": 0.04772962603132569, - "9": 0.02884167541979843, - "10": 0.015909235519620757, - "11": 0.003844712139748019, - "12": -0.058498820185680755, - "13": -0.13720419745170287, - "14": -0.2181830630917866, - "15": -0.2243559774314241 - }, - "test_rmse": { - "1": 4.752414536434205, - "2": 2.450472819906951, - "3": 1.532971382091772, - "4": 1.469471160946046, - "5": 1.5436366576387646, - "6": 1.222557829563998, - "7": 1.1830905687783326, - "8": 1.239949315322786, - "9": 0.9570925772303377, - "10": 0.9645163441554637, - "11": 1.1261129711277975, - "12": 1.060907866385932, - "13": 1.0249289661537424, - "14": 0.9415996416938199, - "15": 0.9931192775349461 - }, - "hyperparameter_posterior_final": { - "mean": [ - 1.3564535683707466, - 0.8623918547679581, - 3.195015639568722 - ], - "median": [ - 1.2861863176137853, - 0.8194868463441703, - 3.0440724074555088 - ], - "std": [ - 0.39761949918655504, - 0.2814676958684645, - 1.1350095402382614 - ] - }, - "column_wilson_converged": true, - "column_surrogate_converged": true, - "column_wilson_stages": [ - { - "stage": 1, - "liquid": 0.22649452787395025, - "vapor": 0.43 - }, - { - "stage": 2, - "liquid": 0.047286250942829886, - "vapor": 0.328247263936976 - }, - { - "stage": 3, - "liquid": 0.021574162657992516, - "vapor": 0.23864312547141986 - }, - { - "stage": 4, - "liquid": 0.01, - "vapor": 0.139891841726997 - } - ], - "column_surrogate_stages": [ - { - "stage": 1, - "liquid": 0.25371911983626494, - "vapor": 0.43 - }, - { - "stage": 2, - "liquid": 0.0432952536085202, - "vapor": 0.34185955991813133 - }, - { - "stage": 3, - "liquid": 0.019708958411344915, - "vapor": 0.23664762680427678 - }, - { - "stage": 4, - "liquid": 0.01, - "vapor": 0.1434869076815258 - } - ] -} \ No newline at end of file diff --git a/paper/phase9_validation/full_run_summary_pre_phase9b_fix.json b/paper/phase9_validation/full_run_summary_pre_phase9b_fix.json deleted file mode 100644 index 3f87236..0000000 --- a/paper/phase9_validation/full_run_summary_pre_phase9b_fix.json +++ /dev/null @@ -1,261 +0,0 @@ -{ - "config": { - "n_init": 10, - "n_test": 10, - "n_iters": 15, - "seed": 10 - }, - "elapsed_hours": 0.4165407505962584, - "final_n_points": 25, - "rhat": { - "1": [ - 1.0036140953257404, - 1.0050428434175722, - 1.006526460705449 - ], - "2": [ - 1.0046231615266175, - 1.005771353983787, - 1.0071164971603306 - ], - "3": [ - 1.0047217352461284, - 1.0059574498192365, - 1.0074044664552049 - ], - "4": [ - 1.0047890797776038, - 1.0063938572409117, - 1.0076646237266043 - ], - "5": [ - 1.0047474666774319, - 1.0065378705940304, - 1.007563513314564 - ], - "6": [ - 1.0047862876662261, - 1.0066441176163075, - 1.007573763207251 - ], - "7": [ - 1.0048650408877222, - 1.0067026725516341, - 1.0077941066760632 - ], - "8": [ - 1.0047428856413676, - 1.0067513991435542, - 1.0080458649468413 - ], - "9": [ - 1.0048227679384483, - 1.006854434926923, - 1.0082155273821978 - ], - "10": [ - 1.004927750960709, - 1.0069726234070515, - 1.008380251188962 - ], - "11": [ - 1.0050405058725085, - 1.0070359599434289, - 1.0084159940708202 - ], - "12": [ - 1.0049926736237276, - 1.0070409590795677, - 1.0085256282119284 - ], - "13": [ - 1.0050708004812319, - 1.0071454712457015, - 1.0086826869792733 - ], - "14": [ - 1.0051262418198017, - 1.007237618927919, - 1.0087852371546187 - ], - "15": [ - 1.0052293412023583, - 1.0072992371956864, - 1.008792527111599 - ] - }, - "ess": { - "1": [ - 1652.2607758751883, - 2880.944135885453, - 735.5361631277776 - ], - "2": [ - 1057.978462033747, - 2734.3433278451816, - 729.9962584307982 - ], - "3": [ - 1445.597718344832, - 2603.6971765532726, - 713.4010282239996 - ], - "4": [ - 1474.5368112305632, - 2569.4887345210623, - 700.2304181588405 - ], - "5": [ - 1481.2259305458538, - 2555.2786310414713, - 694.852884242609 - ], - "6": [ - 1484.3239529572875, - 2537.4204722546365, - 687.2051757438751 - ], - "7": [ - 1481.0455710471356, - 2540.8985748812984, - 678.9936778477942 - ], - "8": [ - 1506.0123091965092, - 2539.679731612028, - 677.0356779877467 - ], - "9": [ - 1498.3246470261593, - 2540.329496411773, - 673.1378008821556 - ], - "10": [ - 1473.8113779222606, - 2488.201686085358, - 669.8741841564987 - ], - "11": [ - 1470.4018882434568, - 2477.317556932017, - 666.8990945263519 - ], - "12": [ - 1491.7915237236152, - 2479.7030417379724, - 661.6058414058225 - ], - "13": [ - 1485.4964733491674, - 2472.948650202977, - 657.3969296908776 - ], - "14": [ - 1482.6335354442174, - 2455.184231118362, - 654.0276615458931 - ], - "15": [ - 1468.2900559280745, - 2428.0938804768984, - 653.1451900064535 - ] - }, - "max_entropy": { - "1": 1.4576908211674722, - "2": 1.2470067603579276, - "3": 0.6695318507453859, - "4": 0.3737886196772696, - "5": 0.2735756995129844, - "6": 0.18311675762785679, - "7": 0.11596694633588284, - "8": 0.04772962603132569, - "9": 0.02884167541979843, - "10": 0.015909235519620757, - "11": 0.003844712139748019, - "12": -0.058498820185680755, - "13": -0.13720419745170287, - "14": -0.2181830630917866, - "15": -0.2243559774314241 - }, - "test_rmse": { - "1": 4.3365979407929265, - "2": 2.480115588840578, - "3": 1.7566801151333729, - "4": 1.4882478796019056, - "5": 1.314304949940719, - "6": 1.1773933074892475, - "7": 1.0748144961007042, - "8": 1.229776327155701, - "9": 1.0757636306596277, - "10": 0.9285915062278457, - "11": 0.9204235514944763, - "12": 1.0730573622067199, - "13": 0.9549790130200277, - "14": 0.9168553362113754, - "15": 0.8870659763377472 - }, - "hyperparameter_posterior_final": { - "mean": [ - 1.3564535683707466, - 0.8623918547679581, - 3.195015639568722 - ], - "median": [ - 1.2861863176137853, - 0.8194868463441703, - 3.0440724074555088 - ], - "std": [ - 0.39761949918655504, - 0.2814676958684645, - 1.1350095402382614 - ] - }, - "column_wilson_converged": true, - "column_surrogate_converged": false, - "column_wilson_stages": [ - { - "stage": 1, - "liquid": 0.22649452787395025, - "vapor": 0.43 - }, - { - "stage": 2, - "liquid": 0.047286250942829886, - "vapor": 0.328247263936976 - }, - { - "stage": 3, - "liquid": 0.021574162657992516, - "vapor": 0.23864312547141986 - }, - { - "stage": 4, - "liquid": 0.01, - "vapor": 0.139891841726997 - } - ], - "column_surrogate_stages": [ - { - "stage": 1, - "liquid": 0.2671311288617889, - "vapor": 0.4567416616275945 - }, - { - "stage": 2, - "liquid": 1.725805727168204, - "vapor": 0.37756564604512194 - }, - { - "stage": 3, - "liquid": 0.16648026026257204, - "vapor": 1.1491305852571772 - }, - { - "stage": 4, - "liquid": 0.017465321298021103, - "vapor": 0.19701226994025983 - } - ] -} \ No newline at end of file diff --git a/paper/phase9_validation/phase9b_curve_comparison.png b/paper/phase9_validation/phase9b_curve_comparison.png deleted file mode 100644 index cc843cb..0000000 Binary files a/paper/phase9_validation/phase9b_curve_comparison.png and /dev/null differ diff --git a/paper/phase9_validation/phase_diagram_fresh.png b/paper/phase9_validation/phase_diagram_fresh.png deleted file mode 100644 index 65f6707..0000000 Binary files a/paper/phase9_validation/phase_diagram_fresh.png and /dev/null differ diff --git a/paper/phase9_validation/rmse_and_entropy.png b/paper/phase9_validation/rmse_and_entropy.png deleted file mode 100644 index abdbf0e..0000000 Binary files a/paper/phase9_validation/rmse_and_entropy.png and /dev/null differ diff --git a/paper/reference/README.md b/paper/reference/README.md index b302bdd..a8a81e1 100644 --- a/paper/reference/README.md +++ b/paper/reference/README.md @@ -3,22 +3,21 @@ Small, version-controlled JSON snapshots of the paper's key scalar results, extracted from the **archived published run** (`results/less_x_new_manuscript_revisions`, **iteration 15** — the run whose R-hat/ESS match paper Fig 10 exactly). The bulk -`results/` archive is *not* committed (it goes to Zenodo, see `paper/DATA.md`); only -these scalars live in the repo. +`results/` archive is *not* committed (no Zenodo deposit either — see `paper/DATA.md`); +only these scalars live in the repo. They serve two purposes: 1. **Characterization / regression** — `tests/regression/` reads these files and pins them against the values reported in the paper (locks them against corruption). -2. **Reproduction diff** — Phase 7's `paper/reproduce.py` regenerates the figures - through the clean API and diffs against these reference targets (with the stated - tolerances). +2. **Reproduction diff** — `paper/reproduce.py` regenerates the figures through the + clean API and diffs against these reference targets (with the stated tolerances). | File | Paper ref | Contents | Source | |---|---|---|---| | `hmc_diagnostics.json` | Fig 10 | HMC R-hat & ESS for the 3 hyperparameters | `rhat_value_15.txt`, `ess_value_15.txt` | | `fig5_error_metrics.json` | Fig 5 | Train/test RMSE & MAE distributions at iters 1 & 15 | `gp_predict_{train,test}_{1,15}` vs `activity_data_1`/`activity_test_points` | | `hyperparameter_posterior.json` | Fig 10 | Posterior mean/median/std/quantiles of the kernel hyperparameters | `param_posterior_samples_15` | -| `mccabe_thiele_stages.json` | Fig 9c | Distillation stage liquid/vapor mole fractions (Wilson vs surrogate) | paper Fig 9c / `run_example.py` | +| `mccabe_thiele_stages.json` | Fig 9c | Distillation stage liquid/vapor mole fractions (Wilson vs surrogate) | transcribed from paper Fig 9c; recomputed by `paper/figures/fig09_mccabe_thiele.py` | **Regenerating** (needs read access to the old-repo archive; pure NumPy): @@ -27,4 +26,4 @@ python paper/extract_reference.py ``` produces all but the McCabe-Thiele table (which is transcribed from the paper — its -recompute needs the Julia VLE backend ported in Phase 6). +recompute needs the Julia VLE backend). diff --git a/paper/reproduce.py b/paper/reproduce.py index e93310a..f3c8610 100644 --- a/paper/reproduce.py +++ b/paper/reproduce.py @@ -7,18 +7,17 @@ expensive, and not what "reproduce the figures" requires) -- it loads what that loop already produced. -As of Phase 9, the default data source is the curated, committed ``paper/data/`` -subset (~16 MB -- exactly the files the figures read, copied from the private -archive; see ``paper/data/README.md``), so this runs with **no private-archive -access** out of the box. Pass ``--archive`` (or set ``BFG_ARCHIVE_DIR``) to point at -the full private archive instead -- e.g. to regenerate Fig 3/4 with more than the 6 -committed early-iteration panels, or Fig 6/7 at iterations other than 1/15. +The default data source is the committed ``paper/data/`` directory (~16 MB -- exactly +the files the figures read; see ``paper/data/README.md``), so this runs from a fresh clone +with no extra downloads. Pass ``--archive`` (or set ``BFG_ARCHIVE_DIR``) to point at +another directory of run artifacts with the same layout instead -- for example the output +of ``paper/full_reproduction.py``, to plot your own run rather than the published one. Usage:: export PYTHON_JULIACALL_HANDLE_SIGNALS=yes # macOS, only needed for Fig 8/9 python paper/reproduce.py # uses the committed paper/data/ - python paper/reproduce.py --archive /path/to/less_x_new_manuscript_revisions + python paper/reproduce.py --archive /path/to/other/run/artifacts python paper/reproduce.py --figures 5 8 9 10 # only regenerate a subset Output goes to ``--out-dir`` (default: ``results_remaked/``, already gitignored) -- @@ -36,8 +35,8 @@ if EXAMPLES_DIR not in sys.path: sys.path.insert(0, EXAMPLES_DIR) -# Phase 9: default to the curated, committed subset (no private-archive access -# needed); override with --archive/$BFG_ARCHIVE_DIR to use the full private archive. +# Default to the committed paper/data/ directory; override with +# --archive/$BFG_ARCHIVE_DIR to plot a different run's artifacts. DEFAULT_ARCHIVE = os.path.join(REPO_ROOT, "paper", "data") DEFAULT_OUT_DIR = os.path.join(REPO_ROOT, "results_remaked") @@ -65,8 +64,8 @@ def main(argv=None): "--archive", default=os.environ.get("BFG_ARCHIVE_DIR", DEFAULT_ARCHIVE), help="Path to the plot-input data (default: $BFG_ARCHIVE_DIR or the " - f"committed paper/data/ subset, {DEFAULT_ARCHIVE}); point this at the " - "full private archive for iterations/figures beyond the curated subset", + f"committed paper/data/ directory, {DEFAULT_ARCHIVE}); point this at another " + "run's artifacts (same layout) to plot that run instead", ) parser.add_argument( "--out-dir", default=DEFAULT_OUT_DIR, help=f"Output directory (default: {DEFAULT_OUT_DIR})" diff --git a/pyproject.toml b/pyproject.toml index a79e04d..5c9d839 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,8 +28,9 @@ classifiers = [ "Operating System :: OS Independent", ] -# Frozen reproduction baseline (see REFACTOR_PLAN.md decision 4). -# tensorflow-macos on Apple Silicon; plain tensorflow elsewhere. +# Frozen reproduction baseline: an exact, verified-working stack rather than floating +# version ranges, since GPflow's TensorFlow dependency makes casual version bumps +# risky. tensorflow-macos on Apple Silicon; plain tensorflow elsewhere. dependencies = [ "numpy>=1.26,<2", "scipy>=1.13,<1.14", @@ -46,7 +47,7 @@ dependencies = [ # package whose TF-backed modules (kernels/means/gp/mixture/acquisition/sampler) all # fail with `ModuleNotFoundError: No module named 'pkg_resources'` -- setuptools' # own deprecation warning recommends exactly this pin. Revisit when the frozen - # GPflow/TF stack is modernized (REFACTOR_PLAN.md section 6). + # GPflow/TF stack is modernized. "setuptools<81", ] @@ -75,7 +76,7 @@ packages = ["src/bits_for_gaps"] [tool.hatch.build.targets.sdist] # hatchling's default sdist includes the whole repo (paper/data's ~16 MB, tests/, # examples/, docs/) -- restrict it to the package + the files a source release needs -# (REFACTOR_PLAN.md §7 decision 3: examples/paper/docs stay repo-only, not shipped). +# (examples/paper/docs stay repo-only, not shipped -- see docs/installation.md). # `only-include` (not `include`, which ADDS to hatchling's default vcs-based file set # rather than replacing it) is what actually makes this an allowlist. only-include = [ @@ -109,10 +110,9 @@ ignore = ["B023"] # on macOS or juliacall SIGBUSes. This is a function-local import ruff's isort # wouldn't touch, but ignore E402 here defensively in case rule selection changes. "examples/vle_distillation/activity_model.py" = ["I001", "E402"] -# - full_reproduction.py sets TF's single-threaded config (a >10x wall-clock fix -# found during Phase 9b's investigation) before importing bits_for_gaps/ -# vle_distillation, so those imports don't trigger TF's default multi-threaded -# init first. +# - full_reproduction.py sets TF's single-threaded config (a >10x wall-clock fix for +# this script's many small ops) before importing bits_for_gaps/vle_distillation, so +# those imports don't trigger TF's default multi-threaded init first. "paper/full_reproduction.py" = ["I001", "E402"] [tool.ruff.format] diff --git a/src/bits_for_gaps/__init__.py b/src/bits_for_gaps/__init__.py index e52a5e9..632fe55 100644 --- a/src/bits_for_gaps/__init__.py +++ b/src/bits_for_gaps/__init__.py @@ -30,7 +30,7 @@ "AnisotropicSE": ("kernels", "AnisotropicSE"), "FixedInverseMean": ("means", "FixedInverseMean"), "adaptiveEntropy": ("sampler", "adaptiveEntropy"), - "BitsForGaps": ("sampler", "BitsForGaps"), # target public API (Phase 4) + "BitsForGaps": ("sampler", "BitsForGaps"), # friendly-named public API facade } diff --git a/src/bits_for_gaps/_util.py b/src/bits_for_gaps/_util.py index 3139b6d..a0463f4 100644 --- a/src/bits_for_gaps/_util.py +++ b/src/bits_for_gaps/_util.py @@ -1,4 +1,4 @@ -"""Small array helpers. Moved from the paper code's ``fxns/util.py``.""" +"""Small array helpers.""" from __future__ import annotations diff --git a/src/bits_for_gaps/acquisition.py b/src/bits_for_gaps/acquisition.py index e0fefc8..6dfb122 100644 --- a/src/bits_for_gaps/acquisition.py +++ b/src/bits_for_gaps/acquisition.py @@ -1,36 +1,29 @@ """Entropy-maximization acquisition function (2nd-order Taylor estimator). -Moved from the paper code's ``driver_new.py`` (``adaptiveEntropy.entropy_objective`` / -``gen_entropy_surface_data_2D`` / ``optimize_2D``). Pure functions over explicit -arguments -- no disk I/O; ``entropy.py`` provides the underlying mixture-entropy math. - -Phase 5: ``entropy_objective`` and ``optimize`` (was ``optimize_2D``) are dimension- -general -- kernel hyperparameters are assigned via ``kernels.assign_hyperparameters`` -(the canonical-order contract, not hardcoded names), and ``optimize``'s Sobol restarts -and bounds scale with ``len(x_bounds)``. ``entropy_surface_2D`` stays 2-D-only: a dense -grid is exponential in d, and it is a visualization diagnostic that does not feed the -acquisition -- it raises a clear error for d != 2 rather than silently degrading. - -Phase 9c: ``entropy_objective`` is called many times per ``optimize``/ -``entropy_surface_2D`` call (once per Sobol restart x scipy.optimize.minimize -iteration, or once per grid point), each time reassigning ``GPmodel.kernel``'s -hyperparameters -- it used to leave the kernel at whichever sample the *last* such call -happened to use, not any meaningful state. Since ``sampler.py``'s ``run()`` calls -``optimize``/``entropy_surface_2D`` on the same ``GPmodel`` it then stores in -``IterationRecord.GPmodel`` (and optionally checkpoints), every iteration's returned -model used to carry this arbitrary leftover state -- the same class of footgun Phase 9b -found in ``mixture.sample_gp_posterior_mixture`` (see -``paper/PHASE9B_INVESTIGATION.md``). Now saves/restores the kernel around each call, so -``GPmodel`` is unchanged after ``entropy_objective`` returns -- behavior-preserving for -the entropy value itself (computed before the restore). - -Phase 9d: the paper derives TWO entropy estimators -- the 2nd-order Taylor -approximation (``entropy.second_order_entropy``, the one that actually drove -acquisition in the paper and remains the default here) and a closed-form lower bound -(``entropy.entropy_lower_bound``, paper Theorem/SI-2), implemented and unit-tested -since Phase 2 but never wired up as a *usable* acquisition objective. ``objective`` -selects between them (default ``"taylor"``, so existing behavior/baselines/reference -values are unchanged unless a caller explicitly opts into ``"lower_bound"``). +Pure functions over explicit arguments -- no disk I/O; ``entropy.py`` provides the +underlying mixture-entropy math. + +``entropy_objective`` and ``optimize`` are dimension-general -- kernel hyperparameters +are assigned via ``kernels.assign_hyperparameters`` (the canonical-order contract, not +hardcoded names), and ``optimize``'s Sobol restarts and bounds scale with +``len(x_bounds)``. ``entropy_surface_2D`` stays 2-D-only: a dense grid is exponential in +d, and it is a visualization diagnostic that does not feed the acquisition -- it raises +a clear error for d != 2 rather than silently degrading. + +CAUTION -- kernel mutation: ``entropy_objective`` is called many times per +``optimize``/``entropy_surface_2D`` call (once per Sobol restart x +``scipy.optimize.minimize`` iteration, or once per grid point), each time reassigning +``GPmodel.kernel``'s hyperparameters. It saves the kernel's hyperparameters before each +call and restores them in a ``finally``, so ``GPmodel`` is unchanged after +``entropy_objective`` returns regardless of how it exits -- ``sampler.py``'s ``run()`` +relies on this, since it calls ``optimize``/``entropy_surface_2D`` on the same +``GPmodel`` it then stores in ``IterationRecord.GPmodel`` (and optionally checkpoints). + +The paper derives TWO entropy estimators -- the 2nd-order Taylor approximation +(``entropy.second_order_entropy``, the one that actually drove acquisition in the paper +and remains the default here) and a closed-form lower bound +(``entropy.entropy_lower_bound``, paper Theorem/SI-2). ``objective`` selects between +them (default ``"taylor"``). """ from __future__ import annotations @@ -45,8 +38,8 @@ from . import entropy as max_ent_design from . import kernels -# Phase 9d: the two paper-derived entropy estimators, selectable via `objective` below. -# Both accept (weights, means, covs_or_variances) positionally -- entropy_lower_bound's +# The two paper-derived entropy estimators, selectable via `objective` below. Both +# accept (weights, means, covs_or_variances) positionally -- entropy_lower_bound's # univariate-GMM assumption holds for every call site below, since entropy_objective # always evaluates the entropy of the GP's (scalar) OUTPUT at one point, never a # multi-point joint, so `means`/`variances` are always 1-D scalar arrays regardless of @@ -121,11 +114,9 @@ def entropy_surface_2D( ) -> np.ndarray: """Entropy field over a 2-D grid spanning ``x_bounds`` at ``mesh`` points per dim. - Moved from ``adaptiveEntropy.gen_entropy_surface_data_2D``. - - 2-D-ONLY VISUALIZATION DIAGNOSTIC (Phase 5): a dense grid is exponential in the - input dimension d, so this is intentionally not generalized to N-D. It does not - feed the acquisition -- see ``optimize`` for the N-D-general acquisition path. + 2-D-ONLY VISUALIZATION DIAGNOSTIC: a dense grid is exponential in the input + dimension d, so this is intentionally not generalized to N-D. It does not feed the + acquisition -- see ``optimize`` for the N-D-general acquisition path. Returns ------- @@ -167,17 +158,16 @@ def optimize( ) -> OptimizeResult: """Multistart optimization of the entropy objective over ``x_bounds`` (N-D). - Moved from ``adaptiveEntropy.optimize_2D``; generalized (Phase 5) to arbitrary - input dimension ``d = len(x_bounds)`` -- this is the acquisition path an N-D run - actually depends on (contrast ``entropy_surface_2D``, a 2-D-only diagnostic). - Sobol-scrambled restarts drawn in the physical space, optimized in GP (transformed) - space. - - For ``d == 2`` this reproduces the pre-Phase-5 ``optimize_2D`` bit-for-bit: the - vectorized bound-scaling below (``lo + x0 * (hi - lo)``) is the same floating-point - operation order as the original per-dimension expression (IEEE 754 addition is - commutative, so ``lo[j] + x0[j] * (hi[j] - lo[j])`` and the original - ``x0[j] * (hi[j] - lo[j]) + lo[j]`` round identically). + Works at arbitrary input dimension ``d = len(x_bounds)`` -- this is the acquisition + path an N-D run actually depends on (contrast ``entropy_surface_2D``, a 2-D-only + diagnostic). Sobol-scrambled restarts drawn in the physical space, optimized in GP + (transformed) space. + + The vectorized bound-scaling below (``lo + x0 * (hi - lo)``) is written in this + operand order deliberately: it matters for bit-exactness against + ``tests/integration/data/synthetic_baseline.json``'s pinned values at ``d == 2`` + (IEEE 754 addition is commutative, so ``x0[j] * (hi[j] - lo[j]) + lo[j]`` would be + equally correct mathematically, but rounds differently). Parameters ---------- diff --git a/src/bits_for_gaps/design.py b/src/bits_for_gaps/design.py index 3b9ed92..84e86bd 100644 --- a/src/bits_for_gaps/design.py +++ b/src/bits_for_gaps/design.py @@ -1,8 +1,6 @@ """Space-filling initial designs over a bounded input space. -Extracted and generalized from the paper code's ``proh_water_class`` (which mixed -design generation with Julia activity-coefficient calls and disk I/O). These are pure, -N-dimensional, and return arrays -- no files, no thermodynamics. +Pure, N-dimensional, and return arrays -- no files, no thermodynamics. """ from __future__ import annotations diff --git a/src/bits_for_gaps/diagnostics.py b/src/bits_for_gaps/diagnostics.py index 40eb0c4..02471cf 100644 --- a/src/bits_for_gaps/diagnostics.py +++ b/src/bits_for_gaps/diagnostics.py @@ -1,8 +1,8 @@ """HMC convergence diagnostics. -Thin wrappers over ``tensorflow_probability.mcmc``, split out of ``driver_new.py``'s -``adaptiveEntropy.run_mcmc`` so ``gp.py``'s HMC driver doesn't need to know the -diagnostics API directly (and so these two lines are independently testable/reusable). +Thin wrappers over ``tensorflow_probability.mcmc``, split out so ``gp.py``'s HMC driver +doesn't need to know the diagnostics API directly (and so these two functions are +independently testable/reusable). """ from __future__ import annotations diff --git a/src/bits_for_gaps/entropy.py b/src/bits_for_gaps/entropy.py index d9514e7..8aeb646 100644 --- a/src/bits_for_gaps/entropy.py +++ b/src/bits_for_gaps/entropy.py @@ -15,8 +15,7 @@ analytic reference for the true mixture entropy. Everything here is pure NumPy/SciPy and dimension-agnostic (univariate or multivariate -components), so it is straightforward to unit-test. Moved verbatim-in-spirit from the -paper code's ``fxns/max_ent_design.py`` (dead commented-out variants removed). +components), so it is straightforward to unit-test. """ from __future__ import annotations @@ -99,10 +98,10 @@ def first_order_entropy_approx( ValueError If the mixture density at some component mean is not positive (e.g. floating- point underflow for very small covariances / far-apart means) -- a data- - dependent runtime condition, not a static invariant, so Phase 9c makes this an - explicit exception rather than a bare ``assert`` (asserts are silently - stripped under ``python -O``, which would let a NaN/garbage entropy value - propagate instead of failing clearly). + dependent runtime condition, not a static invariant, hence an explicit + exception rather than a bare ``assert`` (asserts are silently stripped under + ``python -O``, which would let a NaN/garbage entropy value propagate instead of + failing clearly). """ H = 0 for loopA, wA in enumerate(weights): diff --git a/src/bits_for_gaps/gp.py b/src/bits_for_gaps/gp.py index c8f916e..0ee55da 100644 --- a/src/bits_for_gaps/gp.py +++ b/src/bits_for_gaps/gp.py @@ -1,17 +1,13 @@ """GP construction, log-marginal-likelihood optimization, and HMC sampling. -Moved from the paper code's ``driver_new.py`` (``adaptiveEntropy.build_gp`` / -``maximize_lml`` / ``run_mcmc``). These are pure functions over explicit arguments -- -no disk I/O; the orchestrator (``sampler.py``) decides whether/where to checkpoint. +These are pure functions over explicit arguments -- no disk I/O; the orchestrator +(``sampler.py``) decides whether/where to checkpoint. -Phase 5: ``run_mcmc`` no longer indexes ``GPmodel.trainable_parameters`` by hardcoded -position. It uses ``GPmodel.kernel.hyperparameters`` (see +``run_mcmc`` uses ``GPmodel.kernel.hyperparameters`` (see ``kernels.AnisotropicSE.hyperparameters``) directly as the HMC state -- any kernel that exposes a ``.hyperparameters`` property (a list of ``gpflow.Parameter``, in whatever order that kernel defines as canonical) works here, not just the paper's 3-hyperparameter -2-D kernel. For the paper's kernel this is verified identity-equal to the old -``[trainable_parameters[2], [0], [1]]`` indexing (same Parameter objects, same order), -so d=2 runs are bit-exact with the pre-Phase-5 code. +2-D kernel. """ from __future__ import annotations diff --git a/src/bits_for_gaps/kernels.py b/src/bits_for_gaps/kernels.py index 6c06fc4..518dbca 100644 --- a/src/bits_for_gaps/kernels.py +++ b/src/bits_for_gaps/kernels.py @@ -1,12 +1,11 @@ """GP covariance kernels with hierarchical (prior-bearing) hyperparameters. -Moved from the paper code's ``fxns/my_kermel_fxn.py`` (commented-out variants and the -misspelled module name dropped). The kernel's hyperparameters carry the tfp priors that -make the GP *hierarchical* -- HMC samples these to form the mixture predictive posterior. +The kernel's hyperparameters carry the tfp priors that make the GP *hierarchical* -- HMC +samples these to form the mixture predictive posterior. -Phase 5: generalized to N input dimensions. DESIGN DECISION -- each lengthscale (and the -kernel variance) is its own scalar ``gpflow.Parameter``, not a single vector-valued -Parameter. This is deliberate: the paper's method hinges on per-dimension PRIOR FAMILIES, +Works at N input dimensions. DESIGN DECISION -- each lengthscale (and the kernel +variance) is its own scalar ``gpflow.Parameter``, not a single vector-valued Parameter. +This is deliberate: the paper's method hinges on per-dimension PRIOR FAMILIES, not just per-dimension prior parameters -- ``std_dev`` ~ LogNormal, ``lengthscale_1`` ~ LogNormal, ``lengthscale_2`` ~ Gamma *and unconstrained* (no positivity bijector). A single vector Parameter carries exactly one prior distribution and one transform for the @@ -42,7 +41,7 @@ class AnisotropicSE(gpflow.kernels.Kernel): """Anisotropic squared-exponential kernel with per-dimension lengthscales. - Generalized (Phase 5) to N input dimensions. Each lengthscale -- and the kernel + Works at N input dimensions. Each lengthscale -- and the kernel variance -- is its own ``gpflow.Parameter``, carrying its own prior distribution and (optionally) its own bijector, so per-dimension prior FAMILIES are supported (see the module docstring for why this matters). @@ -196,20 +195,15 @@ def assign_hyperparameters(kernel: gpflow.kernels.Kernel, values: Sequence[float ``.hyperparameters`` property (a list of ``gpflow.Parameter``), not just ``AnisotropicSE``. - Phase 9d: this is called deep inside ``mixture.py``/``acquisition.py``'s hot loops - to replay one posterior/mixture-component sample at a time. An extreme outlier - sample -- most plausibly from ``lengthscale_2``, deliberately left unconstrained - (no positivity bijector; see the module docstring) so nothing bounds how far an - HMC leapfrog step can push it -- can round-trip through a bijector's inverse to a - non-finite unconstrained value, which gpflow's own ``Parameter.assign`` rejects - with a low-level ``InvalidArgumentError`` (``Tensor had NaN/Inf values - [Op:CheckNumerics]``) that doesn't say *which* value or parameter caused it. Not a - hypothetical: this is the exact error hit mid-investigation in Phase 9b/9c (from a - genuinely out-of-range value, in that case an unrelated script bug, not a posterior - sample) -- see ``paper/PHASE9B_INVESTIGATION.md``. Re-raised here with the - parameter name and value attached; behavior-preserving for every value that was - already assignable (which is every value seen in this codebase's tests, reference - regressions, and the from-scratch stochastic reproduction runs). + This is called deep inside ``mixture.py``/``acquisition.py``'s hot loops to replay + one posterior/mixture-component sample at a time. An extreme outlier sample -- most + plausibly from ``lengthscale_2``, deliberately left unconstrained (no positivity + bijector; see the module docstring) so nothing bounds how far an HMC leapfrog step + can push it -- can round-trip through a bijector's inverse to a non-finite + unconstrained value, which gpflow's own ``Parameter.assign`` rejects with a + low-level ``InvalidArgumentError`` (``Tensor had NaN/Inf values + [Op:CheckNumerics]``) that doesn't say *which* value or parameter caused it. + Re-raised here with the parameter name and value attached. """ for param, value in zip(kernel.hyperparameters, values): try: @@ -229,8 +223,8 @@ def assign_hyperparameters(kernel: gpflow.kernels.Kernel, values: Sequence[float def save_hyperparameters(kernel: gpflow.kernels.Kernel) -> List[float]: """Snapshot ``kernel.hyperparameters``' current (constrained) values. - Phase 9c: pairs with :func:`assign_hyperparameters` to save/restore a kernel's - state around code that reassigns it in a loop (``mixture.sample_gp_posterior_mixture``, + Pairs with :func:`assign_hyperparameters` to save/restore a kernel's state around + code that reassigns it in a loop (``mixture.sample_gp_posterior_mixture``, ``acquisition.entropy_objective``) -- see their docstrings. Returns a plain list of floats, not live references, so later mutating the kernel cannot change the snapshot. """ diff --git a/src/bits_for_gaps/means.py b/src/bits_for_gaps/means.py index b95110c..8c795f4 100644 --- a/src/bits_for_gaps/means.py +++ b/src/bits_for_gaps/means.py @@ -1,12 +1,11 @@ """GP mean functions. -Moved from the paper code's ``fxns/my_mean_fxn.py``. The paper's VLE study uses a -zero mean over the log-activity-coefficient output (encoding ideal mixing, gamma -> 1, -in the absence of data). ``FixedInverseMean`` is an alternative physics-informed mean -retained for reference. +The paper's VLE study uses a zero mean over the log-activity-coefficient output +(encoding ideal mixing, gamma -> 1, in the absence of data). ``FixedInverseMean`` is an +alternative physics-informed mean retained for reference. -TODO(Phase 5): ``FixedInverseMean`` assumes the mole fraction is input column 0; -generalize the input-column convention when the kernel goes N-D. +TODO: ``FixedInverseMean`` assumes the mole fraction is input column 0; generalize the +input-column convention for kernels with more input dimensions. """ from __future__ import annotations diff --git a/src/bits_for_gaps/mixture.py b/src/bits_for_gaps/mixture.py index b71ba6c..93b8a89 100644 --- a/src/bits_for_gaps/mixture.py +++ b/src/bits_for_gaps/mixture.py @@ -1,36 +1,28 @@ """Gaussian-mixture predictive posterior from hyperparameter-posterior draws. -Moved from the paper code's ``driver_new.py`` (``adaptiveEntropy.sample_gp_posterior_mixture`` -/ ``gp_predict_2D``). Each mixture component corresponds to one HMC posterior draw of the -GP kernel's hyperparameters (in that kernel's canonical order -- see -``kernels.AnisotropicSE.hyperparameters``); sampling reassigns the given GP model's kernel -parameters, once per draw, to walk through the mixture components, matching the paper code. - -Phase 5: kernel hyperparameters are no longer assigned by hardcoded attribute name -- -``kernels.assign_hyperparameters`` maps trace columns to a kernel's ``.hyperparameters`` -generically, so this works for any dimension/kernel exposing that contract. - -Phase 9c: :func:`sample_gp_posterior_mixture` used to leave ``GPmodel.kernel`` at -whatever the *last* draw's hyperparameters happened to be -- harmless for the values -this module itself returns (computed before the mutation matters), but a footgun for -any caller that reuses the same ``GPmodel`` object afterward: Phase 9b traced a -from-scratch validation script's spurious McCabe-Thiele non-convergence to exactly -this (see ``paper/PHASE9B_INVESTIGATION.md``) -- a test-RMSE step left the kernel at an -arbitrary leftover state before a later step reused the same model for the phase -diagram. Now saves the kernel's hyperparameters before sampling and restores them in a -``finally``, so the caller's model is unchanged after the call regardless of how it -exits. This changes only the model's POST-call state, not any value this module -computes and returns -- behavior-preserving for every existing test/baseline. +Each mixture component corresponds to one HMC posterior draw of the GP kernel's +hyperparameters (Eq 7), in that kernel's canonical order (see +``kernels.AnisotropicSE.hyperparameters``). Sampling reassigns the given GP model's +kernel parameters once per draw to walk through the mixture components, via +``kernels.assign_hyperparameters`` -- which maps trace columns to a kernel's +``.hyperparameters`` generically, so this works for any dimension/kernel exposing that +contract. + +CAUTION -- kernel mutation: :func:`sample_gp_posterior_mixture` reassigns +``GPmodel.kernel``'s hyperparameters once per posterior draw while it runs. It saves the +kernel's hyperparameters beforehand and restores them in a ``finally``, so the caller's +model is unchanged after the call returns regardless of how it exits. A caller that +reuses the same ``GPmodel`` object for a second purpose (e.g. building a phase diagram) +must do so *after* this function returns, not concurrently with it -- reusing it from +code that runs interleaved with this function would see the kernel at an arbitrary +intermediate draw's hyperparameters, not the model's real state. NOTE: ``GPmodel.predict_f_samples`` draws from GPflow/TensorFlow's *ambient* default random generator, not ``numpy``'s -- the ``np.random.seed`` call below seeds which -posterior-sample components are selected, but not the draws themselves. This means -these two functions were never bitwise-reproducible in the original paper code either -(confirmed: two successive calls to ``predict_f_samples`` on the same inputs, same -process, differ). That is why the Phase 2 integration test deliberately excludes this -plotting-only step from its determinism/baseline pins. Phase 9c adds an OPTIONAL -``tf_seed`` to make a single call reproducible on request (seeds TF's global RNG right -before drawing) -- default ``None`` leaves the ambient-RNG behavior unchanged. +posterior-sample components are selected, but not the draws themselves. Two successive +calls to ``predict_f_samples`` on the same inputs, in the same process, will differ. Pass +``tf_seed`` to make a single call's draws reproducible (seeds TF's global RNG right +before drawing) -- the default ``None`` leaves the ambient-RNG behavior unchanged. """ from __future__ import annotations @@ -114,13 +106,13 @@ def predict_grid_2D( ) -> np.ndarray: """Full-grid GP posterior-predictive samples, for 2-D plotting diagnostics only. - Moved from ``adaptiveEntropy.gp_predict_2D``. This does **not** feed the acquisition - (entropy/next-point selection) -- it is an expensive (``size`` full-covariance draws - over an ``n_grid x n_grid`` grid), disk-write-oriented diagnostic kept only for - parity with the paper's plotting pipeline. Callers should treat it as opt-in. + This does **not** feed the acquisition (entropy/next-point selection) -- it is an + expensive (``size`` full-covariance draws over an ``n_grid x n_grid`` grid), + disk-write-oriented diagnostic kept only for parity with the paper's plotting + pipeline. Callers should treat it as opt-in. - 2-D-ONLY (Phase 5): a dense grid is exponential in the input dimension d, so this - is intentionally not generalized to N-D -- see ``acquisition.optimize`` for the + 2-D-ONLY: a dense grid is exponential in the input dimension d, so this is + intentionally not generalized to N-D -- see ``acquisition.optimize`` for the N-D-general acquisition path this diagnostic does not feed. Parameters diff --git a/src/bits_for_gaps/sampler.py b/src/bits_for_gaps/sampler.py index d703ae9..3e973b5 100644 --- a/src/bits_for_gaps/sampler.py +++ b/src/bits_for_gaps/sampler.py @@ -1,34 +1,31 @@ """The BITS for GAPS sequential-design engine. -Phase 4: ``adaptiveEntropy`` is now an *orchestrator* over the decomposed modules -- -:mod:`bits_for_gaps.gp` (GP construction + HMC + R-hat/ESS), :mod:`bits_for_gaps.mixture` -(GMM predictive posterior), :mod:`bits_for_gaps.acquisition` (entropy objective + N-D -optimizer), :mod:`bits_for_gaps.transforms` (per-dimension input/output transforms), and -:mod:`bits_for_gaps.state` (in-memory run history). It no longer contains the algorithm +``adaptiveEntropy`` is an *orchestrator* over decomposed modules -- :mod:`bits_for_gaps.gp` +(GP construction + HMC + R-hat/ESS), :mod:`bits_for_gaps.mixture` (GMM predictive +posterior), :mod:`bits_for_gaps.acquisition` (entropy objective + N-D optimizer), +:mod:`bits_for_gaps.transforms` (per-dimension input/output transforms), and +:mod:`bits_for_gaps.state` (in-memory run history). It does not contain the algorithm math itself -- that lives in the modules above, as pure functions over explicit arguments, independently testable and reusable. -Disk-as-state is retired: :meth:`adaptiveEntropy.run` takes the initial design directly -(in memory) and returns a :class:`bits_for_gaps.state.RunHistory` -- a full run executes -with zero disk writes by default. File output (mirroring the paper code's per-iteration -``np.savetxt``/``pickle`` dump under ``results/{exp_name}/``) is available but opt-in via -``checkpoint_dir``. :meth:`run_model` is kept as a deprecated, disk-based shim for -scripts still relying on the original zero-argument convention. - -The decomposition is behavior-preserving: ``tests/integration/data/synthetic_baseline.json`` -pins the tiny seeded synthetic run's exact outputs from before this decomposition, and -the acquisition/gp/mixture modules were verified bit-exact (atol=1e-12) against the -pre-decomposition methods before this rewrite. - -Phase 5: generalized to N input dimensions. ``optimize`` (was ``optimize_2D``) -- the -acquisition path a run actually depends on -- is dimension-general. ``predict_grid_2D`` -and ``entropy_surface_2D`` stay 2-D-only (dense grids are exponential in d; they are -visualization diagnostics that don't feed the acquisition): :meth:`run` only calls -``entropy_surface_2D`` when ``len(self.XBnds) == 2``, leaving ``entropy_field=None`` -otherwise, and both raise a clear ``ValueError`` if called directly for d != 2. -``call_model`` calls the injected black box as ``FwdModel(*FwdModelArgs, *xStar)`` -- -``xStar``'s components in natural dimension order (was a 2-D-specific, reversed -``FwdModel(*args, x2, x1)`` convention inherited from the VLE example's Julia call). +State is in memory, not on disk: :meth:`adaptiveEntropy.run` takes the initial design +directly and returns a :class:`bits_for_gaps.state.RunHistory` -- a full run executes +with zero disk writes by default. Per-iteration file output (``np.savetxt``/``pickle`` +under ``results/{exp_name}/``) is available but opt-in via ``checkpoint_dir``. +:meth:`run_model` is kept as a deprecated, disk-based shim for scripts relying on the +original zero-argument convention. + +``tests/integration/data/synthetic_baseline.json`` pins a tiny seeded synthetic run's +exact outputs as a regression baseline (atol=1e-10). + +The pipeline is dimension-general: ``optimize`` -- the acquisition path a run actually +depends on -- works at any input dimension. ``predict_grid_2D`` and ``entropy_surface_2D`` +stay 2-D-only (dense grids are exponential in d; they are visualization diagnostics that +don't feed the acquisition): :meth:`run` only calls ``entropy_surface_2D`` when +``len(self.XBnds) == 2``, leaving ``entropy_field=None`` otherwise, and both raise a clear +``ValueError`` if called directly for d != 2. ``call_model`` calls the injected black box +as ``FwdModel(*FwdModelArgs, *xStar)`` -- ``xStar``'s components in natural dimension +order. """ from __future__ import annotations @@ -51,7 +48,7 @@ def _validate_bounds(x_bounds: Bounds) -> None: - """Each entry must be a ``(lo, hi)`` pair with ``lo < hi`` -- Phase 9c.""" + """Each entry must be a ``(lo, hi)`` pair with ``lo < hi``.""" for i, b in enumerate(x_bounds): if len(b) != 2: raise ValueError(f"x_bounds[{i}] must be a (lo, hi) pair, got {b!r}") @@ -78,10 +75,10 @@ def __init__( exp_name, iters, x_bounds, likelihood_var, mean_fxn, kernel_fxn, fwd_model, fwd_model_args - Phase 9c: validates ``x_bounds`` (each ``lo < hi``) and, if ``kernel_fxn`` - exposes an ``.ndim`` (e.g. ``kernels.AnisotropicSE``), that it matches - ``len(x_bounds)`` -- a mismatch here previously surfaced as a cryptic shape - error deep inside GPflow/TensorFlow the first time the kernel was evaluated. + Validates ``x_bounds`` (each ``lo < hi``) and, if ``kernel_fxn`` exposes an + ``.ndim`` (e.g. ``kernels.AnisotropicSE``), that it matches ``len(x_bounds)`` -- + a mismatch otherwise surfaces as a cryptic shape error deep inside + GPflow/TensorFlow the first time the kernel is evaluated. """ _validate_bounds(x_bounds) kernel_ndim = getattr(kernel_fxn, "ndim", None) @@ -103,13 +100,12 @@ def __init__( ## Sequential design self.noIters = iters self.startIter = 0 # iteration offset for resuming a prior run - # (was a hardcoded ``i += 50`` in the manuscript run) self.noRestarts = 10 self.noGaussians = 25 self.entropyMesh = [10 for _ in self.XBnds] - # Phase 9d: "taylor" (default, the paper's 2nd-order Taylor estimator, matches - # all pre-Phase-9d behavior/baselines) or "lower_bound" (the paper's closed- - # form cross-overlap lower bound, Theorem/SI-2) -- see acquisition.py. + # "taylor" (default: the paper's 2nd-order Taylor estimator) or "lower_bound" + # (the paper's closed-form cross-overlap lower bound, Theorem/SI-2) -- see + # acquisition.py. self.acquisitionObjective = "taylor" self.optMethod = None self.optOptions = None @@ -268,9 +264,8 @@ def call_model( """Evaluate the injected black box at ``xStar`` and append it to the design. Calls ``self.FwdModel(*self.FwdModelArgs, *xStar)`` -- ``xStar``'s components - in natural dimension order (Phase 5: was a 2-D-specific, reversed - ``FwdModel(*args, x2, x1)`` convention inherited from the VLE example's Julia - call). Returns the extended ``(XData, yData)`` -- no disk write. + in natural dimension order. Returns the extended ``(XData, yData)`` -- no disk + write. """ xStar = np.asarray(xStar, dtype=float).reshape(-1) XData = np.atleast_2d(np.asarray(XData, dtype=float)) @@ -291,7 +286,7 @@ def call_model( ## ------------------------------------------------------------------ def _validate_config(self) -> None: - """Positive/range checks on the HMC + acquisition config -- Phase 9c. + """Positive/range checks on the HMC + acquisition config. Checked here (not in ``__init__``) because every caller sets these via attribute assignment *after* construction (e.g. ``bfg.noSamples = 5000``) -- @@ -444,11 +439,11 @@ def _write_checkpoint( ) -> None: """Persist one iteration's artifacts to disk (opt-in; off by default). - A Phase-4, best-effort equivalent of the paper code's per-iteration file dump -- - not guaranteed byte-identical to the original file layout (e.g. the original - also wrote an intermediate ``gp_training_data_`` file); intended for users who - want on-disk artifacts, not as the mechanism for state hand-off between - iterations (see module docstring). + A best-effort equivalent of the paper code's per-iteration file dump -- not + guaranteed byte-identical to the original file layout (e.g. the original also + wrote an intermediate ``gp_training_data_`` file); intended for users who want + on-disk artifacts, not as the mechanism for state hand-off between iterations + (see module docstring). """ os.makedirs(checkpoint_dir, exist_ok=True) it = record.iteration @@ -479,16 +474,16 @@ def _write_checkpoint( class BitsForGaps(adaptiveEntropy): """Public-API-friendly constructor for the BITS-for-GAPS sequential-design engine. - A thin wrapper over :class:`adaptiveEntropy` (kept for backward compatibility) - using the target public kwarg names from REFACTOR_PLAN.md Sec 4. Numerically - identical to ``adaptiveEntropy`` -- no new computation, just friendlier - constructor names; all methods (including :meth:`run`) are inherited unchanged. - Advanced/legacy configuration (HMC tuning, restarts, mesh density, ...) is still set - via the same instance attributes as ``adaptiveEntropy`` (e.g. ``.noSamples``). + A thin wrapper over :class:`adaptiveEntropy` (kept for backward compatibility) with + friendlier public kwarg names. Numerically identical to ``adaptiveEntropy`` -- no new + computation, just friendlier constructor names; all methods (including :meth:`run`) + are inherited unchanged. Advanced/legacy configuration (HMC tuning, restarts, mesh + density, ...) is still set via the same instance attributes as ``adaptiveEntropy`` + (e.g. ``.noSamples``). - TODO(Phase 6): once the VLE example is ported onto this API, this can grow an - ``mcmc=MCMCConfig(...)``-style kwarg for HMC tuning, replacing the passthrough - instance attributes inherited from ``adaptiveEntropy`` (e.g. ``.noSamples``). + TODO: an ``mcmc=MCMCConfig(...)``-style kwarg for HMC tuning would be a cleaner + alternative to the passthrough instance attributes inherited from + ``adaptiveEntropy`` (e.g. ``.noSamples``). """ def __init__( diff --git a/src/bits_for_gaps/state.py b/src/bits_for_gaps/state.py index a7ff459..b764f8c 100644 --- a/src/bits_for_gaps/state.py +++ b/src/bits_for_gaps/state.py @@ -1,11 +1,10 @@ """In-memory run state for the sequential-design loop. -Replaces the paper code's disk-as-state convention (``np.savetxt``/``pickle`` under -``results/{exp_name}/``, read back on the next iteration -- see ``driver_new.py``'s -``adaptiveEntropy.run_model``) with plain in-memory records. Disk output is still -available, but as an *opt-in* checkpoint (``sampler.adaptiveEntropy.run``'s -``checkpoint_dir`` argument), not the mechanism by which state is threaded across -iterations. +State is threaded across iterations via plain in-memory records, not disk. Disk output +is still available, but as an *opt-in* checkpoint (``sampler.adaptiveEntropy.run``'s +``checkpoint_dir`` argument, mirroring the paper code's per-iteration ``np.savetxt``/ +``pickle`` dump under ``results/{exp_name}/``), not the mechanism by which state is +threaded across iterations. """ from __future__ import annotations @@ -33,7 +32,7 @@ class IterationRecord: chains_states: np.ndarray # all chains, unconstrained rhat: np.ndarray ess: np.ndarray - entropy_field: Optional[np.ndarray] = None # 2-D only (Phase 5 generalizes) + entropy_field: Optional[np.ndarray] = None # 2-D only; None for N-D runs xStar: Optional[np.ndarray] = None max_entropy: Optional[float] = None lml_result: Optional[Any] = None # scipy OptimizeResult, if maximize_lml ran diff --git a/src/bits_for_gaps/transforms.py b/src/bits_for_gaps/transforms.py index ae0d5dc..f6ac618 100644 --- a/src/bits_for_gaps/transforms.py +++ b/src/bits_for_gaps/transforms.py @@ -1,15 +1,14 @@ """Per-dimension input/output transforms for the GP surrogate. -Lifts the paper code's list-of-lambdas convention (``driver_new.py``'s -``adaptiveEntropy.XTrsfFwd``/``XTrsfBkwd``/``yTrsfFwd``/``yTrsfBkwd``, one identity -lambda per input dimension plus a scalar output lambda) into small, testable classes. -Identity by default -- matches the paper's VLE study, which trains the GP directly on -mole fraction / temperature / log-activity-coefficient with no rescaling. +One elementwise callable per input dimension plus a scalar output callable, lifted into +small, testable classes. Identity by default -- matches the paper's VLE study, which +trains the GP directly on mole fraction / temperature / log-activity-coefficient with no +rescaling. ``forward_fns``/``backward_fns`` (lists of per-dimension callables) are exposed alongside the array-oriented ``forward``/``backward`` methods because the acquisition/ mixture code applies them element-wise to scalars (e.g. a single bound) as well as to -whole columns -- exactly the calling convention the paper code used. +whole columns. """ from __future__ import annotations diff --git a/tests/conftest.py b/tests/conftest.py index e123562..c54b0e3 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -5,7 +5,7 @@ from the archived published run (iteration 15); see ``paper/reference/README.md``. Also puts the repo's ``examples/`` directory AND the repo root on ``sys.path`` so -``examples/`` and ``paper/`` (neither pip-installed -- see REFACTOR_PLAN.md §7.3) are +``examples/`` and ``paper/`` (neither pip-installed -- both are repo-only) are importable in dev/CI as top-level packages, e.g. ``import vle_distillation. activity_model`` and ``import paper.figures.fig10_traces``. This also makes ``examples/vle_distillation/juliapkg.json`` discoverable by juliapkg (which scans diff --git a/tests/integration/test_bits_for_gaps_facade.py b/tests/integration/test_bits_for_gaps_facade.py index d6206e7..071666a 100644 --- a/tests/integration/test_bits_for_gaps_facade.py +++ b/tests/integration/test_bits_for_gaps_facade.py @@ -1,4 +1,4 @@ -"""Parity test for the ``BitsForGaps`` public-API facade (REFACTOR_PLAN.md Sec 4). +"""Parity test for the ``BitsForGaps`` public-API facade. ``BitsForGaps`` is a thin, renamed-kwarg subclass of ``adaptiveEntropy`` -- it adds no new computation, so a run through it must reproduce the exact same pinned baseline as diff --git a/tests/integration/test_end_to_end.py b/tests/integration/test_end_to_end.py index 50d045e..a329c06 100644 --- a/tests/integration/test_end_to_end.py +++ b/tests/integration/test_end_to_end.py @@ -1,12 +1,11 @@ """Seeded end-to-end integration test for the BITS-for-GAPS sampler. Runs the full sequential-design decision pipeline of ``adaptiveEntropy`` on a *synthetic, -pure-Python* black box (a smooth 2-D function -- NO Julia), via the Phase 4 in-memory -``run()`` API. The run must complete and its outputs (the selected next point, the -R-hat/ESS shapes, and the entropy field) must be stable across two runs with the same -seed, and reproduce ``tests/integration/data/synthetic_baseline.json`` -- a hard pin of -this run's exact outputs captured from the pre-Phase-4 (monolithic, disk-based) -``sampler.py``. +pure-Python* black box (a smooth 2-D function -- NO Julia), via the in-memory ``run()`` +API. The run must complete and its outputs (the selected next point, the R-hat/ESS +shapes, and the entropy field) must be stable across two runs with the same seed, and +reproduce ``tests/integration/data/synthetic_baseline.json`` -- a hard pin of this run's +exact outputs. We deliberately don't pass ``predict_grid=True``: it re-pickles the model and computes the full-grid posterior-sample array used only for figures, takes ~20 s (100 full- @@ -38,8 +37,7 @@ def _true_f(x1, x2): def _fwd_model(x1, x2): - # Phase 5: the sampler calls FwdModel(*args, *xStar) -- natural dimension order - # (was a reversed, 2-D-specific FwdModel(*args, x2, x1) convention pre-Phase-5). + # The sampler calls FwdModel(*args, *xStar) in natural dimension order. return [float(_true_f(x1, x2))] @@ -137,7 +135,7 @@ def test_next_point_appended_via_injected_fwd_model(run_a): @pytest.mark.slow def test_stable_across_two_runs_with_same_seed(run_a, run_b): # Same seed, same process => the sampler must be deterministic. This is the guard - # against nondeterminism sneaking in during the Phase 4 decomposition. + # against nondeterminism sneaking into the sampler's internals. a, b = run_a, run_b np.testing.assert_allclose(a["rhat"], b["rhat"], atol=1e-10) np.testing.assert_allclose(a["ess"], b["ess"], atol=1e-10) @@ -147,10 +145,10 @@ def test_stable_across_two_runs_with_same_seed(run_a, run_b): @pytest.mark.slow -def test_matches_pre_phase4_baseline(run_a): - # Hard pin against tests/integration/data/synthetic_baseline.json, captured from the - # monolithic (pre-decomposition) sampler.py. The Phase 4 module split -- and the - # disk-as-state removal -- must reproduce these exact numbers, not merely match itself. +def test_matches_synthetic_baseline(run_a): + # Hard pin against tests/integration/data/synthetic_baseline.json -- captured once + # and never regenerated, so this must reproduce these exact numbers, not merely + # match itself. with open(BASELINE_PATH) as f: base = json.load(f) r = run_a @@ -166,8 +164,8 @@ def test_matches_pre_phase4_baseline(run_a): @pytest.mark.slow def test_run_writes_no_files_by_default(tmp_path, monkeypatch): - # Phase 4 retires disk-as-state: a full run must execute with zero disk writes - # unless the caller explicitly opts in via checkpoint_dir. + # A full run must execute with zero disk writes unless the caller explicitly opts + # in via checkpoint_dir. monkeypatch.chdir(tmp_path) X_init, y_init = _initial_design(n=12) s = _build_sampler() @@ -217,12 +215,11 @@ def test_run_with_initial_lml_maximization(): @pytest.mark.slow def test_run_with_lower_bound_acquisition_objective_completes(): - # Phase 9d: acquisitionObjective="lower_bound" selects the paper's closed-form - # entropy lower bound (Theorem/SI-2) instead of the default 2nd-order Taylor - # estimator -- implemented and unit-tested since Phase 2 but never wired up as a - # usable acquisition objective before now. Default ("taylor") behavior/baselines - # are covered by every other test in this file; this just confirms the - # alternative objective runs a full iteration end-to-end without error. + # acquisitionObjective="lower_bound" selects the paper's closed-form entropy lower + # bound (Theorem/SI-2) instead of the default 2nd-order Taylor estimator. Default + # ("taylor") behavior/baselines are covered by every other test in this file; this + # just confirms the alternative objective runs a full iteration end-to-end without + # error. X_init, y_init = _initial_design(n=12) s = _build_sampler() s.acquisitionObjective = "lower_bound" @@ -244,14 +241,11 @@ def test_run_rejects_unknown_acquisition_objective(): @pytest.mark.slow def test_iteration_record_gpmodel_survives_downstream_mixture_sampling(): - # Phase 9c regression guard for the exact real-world bug Phase 9b found (see - # paper/PHASE9B_INVESTIGATION.md): a caller that reuses history.last.GPmodel for a - # second purpose (e.g. building a surrogate phase diagram) AFTER something else + # Regression guard: a caller that reuses history.last.GPmodel for a second purpose + # (e.g. building a surrogate phase diagram) AFTER something else # (mixture.sample_gp_posterior_mixture, e.g. for a test-RMSE metric) has already # been called on that same object must not find it left at an arbitrary leftover - # hyperparameter state. run()'s own internal acquisition search (entropy_objective, - # via optimize()/entropy_surface_2D()) used to leave exactly this kind of leftover - # state on record.GPmodel even before any downstream caller touched it. + # hyperparameter state. from bits_for_gaps import mixture X_init, y_init = _initial_design(n=12) diff --git a/tests/integration/test_nd_synthetic.py b/tests/integration/test_nd_synthetic.py index 299ff7e..0a30d9e 100644 --- a/tests/integration/test_nd_synthetic.py +++ b/tests/integration/test_nd_synthetic.py @@ -1,9 +1,9 @@ -"""Phase 5: N-D synthetic end-to-end tests via the ``BitsForGaps`` public API. +"""N-D synthetic end-to-end tests via the ``BitsForGaps`` public API. Runs a tiny, fully-seeded ``BitsForGaps.run(...)`` on pure-Python synthetic black boxes -(no Julia) at d=1 and d=3, proving the "already general" claim in code rather than just -at d=2. Mirrors ``test_end_to_end.py``'s tiny configuration and assertion style, but -parametrized over dimension. +(no Julia) at d=1 and d=3, proving the pipeline is dimension-general in code rather than +just at the paper's d=2. Mirrors ``test_end_to_end.py``'s tiny configuration and +assertion style, but parametrized over dimension. The 3-D kernel deliberately mixes prior families across dimensions (LogNormal-positive, Gamma-unconstrained, LogNormal-positive) -- the same per-dimension-Parameter design diff --git a/tests/regression/test_hmc_diagnostics.py b/tests/regression/test_hmc_diagnostics.py index ec160c5..849d91f 100644 --- a/tests/regression/test_hmc_diagnostics.py +++ b/tests/regression/test_hmc_diagnostics.py @@ -3,7 +3,7 @@ Reads the committed reference snapshot (``paper/reference/hmc_diagnostics.json``) extracted from the archived iteration-15 run and pins it against the values reported in the paper. This locks the reference file against corruption and documents the reproduction tolerances -that ``paper/reproduce.py`` will diff against in Phase 7. +that ``paper/reproduce.py`` diffs against. Pure (JSON only) -- no TensorFlow, no Julia. """ diff --git a/tests/regression/test_hyperparameter_posterior.py b/tests/regression/test_hyperparameter_posterior.py index 9343480..8d32cce 100644 --- a/tests/regression/test_hyperparameter_posterior.py +++ b/tests/regression/test_hyperparameter_posterior.py @@ -3,7 +3,7 @@ Reads ``paper/reference/hyperparameter_posterior.json`` -- the posterior mean/median/std and 5%/95% quantiles of the three kernel hyperparameters at the published iteration -- and checks the summary is well-formed and physically sensible. This pins the target the -mixture-posterior code must reproduce after the Phase 4/5 refactor. +mixture-posterior code must reproduce. Pure (JSON only) -- no TensorFlow, no Julia. """ diff --git a/tests/regression/test_mccabe_thiele.py b/tests/regression/test_mccabe_thiele.py index 3ab734b..57e8506 100644 --- a/tests/regression/test_mccabe_thiele.py +++ b/tests/regression/test_mccabe_thiele.py @@ -3,14 +3,14 @@ The stage table is the paper's end-to-end validation: the distillation column designed with the BITS-for-GAPS GP *surrogate* reproduces the column designed with the *Wilson* ground-truth activity model. Reproducing it requires the Julia/Clapeyron VLE distillation -backend, ported in Phase 6 (``examples/vle_distillation/``); the recompute test below is -gated behind ``@pytest.mark.vle`` (needs Julia -- deselected by default). +backend (``examples/vle_distillation/``); the recompute test below is gated behind +``@pytest.mark.vle`` (needs Julia -- deselected by default). The consistency checks below read only the committed reference JSON and pin the paper's claim (surrogate == Wilson within tolerance), so they run in the default suite. The recompute itself (``wilson_column``/``surrogate_column``) lives in -``paper.figures.fig09_mccabe_thiele`` (Phase 7) -- this test imports it rather than +``paper.figures.fig09_mccabe_thiele`` -- this test imports it rather than reimplementing it, since Fig 9's plotting script needs the exact same columns. """ @@ -75,9 +75,9 @@ def test_surrogate_matches_wilson(reference): @pytest.mark.vle def test_reproduce_stage_table_with_distillation_backend(reference): - # Phase 6/7: recompute the stage table through the ported VLE distillation backend - # and diff against the reference (transcribed from paper Fig 9c, hence the looser- - # than-report atol below -- see paper/reference/README.md and HANDOFF.md). + # Recompute the stage table through the VLE distillation backend and diff against + # the reference (transcribed from paper Fig 9c, hence the looser-than-report atol + # below -- see paper/reference/README.md). from paper.figures.fig09_mccabe_thiele import surrogate_column, wilson_column g = reference("mccabe_thiele_stages.json") diff --git a/tests/regression/test_paper_figures.py b/tests/regression/test_paper_figures.py index 6dd2fcc..9a38555 100644 --- a/tests/regression/test_paper_figures.py +++ b/tests/regression/test_paper_figures.py @@ -1,16 +1,16 @@ """Regression: recompute the paper's quantitative figures and diff against ``paper/reference/*``. -As of Phase 9, the data these tests read comes from the curated, COMMITTED -``paper/data/`` subset (see ``paper/data/README.md``) -- no private-archive access -needed, so the three read-only tests below (Fig 10, Fig 5, the hyperparameter -posterior behind Fig 11) run in the DEFAULT suite. Only +The data these tests read comes from the curated, COMMITTED ``paper/data/`` subset +(see ``paper/data/README.md``) -- no private-archive access needed, so the three +read-only tests below (Fig 10, Fig 5, the hyperparameter posterior behind Fig 11) run +in the DEFAULT suite. Only ``test_fig08_wilson_curve_matches_archived_ground_truth`` stays gated behind ``@pytest.mark.vle``: it *recomputes* the Wilson curve via live Clapeyron calls (needs Julia), unlike the others, which only read committed text files. -Point ``BFG_ARCHIVE_DIR`` at the full private archive instead of ``paper/data/`` if -you want to check iterations beyond the curated subset. +Point ``BFG_ARCHIVE_DIR`` at another directory of run artifacts (same layout) instead +of ``paper/data/`` to check a different run. Fig 9's stage table already has its own gated test (``test_mccabe_thiele.py``) -- not duplicated here. Fig 8 (phase diagram) has no dedicated reference file (it's a diff --git a/tests/unit/test_acquisition.py b/tests/unit/test_acquisition.py index 356e027..a0f57d2 100644 --- a/tests/unit/test_acquisition.py +++ b/tests/unit/test_acquisition.py @@ -1,12 +1,10 @@ """Unit tests for ``acquisition.py``'s entropy-maximization objective/optimizer. -Phase 9c: the headline behavior under test is that ``entropy_objective`` (and its -callers ``optimize``/``entropy_surface_2D``) leave the caller's ``GPmodel`` kernel -unchanged -- the same class of state-mutation footgun Phase 9b found in -``mixture.sample_gp_posterior_mixture`` (see ``paper/PHASE9B_INVESTIGATION.md``): since -``sampler.py``'s ``run()`` calls these on the same ``GPmodel`` object it later stores in -``IterationRecord.GPmodel``, an unrestored kernel would silently corrupt every -iteration's returned model. +The headline behavior under test is that ``entropy_objective`` (and its callers +``optimize``/``entropy_surface_2D``) leave the caller's ``GPmodel`` kernel unchanged -- +since ``sampler.py``'s ``run()`` calls these on the same ``GPmodel`` object it later +stores in ``IterationRecord.GPmodel``, an unrestored kernel would silently corrupt +every iteration's returned model. """ import gpflow @@ -101,7 +99,7 @@ def test_entropy_surface_2d_restores_kernel_state(gp_model, trace): ## --------------------------------------------------------------------------- -## Phase 9d: selectable acquisition objective (entropy_lower_bound wired up). +## Selectable acquisition objective (entropy_lower_bound as an alternative to Taylor). ## --------------------------------------------------------------------------- @@ -156,8 +154,8 @@ def test_entropy_estimators_registry_has_both_paper_estimators(): def test_entropy_surface_2d_rejects_non_2d_bounds(gp_model, trace): - # entropy_surface_2D is a 2-D-only visualization diagnostic (Phase 5) -- it must - # raise a clear error for other dimensions rather than silently misbehaving. + # entropy_surface_2D is a 2-D-only visualization diagnostic -- it must raise a + # clear error for other dimensions rather than silently misbehaving. with pytest.raises(ValueError, match="2-D-only"): entropy_surface_2D( trace, diff --git a/tests/unit/test_design.py b/tests/unit/test_design.py index 6981625..2982d17 100644 --- a/tests/unit/test_design.py +++ b/tests/unit/test_design.py @@ -49,10 +49,6 @@ def test_full_factorial_exact_grid_size_needs_no_trimming(): ## NOTE: `full_factorial_design`'s "grid too small" `ValueError` (design.py:74-75) is -## unreachable via any (bounds, n_train, n_test) combination: `levels = -## ceil(n_total ** (1/d))` guarantees `levels ** d >= n_total` for every positive -## integer `n_total` and `d` (verified by brute-force search over n_total <= 2000, -## d <= 4 -- no counterexample). It is defensive dead code, not a bug (every reachable -## call returns a correctly-shaped grid); left as-is since this pass is -## behavior-preserving and the guard is harmless. Not given a test here since there is -## no real input that exercises it. +## unreachable for any (bounds, n_train, n_test): `levels = ceil(n_total ** (1/d))` +## guarantees `levels ** d >= n_total`. Defensive dead code, not a bug -- not exercised +## here since no input triggers it. diff --git a/tests/unit/test_distillation.py b/tests/unit/test_distillation.py index 0c6ebf8..e83376d 100644 --- a/tests/unit/test_distillation.py +++ b/tests/unit/test_distillation.py @@ -94,7 +94,7 @@ def test_resolve_fixed_indices_rejects_unknown_name(): ## --------------------------------------------------------------------------- -## Phase 9c: solve_column's retry-on-non-convergence orchestration. +## solve_column's retry-on-non-convergence orchestration. ## ## Stubs _try_solve_column rather than hunting for a real equilibrium curve that ## reproducibly fails to converge -- fsolve's behavior on synthetic curves is diff --git a/tests/unit/test_entropy.py b/tests/unit/test_entropy.py index fc7aa5e..5cfca57 100644 --- a/tests/unit/test_entropy.py +++ b/tests/unit/test_entropy.py @@ -2,7 +2,7 @@ These are pure NumPy/SciPy (no TensorFlow), so they run fast and in CI without the GP stack. They combine closed-form correctness checks (no magic numbers) with one -regression pin against the paper code's 5-component mixture example (huber_et_al.py). +regression pin against the paper's own 5-component mixture example. """ import numpy as np @@ -59,8 +59,8 @@ def test_lower_bound_below_second_order_for_mixture(): def test_huber_5d_mixture_regression(): - # Regression pin against the paper code (fxns/max_ent_design.second_order_entropy) - # on the 5-component bivariate mixture from huber_et_al.py, with means[4] = [1, 1]. + # Regression pin on the paper's own 5-component bivariate mixture example, + # with means[4] = [1, 1]. means = np.array([[0, 0], [3, 2], [1, -0.5], [2.5, 1.5], [1, 1]], dtype=float) covs = np.array( [ @@ -77,7 +77,7 @@ def test_huber_5d_mixture_regression(): ## --------------------------------------------------------------------------- -## Phase 9c: assert -> explicit exception for the runtime, data-dependent density +## Explicit exception (not a bare assert) for the runtime, data-dependent density ## check (asserts are silently stripped under python -O). ## --------------------------------------------------------------------------- diff --git a/tests/unit/test_entropy_mc_validation.py b/tests/unit/test_entropy_mc_validation.py index 3587b3e..953468c 100644 --- a/tests/unit/test_entropy_mc_validation.py +++ b/tests/unit/test_entropy_mc_validation.py @@ -1,5 +1,5 @@ """Monte-Carlo validation of the entropy estimators against the true GMM differential -entropy -- Phase 9d. +entropy. ``test_entropy.py`` has closed-form checks (single-component exactness) and a regression pin against the paper code's output, but nothing that validates diff --git a/tests/unit/test_gp.py b/tests/unit/test_gp.py index ec97eb3..f998ece 100644 --- a/tests/unit/test_gp.py +++ b/tests/unit/test_gp.py @@ -1,8 +1,7 @@ """Unit tests for ``gp.py``'s GP construction and log-marginal-likelihood optimization. -Phase 9d: this module had no direct unit tests before (only indirect coverage via the -integration suite's full HMC runs, which never exercise the ``summarize``/``debug_cov`` -diagnostic branches). ``run_mcmc`` itself stays integration-tested (it's inherently a +Covers the ``summarize``/``debug_cov`` diagnostic branches the integration suite's full +HMC runs don't exercise. ``run_mcmc`` itself stays integration-tested (it's inherently a full HMC run, not a cheap unit) -- see ``tests/integration/test_end_to_end.py``. """ diff --git a/tests/unit/test_kernels.py b/tests/unit/test_kernels.py index 35898cc..9274fd7 100644 --- a/tests/unit/test_kernels.py +++ b/tests/unit/test_kernels.py @@ -1,10 +1,10 @@ """Unit tests for the AnisotropicSE covariance kernel. -Pins the 2-D kernel behavior (symmetry, positive semi-definiteness, the (n, m) -cross-covariance shape, the K_diag shortcut, per-dimension lengthscale scaling) that -predates the Phase 5 N-D generalization, plus new tests for the generalized API: the -canonical ``hyperparameters`` order, ``assign_hyperparameters`` round-tripping, and -explicit N-D (1-D and 3-D) construction with per-dimension prior families. +Covers both the 2-D kernel behavior (symmetry, positive semi-definiteness, the (n, m) +cross-covariance shape, the K_diag shortcut, per-dimension lengthscale scaling) and the +generalized N-D API: the canonical ``hyperparameters`` order, ``assign_hyperparameters`` +round-tripping, and explicit N-D (1-D and 3-D) construction with per-dimension prior +families. """ import gpflow @@ -93,7 +93,7 @@ def test_anisotropic_lengthscale_scaling(kernel): ## --------------------------------------------------------------------------- -## Phase 5: N-D generalization -- canonical ordering, generic assignment, N-D construction +## N-D generalization -- canonical ordering, generic assignment, N-D construction ## --------------------------------------------------------------------------- @@ -134,8 +134,9 @@ def test_assign_hyperparameters_round_trips(kernel): ## --------------------------------------------------------------------------- -## Phase 9c: save/restore -- the save_hyperparameters half of the mutation-footgun fix -## (mixture.py/acquisition.py use these together; see their tests for the full loop). +## Kernel hyperparameter save/restore (mixture.py/acquisition.py use these together +## to guard against leaving a GP's kernel mutated mid-loop; see their tests for the +## full loop). ## --------------------------------------------------------------------------- @@ -219,12 +220,9 @@ def test_variance_prior_requires_lengthscale_priors(): ## --------------------------------------------------------------------------- -## Phase 9d: assign_hyperparameters raises a clear, specific error instead of a -## low-level gpflow/TF traceback for a value that can't round-trip through a -## parameter's transform -- the exact error class hit mid-investigation in Phase 9b/9c -## (see paper/PHASE9B_INVESTIGATION.md and kernels.py's assign_hyperparameters -## docstring). Behavior-preserving for every assignable value (every value seen in -## this codebase's tests/reference regressions/from-scratch reproduction runs). +## assign_hyperparameters raises a clear, specific error instead of a low-level +## gpflow/TF traceback for a value that can't round-trip through a parameter's +## transform (see kernels.py's assign_hyperparameters docstring). ## --------------------------------------------------------------------------- diff --git a/tests/unit/test_mixture.py b/tests/unit/test_mixture.py index 6fea968..d368c23 100644 --- a/tests/unit/test_mixture.py +++ b/tests/unit/test_mixture.py @@ -1,9 +1,10 @@ """Unit tests for ``mixture.py``'s posterior-mixture sampling. -Phase 9c: the headline behavior under test is that ``sample_gp_posterior_mixture`` -(and ``predict_grid_2D``, which calls it) leave the caller's ``GPmodel`` unchanged -- -Phase 9b traced a real bug to this function leaving the kernel at an arbitrary leftover -hyperparameter state (see ``paper/PHASE9B_INVESTIGATION.md``). +The headline behavior under test is that ``sample_gp_posterior_mixture`` (and +``predict_grid_2D``, which calls it) leave the caller's ``GPmodel`` unchanged -- both +mutate the kernel's hyperparameters once per posterior draw while running, so a caller +reusing the same model object afterward must not find it left at an arbitrary leftover +state. """ import gpflow @@ -93,8 +94,8 @@ def test_sample_gp_posterior_mixture_default_tf_seed_is_unset(gp_model, trace): def test_predict_grid_2d_rejects_non_2d_bounds(gp_model, trace): - # predict_grid_2D is a 2-D-only visualization diagnostic (Phase 5) -- it must - # raise a clear error for other dimensions rather than silently misbehaving. + # predict_grid_2D is a 2-D-only visualization diagnostic -- it must raise a clear + # error for other dimensions rather than silently misbehaving. with pytest.raises(ValueError, match="2-D-only"): predict_grid_2D( trace, diff --git a/tests/unit/test_sampler_validation.py b/tests/unit/test_sampler_validation.py index 4eb52f2..341ce94 100644 --- a/tests/unit/test_sampler_validation.py +++ b/tests/unit/test_sampler_validation.py @@ -1,6 +1,6 @@ -"""Unit tests for Phase 9c's public-API input validation on ``adaptiveEntropy``/ -``BitsForGaps``: clear, early errors instead of a cryptic failure deep inside -GPflow/TensorFlow or a confusing downstream shape mismatch. +"""Unit tests for ``adaptiveEntropy``/``BitsForGaps``'s public-API input validation: +clear, early errors instead of a cryptic failure deep inside GPflow/TensorFlow or a +confusing downstream shape mismatch. """ import gpflow