mkdir -p ~/.kaggle
python3 -c "import json;print(json.load(open('kaggle.json'))['key'],end='')" > ~/.kaggle/access_token
chmod 600 ~/.kaggle/access_token

%matplotlib inline
import time, warnings; warnings.filterwarnings('ignore')   # keep output clean
import numpy as np, pandas as pd                            # numerics + dataframes
import matplotlib.pyplot as plt                             # static plots (embed in HTML+PDF)
plt.rcParams['figure.dpi'] = 120                            # crisp figures
RANDOM_STATE = 0                                            # single seed used everywhere
np.random.seed(RANDOM_STATE)                                # reproducible sampling
NEG_WORD, POS_WORD = 'benign', 'attack'                      # class names (overridden by some loaders)

import os, glob
# CIRA-CIC-DoHBrw-2020 (MontazeriShatoori et al., 2020): flow-statistics of DNS-over-HTTPS traffic,
# labelled benign vs malicious (DoH tunneling via dns2tcp / DNSCat2 / Iodine). Self-contained.
os.environ.setdefault('KAGGLE_KEY', open(os.path.expanduser('~/.kaggle/access_token')).read().strip())
DEST = '/tmp/kg_doh2020'; os.makedirs(DEST, exist_ok=True)
if not glob.glob(DEST + '/**/*.csv', recursive=True):
    import kaggle; kaggle.api.authenticate()
    print('downloading CIRA-CIC-DoHBrw-2020 (one-time)...')
    kaggle.api.dataset_download_files('bcccdatasets/bccc-cira-cic-dohbrw-2020', path=DEST, unzip=True, quiet=True)
f = sorted(glob.glob(DEST + '/**/*.csv', recursive=True), key=os.path.getsize, reverse=True)[0]
df = pd.read_csv(f, low_memory=False); df.columns = [str(c).strip() for c in df.columns]
assert len(df) >= 400_000, f'floor not met: {len(df):,}'
LABEL = [c for c in df.columns if c.lower() == 'label'][0]
df['y'] = (df[LABEL].astype(str).str.strip().str.lower() != 'benign').astype(int)
df['family'] = df[LABEL].astype(str).str.strip()               # Benign / Malicious
# Drop the label + any IP/timestamp identifiers so the model uses only flow-timing statistics.
DROP = [LABEL, 'y', 'family'] + [c for c in df.columns if c.lower() in
        ('sourceip','destinationip','source_ip','destination_ip','timestamp','flowid','flow_id')]
feat = [c for c in df.columns if c not in DROP]
from sklearn.preprocessing import LabelEncoder
X = df[feat].copy()
idlike = [c for c in X.select_dtypes(include='object').columns if X[c].nunique() > 0.5*len(X)]
X = X.drop(columns=idlike)                                     # drop id/timestamp-like leaky columns
for c in X.select_dtypes(include='object').columns:
    X[c] = LabelEncoder().fit_transform(X[c].astype(str))
X = X.apply(pd.to_numeric, errors='coerce').replace([np.inf,-np.inf],np.nan).fillna(0.0)
X = X.clip(-1e15, 1e15); X = X.loc[:, X.nunique() > 1]         # float32-safe; drop constants
import re
_seen, _cols = {}, []
for _c in X.columns:                                           # unique LightGBM-safe names
    _c = re.sub(r'[^0-9A-Za-z_]+', '_', str(_c)).strip('_') or 'f'
    _seen[_c] = _seen.get(_c, -1) + 1
    _cols.append(_c if _seen[_c] == 0 else f'{_c}_{_seen[_c]}')
X.columns = _cols; feat = list(X.columns)
y = df['y'].to_numpy(); family = df['family'].to_numpy()
print(f'loaded {len(df):,} DoH flows x {len(feat)} timing features; malicious rate {y.mean():.4f}')

loaded 499,106 DoH flows x 28 timing features; malicious rate 0.5000

# --- EDA 1: class balance and the attack-family mix ---
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
df['y'].map({0:NEG_WORD,1:POS_WORD}).value_counts().plot.bar(               # counts per class
    ax=ax[0], color=['#2a9d8f','#e76f51']); ax[0].set_yscale('log')
ax[0].set_title(f'Class balance ({NEG_WORD} vs {POS_WORD})'); ax[0].set_ylabel('records (log)')
df.loc[df.y==1,'family'].value_counts().head(8).plot.barh(                  # top attack families
    ax=ax[1], color='#e76f51'); ax[1].invert_yaxis(); ax[1].set_title('Top attack families')
