verify_external_dataset.py

from pathlib import Path
# -*- coding: utf-8 -*-
"""External-dataset validation: 5035 conductivity experiments (Zenodo 7244939).

Track A : formulation-level @ 20 degC  (504 formulations, 4 comps -> ln sigma)
Track B : all rows (5035, 4 comps + temperature -> ln sigma)

Experiments
  Exp-1 : 5-fold CV ceiling R^2 (UWE surrogate vs HGB vs GP) -> noise floor
  Exp-2 : surrogate R^2 vs n in {47,100,200,504} (LOOCV/10-fold) -> small-sample info loss
  Exp-3 : active-learning retrospective @ n=47 pools (electrolyte scale): 6 pools x 12 seeds
          strategies: ucb_ensemble / gp_ucb / gp_oob / pure_unc / random
  Exp-4 : AL retrospective on full 5035 pool (n0=20, budget=25, 6 seeds)
"""
import os, time, warnings, json
import numpy as np
import pandas as pd
warnings.filterwarnings("ignore")
os.environ.setdefault("LOKY_MAX_CPU_COUNT", "1")

from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel
from sklearn.model_selection import KFold, LeaveOneOut
from sklearn.metrics import r2_score

from smallmatprep.modeling.uwe import make_default_models, uwe_fit_predict
from smallmatprep.active.acquirers import ACQUISITION_REGISTRY

RAW = str(Path(__file__).resolve().parent / "data" / "conductivity_5035.csv")
RESULTS = str(Path(__file__).resolve().parent / "results")
COMP = ["PC", "EC", "EMC", "LiPF_6"]

def load_clean():
    df = pd.read_csv(RAW, sep=";")
    df = df.iloc[2:].copy()  # drop symbol/unit rows
    num = ["temperature", "PC", "EC", "EMC", "LiPF_6", "EIS_conductivity"]
    for c in num:
        df[c] = pd.to_numeric(df[c], errors="coerce")
    df = df.dropna(subset=num).reset_index(drop=True)
    df["ln_cond"] = np.log(df["EIS_conductivity"])
    return df

def track_a(df):
    d = df[df["temperature"] == 20.0].reset_index(drop=True)
    X = d[COMP].values.astype(float)
    y = d["ln_cond"].values.astype(float)
    return X, y

def track_b(df):
    X = df[COMP + ["temperature"]].values.astype(float)
    y = df["ln_cond"].values.astype(float)
    return X, y

def cv_uwe(X, y, cv, seed=0):
    """UWE ensemble evaluated with a CV splitter (no leakage)."""
    models = make_default_models(seed)
    preds, truths = [], []
    for tr, va in cv.split(X):
        uw = (len(X[tr]) <= 500)  # large-n: skip expensive LOO weight updates
        mu, _, _, _ = uwe_fit_predict(models, X[tr], y[tr], X[va], update_weights=uw)
        preds.append(mu); truths.append(y[va])
    preds = np.concatenate(preds); truths = np.concatenate(truths)
    return r2_score(truths, preds), float(np.sqrt(np.mean((truths - preds) ** 2)))

def cv_hgb(X, y, cv, seed=0):
    m = HistGradientBoostingRegressor(max_iter=300, learning_rate=0.05, max_depth=3,
                                      random_state=seed, early_stopping=False)
    preds, truths = [], []
    for tr, va in cv.split(X):
        m.fit(X[tr], y[tr]); preds.append(m.predict(X[va])); truths.append(y[va])
    preds = np.concatenate(preds); truths = np.concatenate(truths)
    return r2_score(truths, preds), float(np.sqrt(np.mean((truths - preds) ** 2)))

def cv_gp(X, y, cv, seed=0):
    kernel = ConstantKernel(1.0, (1e-3, 1e3)) * RBF(np.ones(X.shape[1]), (1e-2, 1e2)) + WhiteKernel(0.1, (1e-3, 1e2))
    preds, truths = [], []
    for tr, va in cv.split(X):
        gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True, alpha=1e-6, random_state=seed)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            gp.fit(X[tr], y[tr])
        preds.append(gp.predict(X[va])); truths.append(y[va])
    preds = np.concatenate(preds); truths = np.concatenate(truths)
    return r2_score(truths, preds), float(np.sqrt(np.mean((truths - preds) ** 2)))

