Day 166 — Playground Series 実戦投入④ — 性質の違うモデルを混ぜる。log_Dur は Ridge に −0.047 と桁違いに効いたが、Ridge の重みは 0。OOF で見えた混ぜる改善 −0.00028 は、提出する形（5モデル平均）では −0.00006 に縮んだ

oof = np.zeros((4, n_models, N_TRAIN)); tp = np.zeros((4, n_models, N_TEST))
for s in range(4):                                   # 分割seed 0〜3
    for k, (tr, va) in enumerate(KFold(5, shuffle=True, random_state=s).split(train)):
        for j, name in enumerate(MODELS):
            pv, pt = fit_predict(name, tr_df, y_log, [train.iloc[va], test])
            oof[s, j, va] = pv; tp[s, j] += pt / 5   # tp = 5モデル平均
w = best_w(y, oof[0, ids], grid(3))                  # 重みは分割seed 0 の OOF だけで決める

"""Day166: Playground S5E5 縮小シミュレータで「性質の違うモデルを混ぜる」を測る。
モデル3つ（HistGB / Ridge+log特徴量 / k近傍）を、分割seed 0〜3 の5分割で OOF 予測し、
(1) Ridge に log_Dur を足すと効くか (2) 混ぜる重みを OOF で決めたときの楽観 (3) 重みの安定性 を数える。
使い方: OMP_NUM_THREADS=2 python day166.py <world>   -> raw_w<world>.npz
"""
import sys, time
import numpy as np, pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.neighbors import KNeighborsRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold

T0 = time.time()
W = int(sys.argv[1])
N_ORIG, N_TRAIN, N_TEST = 15000, 30000, 60000


def keytel_rate(sex, hr, w, age):
    m = (-55.0969 + 0.6309 * hr + 0.1988 * w + 0.2017 * age) / 4.184
    f = (-20.4022 + 0.4472 * hr - 0.1263 * w + 0.074 * age) / 4.184
    return np.where(sex == 1, m, f)


def make_original(n, rng):  # Day163-165 と同じ
    sex = rng.integers(0, 2, n); age = rng.integers(20, 80, n)
    h = np.where(sex == 1, rng.normal(175, 7, n), rng.normal(162, 6.5, n)).round()
    w = (rng.normal(24.5, 2.5, n) * (h / 100) ** 2).round()
    dur = rng.integers(1, 31, n).astype(float)
    hr = np.clip(80 + 1.2 * dur + rng.normal(0, 5, n), 67, 128).round()
    temp = (37 + 0.35 * dur ** 0.5 + rng.normal(0, 0.2, n)).round(1)
    rate = np.clip(keytel_rate(sex, hr, w, age), 0.3, None)
    cal = np.clip((rate * dur * np.exp(rng.normal(0, 0.05, n))).round(), 1, None)
    return pd.DataFrame(dict(Sex=sex, Age=age, Height=h, Weight=w, Duration=dur,
                             Heart_Rate=hr, Body_Temp=temp, Calories=cal))


def synth_from(orig, n, rng):  # Day163-165 と同じ
    s = orig.sample(n, replace=True, random_state=int(rng.integers(1 << 30))).reset_index(drop=True)
    s["Age"] = np.clip((s.Age + rng.normal(0, 1.5, n)).round(), 20, 79)
    s["Height"] = (s.Height + rng.normal(0, 2, n)).round()
    s["Weight"] = (s.Weight + rng.normal(0, 2, n)).round()
    s["Duration"] = np.clip((s.Duration + rng.normal(0, 1.0, n)).round(), 1, 30)
    s["Heart_Rate"] = (s.Heart_Rate + rng.normal(0, 2, n)).round()
    s["Body_Temp"] = (s.Body_Temp + rng.normal(0, 0.1, n)).round(1)
    s["Calories"] = np.clip((s.Calories * np.exp(rng.normal(0, 0.03, n))).round(), 1, None)
    return s


