verify_optimizations.py

# -*- coding: utf-8 -*-
"""Verify four information-utilization optimizations for the unified loop.

Exp1 (L1): profile Bayesian shrinkage (median/mean profile -> shrink to uniform)
Exp2 (L2): UWE weight shrinkage (ML inverse-variance -> blend to uniform)
Exp3 (L3): EI vs UCB vs pure-unc acquisition on the synthetic 300-pool benchmark
Exp4 (L3): multi-target sharing on the electrolyte dataset P2 (chained auxiliary surrogates)
"""
import os, sys, time, warnings
import numpy as np
import pandas as pd
from scipy.stats import norm
warnings.filterwarnings("ignore")
os.environ.setdefault("LOKY_MAX_CPU_COUNT", "1")

from sklearn.linear_model import Ridge
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import LeaveOneOut
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C, WhiteKernel

T = 100.0
DATA = str(Path(__file__).resolve().parent / "data" / "electrolyte_data.csv")
COMP = list("ABCDEFGHIJKLMNOPQRSTUVWX")
AUX_COLS = ["P9", "P15", "P10", "P5", "P17", "P3", "P14"]
OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results")

def make_models(seed=0):
    return [
        ("Ridge10", Ridge(alpha=10.0)),
        ("Ridge1", Ridge(alpha=1.0)),
        ("KNN3", KNeighborsRegressor(n_neighbors=3)),
        ("KNN5", KNeighborsRegressor(n_neighbors=5)),
        ("RF", RandomForestRegressor(n_estimators=40, max_depth=3, random_state=seed, n_jobs=1)),
    ]

# ---------------------------------------------------------------- Exp 1
def exp1_profile_shrinkage(rng):
    d, n_repeat, n_eval, miss_rate = 6, 10, 100, 0.20
    n_c_vals = [4, 8, 16, 32]
    kappas = [0.0, 1.0, 3.0, 10.0, 30.0]
    rows = []
    for rep in range(n_repeat):
        alpha = rng.uniform(0.5, 5.0, d)
        for n_c in n_c_vals:
            Xc = rng.dirichlet(alpha, n_c) * T
            Xe = rng.dirichlet(alpha, n_eval) * T
            mask = rng.random(Xe.shape) < miss_rate
            Xm = Xe.copy(); Xm[mask] = np.nan
            med = np.median(Xc / T, axis=0)
            cnt = (Xc / T).mean(axis=0)
            oracle = alpha / alpha.sum()
            for kap in kappas:
                p_med = (n_c * med + kap / d) / (n_c + kap)
                p_cnt = (n_c * cnt + kap / d) / (n_c + kap)
                mae_med = _impute_mae(Xe, Xm, mask, p_med)
                mae_cnt = _impute_mae(Xe, Xm, mask, p_cnt)
                rows.append(dict(rep=rep, n_c=n_c, kappa=kap, profile="median", mae=mae_med))
                rows.append(dict(rep=rep, n_c=n_c, kappa=kap, profile="count", mae=mae_cnt))
            rows.append(dict(rep=rep, n_c=n_c, kappa=np.nan, profile="oracle", mae=_impute_mae(Xe, Xm, mask, oracle)))
    df = pd.DataFrame(rows)
    tbl = df.groupby(["n_c", "kappa", "profile"])["mae"].mean().unstack(level=[1, 2])
    return df, tbl

def _impute_mae(Xe, Xm, mask, p):
    Xf = Xm.copy()
    errs = []
    for i in range(Xm.shape[0]):
        miss = np.isnan(Xm[i])
        if not miss.any():
            continue
        known_sum = np.nansum(Xm[i])
        R = T - known_sum
        ps = p[miss]
        if ps.sum() <= 0:
            continue
        Xf[i, miss] = ps / ps.sum() * R
        errs.append(np.abs(Xf[i, miss] - Xe[i, miss]))
    if not errs:
        return np.nan
    return float(np.mean(np.concatenate(errs)))

