Update baseline evaluation script: add multi-seed, 80/20 train/calib split, freeze tau, and compute exact matched coverage

TasfinMahmud · TasfinMahmud · commit edae1bcc2559 · 2026-06-30T19:32:45.000+06:00
diff --git a/evaluate_baselines.py b/evaluate_baselines.py
@@ -1,24 +1,40 @@
 import sys
 import traceback
+import os
+import numpy as np
+import scipy.stats as stats
 
 try:
+    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
     import os
+    import site
+    import sys
+
+    # Inject pip-installed NVIDIA DLLs into Windows PATH so TF 2.10 can use the GPU without a system CUDA install
+    try:
+        site_packages = site.getsitepackages()[1] if len(site.getsitepackages()) > 1 else site.getsitepackages()[0]
+        nvidia_base = os.path.join(site_packages, "nvidia")
+        dll_paths = [
+            os.path.join(nvidia_base, "cudnn", "bin"),
+            os.path.join(nvidia_base, "cublas", "bin"),
+            os.path.join(nvidia_base, "cuda_runtime", "bin"),
+        ]
+        os.environ["PATH"] = os.pathsep.join(dll_paths) + os.pathsep + os.environ.get("PATH", "")
+    except Exception:
+        pass
 
-    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
-    import numpy as np
-    import scipy.stats as stats
     import tensorflow as tf
 
     from cnsd import Dataset
     from cnsd.diagnosis.system import CNSD
-    from cnsd.perception.cnn import _train_cnn  # <--- WE IMPORT JUST THE CNN TRAINER
+    from cnsd.perception.cnn import _train_cnn
     from cnsd.physics import PhysicsConfig
     from validate_pu import load_pu_domain_split
 
     print('Loading Authentic PU dataset (Cross-Domain RPM Split)...')
-    (X_train, y_train, cond_train), (X_target, y_target, cond_target) = load_pu_domain_split()
+    (X_train_full, y_train_full, cond_train_full), (X_target, y_target, cond_target) = load_pu_domain_split()
 
-    unique_rpm = set(cond_train).union(set(cond_target))
+    unique_rpm = set(cond_train_full).union(set(cond_target))
     rpm_map = {float(r): float(r) for r in unique_rpm}
     pu_physics = PhysicsConfig(
         bearing={'n_balls': 8, 'd_ball': 6.75, 'd_pitch': 28.5, 'contact_angle': 0.0},
@@ -32,96 +48,22 @@
         2: ('Inner Race', 'High'),
     }
 