def add_feats(d):
    d = d.copy()
    d["Dur_x_HR"] = d.Duration * d.Heart_Rate                       # 決定6
    d["Keytel_kcal"] = np.clip(keytel_rate(d.Sex.values, d.Heart_Rate.values, d.Weight.values, d.Age.values),
                               0.3, None) * d.Duration              # 決定6
    d["log_Dur"] = np.log(d.Duration)
    d["log_HR"] = np.log(d.Heart_Rate)
    d["log_W"] = np.log(d.Weight)
    d["log_Keytel"] = np.log(d.Keytel_kcal)
    return d


BASE = ["Sex", "Age", "Height", "Weight", "Duration", "Heart_Rate", "Body_Temp", "Dur_x_HR", "Keytel_kcal"]
SETS = {"hgb": BASE,
        "ridge_raw": BASE,
        "ridge_logDur": BASE + ["log_Dur"],                          # 1列だけ足す（Day164 の問いを線形モデルで）
        "ridge_log": BASE + ["log_Dur", "log_HR", "log_W", "log_Keytel"],
        "knn": ["log_Dur", "Heart_Rate", "Weight", "Age", "Sex"]}
MODELS = list(SETS)


def make(name):
    if name == "hgb":  # 決定7の既定値、モデル側の乱数なし
        return HistGradientBoostingRegressor(early_stopping=False, random_state=0, learning_rate=0.1,
                                             max_iter=300, max_leaf_nodes=31, min_samples_leaf=20)
    if name.startswith("ridge"):
        return Ridge(alpha=1.0)
    return KNeighborsRegressor(n_neighbors=15, algorithm="kd_tree")


def fit_predict(name, tr_df, ytr_log, preds_dfs):
    cols = SETS[name]
    if name == "hgb":
        m = make(name).fit(tr_df[cols], ytr_log)
        return [m.predict(d[cols]) for d in preds_dfs]
    sc = StandardScaler().fit(tr_df[cols])
    m = make(name).fit(sc.transform(tr_df[cols]), ytr_log)
    return [m.predict(sc.transform(d[cols])) for d in preds_dfs]


rng = np.random.default_rng(2000 + W)        # Day164/165 のワールド W と同じデータ
orig = make_original(N_ORIG, rng)
pool = synth_from(orig, N_TRAIN + N_TEST, rng)
orig, pool = add_feats(orig), add_feats(pool)
train, test = pool.iloc[:N_TRAIN].reset_index(drop=True), pool.iloc[N_TRAIN:].reset_index(drop=True)

oof = np.zeros((4, len(MODELS), N_TRAIN))      # [分割seed, モデル, 行]  すべて log1p 空間
tp = np.zeros((4, len(MODELS), N_TEST))
tpf = np.zeros((5, len(MODELS), N_TEST), dtype=np.float32)   # 分割seed1 の fold ごとの test 予測（1モデルずつ）
for s in range(4):
    for k, (tr, va) in enumerate(KFold(5, shuffle=True, random_state=s).split(train)):
        tr_df = pd.concat([train.iloc[tr], orig], ignore_index=True)
        y = np.log1p(tr_df.Calories.values)
        for j, name in enumerate(MODELS):
            pv, pt = fit_predict(name, tr_df, y, [train.iloc[va], test])
            oof[s, j, va] = pv
            tp[s, j] += pt / 5
            if s == 1:
                tpf[k, j] = pt
    print(f"world {W} split {s} done {time.time()-T0:.0f}s", flush=True)
np.savez_compressed(f"raw_w{W}.npz", oof=oof, tp=tp, tpf=tpf, ytr=train.Calories.values, yt=test.Calories.values,
                    models=np.array(MODELS))
print(f"world {W} wall {time.time()-T0:.0f}s")

"""Day166 集計: raw_w*.npz を読み、表1〜5 を出す。使い方: python analyze166.py"""
import glob, itertools
import numpy as np, pandas as pd
from sklearn.model_selection import KFold

pd.set_option("display.width", 250)
files = sorted(glob.glob("raw_w*.npz"))
M = ["hgb", "ridge_raw", "ridge_logDur", "ridge_log", "knn"]
ix = {m: i for i, m in enumerate(M)}


