diff --git a/README.md b/README.md index cec4852..24886d5 100644 --- a/README.md +++ b/README.md @@ -109,3 +109,23 @@ while updating observations/pending files, then resume with new request IDs. Preview mode can still display already-reserved conditions because it does not read the ledger. Store the ledger alongside its recovery journal files outside Git and use SQLite-aware backups while it is active. + +## Acquisition bake-off (v0.3) + +`medialoop-bakeoff` compares the current planner (expected improvement) with +random search, fixed-kappa UCB, and one-draw Thompson sampling on a **fixed +grid** inside the published windows. A fifth policy, repeat expected +improvement, re-proposes the same candidate inside a batch so violations are +visible. Guarded policies fantasize at the posterior mean within a batch +(kriging believer) and call the cartoon only once per chosen grid point. + +```bash +PYTHONPATH=src python3 -m medialoop.bakeoff --objective cardiac --seeds 6 +``` + +Report simple regret against the best point **on that grid**, cumulative +regret, and batch-uniqueness violations. This is still the synthetic surface. +It is not a dose, and it is not a reason to plate a well. A neural policy is +intentionally absent: a course project can add one behind the same regret +and violation columns without replacing the expected-improvement baseline. + diff --git a/pyproject.toml b/pyproject.toml index 6296efa..3a156dc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "medialoop" -version = "0.2.0" +version = "0.3.0" description = "Closed-loop in-vitro media search against a synthetic surface." readme = "README.md" requires-python = ">=3.10" @@ -17,6 +17,7 @@ dependencies = [ [project.scripts] medialoop = "medialoop.cli:main" medialoop-plan = "medialoop.planner:main" +medialoop-bakeoff = "medialoop.bakeoff:main" [tool.setuptools.packages.find] where = ["src"] diff --git a/src/medialoop/__init__.py b/src/medialoop/__init__.py index 9df9b6f..85a4a17 100644 --- a/src/medialoop/__init__.py +++ b/src/medialoop/__init__.py @@ -1,3 +1,3 @@ """Closed-loop proposals for in-vitro differentiation media screens.""" -__version__ = "0.2.0" +__version__ = "0.3.0" diff --git a/src/medialoop/bakeoff.py b/src/medialoop/bakeoff.py new file mode 100644 index 0000000..1617bb3 --- /dev/null +++ b/src/medialoop/bakeoff.py @@ -0,0 +1,188 @@ +"""Fixed-grid acquisition bake-off on the synthetic media surface. + +The current planner is expected improvement. This module compares it with +random search, fixed-kappa UCB, and one-draw Thompson sampling. A repeat +policy deliberately re-proposes the same point inside a batch so constraint +violations are countable. None of these scores are cell readouts. +""" + +from __future__ import annotations + +import argparse +import json +import sys + +import numpy as np + +from medialoop.gp import expected_improvement, predict, thompson_draw, upper_confidence_bound +from medialoop.loop import _row, inside_bounds +from medialoop.surfaces import OBJECTIVES + +POLICIES = ( + "random", + "expected_improvement", + "ucb", + "thompson", + "repeat_expected_improvement", +) + + +def candidate_grid() -> np.ndarray: + """Published-window lattice. The cardiac cartoon peak (8 µM, 4 µM, 0, 0) is on it.""" + chir = (0.0, 3.0, 6.0, 8.0, 12.0) + iwp = (0.0, 2.0, 4.0, 5.0) + sb = (0.0, 5.0, 10.0) + ldn = (0.0, 100.0, 250.0) + rows = [(a, b, c, d) for a in chir for b in iwp for c in sb for d in ldn] + return np.asarray(rows, dtype=float) + + +def _scores(policy, mu, sigma, fantasy_y, rng, kappa): + if policy in ("expected_improvement", "repeat_expected_improvement"): + return expected_improvement(mu, sigma, float(np.max(fantasy_y))) + if policy == "ucb": + return upper_confidence_bound(mu, sigma, kappa) + if policy == "thompson": + return thompson_draw(mu, sigma, rng) + if policy == "random": + return rng.random(len(mu)) + raise ValueError(f"unknown policy: {policy}") + + +def campaign( + objective, + policy: str, + *, + seed: int, + n_init: int = 4, + rounds: int = 4, + batch_size: int = 2, + noise: float = 0.02, + kappa: float = 1.5, +) -> dict: + if policy not in POLICIES: + raise ValueError(f"unknown policy: {policy}") + if n_init < 2 or rounds < 1 or batch_size < 1: + raise ValueError("need n_init >= 2, rounds >= 1, batch_size >= 1") + if not np.isfinite(noise) or noise < 0 or not np.isfinite(kappa) or kappa < 0: + raise ValueError("noise and kappa must be finite and nonnegative") + grid = candidate_grid() + if len(grid) <= n_init + rounds: + raise ValueError("grid is smaller than the evaluation budget") + truth = np.array([objective(_row(row)) for row in grid]) + oracle = float(truth.max()) + rng = np.random.default_rng(seed) + observed = list(map(int, rng.choice(len(grid), size=n_init, replace=False))) + y_obs = truth[observed] + rng.normal(0.0, noise, size=n_init) + violations = 0 + regrets = [] + + def mark() -> None: + best = float(truth[list(dict.fromkeys(observed))].max()) + regrets.append(oracle - best) + + mark() + for _ in range(rounds): + picks: list[int] = [] + fantasy_x = grid[observed] + fantasy_y = np.array(y_obs, dtype=float) + for _step in range(batch_size): + if policy == "repeat_expected_improvement": + pool = [i for i in range(len(grid)) if i not in observed] + else: + pool = [i for i in range(len(grid)) if i not in observed and i not in picks] + if not pool: + break + if policy == "random": + choice = int(rng.choice(pool)) + else: + mu, sigma = predict(fantasy_x, fantasy_y, grid[pool], noise=noise) + score = _scores(policy, mu, sigma, fantasy_y, rng, kappa) + tied = np.flatnonzero(np.isclose(score, score.max(), rtol=1e-10, atol=1e-15)) + # The repeat policy always re-selects the lowest index so a violation is forced. + local = int(tied[0] if policy == "repeat_expected_improvement" else rng.choice(tied)) + choice = pool[local] + if choice in observed or choice in picks: + violations += 1 + picks.append(choice) + if policy != "repeat_expected_improvement": + mu_c, _sigma_c = predict(fantasy_x, fantasy_y, grid[[choice]], noise=noise) + fantasy_y = np.append(fantasy_y, float(mu_c[0])) + fantasy_x = np.vstack([fantasy_x, grid[choice]]) + for choice in picks: + y_obs = np.append(y_obs, truth[choice] + float(rng.normal(0.0, noise))) + observed.append(choice) + mark() + chosen = grid[observed] + if not inside_bounds(chosen): + violations += int(len(observed)) + unique = list(dict.fromkeys(observed)) + return { + "policy": policy, + "best_true": float(truth[unique].max()), + "oracle_max_on_grid": oracle, + "final_simple_regret": regrets[-1], + "cumulative_regret": float(sum(regrets)), + "simple_regret": [float(value) for value in regrets], + "violations": int(violations), + "n_evaluations": len(observed), + "n_unique": len(unique), + "seed": seed, + "batch_size": batch_size, + "rounds": rounds, + "n_init": n_init, + "note": ( + "Synthetic surface. Simple regret is the gap to the best point on this fixed grid, " + "not a biological optimum. Repeat expected improvement is a negative control that " + "re-proposes the same candidate inside a batch." + ), + } + + +def summarize(objective: str, seeds: int, **kwargs) -> dict: + fn = OBJECTIVES[objective] + by_policy = {} + for policy in POLICIES: + runs = [campaign(fn, policy, seed=seed, **kwargs) for seed in range(seeds)] + by_policy[policy] = { + "mean_best_true": round(sum(item["best_true"] for item in runs) / seeds, 4), + "mean_final_simple_regret": round(sum(item["final_simple_regret"] for item in runs) / seeds, 4), + "mean_cumulative_regret": round(sum(item["cumulative_regret"] for item in runs) / seeds, 4), + "total_violations": int(sum(item["violations"] for item in runs)), + "mean_unique": round(sum(item["n_unique"] for item in runs) / seeds, 2), + } + return { + "objective": objective, + "seeds": seeds, + "grid_size": int(len(candidate_grid())), + "policies": by_policy, + "current_planner": "expected_improvement", + "note": ( + "Scores are on a synthetic response surface, not cells. " + "Bounds are published in-vitro windows, not doses for a person. " + "Do not plate a proposal from this bake-off." + ), + } + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description="Acquisition bake-off on the synthetic media surface") + parser.add_argument("--objective", choices=sorted(OBJECTIVES), default="cardiac") + parser.add_argument("--seeds", type=int, default=6) + parser.add_argument("--rounds", type=int, default=4) + parser.add_argument("--batch", type=int, default=2) + parser.add_argument("--init", type=int, default=4) + args = parser.parse_args(argv) + if args.seeds < 1: + parser.error("seeds must be positive") + json.dump( + summarize(args.objective, args.seeds, n_init=args.init, rounds=args.rounds, batch_size=args.batch), + sys.stdout, + indent=2, + ) + sys.stdout.write("\n") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/medialoop/gp.py b/src/medialoop/gp.py index 5a7244a..f400a23 100644 --- a/src/medialoop/gp.py +++ b/src/medialoop/gp.py @@ -61,3 +61,17 @@ def expected_improvement(mu: np.ndarray, sigma: np.ndarray, best: float, xi: flo improve = mu - best - xi z = improve / sigma return improve * _Phi(z) + sigma * _phi(z) + + +def upper_confidence_bound(mu: np.ndarray, sigma: np.ndarray, kappa: float = 1.5) -> np.ndarray: + """Fixed-kappa UCB. kappa is a baseline, not a calibrated confidence level.""" + if not np.isfinite(kappa) or kappa < 0: + raise ValueError("kappa must be finite and nonnegative") + return np.asarray(mu, dtype=float) + kappa * np.maximum(np.asarray(sigma, dtype=float), 0.0) + + +def thompson_draw(mu: np.ndarray, sigma: np.ndarray, rng: np.random.Generator) -> np.ndarray: + """One posterior draw per candidate. Not a multi-sample Thompson average.""" + sigma = np.maximum(np.asarray(sigma, dtype=float), 1e-9) + return rng.normal(np.asarray(mu, dtype=float), sigma) + diff --git a/tests/test_bakeoff.py b/tests/test_bakeoff.py new file mode 100644 index 0000000..816b5f3 --- /dev/null +++ b/tests/test_bakeoff.py @@ -0,0 +1,39 @@ +import unittest + +from medialoop.bakeoff import POLICIES, campaign, candidate_grid +from medialoop.loop import inside_bounds +from medialoop.surfaces import cardiac + + +class BakeoffTests(unittest.TestCase): + def test_grid_contains_the_cardiac_cartoon_peak_and_stays_in_bounds(self): + grid = candidate_grid() + self.assertTrue(inside_bounds(grid)) + self.assertIn((8.0, 4.0, 0.0, 0.0), {tuple(row) for row in grid.tolist()}) + + def test_guarded_policies_do_not_repeat_or_leave_the_box(self): + for policy in ("random", "expected_improvement", "ucb", "thompson"): + result = campaign(cardiac, policy, seed=1, n_init=4, rounds=3, batch_size=2) + self.assertEqual(result["violations"], 0, policy) + self.assertEqual(result["n_unique"], result["n_evaluations"], policy) + self.assertGreaterEqual(result["final_simple_regret"], -1e-9) + + def test_repeat_policy_records_batch_violations(self): + rounds, batch = 3, 2 + result = campaign(cardiac, "repeat_expected_improvement", seed=0, n_init=4, rounds=rounds, batch_size=batch) + self.assertEqual(result["violations"], (batch - 1) * rounds) + self.assertLess(result["n_unique"], result["n_evaluations"]) + + def test_expected_improvement_regrets_less_than_random_on_the_cartoon(self): + ei = [campaign(cardiac, "expected_improvement", seed=s, rounds=5, batch_size=2)["final_simple_regret"] for s in range(6)] + rnd = [campaign(cardiac, "random", seed=s, rounds=5, batch_size=2)["final_simple_regret"] for s in range(6)] + self.assertLess(sum(ei) / len(ei), sum(rnd) / len(rnd) - 0.05) + + def test_unknown_policy_is_rejected(self): + self.assertNotIn("dqn", POLICIES) + with self.assertRaises(ValueError): + campaign(cardiac, "dqn", seed=0) + + +if __name__ == "__main__": + unittest.main()