smallmatprep/verify_quantile_smooth.py at master · wenyu2026/smallmatprep · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
# -*- coding: utf-8 -*-
"""验证：大池默认改 quantile_ucb（无多样性）在平滑合成池上是否回退
（复刻 verify_acquisition_ucb.py 的池生成与 loop_once，新增 quantile 分支）
"""
import os, warnings, time
import numpy as np
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

rng = np.random.default_rng(20260804)
T = 100.0

def set_seed(s): np.random.seed(s)

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

def make_models():
    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=0, n_jobs=1))]

def loop_once(strategy, pool, y_pool, n0=12, budget=20, seed=0, sigma=None):
    set_seed(seed)
    idx = np.arange(len(pool))
    observed = list(rng.choice(idx, size=n0, replace=False))
    if sigma is None: sigma = 0.08 * np.std(y_pool)
    Xc, yc = pool[observed].copy(), y_pool[observed] + rng.normal(0, sigma, len(observed))
    best = yc.max()
    models = make_models(); 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 [m.fit(Xc, yc).predict(pool) for _, m in models]
    for t in range(budget):
        preds = fit(np.asarray(Xc), np.asarray(yc), do_loo=(t % 5 == 0))
        mu = sum(w[k] * preds[k] for k in range(len(models)))
        unc = np.sqrt(sum(w[k] * (preds[k] - mu) ** 2 for k in range(len(models))))
        mask = np.isin(idx, observed)
        if strategy == "random":
            pick = int(rng.choice(idx[~mask]))
        elif strategy == "ucb_div":   # 现状：0.5μ+0.3σ+0.2D（多样性）
            mu_n = (mu - mu.min()) / (mu.max() - mu.min() + 1e-12)
            unc_n = unc / (unc.max() + 1e-12)
            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; score[mask] = -np.inf
            pick = int(np.argmax(score))
        elif strategy == "quantile":  # 新默认：μ+1.645σ，无多样性
            score = mu + 1.645 * unc; score[mask] = -np.inf
            pick = int(np.argmax(score))
        else:  # gp_bo
            from sklearn.gaussian_process import GaussianProcessRegressor
            from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C, WhiteKernel
            gp = GaussianProcessRegressor(kernel=C(1.0, (1e-3, 1e3)) * RBF(np.ones(pool.shape[1]), (1e-2, 1e2)) + WhiteKernel(0.1, (1e-3, 1e2)),
                                          normalize_y=True, alpha=1e-6, random_state=0)
            gp.fit(np.asarray(Xc), np.asarray(yc))
            m, s = gp.predict(pool, return_std=True)
            score = m + 2.0 * s; 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

if __name__ == "__main__":
    t0 = time.time()
    pool_size, d, n0, budget, n_seed = 300, 6, 12, 20, 10
    set_seed(77)
    pool = rng.dirichlet(np.ones(d), size=pool_size) * T
    y_pool = true_f(pool); ymax = y_pool.max()
    results = {}
    for s in ["random", "ucb_div", "quantile", "gp_bo"]:
        vals = [loop_once(s, pool, y_pool, n0=n0, budget=budget, seed=rep) for rep in range(n_seed)]
        arr = np.array(vals); results[s] = arr
        print(f"{s:10s} mean={arr.mean():.4f} success90={(arr >= 0.90*ymax).mean()*100:.0f}% "
              f"vs_rand={(arr.mean()/results['random'].mean()-1)*100:+.1f}%")
    print(f"pool max = {ymax:.4f}, 耗时 {time.time()-t0:.1f}s")