def al_loop(X, y_ln, sigma_pool, strategy, n0, budget, seed):
    """Retrospective AL loop using framework acquirers. Thresholds on raw sigma scale."""
    rng = np.random.default_rng(seed)
    n = len(y_ln)
    idx = np.arange(n)
    obs = list(rng.choice(idx, size=n0, replace=False))
    Xo = np.asarray(X[obs], dtype=float)
    yo = np.asarray(y_ln[obs], dtype=float)
    sig_max = float(sigma_pool.max())
    t90, t95 = np.log(0.90 * sig_max), np.log(0.95 * sig_max)
    best_ln = float(yo.max())
    succ90 = bool(best_ln >= t90)
    steps95 = None
    curve = []
    for t in range(min(budget, n - n0)):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            scores = np.asarray(ACQUISITION_REGISTRY[strategy](Xo, yo, X, observed_idx=obs, seed=seed + t), dtype=float)
        scores = np.nan_to_num(scores, nan=-np.inf, posinf=1e12, neginf=-np.inf)
        scores[obs] = -np.inf
        pick = int(np.argmax(scores))
        obs.append(pick)
        Xo = np.vstack([Xo, X[pick]])
        yo = np.append(yo, y_ln[pick])
        best_ln = max(best_ln, float(y_ln[pick]))
        if not succ90 and best_ln >= t90:
            succ90 = True
        if steps95 is None and best_ln >= t95:
            steps95 = t + 1
        curve.append(best_ln)
    return {"best_ln": best_ln, "succ90": succ90, "steps95": steps95,
            "best_pct": float(np.exp(best_ln) / sig_max * 100.0), "curve": np.array(curve)}