-    # Shuffle train data exactly like validate_pu.py to ensure parity
-    train_indices = np.arange(len(y_train))
-    np.random.seed(42)
-    np.random.shuffle(train_indices)
-    X_train, y_train, cond_train = (
-        X_train[train_indices],
-        y_train[train_indices],
-        cond_train[train_indices],
-    )
-
-    train_ds = Dataset.from_arrays(
-        X_train,
-        y_train,
-        cond_train,
-        fs=64000,
-        physics=pu_physics,
-        taxonomy=pu_taxonomy,
-        name='PU_Train',
-    )
-    test_ds = Dataset.from_arrays(
-        X_target,
-        y_target,
-        cond_target,
-        fs=64000,
-        physics=pu_physics,
-        taxonomy=pu_taxonomy,
-        name='PU_Test',
-    )
-
-    # Train Primary Model
-    print('\n[A] Training Primary CNSD Model (Seed=42)...')
-    np.random.seed(42)
-    tf.random.set_seed(42)
-    model = CNSD()
-    model.fit(train_ds, epochs=20)
-
-    # Train Ensemble Models
-    # FIX: We only need the CNN part for the ensemble, not the full Physics+DoWhy SCM pipeline!
-    # Training the full CNSD() 4 times was causing the joblib multiprocessing memory leak.
-    ens = []
-    nc = int(train_ds.y.max()) + 1
-    for s in [100, 101, 102]:
-        print(f'\n[B] Training Ensemble Member CNN only (Seed={s})...')
-        np.random.seed(s)
-        tf.random.set_seed(s)
-        # Train just the Keras CNN directly, skipping the heavy memory SCM
-        m_cnn = _train_cnn(train_ds.X, train_ds.y, num_classes=nc, epochs=20, seed=s)
-        ens.append(m_cnn)
-
-    # ---- shared: test signals, labels, condition ----
-    sig = np.stack([test_ds.X[i].reshape(-1) for i in range(len(test_ds.X))]).astype(np.float32)
-    yte = test_ds.y
-    cond_te = test_ds.cond
-    Xin = sig[..., None]  # (n, L, 1) for the conv nets
-
-    def gap_by_score(score, correct, thr):
-        hi = score >= thr
-        lo = ~hi
-        ah = correct[hi].mean() if hi.any() else float('nan')
-        al = correct[lo].mean() if lo.any() else float('nan')
-        return ah - al, ah, al, int(hi.sum()), int(lo.sum())
-
-    print('=' * 60, '\n1. PHYSICS gap across tau\n', '=' * 60)
-    probs = model.cnn.predict(Xin, batch_size=128, verbose=0)
-    pred = probs.argmax(1)
-    correct = pred == yte
-    print(f'baseline CNN acc on test: {correct.mean():.3f}')
-    for tau in [1.0, 2.0, 3.0]:
-        model.symbolic.tau = float(tau)
-        verds = np.array(
-            [
-                model.symbolic.diagnose(sig[i], pred[i], cond_te[i])['verdict']
-                for i in range(len(sig))
-            ]
-        )
-        conf = verds == 'CONFIRMED'
-        cnfl = verds == 'CONFLICT'
-        ca = correct[conf].mean() if conf.any() else float('nan')
-        fa = correct[cnfl].mean() if cnfl.any() else float('nan')
-        print(
-            f'  tau={tau:.1f}  GAP={ca - fa:+.3f}  CONF={ca:.3f}(n={conf.sum()})  CNFL={fa:.3f}(n={cnfl.sum()})'
-        )
-
-    print('\n', '=' * 60, '\n2. SOFTMAX confidence\n', '=' * 60)
-    softmax_conf = probs.max(1)
-    for thr in [0.5, 0.7, 0.9, 0.95, 0.99]:
-        g, ah, al, nh, nl = gap_by_score(softmax_conf, correct, thr)
-        print(f'  conf>={thr:.2f}  GAP={g:+.3f}  HI={ah:.3f}(n={nh})  LO={al:.3f}(n={nl})')
-
-    print('\n', '=' * 60, '\n3. MC-DROPOUT (30 passes)\n', '=' * 60)
+    def get_matched_coverage_gap(score, correct, target_n):
+        if target_n == 0:
+            return float('nan')
+        # Score is higher for MORE confident
+        # Sort descending by score
+        sorted_indices = np.argsort(score)[::-1]
+        hi_indices = sorted_indices[:target_n]
+        lo_indices = sorted_indices[target_n:]
+        
+        hi_mask = np.zeros(len(score), dtype=bool)
+        hi_mask[hi_indices] = True
+        lo_mask = ~hi_mask
+        
+        ah = correct[hi_mask].mean() if hi_mask.any() else float('nan')
+        al = correct[lo_mask].mean() if lo_mask.any() else float('nan')
+        return ah - al
 
     class AlwaysDropout(tf.keras.layers.Dropout):
         def call(self, inputs, training=None):
@@ -132,69 +74,175 @@ def clone_for_mc(layer):
             return AlwaysDropout(layer.rate)
         return layer.__class__.from_config(layer.get_config())
 
