培地成分を対象としたベイズ最適化の活用入門（バッチ最適化編）

!pip -q install imageio openpyxl scikit-optimize

import pandas as pd
from io import StringIO

csv_text = """Glucose (g),NH4Cl (g),0D600
5.0,5.0,1.079207
5.0,95.0,0.000000
95.0,5.0,0.000000
95.0,95.0,0.000000
50.0,50.0,0.510030
"""

# 文字列CSV → DataFrame
df_raw = pd.read_csv(StringIO(csv_text))

# 列名を標準化：x=glucose, y=NH4Cl, od=OD600
seed_df = df_raw.rename(columns={"Glucose (g)": "x", "NH4Cl (g)": "y", "0D600": "od"})

# 確認
seed_df.round(3)

from od_objective_2d import ODObjective2D, BOUNDS2D
from skopt import Optimizer
from skopt.space import Real
import numpy as np

# 目的関数（細胞くん）の生成
obj = ODObjective2D()

# 探索範囲（x=glucose, y=NH4Cl の下限/上限）を定義
space = [Real(BOUNDS2D["x"][0], BOUNDS2D["x"][1], name="x"),
         Real(BOUNDS2D["y"][0], BOUNDS2D["y"][1], name="y")]

# ベイズ最適化器の生成
opt = Optimizer(
    dimensions=space,
    base_estimator="GP",  # 内部の代理モデルにガウス過程 (Gaussian Process) を使用
    acq_func="EI",  # 獲得価値関数：改善期待値 (Expected Improvement) を使って「次に測るべき点」を決めます。
    random_state=7
)

best_od = -np.inf
best_xy = (None, None)
rows = []
for i, r in seed_df.reset_index(drop=True).iterrows():
    xi, yi, odi = float(r["x"]), float(r["y"]), float(r["od"])
    opt.tell([xi, yi], -odi)
    if odi > best_od:
        best_od, best_xy = odi, (xi, yi)
    rows.append({"iter": i+1, "phase": "seed", "x": xi, "y": yi, "od": odi,
                 "best_so_far": best_od, "best_x": best_xy[0], "best_y": best_xy[1]})

n_rounds = 10      # ラウンド数（各ラウンドで batch_size 個の条件を提案）
batch_size = 2     # 1ラウンドで何条件提案するか

print("\n=== ベイズ最適化（バッチ提案）を開始 ===")
for r in range(1, n_rounds + 1):
    # 1) ask: 次に測るべき点を batch_size 個まとめて提案
    X_next = opt.ask(n_points=batch_size, strategy="cl_min")
    print(f"[BO round {r}] 提案点:")
    for j, (x_next, y_next) in enumerate(X_next, start=1):
        print(f"  cond{j}: Glucose={x_next:.1f}, NH4Cl={y_next:.1f}")

    # 2) 実験（デモ: 目的関数で代用。実運用では実測ODに置換）
    od_list = [float(obj(xn, yn)) for (xn, yn) in X_next]
    for j, odm in enumerate(od_list, start=1):
        print(f"    → cond{j} 測定 OD= {odm:.3f}")

    # 3) tell: まとめて BO・モデル に渡す
    opt.tell(X_next, [-od for od in od_list])

    # 4) ベスト更新チェック & ログ追加
    for j, ((xn, yn), odm) in enumerate(zip(X_next, od_list), start=1):
        improved = ""
        if odm > best_od:
            best_od, best_xy = odm, (xn, yn)
            improved = "  ← ★ベスト更新！"

        rows.append({
            "iter": len(rows) + 1,
            "phase": "bo",
            "round": r,
            "cond_in_round": j,
            "x": xn,
            "y": yn,
            "od": odm,
            "best_so_far": best_od,
            "best_x": best_xy[0],
            "best_y": best_xy[1],
        })
        print(f"    cond{j}: 現在のbest = {best_od:.3f} @ ({best_xy[0]:.1f}, {best_xy[1]:.1f}){improved}")

# skopt 内部での最良点（-OD を最小化した結果）を確認
best_idx = int(np.argmin(opt.yi))
skopt_best_xy = tuple(opt.Xi[best_idx])
skopt_best_od = -opt.yi[best_idx]

print("\n=== 最終サマリ ===")
print(f"ログ上の best: OD={best_od:.3f} @ ({best_xy[0]:.1f}, {best_xy[1]:.1f})")
print(f"skopt の best: OD={skopt_best_od:.3f} @ ({skopt_best_xy[0]:.1f}, {skopt_best_xy[1]:.1f})")

n_rounds = 10      # ラウンド数（各ラウンドで batch_size 個の条件を提案）
batch_size = 2     # 1ラウンドで何条件提案するか