def main():
    t0 = time.time()
    df = load_clean()
    Xa, ya = track_a(df)
    Xb, yb = track_b(df)
    print(f"Track A: {len(ya)} formulations @20C, d={Xa.shape[1]}, ln_cond [{ya.min():.3f},{ya.max():.3f}]")
    print(f"Track B: {len(yb)} rows, d={Xb.shape[1]}, ln_cond [{yb.min():.3f},{yb.max():.3f}]")
    out = {}
    summary = {"track_a_n": int(len(ya)), "track_b_n": int(len(yb))}

    # ---------------- Exp-1: CV ceiling ----------------
    print("\n=== Exp-1  CV ceiling ===")
    kf5 = KFold(n_splits=5, shuffle=True, random_state=42)
    exp1_rows = []
    for name, X, y in [("A_20C", Xa, ya), ("B_allT", Xb, yb)]:
        r2u, rmse_u = cv_uwe(X, y, kf5)
        r2h, rmse_h = cv_hgb(X, y, kf5)
        row = {"track": name, "n": len(y), "d": X.shape[1],
               "uwe_r2": r2u, "uwe_rmse": rmse_u, "hgb_r2": r2h, "hgb_rmse": rmse_h}
        if name == "A_20C":
            r2g, rmse_g = cv_gp(X, y, kf5)
            row["gp_r2"], row["gp_rmse"] = r2g, rmse_g
        exp1_rows.append(row)
        best = min(row[k] for k in ["uwe_rmse", "hgb_rmse"] if k in row)
        r2_best = max(row[k] for k in ["uwe_r2", "hgb_r2", "gp_r2"] if k in row)
        var_y = float(y.var())
        row["r2_ceiling"] = r2_best
        row["noise_floor_rmse"] = best
        row["r2_noise_bound"] = 1.0 - best ** 2 / var_y
        print(f"  {name}: UWE r2={r2u:.3f} HGB r2={r2h:.3f}" + (f" GP r2={r2g:.3f}" if name == "A_20C" else "") +
              f" | ceiling r2={r2_best:.3f} noise-bound r2={row['r2_noise_bound']:.3f}")
    exp1 = pd.DataFrame(exp1_rows)
    exp1.to_csv(os.path.join(RESULTS, "external_exp1_ceiling.csv"), index=False, encoding="utf-8-sig")
    out["exp1"] = exp1_rows

    # ---------------- Exp-2: R2 vs n (small-sample info loss) ----------------
    print("\n=== Exp-2  surrogate R2 vs n ===")
    ns = [47, 100, 200, 504]
    exp2_rows = []
    for n in ns:
        if n < len(ya):
            rng = np.random.default_rng(123)
            sub = rng.choice(len(ya), size=n, replace=False)
            Xs, ys = Xa[sub], ya[sub]
        else:
            Xs, ys = Xa, ya
        if n <= 47:
            cv = LeaveOneOut()
        else:
            cv = KFold(n_splits=10, shuffle=True, random_state=42)
        r2u, rmse_u = cv_uwe(Xs, ys, cv)
        exp2_rows.append({"n": n, "r2_uwe": r2u, "rmse_uwe": rmse_u})
        print(f"  n={n}: UWE r2={r2u:.3f} rmse={rmse_u:.4f}")
    exp2 = pd.DataFrame(exp2_rows)
    exp2.to_csv(os.path.join(RESULTS, "external_exp2_ncurve.csv"), index=False, encoding="utf-8-sig")
    out["exp2"] = exp2_rows

    # ---------------- Exp-3: AL retrospective @ n=47 pools ----------------
    print("\n=== Exp-3  AL loop @ n=47 pools (electrolyte scale) ===")
    strategies = ["ucb_ensemble", "gp_ucb", "gp_oob", "pure_unc", "random"]
    n0, budget, n_seed, n_pool = 8, 12, 12, 6
    rec = {s: [] for s in strategies}
    for pi in range(n_pool):
        rng = np.random.default_rng(1000 + pi)
        sub = rng.choice(len(ya), size=47, replace=False)
        Xp, yp, sp = Xa[sub], ya[sub], np.exp(ya[sub])
        for s in strategies:
            for seed in range(n_seed):
                r = al_loop(Xp, yp, sp, s, n0, budget, seed)
                rec[s].append(r)
    exp3_rows = []
    for s in strategies:
        bp = np.array([r["best_pct"] for r in rec[s]])
        s90 = np.array([r["succ90"] for r in rec[s]])
        st = np.array([r["steps95"] if r["steps95"] is not None else np.inf for r in rec[s]])
        st_f = st[np.isfinite(st)]
        exp3_rows.append({
            "strategy": s,
            "mean_best_pct": float(bp.mean()), "std_best_pct": float(bp.std()),
            "succ90_pct": float(s90.mean() * 100),
            "steps95_med": float(np.median(st_f)) if len(st_f) else float("nan"),
            "steps95_reach": int(len(st_f)), "runs": int(len(bp)),
        })
        print(f"  {s:14s} mean_best={bp.mean():6.2f}%  succ90={s90.mean()*100:3.0f}%  "
              f"steps95 med={np.median(st_f) if len(st_f) else float('nan'):5.1f} ({len(st_f)}/{len(bp)})")
    exp3 = pd.DataFrame(exp3_rows)
    exp3.to_csv(os.path.join(RESULTS, "external_exp3_al47.csv"), index=False, encoding="utf-8-sig")
    out["exp3"] = exp3_rows

    # ---------------- Exp-4: AL on full 5035 pool ----------------
    print("\n=== Exp-4  AL loop @ full 5035 pool ===")
    strat4 = ["ucb_ensemble", "gp_ucb", "random"]
    n0b, budget_b, n_seed_b = 20, 25, 6
    sp_b = np.exp(yb)
    rec4 = {s: [] for s in strat4}
    for s in strat4:
        for seed in range(n_seed_b):
            r = al_loop(Xb, yb, sp_b, s, n0b, budget_b, seed)
            rec4[s].append(r)
    exp4_rows = []
    for s in strat4:
        bp = np.array([r["best_pct"] for r in rec4[s]])
        s90 = np.array([r["succ90"] for r in rec4[s]])
        st = np.array([r["steps95"] if r["steps95"] is not None else np.inf for r in rec4[s]])
        st_f = st[np.isfinite(st)]
        exp4_rows.append({
            "strategy": s,
            "mean_best_pct": float(bp.mean()), "std_best_pct": float(bp.std()),
            "succ90_pct": float(s90.mean() * 100),
            "steps95_med": float(np.median(st_f)) if len(st_f) else float("nan"),
            "steps95_reach": int(len(st_f)), "runs": int(len(bp)),
        })
        print(f"  {s:14s} mean_best={bp.mean():6.2f}%  succ90={s90.mean()*100:3.0f}%  "
              f"steps95 med={np.median(st_f) if len(st_f) else float('nan'):5.1f} ({len(st_f)}/{len(bp)})")
    exp4 = pd.DataFrame(exp4_rows)
    exp4.to_csv(os.path.join(RESULTS, "external_exp4_al5035.csv"), index=False, encoding="utf-8-sig")
    out["exp4"] = exp4_rows

    # ---------------- save summary ----------------
    summary["elapsed_s"] = round(time.time() - t0, 1)
    with open(os.path.join(RESULTS, "external_validation_summary.json"), "w", encoding="utf-8") as f:
        json.dump(summary, f, ensure_ascii=False, indent=2, default=float)
    print(f"\n总耗时 {time.time()-t0:.1f}s")

if __name__ == "__main__":
    main()