# ---------------------------------------------------------------- Exp 2
def true_f(X):
    x = X / T
    return 1.5 * x[:, 0] + 20.0 * x[:, 1] * x[:, 2] + 0.8 * np.sin(6 * np.pi * x[:, 3]) + 0.5 * x[:, 4]

def exp2_uwe_shrinkage(rng):
    d, n_train, n_test, n_repeat = 6, 30, 200, 40
    sigma_fracs = [0.1, 0.3]
    lambdas = [0.0, 0.25, 0.5, 0.75, 1.0]
    rows = []
    Xt = rng.dirichlet(np.ones(d), n_test) * T
    yt = true_f(Xt)
    for sf in sigma_fracs:
        for rep in range(n_repeat):
            Xtr = rng.dirichlet(np.ones(d), n_train) * T
            ytr0 = true_f(Xtr)
            sigma = sf * np.std(ytr0)
            ytr = ytr0 + rng.normal(0, sigma, n_train)
            preds = np.column_stack([m.fit(Xtr, ytr).predict(Xt) for _, m in make_models(rep)])
            # LOOCV variances -> ML weights
            vs = []
            for _, m in make_models(rep):
                e = []
                for tr, va in LeaveOneOut().split(Xtr):
                    m.fit(Xtr[tr], ytr[tr])
                    e.append((ytr[va] - m.predict(Xtr[va]))[0])
                vs.append(np.mean(np.array(e) ** 2))
            vs = np.array(vs) + 1e-12
            w_ml = 1.0 / vs; w_ml = w_ml / w_ml.sum()
            for lam in lambdas:
                w = lam * w_ml + (1.0 - lam) * np.ones(len(vs)) / len(vs)
                mu = preds @ w
                rows.append(dict(sigma_frac=sf, rep=rep, lam=lam, rmse=float(np.sqrt(np.mean((mu - yt) ** 2)))))
    df = pd.DataFrame(rows)
    tbl = df.groupby(["sigma_frac", "lam"])["rmse"].agg(["mean", "std"]).round(5)
    return df, tbl

# ---------------------------------------------------------------- Exp 3
def _gp_fit(Xc, yc, seed, d):
    gp = GaussianProcessRegressor(
        kernel=C(1.0, (1e-3, 1e3)) * RBF(np.ones(d), (1e-2, 1e2)) + WhiteKernel(0.1, (1e-3, 1e2)),
        normalize_y=True, alpha=1e-6, random_state=seed)
    gp.fit(np.asarray(Xc, dtype=float), np.asarray(yc, dtype=float))
    return gp

def exp3_ei_vs_ucb(rng):
    pool_size, d, n0, budget, n_seed = 300, 6, 12, 20, 10
    pool = rng.dirichlet(np.ones(d), size=pool_size) * T
    y_pool = true_f(pool)
    ymax = y_pool.max()
    strategies = ["gp_ei", "gp_ucb15", "gp_ucb20", "ucb_ensemble", "pure_unc", "random"]
    results = {s: [] for s in strategies}
    for s in strategies:
        for rep in range(n_seed):
            results[s].append(_loop_once_e3b(s, pool, y_pool, d, n0, budget, rep))
    out = {}
    for s in strategies:
        arr = np.array(results[s])
        out[s] = dict(mean=arr.mean(), std=arr.std(),
                      succ90=(arr >= 0.90 * ymax).mean() * 100,
                      vs_random=(arr.mean() / np.mean(results["random"]) - 1) * 100)
    return out, ymax

