1"""Diagnostic baseline: classical regressors on simple spectral features, with
2the same held-out split as train.py, to show what the spectrum identifies in
3principle independently of the network.
5 python3 ml/baseline.py --data DATASET_DIR [--val-every 8]
7Features per clip, from the cached spectrogram: the mean log spectrum over
8the whole clip, over the first 0.5 s, and over the last 1.5 s (1024 bins each,
9pooled 4:1 to 256), and the frequencies of the first 8 spectral peaks of the
10whole-clip mean. Models: ridge regression (alpha by cross-validation) and
11gradient boosting, one per target, on standardized log targets.
12"""
14import argparse
15import math
16import os
17import sys
19import numpy as np
20from sklearn.ensemble import HistGradientBoostingRegressor
21from sklearn.linear_model import RidgeCV
22from sklearn.pipeline import make_pipeline
23from sklearn.preprocessing import StandardScaler
25sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
26from features import FS, N_FFT, HOP, build_cache # noqa: E402
27from train import TARGETS, DERIVED, metrics # noqa: E402
29DERIVED = dict(DERIVED, **{"kappa/area": lambda d: d["h"] + 0.5 * (d["E"] - d["rho"]) - d["Lx"] - d["Ly"]})
32def peak_freqs(spec_db, n=8, min_prom=6.0):
33 """Frequencies [Hz] of the first n peaks of a mean spectrum in dB, by increasing frequency; missing ones repeat the last."""
34 from scipy.signal import find_peaks
35 pk, _ = find_peaks(spec_db, prominence=min_prom)
36 f = pk * FS / N_FFT
37 f = list(f[:n])
38 while len(f) < n:
39 f.append(f[-1] if f else 0.0)
40 return np.log(np.array(f) + 1.0)
43def features_of(x):
44 """x: (1024, T) in [0, 1] (dB / 80 + 1)."""
45 db = (x.astype(np.float32) - 1) * 80
46 fps = FS / HOP
47 whole = db.mean(1)
48 early = db[:, : int(0.5 * fps)].mean(1)
49 late = db[:, int(1.5 * fps):].mean(1)
50 pool = lambda v: v.reshape(-1, 4).mean(1)
51 return np.concatenate([pool(whole), pool(early), pool(late), pool(early - late), peak_freqs(whole)])
54def main():
55 ap = argparse.ArgumentParser()
56 ap.add_argument("--data", required=True)
57 ap.add_argument("--val-every", type=int, default=8)
58 ap.add_argument("--max-p2p", type=float, default=None, help="keep only clips whose initial shape has a peak-to-peak height below this [m]")
59 args = ap.parse_args()
60 rows, feats = build_cache(args.data, os.path.join(args.data, "ml", "cache"))
61 if args.max_p2p is not None:
62 keep = [i for i, r in enumerate(rows) if r["peakToPeak"] < args.max_p2p]
63 rows, feats = [rows[i] for i in keep], feats[keep]
64 X = np.stack([features_of(f) for f in feats])
65 y = np.stack([[math.log(r[t]) for t in TARGETS] for r in rows]).astype(np.float64)
66 idx = np.array([r["index"] for r in rows])
67 val = idx % args.val_every == 0
68 tr = ~val
69 print(f"{len(rows)} clips: {tr.sum()} train, {val.sum()} val; {X.shape[1]} features")
70 mu = y[tr].mean(0)
71 for name, make in [("ridge", lambda: make_pipeline(StandardScaler(), RidgeCV(alphas=np.logspace(-2, 4, 25)))),
72 ("gbdt", lambda: HistGradientBoostingRegressor(max_iter=400, learning_rate=0.05, max_leaf_nodes=15, l2_regularization=1.0))]:
73 pred = np.zeros((val.sum(), len(TARGETS)))
74 for i, t in enumerate(TARGETS):
75 m = make().fit(X[tr], y[tr, i])
76 pred[:, i] = m.predict(X[val])
77 d_pred = {t: pred[:, i] for i, t in enumerate(TARGETS)}
78 d_true = {t: y[val, i] for i, t in enumerate(TARGETS)}
79 d_mean = {t: np.full(val.sum(), mu[i]) for i, t in enumerate(TARGETS)}
80 print(f"--- {name}")
81 for t in TARGETS:
82 r = metrics(d_pred[t], d_true[t], d_mean[t])
83 print(f"{t:>10}: rel RMS {r['rel_rms_pct']:5.1f}% (baseline {r['baseline_rel_rms_pct']:5.1f}%) R^2 {r['r2']:5.2f}")
84 for t, f in DERIVED.items():
85 r = metrics(f(d_pred), f(d_true), f(d_mean))
86 print(f"{t:>10}: rel RMS {r['rel_rms_pct']:5.1f}% (baseline {r['baseline_rel_rms_pct']:5.1f}%) R^2 {r['r2']:5.2f}")
89if __name__ == "__main__":
90 main()