-    mc_model = tf.keras.models.clone_model(model.cnn, clone_function=clone_for_mc)
-    mc_model.set_weights(model.cnn.get_weights())
-
-    T = 30
-    mc = []
-    for _ in range(T):
-        mc.append(mc_model.predict(Xin, batch_size=128, verbose=0))
-    mc = np.stack(mc)  # (T,n,classes)
-    mc_mean = mc.mean(0)
-    mc_pred = mc_mean.argmax(1)
-    mc_correct = mc_pred == yte
-    eps = 1e-12
-    certainty = (mc_mean * np.log(mc_mean + eps)).sum(1)  # = -entropy; higher = more certain
-    print(f'MC mean acc: {mc_correct.mean():.3f}')
-    for q in [50, 30, 20, 10]:  # keep top (100-q)% most certain as 'reliable'
-        thr = np.percentile(certainty, q)
-        g, ah, al, nh, nl = gap_by_score(certainty, mc_correct, thr)
-        print(f'  top {100 - q}% certain  GAP={g:+.3f}  HI={ah:.3f}(n={nh})  LO={al:.3f}(n={nl})')
-
-    print('\n', '=' * 60, '\n4. ENSEMBLE disagreement\n', '=' * 60)
-    ens_preds = np.stack([m.predict(Xin, batch_size=128, verbose=0).argmax(1) for m in ens])
-    vote = stats.mode(ens_preds, axis=0, keepdims=False).mode
-    unanimous = (ens_preds == vote).all(0)
-    ens_correct = vote == yte
-    print(f'ensemble vote acc: {ens_correct.mean():.3f} | unanimous: {unanimous.mean():.1%}')
-    g, ah, al, nh, nl = gap_by_score(
-        unanimous.astype(float), ens_correct, 1.0
-    )  # unanimous=reliable
-    print(f'  unanimous=reliable  GAP={g:+.3f}  HI={ah:.3f}(n={nh})  LO={al:.3f}(n={nl})')
-
-    print('\n', '=' * 60, '\n5. NOISE TEST: physics CONFLICT on unanimous-but-wrong\n', '=' * 60)
-    rng = np.random.RandomState(42)
-    sig_power = (sig**2).mean()
-    model.symbolic.tau = 3.0
-    print(
-        f'{"noise":>6} | {"ens_acc":>7} | {"unanim-wrong":>12} | {"phys catches":>12} | {"base rate":>9}'
+    seeds = [42, 43, 44, 45, 46]
+    
+    results = {
+        'phys_gap': [],
+        'soft_gap': [],
+        'mc_gap': [],
+        'ens_gap': [],
+        'noise_catch': {db: [] for db in [np.inf, 20, 10, 5, 0]}
+    }
+    
+    # Test data processing
+    test_ds = Dataset.from_arrays(
+        X_target, y_target, cond_target, fs=64000, physics=pu_physics, taxonomy=pu_taxonomy, name='PU_Test'
     )
-    print('-' * 60)
+    sig_te = np.stack([test_ds.X[i].reshape(-1) for i in range(len(test_ds.X))]).astype(np.float32)
+    yte = test_ds.y
+    cond_te = test_ds.cond
+    Xin_te = sig_te[..., None]
+    
+    for seed in seeds:
+        print(f"\n{'='*80}\n=== RUNNING SEED {seed} ===\n{'='*80}")
+        tf.keras.backend.clear_session()
+        np.random.seed(seed)
+        tf.random.set_seed(seed)
+        
+        # 1. 80/20 Split for Calibration
+        indices = np.arange(len(y_train_full))
+        np.random.shuffle(indices)
+        
+        split_idx = int(0.8 * len(indices))
+        train_idx = indices[:split_idx]
+        calib_idx = indices[split_idx:]
+        
+        X_tr, y_tr, cond_tr = X_train_full[train_idx], y_train_full[train_idx], cond_train_full[train_idx]
+        X_ca, y_ca, cond_ca = X_train_full[calib_idx], y_train_full[calib_idx], cond_train_full[calib_idx]
+        
+        train_ds = Dataset.from_arrays(X_tr, y_tr, cond_tr, fs=64000, physics=pu_physics, taxonomy=pu_taxonomy, name='PU_Train')
+        calib_ds = Dataset.from_arrays(X_ca, y_ca, cond_ca, fs=64000, physics=pu_physics, taxonomy=pu_taxonomy, name='PU_Calib')
+        
+        # 2. Train Primary Model (Bypass SCM to prevent multiprocess deadlocks)
+        model = CNSD()
+        nc = int(train_ds.y.max()) + 1
+        model.cnn = _train_cnn(train_ds.X, train_ds.y, num_classes=nc, epochs=20, seed=seed)
+        model.symbolic = model._build_symbolic(train_ds)
+        model._fitted = True
+        
+        # 3. Train Ensemble Models
+        ens = []
+        nc = int(train_ds.y.max()) + 1
+        for s in [seed*10, seed*10+1, seed*10+2]:
+            np.random.seed(s)
+            tf.random.set_seed(s)
+            m_cnn = _train_cnn(train_ds.X, train_ds.y, num_classes=nc, epochs=20, seed=s)
+            ens.append(m_cnn)
+            
+        # 4. Calibrate Tau
+        sig_ca = np.stack([calib_ds.X[i].reshape(-1) for i in range(len(calib_ds.X))]).astype(np.float32)
+        Xin_ca = sig_ca[..., None]
+        probs_ca = model.cnn.predict(Xin_ca, batch_size=128, verbose=0)
+        pred_ca = probs_ca.argmax(1)
+        correct_ca = (pred_ca == calib_ds.y)
+        
+        best_tau = 1.0
+        best_gap = -100.0
+        for tau in [1.0, 1.5, 2.0, 2.5, 3.0]:
+            model.symbolic.tau = float(tau)
+            verds = np.array([model.symbolic.diagnose(sig_ca[i], pred_ca[i], calib_ds.cond[i])['verdict'] for i in range(len(sig_ca))])
+            conf = (verds == 'CONFIRMED')
+            cnfl = (verds == 'CONFLICT')
+            ca = correct_ca[conf].mean() if conf.any() else 0.0
+            fa = correct_ca[cnfl].mean() if cnfl.any() else 1.0
+            gap = ca - fa
+            if gap > best_gap:
+                best_gap = gap
+                best_tau = tau
+                
+        print(f"Calibrated best tau = {best_tau} (Calib GAP = {best_gap:+.3f})")
+        model.symbolic.tau = float(best_tau)
+        
+        # 5. Evaluate on Test Set
+        probs_te = model.cnn.predict(Xin_te, batch_size=128, verbose=0)
+        pred_te = probs_te.argmax(1)
+        correct_te = (pred_te == yte)
+        
+        # Physics evaluation
+        verds = np.array([model.symbolic.diagnose(sig_te[i], pred_te[i], cond_te[i])['verdict'] for i in range(len(sig_te))])
+        conf = (verds == 'CONFIRMED')
+        cnfl = (verds == 'CONFLICT')
+        ca = correct_te[conf].mean() if conf.any() else float('nan')
+        fa = correct_te[cnfl].mean() if cnfl.any() else float('nan')
+        phys_gap = ca - fa
+        target_n = int(conf.sum())
+        results['phys_gap'].append(phys_gap)
+        print(f"Physics GAP={phys_gap:+.3f} (Coverage N={target_n})")
+        
+        # Softmax evaluation at matched coverage
+        softmax_score = probs_te.max(1)
+        soft_gap = get_matched_coverage_gap(softmax_score, correct_te, target_n)
+        results['soft_gap'].append(soft_gap)
+        
+        # MC-Dropout at matched coverage
+        mc_model = tf.keras.models.clone_model(model.cnn, clone_function=clone_for_mc)
+        mc_model.set_weights(model.cnn.get_weights())
+        T = 30
+        mc_preds = []
+        for _ in range(T):
+            mc_preds.append(mc_model.predict(Xin_te, batch_size=128, verbose=0))
+        mc_preds = np.stack(mc_preds)
+        mc_mean = mc_preds.mean(0)
+        mc_pred_class = mc_mean.argmax(1)
+        mc_correct = (mc_pred_class == yte)
+        eps = 1e-12
+        mc_score = (mc_mean * np.log(mc_mean + eps)).sum(1) # Certainty (negative entropy)
+        mc_gap = get_matched_coverage_gap(mc_score, mc_correct, target_n)
+        results['mc_gap'].append(mc_gap)
+        
+        # Ensemble at matched coverage
+        ens_preds_probs = np.stack([m.predict(Xin_te, batch_size=128, verbose=0) for m in ens]) # (3, n, c)
+        ens_mean = ens_preds_probs.mean(0) # (n, c)
+        ens_pred_class = ens_mean.argmax(1)
+        ens_correct = (ens_pred_class == yte)
+        # Score = negative entropy of ensemble mean
+        ens_score = (ens_mean * np.log(ens_mean + eps)).sum(1)
+        ens_gap = get_matched_coverage_gap(ens_score, ens_correct, target_n)
+        results['ens_gap'].append(ens_gap)
+        
+        print(f"Matched Coverage GAPs -> Softmax:{soft_gap:+.3f} | MC-Drop:{mc_gap:+.3f} | Ens:{ens_gap:+.3f}")
+        
+        # 6. Noise Test
+        rng = np.random.RandomState(seed)
+        sig_power = (sig_te ** 2).mean()
+        for snr_db in [np.inf, 20, 10, 5, 0]:
+            if np.isinf(snr_db):
+                sig_n = sig_te
+            else:
+                npow = sig_power / (10 ** (snr_db / 10))
+                sig_n = sig_te + rng.randn(*sig_te.shape).astype(np.float32) * np.sqrt(npow)
+            
+            Xin_n = sig_n[..., None]
+            # Use ensemble mode vote to define 'unanimous' exactly like Abhi's template
+            ep_class = np.stack([m.predict(Xin_n, batch_size=128, verbose=0).argmax(1) for m in ens])
+            v = stats.mode(ep_class, axis=0, keepdims=False).mode
+            unan = (ep_class == v).all(0)
+            ok = (v == yte)
+            
+            pv = np.array([model.symbolic.diagnose(sig_n[i], v[i], cond_te[i])['verdict'] for i in range(len(sig_n))])
+            pc = (pv == 'CONFLICT')
+            uw = unan & (~ok)
+            catch = pc[uw].mean() if uw.sum() > 0 else float('nan')
+            results['noise_catch'][snr_db].append(catch)
+            s = 'clean' if np.isinf(snr_db) else f'{snr_db}dB'
+            print(f"  Noise={s:>5} | catch_rate={catch:.3f}")
+            
+    print("\n" + "="*60 + "\nFINAL AGGREGATED RESULTS (5 Seeds)\n" + "="*60)
+    print(f"Physics GAP: {np.nanmean(results['phys_gap']):+.3f} ± {np.nanstd(results['phys_gap']):.3f}")
+    print(f"Softmax GAP: {np.nanmean(results['soft_gap']):+.3f} ± {np.nanstd(results['soft_gap']):.3f}")
+    print(f"MC-Drop GAP: {np.nanmean(results['mc_gap']):+.3f} ± {np.nanstd(results['mc_gap']):.3f}")
+    print(f"Ensemble GAP: {np.nanmean(results['ens_gap']):+.3f} ± {np.nanstd(results['ens_gap']):.3f}")
+    
+    print("\nNoise Test Catch Rate (Physics catches Ensemble's confident errors):")
     for snr_db in [np.inf, 20, 10, 5, 0]:
-        if np.isinf(snr_db):
-            sig_n = sig
-        else:
-            npow = sig_power / (10 ** (snr_db / 10))
-            sig_n = sig + rng.randn(*sig.shape).astype(np.float32) * np.sqrt(npow)
-        Xin_n = sig_n[..., None]
-        ep = np.stack([m.predict(Xin_n, batch_size=128, verbose=0).argmax(1) for m in ens])
-        v = stats.mode(ep, axis=0, keepdims=False).mode
-        unan = (ep == v).all(0)
-        ok = v == yte
-        pv = np.array(
-            [
-                model.symbolic.diagnose(sig_n[i], v[i], cond_te[i])['verdict']
-                for i in range(len(sig_n))
-            ]
-        )
-        pc = pv == 'CONFLICT'
-        uw = unan & (~ok)
-        catch = pc[uw].mean() if uw.sum() > 0 else float('nan')
         s = 'clean' if np.isinf(snr_db) else f'{snr_db}dB'
-        print(f'{s:>6} | {ok.mean():>7.3f} | {uw.sum():>12} | {catch:>12.3f} | {pc.mean():>9.3f}')
-
-    print('================ DONE ================')
-
+        vals = [v for v in results['noise_catch'][snr_db] if not np.isnan(v)]
+        m = np.nanmean(vals) if len(vals) > 0 else float('nan')
+        std = np.nanstd(vals) if len(vals) > 0 else float('nan')
+        print(f"  {s:>6}: {m:.3f} ± {std:.3f}")
+        
+    print("================ DONE ================")
+    
 except Exception:
     with open('crash_traceback.txt', 'w') as f:
         traceback.print_exc(file=f)