def _loop_once_e3b(strategy, pool, y_pool, d, n0, budget, seed):
    rng = np.random.default_rng(seed)
    idx = np.arange(len(pool))
    observed = list(rng.choice(idx, size=n0, replace=False))
    sigma = 0.08 * np.std(y_pool)
    Xc = pool[observed].copy(); yc = y_pool[observed] + rng.normal(0, sigma, len(observed))
    best = yc.max()
    models = make_models(seed); w = np.ones(len(models)) / len(models)
    def fit(Xc, yc, do_loo):
        nonlocal w
        if do_loo:
            vs = []
            for _, m in models:
                e = []
                for tr, va in LeaveOneOut().split(Xc):
                    m.fit(Xc[tr], yc[tr]); e.append((yc[va] - m.predict(Xc[va]))[0])
                vs.append(np.mean(np.array(e) ** 2))
            vs = np.array(vs) + 1e-12
            w = 1.0 / vs; w = w / w.sum()
        return np.column_stack([m.fit(Xc, yc).predict(pool) for _, m in models])
    for t in range(budget):
        mask = np.isin(idx, observed)
        if strategy in ("gp_ei", "gp_ucb15", "gp_ucb20"):
            gp = _gp_fit(Xc, yc, seed, d)
            mu, s = gp.predict(pool, return_std=True)
            if strategy == "gp_ei":
                z = (mu - best) / (s + 1e-12)
                score = (mu - best) * norm.cdf(z) + s * norm.pdf(z)
            else:
                kap = 1.5 if strategy == "gp_ucb15" else 2.0
                score = mu + kap * s
            score[mask] = -np.inf
            pick = int(np.argmax(score))
        else:
            preds = fit(np.asarray(Xc), np.asarray(yc), do_loo=(t % 5 == 0))
            mu = preds @ w
            unc = np.sqrt(np.maximum(np.sum(w[None, :] * (preds - mu[:, None]) ** 2, axis=1), 1e-12))
            if strategy == "random":
                pick = int(rng.choice(idx[~mask]))
            else:
                mu_n = (mu - mu.min()) / (mu.max() - mu.min() + 1e-12)
                unc_n = unc / (unc.max() + 1e-12)
                if strategy == "ucb_ensemble":
                    D = np.full(len(pool), np.inf)
                    for o in pool[observed]:
                        D = np.minimum(D, np.linalg.norm(pool - o, axis=1))
                    D_n = D / (D.max() + 1e-12)
                    score = 0.5 * mu_n + 0.3 * unc_n + 0.2 * D_n
                else:
                    score = unc_n.copy()
                score[mask] = -np.inf
                pick = int(np.argmax(score))
        observed.append(pick)
        y_new = y_pool[pick] + rng.normal(0, sigma)
        Xc = np.vstack([Xc, pool[pick]]); yc = np.append(yc, y_new)
        best = max(best, y_new)
    return best

# ---------------------------------------------------------------- Exp 4
def exp4_multitarget(rng):
    df = pd.read_csv(DATA, encoding="utf-8-sig")
    X = df[COMP].fillna(0.0).values.astype(float)
    y_main = df["P2"].values.astype(float)
    valid = np.isfinite(y_main)
    X, y_main = X[valid], y_main[valid]
    y_aux = {}
    for c in AUX_COLS:
        col = df[c].values.astype(float)
        y_aux[c] = col[valid]
    n0, budget, n_seed = 8, 12, 12
    arms = ["single_pureunc", "single_gp", "multi_pureunc", "multi_gp"]
    results = {a: [] for a in arms}
    steps = {a: [] for a in arms}
    for a in arms:
        for rep in range(n_seed):
            best, curve, st = _loop_once_electrolyte(a, X, y_main, y_aux, n0, budget, rep)
            results[a].append(best); steps[a].append(st if st is not None else np.nan)
    ymax = y_main.max()
    out = {}
    for a in arms:
        arr = np.array(results[a]); st = np.array(steps[a]); finite = st[np.isfinite(st)]
        out[a] = dict(mean=arr.mean(), succ90=(arr >= 0.90 * ymax).mean() * 100,
                      steps_med=np.median(finite) if len(finite) else np.nan,
                      reach=len(finite))
    return out, ymax