plt.tight_layout(); plt.show()

# --- EDA 2: feature correlation + a 2-D PCA projection ---
from sklearn.preprocessing import StandardScaler                 # scale before PCA
from sklearn.decomposition import PCA
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
topv = X[feat].var().sort_values().tail(12).index                # 12 highest-variance features
im = ax[0].imshow(X[topv].corr(), cmap='coolwarm', vmin=-1, vmax=1)  # correlation heatmap
ax[0].set_xticks(range(len(topv))); ax[0].set_xticklabels(topv, rotation=90, fontsize=7)
ax[0].set_yticks(range(len(topv))); ax[0].set_yticklabels(topv, fontsize=7)
ax[0].set_title('Feature correlation (top-variance)'); fig.colorbar(im, ax=ax[0], shrink=0.7)
samp = X.sample(min(5000, len(X)), random_state=RANDOM_STATE)     # subsample for a fast PCA
pc = PCA(n_components=2).fit_transform(StandardScaler().fit_transform(samp))
ys = y[samp.index]                                               # aligned labels for coloring
for lab,c in [(0,'#2a9d8f'),(1,'#e76f51')]:
    ax[1].scatter(pc[ys==lab,0], pc[ys==lab,1], s=4, alpha=0.4, color=c,
                  label={0:NEG_WORD,1:POS_WORD}[lab])
ax[1].set_title('PCA projection (2 components)'); ax[1].legend(); ax[1].set_xlabel('PC1'); ax[1].set_ylabel('PC2')
plt.tight_layout(); plt.show()

# --- Model comparison: four learners on the same held-out split ---
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, roc_auc_score
import xgboost as xgb, lightgbm as lgb

# Stratified split keeps the class ratio in both halves.
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=RANDOM_STATE, stratify=y)
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
N_MATERIALIZED = len(y)                                          # the full corpus we loaded (see printed count)
# HONEST DISCLOSURE: we do NOT train on all N. We fit on a STRATIFIED subsample (<=120k) because
# these learners saturate long before then on this data. Every headline below is a SUBSAMPLE
# number, not a full-corpus number — saying otherwise would be the fabrication this course forbids.
if len(Xtr) > 120_000:
    Xtr, _, ytr, _ = train_test_split(Xtr, ytr, train_size=120_000, random_state=RANDOM_STATE,
                                      stratify=ytr)               # genuinely stratified, not random
MAJORITY_BASELINE = max(np.mean(yte), 1 - np.mean(yte))          # accuracy of 'always predict majority'
print(f'materialized {N_MATERIALIZED:,} rows | trained on {len(Xtr):,} (stratified subsample) | '
      f'held-out {len(yte):,}')
print(f'MAJORITY-CLASS BASELINE accuracy = {MAJORITY_BASELINE:.4f}  '
      f'(any model must beat THIS, not 0.5, to be interesting)')

models = {                                                        # four standard, diverse learners
    'LogisticRegression': make_pipeline(StandardScaler(), LogisticRegression(max_iter=300)),  # scaled!
    'RandomForest': RandomForestClassifier(n_estimators=60, n_jobs=-1, random_state=RANDOM_STATE),
    'XGBoost': xgb.XGBClassifier(n_estimators=80, max_depth=6, tree_method='hist', n_jobs=-1,
                                 eval_metric='logloss', random_state=RANDOM_STATE),
    'LightGBM': lgb.LGBMClassifier(n_estimators=80, n_jobs=-1, verbose=-1, random_state=RANDOM_STATE),
}
rows, fitted = [], {}
for name, m in models.items():                                    # fit + score each model
    t = time.perf_counter(); m.fit(Xtr, ytr); fitted[name] = m
    p = m.predict_proba(Xte)[:, 1]                                # positive-class probability on held-out
    rows.append({'model': name, 'accuracy': round(accuracy_score(yte, (p>0.5).astype(int)), 6),
                 'roc_auc': round(roc_auc_score(yte, p), 6),      # 6 dp: a 1.000000 is a red flag, not a win
                 'train_s': round(time.perf_counter()-t, 1)})
rows.append({'model': 'MajorityBaseline', 'accuracy': round(MAJORITY_BASELINE, 4),
             'roc_auc': 0.5, 'train_s': 0.0})            # show the baseline IN the ranking table