def sle(y, lp):
    return (np.log1p(np.clip(np.expm1(lp), 1, None)) - np.log1p(y)) ** 2


def rmsle(y, lp):
    return np.sqrt(sle(y, lp).mean())


def folds(y, lp, s):
    return np.array([np.sqrt(sle(y[va], lp[va]).mean()) for _, va in KFold(5, shuffle=True, random_state=s).split(y)])


def grid(k, step=0.05):
    n = int(round(1 / step))
    return [np.array(c) * step for c in itertools.product(range(n + 1), repeat=k) if sum(c) == n]


def best_w(y, P, G):  # P: (k, rows)
    sc = [rmsle(y, np.tensordot(w, P, axes=1)) for w in G]
    return G[int(np.argmin(sc))]


def compare(y, yt, oof, tp, a, b):
    """a が b より良いか（a-b、負=a が良い）。新分割R2（分割seed1〜3の15 fold差）と test の真の効果"""
    d15 = np.concatenate([folds(y, oof[s, a], s) - folds(y, oof[s, b], s) for s in (1, 2, 3)])
    dF, seF = d15.mean(), d15.std(ddof=1) / np.sqrt(15)
    pa, pb = tp[1:, a].mean(0), tp[1:, b].mean(0)
    ea, eb = sle(yt, pa), sle(yt, pb)
    dT = np.sqrt(ea.mean()) - np.sqrt(eb.mean())
    seT = (ea - eb).std(ddof=1) / np.sqrt(len(yt)) / (2 * np.sqrt(eb.mean()))
    return dF, seF, dT, seT


