phase3_folds.py

"""PHASE 3 (refit) — stability selection run INSIDE each training fold only.

Phase 3 as originally run saw all 1567 lots, which is what Phase 4b measured as
+0.146 AUC of selection leakage. This reruns the whole Phase 3 apparatus — frequency
thresholds, penalty sensitivity, permutation null calibration — separately on each
chronological training prefix.

Stricter than phase4.reselect() in one way that matters: the POOL (which columns are
constant / too sparse to use) is recomputed from the fold's rows only. The earlier
version decided that globally, which quietly used test-period rows.
"""
import numpy as np, pandas as pd
from sklearn.linear_model import LogisticRegression
from secom_common import load

NRES, PI, NPERM = 100, .60, 20
CG = {"C=0.03": .03, "C=0.05": .05}
CUTS = (.5, .6, .7, .8)
FULL = ["S59", "S103", "S129", "S21", "S348", "S64"]   # full-record Phase 3 result

X, y, ts = load()
o = np.argsort(ts.values)
X, y, ts = X.iloc[o].reset_index(drop=True), y.iloc[o].reset_index(drop=True), ts.iloc[o].reset_index(drop=True)
yv = y.values
n = len(X)


def fold_pool(tr):
    """Constant / high-missing decided on THIS fold's rows only."""
    sub = X.iloc[tr]
    nun, mr = sub.nunique(dropna=True), sub.isna().mean()
    return list(sub.columns[(nun > 1) & (mr <= .5)])


def stability(V, yy, C, nres, r):
    p = V.shape[1]
    sel, qs = np.zeros(p), []
    for _ in range(nres):
        idx = np.concatenate([r.choice(np.nonzero(yy == k)[0], int(.5*(yy == k).sum()),
                                       replace=False) for k in (0, 1)])
        A = V[idx]
        med = np.nanmedian(A, 0); med = np.where(np.isnan(med), 0., med)
        A = np.where(np.isnan(A), med, A)
        s = A.std(0); s[s == 0] = 1
        A = (A - A.mean(0)) / s
        c = LogisticRegression(penalty="l1", C=C, solver="liblinear",
                               class_weight="balanced", max_iter=2000).fit(A, yy[idx]).coef_[0] != 0
        sel += c; qs.append(c.sum())
    return sel/nres, float(np.mean(qs))


print("=" * 86)
print("PHASE 3 REFIT — STABILITY SELECTION INSIDE EACH TRAINING FOLD")
print("=" * 86)

freqs, chosen = {}, {}
for f in CUTS:
    cut = int(f*n); tr = np.arange(cut)
    pool = fold_pool(tr)
    V, ytr = X.iloc[tr][pool].values, yv[tr]
    print(f"\n{'-'*86}\nFOLD: first {f:.0%} = {cut} lots, {int(ytr.sum())} failures "
          f"({ytr.mean():.2%}), up to {ts[cut-1]}")
    print(f"  pool recomputed on this fold: {len(pool)} signals "
          f"({590-len(pool)} constant or >50% missing WITHIN the fold)")

    F = pd.DataFrame(index=pool)
    qb = {}
    for name, C in CG.items():
        F[name], qb[name] = stability(V, ytr, C, NRES, np.random.default_rng(7))
    F["freq"] = F["C=0.03"]
    freqs[f] = F
    sel = list(F.index[F.freq > PI])
    chosen[f] = sel

    print(f"  q per fit: C=0.03 -> {qb['C=0.03']:.1f}   C=0.05 -> {qb['C=0.05']:.1f}")
    print(f"  frequency thresholds:  " + "  ".join(
        f">{t:.0%}: {int((F.freq > t).sum())}" for t in (.9, .8, .7, .6, .5)))
    print(f"  never selected in {NRES} fits: {int((F.freq == 0).sum())}/{len(pool)}")

    # null calibration inside this fold
    rp = np.random.default_rng(11)
    nullF = np.array([stability(V, rp.permutation(ytr), CG["C=0.03"], NRES, rp)[0]
                      for _ in range(NPERM)])
    print(f"  null calibration ({NPERM} perms x {NRES} fits): "
          f"max freq reached by noise = {nullF.max():.2f}, "
          f"signals >60% under null: mean {(nullF>PI).sum(1).mean():.2f}, max {(nullF>PI).sum(1).max()}")
    obs = len(sel); nul = (nullF > PI).sum(1).mean()
    print(f"  observed >60%: {obs} signals -> empirical FDR at 60% = "
          f"{(nul/obs if obs else float('nan')):.2f}")

    a, b = [set(F.index[F[k] > PI]) for k in CG]
    print(f"  penalty sensitivity: C=0.03 -> {len(a)}, C=0.05 -> {len(b)}, both -> {len(a&b)}")
    print(f"  SELECTED >60%: {sel}")
    print(f"    with frequencies: " + ", ".join(f"{s}={F.freq[s]:.2f}" for s in sel))
    above_noise = [s for s in sel if F.freq[s] > nullF.max()]
    print(f"    above this fold's noise ceiling ({nullF.max():.2f}): {above_noise}")

# --- cross-fold summary --------------------------------------------------
print(f"\n{'='*86}\nCROSS-FOLD SUMMARY\n{'='*86}")
allsig = sorted({s for v in chosen.values() for s in v} | set(FULL))
T = pd.DataFrame(index=allsig)
for f in CUTS:
    T[f"{f:.0%}"] = [freqs[f].freq.get(s, np.nan) for s in allsig]
T["folds>60%"] = [sum(s in chosen[f] for f in CUTS) for s in allsig]
T["full_record"] = [s in FULL for s in allsig]
T = T.sort_values(["folds>60%"] + [f"{f:.0%}" for f in CUTS], ascending=False)
print("\nselection frequency per fold (NaN = not in that fold's pool):")
print(T.to_string(float_format=lambda v: f"{v:.2f}"))

print(f"\n  signals selected in ALL {len(CUTS)} folds : "
      f"{[s for s in allsig if T.loc[s,'folds>60%']==len(CUTS)]}")
print(f"  in 3+ folds                        : {[s for s in allsig if T.loc[s,'folds>60%']>=3]}")
print(f"  in 2+ folds                        : {[s for s in allsig if T.loc[s,'folds>60%']>=2]}")
print(f"  in exactly 1 fold                  : {[s for s in allsig if T.loc[s,'folds>60%']==1]}")

print(f"\n  FULL-RECORD Phase 3 result was: {FULL}")
for s in FULL:
    k = int(T.loc[s, "folds>60%"])
    fr = " ".join(f"{T.loc[s, f'{f:.0%}']:.2f}" if not np.isnan(T.loc[s, f"{f:.0%}"]) else "  - "
                  for f in CUTS)
    print(f"    {s:<6} folds>60%: {k}/4   freqs: {fr}   "
          f"{'reproduces' if k >= 3 else 'partial' if k >= 1 else 'DOES NOT REPRODUCE'}")

T.to_csv("phase3_folds.csv")
print("\n  -> phase3_folds.csv written")