comparison = pd.DataFrame(rows).sort_values('roc_auc', ascending=False).reset_index(drop=True)
_ranked = comparison[comparison.model != 'MajorityBaseline']
best_name = _ranked.iloc[0]['model']; best = fitted[best_name]  # winner by ROC-AUC (excl. baseline)
print('best model:', best_name); comparison

materialized 499,106 rows | trained on 120,000 (stratified subsample) | held-out 124,777
MAJORITY-CLASS BASELINE accuracy = 0.5000  (any model must beat THIS, not 0.5, to be interesting)

best model: XGBoost

# --- Results for the best model: confusion, ROC, PR, importances, per-family recall ---
from sklearn.metrics import confusion_matrix, roc_curve, precision_recall_curve, recall_score
pb = best.predict_proba(Xte)[:, 1]; pred = (pb > 0.5).astype(int)
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
# (1) confusion matrix
cm = confusion_matrix(yte, pred); ax[0].imshow(cm, cmap='Blues')
ax[0].set_title(f'{best_name}: confusion'); ax[0].set_xticks([0,1]); ax[0].set_yticks([0,1])
ax[0].set_xticklabels([NEG_WORD,POS_WORD]); ax[0].set_yticklabels([NEG_WORD,POS_WORD])
for (i,j),v in np.ndenumerate(cm): ax[0].text(j,i,f'{v:,}',ha='center',va='center')
# (2) ROC and PR curves
fpr,tpr,_ = roc_curve(yte, pb); prec,rec,_ = precision_recall_curve(yte, pb)
ax[1].plot(fpr,tpr,color='#264653'); ax[1].plot([0,1],[0,1],'--',c='grey')
ax[1].set_title(f'ROC (AUC={roc_auc_score(yte,pb):.4f})'); ax[1].set_xlabel('FPR'); ax[1].set_ylabel('TPR')
ax[2].plot(rec,prec,color='#e76f51'); ax[2].set_title('Precision-Recall'); ax[2].set_xlabel('recall'); ax[2].set_ylabel('precision')
plt.tight_layout(); plt.show()

# (3) feature importances + (4) per-attack-family recall
fig, ax = plt.subplots(1, 2, figsize=(13, 5))
imp, names = None, feat                                           # importances, robust to the scaled-LR pipeline
if hasattr(best, 'feature_importances_'):                          # tree models
    imp = best.feature_importances_; names = list(getattr(best, 'feature_names_in_', feat))[:len(imp)]
elif hasattr(best, 'named_steps') and 'logisticregression' in getattr(best, 'named_steps', {}):
    imp = np.abs(best.named_steps['logisticregression'].coef_[0]); names = feat  # LR pipeline
elif hasattr(best, 'coef_'):
    imp = np.abs(best.coef_[0]); names = feat
if imp is not None:
    pd.Series(imp, index=names[:len(imp)]).sort_values().tail(12).plot.barh(ax=ax[0], color='#264653')
ax[0].set_title(f'{best_name}: top importances / |coef|')
# Per-family recall, WORST-first so rare, hard classes are visible, not just the dominant floods.
fam_te = df.loc[Xte.index, 'family']
fr = {}
for fam, cnt in fam_te[yte==1].value_counts().items():
    if cnt < 5: continue                                          # need a few positives for a meaningful recall
    mask = (fam_te==fam).to_numpy(); fr[fam] = recall_score(yte[mask], pred[mask], zero_division=0)
srt = pd.Series(fr).sort_values()
show = pd.concat([srt.head(9), srt.tail(3)]) if len(srt) > 12 else srt   # worst 9 + best 3
show = show[~show.index.duplicated()]
show.plot.barh(ax=ax[1], color=['#e76f51' if v < 0.5 else '#2a9d8f' for v in show]); ax[1].set_xlim(0,1)
ax[1].set_title('Per-family recall (worst first; red < 0.5)')
plt.tight_layout(); plt.show()
# Operational numbers, not just figures: false-positive rate and the worst per-family recalls.
tn, fp = int(cm[0,0]), int(cm[0,1])
fpr_op = fp/(fp+tn) if (fp+tn) > 0 else float('nan')             # benign wrongly flagged @0.5
print(f'operational FALSE-POSITIVE RATE @0.5 = {fpr_op:.4f}  ({fp:,} benign flagged of {fp+tn:,})')
print('worst per-family recalls:', {k: round(v, 3) for k, v in srt.head(6).items()})

operational FALSE-POSITIVE RATE @0.5 = 0.0002  (14 benign flagged of 62,389)
worst per-family recalls: {'Malicious': 1.0}