G3, G2 = grid(3), grid(2)
t1, t2, t3, t4, t5 = [], [], [], [], []
for f in files:
    w = int(f[-5])
    z = np.load(f); oof, tp, y, yt = z["oof"], z["tp"], z["ytr"], z["yt"]
    # 表1 単独モデルの点数
    for m in M:
        t1.append(dict(world=w, model=m, oof_s0=rmsle(y, oof[0, ix[m]]),
                       test=rmsle(yt, tp[1:, ix[m]].mean(0))))
    # 表2 線形モデルへの log 特徴量
    for a, b in [("ridge_logDur", "ridge_raw"), ("ridge_log", "ridge_raw"), ("ridge_log", "ridge_logDur")]:
        dF, seF, dT, seT = compare(y, yt, oof, tp, ix[a], ix[b])
        t2.append(dict(world=w, pair=f"{a} - {b}", dF=dF * 1e4, zF=dF / seF, dT=dT * 1e4, zT=dT / seT))
    # 表3 混ぜる: 分割seed0 の OOF で重みを決め、新分割R2 と test で測る
    for tag, names, G in [("hgb+ridge_log", ["hgb", "ridge_log"], G2),
                          ("hgb+ridge_log+knn", ["hgb", "ridge_log", "knn"], G3)]:
        ids = [ix[n] for n in names]
        wts = best_w(y, oof[0, ids], G)
        for label, wv in [("OOF最適", wts), ("等分", np.ones(len(ids)) / len(ids))]:
            blend_oof = np.tensordot(wv, oof[:, ids], axes=([0], [1]))  # (4, rows)
            blend_tp = np.tensordot(wv, tp[:, ids], axes=([0], [1]))
            oo, tt = oof.copy(), tp.copy()
            oo[:, 0] = blend_oof; tt[:, 0] = blend_tp                    # 0番目(hgb)を混ぜた予測で置換
            dS = rmsle(y, blend_oof[0]) - rmsle(y, oof[0, 0])
            fs0 = folds(y, blend_oof[0], 0) - folds(y, oof[0, 0], 0)
            seS = fs0.std(ddof=1) / np.sqrt(5)
            # 比較相手の hgb 単独は元の配列から
            d15 = np.concatenate([folds(y, blend_oof[s], s) - folds(y, oof[s, 0], s) for s in (1, 2, 3)])
            dF, seF = d15.mean(), d15.std(ddof=1) / np.sqrt(15)
            ea, eb = sle(yt, blend_tp[1:].mean(0)), sle(yt, tp[1:, 0].mean(0))
            dT = np.sqrt(ea.mean()) - np.sqrt(eb.mean())
            seT = (ea - eb).std(ddof=1) / np.sqrt(len(yt)) / (2 * np.sqrt(eb.mean()))
            t3.append(dict(world=w, mix=tag, weights=label + str(np.round(wv, 2).tolist()), dS=dS * 1e4, zS=dS / seS,
                           dF=dF * 1e4, zF=dF / seF, dT=dT * 1e4, zT=dT / seT,
                           opt_S=(dT - dS) * 1e4, opt_F=(dT - dF) * 1e4))
    # 表4 重みの安定性: 分割seed 0〜3 それぞれで決めた重み、test で決めた重み
    ids = [ix["hgb"], ix["ridge_log"], ix["knn"]]
    for s in range(4):
        t4.append(dict(world=w, source=f"OOF split{s}", **dict(zip(["hgb", "ridge_log", "knn"], best_w(y, oof[s, ids], G3)))))
    t4.append(dict(world=w, source="test(正解)", **dict(zip(["hgb", "ridge_log", "knn"], best_w(yt, tp[1:, ids].mean(0), G3)))))
    # 表5 重みを決める行と測る行を分ける: 行の半分で重みを決め、残り半分で測る（OOF split0、20回）
    rng = np.random.default_rng(w)
    gaps = []
    for _ in range(20):
        perm = rng.permutation(len(y)); A, B = perm[: len(y) // 2], perm[len(y) // 2:]
        P = oof[0, ids]
        wv = best_w(y[A], P[:, A], G3)
        gA = rmsle(y[A], np.tensordot(wv, P[:, A], axes=1)) - rmsle(y[A], P[0, A])
        gB = rmsle(y[B], np.tensordot(wv, P[:, B], axes=1)) - rmsle(y[B], P[0, B])
        gaps.append((gA, gB))
    gaps = np.array(gaps)
    t5.append(dict(world=w, 決めた半分=gaps[:, 0].mean() * 1e4, 測った半分=gaps[:, 1].mean() * 1e4,
                   楽観=(gaps[:, 1] - gaps[:, 0]).mean() * 1e4, 楽観SE=(gaps[:, 1] - gaps[:, 0]).std(ddof=1) / np.sqrt(20) * 1e4))

T1 = pd.DataFrame(t1).pivot(index="world", columns="model", values=["oof_s0", "test"])
print("=== 表1 単独モデルの RMSLE（oof_s0 = 分割seed0のOOF、test = 6万行）==="); print(T1.round(5).to_string())
print("\n=== 表2 log特徴量は線形モデルに効くか（差 ×1e-4、負 = 前者が良い）==="); print(pd.DataFrame(t2).round(2).to_string(index=False))
print("\n=== 表3 混ぜる: hgb 単独との差 ×1e-4（負 = 混ぜた方が良い）、opt = 真の差 − 見えた差 ==="); T3 = pd.DataFrame(t3); print(T3.round(2).to_string(index=False))
print("\n=== 表4 重みの安定性 ==="); print(pd.DataFrame(t4).round(2).to_string(index=False))
print("\n=== 表5 行を分けて重みを測る（hgb単独との差 ×1e-4、20回平均）==="); print(pd.DataFrame(t5).round(2).to_string(index=False))

T2, T4 = pd.DataFrame(t2), pd.DataFrame(t4)
assert (T2[T2.pair == "ridge_logDur - ridge_raw"].zT < -10).all(), "log_Dur 1列で Ridge は全ワールドで大きく良くなる"
assert (T1["test"]["ridge_log"] > T1["test"]["hgb"]).all(), "log 4列の Ridge でも HistGB には全ワールドで負ける"
assert (T4[T4.source.str.startswith("OOF")].ridge_log == 0).all(), "OOF で決めた重みで Ridge は 16/16 回すべて 0"
print("all asserts passed")

"""Day166 補助: (1) 3モデルの残差相関 (2) test 予測を 5モデル平均(1分割) と 15モデル平均(3分割) で作ったときの混ぜる効果"""
import glob
import numpy as np
M = ["hgb", "ridge_raw", "ridge_logDur", "ridge_log", "knn"]; ix = {m: i for i, m in enumerate(M)}
def sle(y, lp): return (np.log1p(np.clip(np.expm1(lp), 1, None)) - np.log1p(y)) ** 2
def rm(y, lp): return np.sqrt(sle(y, lp).mean())
print("world | 残差相関 hgb-ridge_log / hgb-knn / ridge_log-knn (OOF split0, log空間) | test RMSLE hgb: 1分割(5モデル平均) 3分割(15モデル平均) | 混ぜる(0.9/0/0.1)との差 ×1e-4: OOF split0, test 1分割, test 3分割")
for f in sorted(glob.glob("raw_w*.npz")):
    z = np.load(f); oof, tp, y, yt = z["oof"], z["tp"], z["ytr"], z["yt"]
    ly = np.log1p(y); lyt = np.log1p(yt)
    r = {m: oof[0, ix[m]] - ly for m in ["hgb", "ridge_log", "knn"]}
    c = lambda a, b: np.corrcoef(r[a], r[b])[0, 1]
    h1, h3 = rm(yt, tp[1, 0]), rm(yt, tp[1:, 0].mean(0))
    w = np.array([0.9, 0, 0.1]); ids = [0, 3, 4]
    bl = lambda P: np.tensordot(w, P[ids], axes=1)
    d_oof = rm(y, bl(oof[0])) - rm(y, oof[0, 0])
    d_t1 = rm(yt, bl(tp[1])) - rm(yt, tp[1, 0])
    d_t3 = rm(yt, bl(tp[1:].mean(0))) - rm(yt, tp[1:, 0].mean(0))
    print(f"{f[-5]} | {c('hgb','ridge_log'):.3f} / {c('hgb','knn'):.3f} / {c('ridge_log','knn'):.3f} | {h1:.5f} {h3:.5f} | {d_oof*1e4:+.2f} {d_t1*1e4:+.2f} {d_t3*1e4:+.2f}")

print("\n--- 1モデル(fold単独) と 5モデル平均 の test 差（分割seed1、混ぜる 0.9/0/0.1 − hgb、×1e-4）---")
for f in sorted(glob.glob("raw_w*.npz")):
    z = np.load(f); tpf, tp, yt = z["tpf"].astype(np.float64), z["tp"], z["yt"]
    w = np.array([0.9, 0, 0.1]); ids = [0, 3, 4]
    single = [rm(yt, np.tensordot(w, tpf[k][ids], axes=1)) - rm(yt, tpf[k, 0]) for k in range(5)]
    avg5 = rm(yt, np.tensordot(w, tp[1][ids], axes=1)) - rm(yt, tp[1, 0])
    h_single = np.mean([rm(yt, tpf[k, 0]) for k in range(5)]); h_avg = rm(yt, tp[1, 0])
    print(f"world {f[-5]}: 1モデル平均 {np.mean(single)*1e4:+.2f} (範囲 {min(single)*1e4:+.2f}〜{max(single)*1e4:+.2f}) | 5モデル平均 {avg5*1e4:+.2f} | hgb 単独 test RMSLE 1モデル {h_single:.5f} 5モデル平均 {h_avg:.5f}")

print("\n--- 固定重み 0.9/0/0.1 の本当の効果（test、分割seed1〜3 平均、hgb 単独との差）---")
for f in sorted(glob.glob("raw_w*.npz")):
    z = np.load(f); tp, yt = z["tp"], z["yt"]
    w = np.array([0.9, 0, 0.1]); ids = [0, 3, 4]
    ea, eb = sle(yt, np.tensordot(w, tp[1:][:, ids].mean(0), axes=1)), sle(yt, tp[1:, 0].mean(0))
    d = np.sqrt(ea.mean()) - np.sqrt(eb.mean()); se = (ea - eb).std(ddof=1) / np.sqrt(len(yt)) / (2 * np.sqrt(eb.mean()))
    print(f"world {f[-5]}: dT {d*1e4:+.2f}e-4  zT {d/se:+.2f}")