def _loop_once_electrolyte(strategy, X, y_main, y_aux, n0, budget, seed):
    rng = np.random.default_rng(seed)
    n = len(y_main); idx = np.arange(n)
    observed = list(rng.choice(idx, size=n0, replace=False))
    best = y_main[observed].max()
    ymax = y_main.max(); target = 0.95 * ymax
    steps = None; curve = []
    def aug_features(X_all, obs):
        cols = []
        for c in AUX_COLS:
            yc = y_aux[c]
            m_obs = np.isin(np.arange(len(yc)), obs) & np.isfinite(yc)
            if m_obs.sum() < 5:
                cols.append(np.zeros(len(X_all))); continue
            m = Ridge(alpha=10.0)
            m.fit(X_all[m_obs], yc[m_obs])
            cols.append(m.predict(X_all))
        return np.column_stack(cols)
    X_cur = X if not strategy.startswith("multi") else aug_features(X, observed)
    for t in range(min(budget, n - n0)):
        Xo = X[observed]; yo = y_main[observed]
        if strategy.startswith("multi"):
            Xo_aug = aug_features(X, observed)[observed]
            Xall = aug_features(X, observed)
            Xtr, Xpool = Xo_aug, Xall
        else:
            Xtr, Xpool = Xo, X_cur
        if strategy.endswith("gp"):
            gp = _gp_fit(Xtr, yo, seed, Xtr.shape[1])
            mu, s = gp.predict(Xpool, return_std=True)
            score = mu + 2.0 * s
        else:
            preds = np.column_stack([m.fit(Xtr, yo).predict(Xpool) for _, m in make_models(seed)])
            vs = []
            for _, m in make_models(seed):
                e = []
                for tr, va in LeaveOneOut().split(Xtr):
                    m.fit(Xtr[tr], yo[tr]); e.append((yo[va] - m.predict(Xtr[va]))[0])
                vs.append(np.mean(np.array(e) ** 2))
            w = 1.0 / (np.array(vs) + 1e-12); w = w / w.sum()
            mu = preds @ w
            unc = np.sqrt(np.maximum(np.sum(w[None, :] * (preds - mu[:, None]) ** 2, axis=1), 1e-12))
            score = unc / (unc.max() + 1e-12)
        score = np.asarray(score, dtype=float).copy()
        score[np.isin(idx, observed)] = -np.inf
        pick = int(np.argmax(score))
        observed.append(pick)
        best = max(best, y_main[pick]); curve.append(best)
        if steps is None and best >= target:
            steps = float(t + 1)
    return best, np.array(curve), steps

if __name__ == "__main__":
    t0 = time.time()
    os.makedirs(OUT, exist_ok=True)
    rng = np.random.default_rng(20260804)
    print("=== Exp1: profile shrinkage (L1) — MAE vs kappa, by n_c ===")
    df1, tbl1 = exp1_profile_shrinkage(rng)
    print(tbl1.round(3).to_string())
    df1.to_csv(os.path.join(OUT, "opt_exp1_profile_shrinkage.csv"), index=False)
    print("\n=== Exp2: UWE weight shrinkage (L2) — test RMSE ===")
    df2, tbl2 = exp2_uwe_shrinkage(rng)
    print(tbl2.to_string())
    df2.to_csv(os.path.join(OUT, "opt_exp2_uwe_shrinkage.csv"), index=False)
    print("\n=== Exp3: EI vs UCB (L3) ===")
    out3, ymax3 = exp3_ei_vs_ucb(rng)
    for s, v in out3.items():
        print(f"  {s:14s} mean={v['mean']:.4f} succ90={v['succ90']:.0f}% vs_random={v['vs_random']:+.1f}%")
    print(f"  pool max = {ymax3:.4f}")
    print("\n=== Exp4: multi-target sharing (L3) on the electrolyte dataset P2 ===")
    out4, ymax4 = exp4_multitarget(rng)
    for a, v in out4.items():
        print(f"  {a:16s} mean_best={v['mean']:6.2f} succ90={v['succ90']:.0f}% "
              f"steps95_med={v['steps_med']:.1f} reach={v['reach']}/12")
    print(f"  P2 max = {ymax4:.2f}")
    print(f"\n总耗时 {time.time()-t0:.1f}s")