# --- Validity audit: is the score real detection, or a data shortcut? ---
from sklearn.metrics import roc_auc_score
samp = X.sample(min(60_000, len(X)), random_state=1); ysamp = y[samp.index]
aucs = {}
for c in feat:                                                    # AUC of EACH feature alone
    col = samp[c].to_numpy(float)
    if col.std()==0: continue
    a = roc_auc_score(ysamp, col); aucs[c] = max(a, 1-a)          # direction-agnostic
best_auc = max(aucs.values()); best_col = max(aucs, key=aucs.get)
dup_rate = 1 - X.drop_duplicates().shape[0]/len(X)               # exact-duplicate feature rows (whole set)
# The statistic that actually inflates a held-out score is TRAIN/TEST CONTAMINATION: how many test
# rows are exact duplicates of a training row. Measure it directly on the split used above.
_trkeys = set(map(tuple, np.round(Xtr.to_numpy(), 6)))
_te = np.round(Xte.to_numpy(), 6)[:50_000]
contam = float(np.mean([tuple(r) in _trkeys for r in _te]))      # fraction of test rows seen in train
# Trust grade reflects BOTH failure modes and takes the WORSE of the two: a near-perfect single
# feature (shortcut) OR heavy train/test contamination each independently invalidate the headline.
_ga = 'F' if best_auc>=0.999 else 'D' if best_auc>=0.99 else 'C' if best_auc>=0.95 else 'B' if best_auc>=0.85 else 'A'
_gc = 'F' if contam>=0.5 else 'D' if contam>=0.3 else 'C' if contam>=0.15 else 'B' if contam>=0.05 else 'A'
grade = max(_ga, _gc)                                            # 'max' letter = worse grade (A best, F worst)
print(f'best single-feature AUC = {best_auc:.4f}  (feature: {best_col})')
print(f'   note: a near-1.0 single-feature AUC means this feature is *near-sufficient* (a shortcut),\n'
      f'   which may be legitimate signal OR an artifact — it is NOT the same as target leakage.')
print(f'exact-duplicate row rate (whole corpus) = {dup_rate:.3f}')
print(f'TRAIN/TEST exact-row contamination       = {contam:.3f}  (single-feat grade {_ga}, contam grade {_gc})')
print(f'==> data trust grade: {grade}   (worse of the two; F = shortcut and/or heavy contamination)')
s = pd.Series(aucs).sort_values().tail(15)
fig, ax = plt.subplots(figsize=(8,5))
s.plot.barh(ax=ax, color=['#e76f51' if v>=0.99 else '#457b9d' for v in s]); ax.axvline(0.5,ls='--',c='grey')
ax.set_xlim(0.5,1.0); ax.set_title('Single-feature ROC-AUC (red = near-perfect shortcut)'); ax.set_xlabel('AUC alone')
plt.tight_layout(); plt.show()

best single-feature AUC = 0.7762  (feature: PacketLengthMean)
   note: a near-1.0 single-feature AUC means this feature is *near-sufficient* (a shortcut),
   which may be legitimate signal OR an artifact — it is NOT the same as target leakage.
exact-duplicate row rate (whole corpus) = 0.014
TRAIN/TEST exact-row contamination       = 0.014  (single-feat grade A, contam grade A)
==> data trust grade: A   (worse of the two; F = shortcut and/or heavy contamination)

# --- Ablation: SHOW the inflation empirically, don't just narrate it ---
from sklearn.base import clone
def _retrain_auc(Xa, ya):                                        # re-split, stratified-subsample, refit best family
    xtr, xte, ytr2, yte2 = train_test_split(Xa, ya, test_size=0.25, random_state=RANDOM_STATE, stratify=ya)
    if len(xtr) > 120_000:
        xtr, _, ytr2, _ = train_test_split(xtr, ytr2, train_size=120_000, random_state=RANDOM_STATE, stratify=ytr2)
    m = clone(best); m.fit(xtr, ytr2)
    return roc_auc_score(yte2, m.predict_proba(xte)[:, 1])
base_auc = roc_auc_score(yte, best.predict_proba(Xte)[:, 1])     # (0) the headline held-out AUC
Xdd = X.drop_duplicates(); ydd = y[Xdd.index]                    # (1) de-duplicated corpus
auc_dedup = _retrain_auc(Xdd, ydd)
auc_noshort = _retrain_auc(X.drop(columns=[best_col]), y) if best_col in X.columns else base_auc  # (2) drop shortcut
ablation = pd.DataFrame([
    {'setting': 'headline (as-is)',              'held_out_auc': round(base_auc, 6)},
    {'setting': f'de-duplicated ({1-len(Xdd)/len(X):.0%} rows removed)', 'held_out_auc': round(auc_dedup, 6)},
    {'setting': f'shortcut feature dropped ({best_col})', 'held_out_auc': round(auc_noshort, 6)},
])
print('Ablation — how much of the headline survives once each artifact is removed:')
ablation

Ablation — how much of the headline survives once each artifact is removed:

# --- Reproducibility & robustness ---
import sklearn
from sklearn.model_selection import StratifiedKFold, cross_val_score
print(f'seed={RANDOM_STATE} | numpy {np.__version__} | sklearn {sklearn.__version__} | '
      f'xgboost {xgb.__version__} | lightgbm {lgb.__version__}')
# 3-fold cross-validated ROC-AUC of the winning model (fresh clone, bounded subsample) -> mean +/- std.
from sklearn.base import clone
cvX, cvy = Xtr.iloc[:40_000], ytr[:40_000]
def _auc_scorer(est, Xv, yv):                                   # robust to xgboost's 2-col predict_proba
    p = est.predict_proba(Xv)
    p = p[:, 1] if getattr(p, 'ndim', 1) == 2 else p
    return roc_auc_score(yv, p)
try:
    cv = cross_val_score(clone(best), cvX, cvy,
                         cv=StratifiedKFold(3, shuffle=True, random_state=RANDOM_STATE),
                         scoring=_auc_scorer, error_score='raise')
    assert np.all(np.isfinite(cv)), 'non-finite CV folds'   # FAIL CLOSED: never narrate a NaN as evidence
    print(f'{best_name} 3-fold CV ROC-AUC = {cv.mean():.4f} +/- {cv.std():.4f}  '
          f'(mean +/- std across 3 stratified folds; a small std means a stable estimate on this split)')
except Exception as e:
    print(f'CV UNAVAILABLE ({type(e).__name__}: {str(e)[:60]}); rely on the single held-out AUC above — '
          f'we do NOT report a CV number we could not compute')

seed=0 | numpy 2.3.5 | sklearn 1.9.0 | xgboost 1.6.2 | lightgbm 4.7.0

XGBoost 3-fold CV ROC-AUC = 1.0000 +/- 0.0000  (mean +/- std across 3 stratified folds; a small std means a stable estimate on this split)

Reported approach	Known caveat
MontazeriShatoori et al. (2020) — time-series classifiers	few tunneling tools; tool-specific timing can be a shortcut
Encrypted-traffic side-channel detectors	sensitive to padding/regime shift

Property	Value
Source	Kaggle `bcccdatasets/bccc-cira-cic-dohbrw-2020`
Rows	499,106 DoH flows. The exact 50/50 balance is constructed: BCCC re-balanced the natively skewed CIRA-CIC-DoHBrw-2020 capture, which is why the majority baseline prints as exactly 0.5000. Do not read that balance as the field rate of malicious DoH
Label	`Label` benign vs malicious (tunneled)
Access	Kaggle API token required

	setting	held_out_auc
0	headline (as-is)	0.999999
1	de-duplicated (1% rows removed)	0.999999
2	shortcut feature dropped (PacketLengthMean)	0.999999

Goal of this notebook¶

What you will learn¶

Where this connects to the course text¶

Encrypted DNS-over-HTTPS Tunnel Detection (CIRA-CIC-DoHBrw-2020)¶

Model comparison + validity audit on ~499k DoH flows (benign vs tunneled)¶

1. Research problem¶

2. Literature review¶

3. Dataset provenance & honesty caveats¶

Before you run this: getting the data¶

4. Solution design¶

5. Implementation architecture¶

6. Data acquisition & preparation¶

7. Exploratory data analysis¶

8. Model comparison¶

9. Results¶

10. Validity audit — is the score real?¶

11. Ablation — does the headline survive removing the artifacts?¶

12. Reproducibility & robustness¶

13. Scientific conclusion¶

References¶

	model	accuracy	roc_auc	train_s
0	XGBoost	0.999800	0.999999	0.3
1	RandomForest	0.999687	0.999998	1.6
2	LightGBM	0.999719	0.999953	1.1
3	LogisticRegression	0.905359	0.960478	0.3
4	MajorityBaseline	0.500000	0.500000	0.0