Add dataset generation scripts, ML fitting code and reference recordings
26 changed files+2410−0
.gitignoremodified+2−0View file
@@ -1,2 +1,4 @@
11 node_modules/
22 dist/
3+dataset/
4+__pycache__/
hit1.wavadded+0−0View file
Binary file not shown.
hit1_aligned.wavadded+0−0View file
Binary file not shown.
ml/README.mdadded+24−0View file
@@ -0,0 +1,24 @@
1+# Predicting the sheet parameters from the sound
2+
3+`ml/` holds a first inverse model: a convolutional network that maps the spectrogram of a single strike to the logarithms of E, rho, h, Lx, and Ly. It is trained on the clips of `scripts/generate.ts` (see `scripts/README-dataset.md`) and meant to be retrained as the clips accumulate, so that the error can be followed against the size of the dataset.
4+
5+ python3 ml/train.py --data DATASET_DIR --out RUN_DIR [--epochs 60] [--val-every 8] [--device cuda]
6+ ml/retrain-loop.sh DATASET_DIR [interval s] [epochs] # retrain periodically
7+
8+Requires numpy, torch, and matplotlib (for the scatter plots). On the office machine, `module load python`; its PyTorch cannot find a cuDNN engine for the convolution backward pass, so `train.py` disables cuDNN unless `CUDNN=1` is set.
9+
10+## Input
11+
12+`ml/features.py` turns each WAV into a log-magnitude spectrogram: 8192-point Hann windows (171 ms, 5.9 Hz bins) every 1024 samples (21 ms), 0 to 6 kHz, in dB relative to the clip's loudest bin, floored at -80 dB and scaled to [0, 1]; 1024 by 141 values per clip, cached as float16 in `DATASET_DIR/ml/cache/`. The long window is chosen because the sheet's low partials lie tens of hertz apart and ring for seconds, so frequency resolution matters more than time resolution. The clip's level does not enter: the model is meant to transfer to recordings without a calibrated level.
13+
14+## Model and training
15+
16+Six blocks of two 3x3 convolutions with batch normalization and a 2x2 max pool (16 to 128 channels), global average pooling, and a two-layer head. A second input channel carries the frequency coordinate of each row, so that absolute frequency, which is what the physical scale of the sheet sets, survives the pooling. The targets are standardized logarithms and the loss is their mean squared error; AdamW with a one-cycle schedule; the epoch with the lowest validation loss is kept. Every eighth clip by index is held out.
17+
18+## Report
19+
20+`RUN_DIR/metrics.json` and `RUN_DIR/report.html` (also copied to `DATASET_DIR/ml/report.html`) give, for each target and for E/rho, kappa = h sqrt(E / (12 rho (1 - nu^2))), the area, and the aspect ratio, the relative RMS error and the relative MAE on the held-out clips, the error of predicting the training mean, and R^2, with scatter plots. `DATASET_DIR/ml/history.jsonl` collects one line per run, so the report's last table shows the error against the number of clips.
21+
22+## What to expect
23+
24+In the linear regime the motion depends on E and rho only through E/rho, and the bending frequencies scale as kappa / L^2; the level-invariant spectrogram therefore determines E/rho, the aspect ratio, and the combination kappa / (Lx Ly) well, while E and rho separately, and h against the size, rest on the nonlinearity (amplitude-dependent frequency shifts at the strike speeds used) and on the contact dynamics. The derived rows of the report show which combinations are identified.
ml/baseline.pyadded+90−0View file
@@ -0,0 +1,90 @@
1+"""Diagnostic baseline: classical regressors on simple spectral features, with
2+the same held-out split as train.py, to show what the spectrum identifies in
3+principle independently of the network.
4+
5+ python3 ml/baseline.py --data DATASET_DIR [--val-every 8]
6+
7+Features per clip, from the cached spectrogram: the mean log spectrum over
8+the whole clip, over the first 0.5 s, and over the last 1.5 s (1024 bins each,
9+pooled 4:1 to 256), and the frequencies of the first 8 spectral peaks of the
10+whole-clip mean. Models: ridge regression (alpha by cross-validation) and
11+gradient boosting, one per target, on standardized log targets.
12+"""
13+
14+import argparse
15+import math
16+import os
17+import sys
18+
19+import numpy as np
20+from sklearn.ensemble import HistGradientBoostingRegressor
21+from sklearn.linear_model import RidgeCV
22+from sklearn.pipeline import make_pipeline
23+from sklearn.preprocessing import StandardScaler
24+
25+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
26+from features import FS, N_FFT, HOP, build_cache # noqa: E402
27+from train import TARGETS, DERIVED, metrics # noqa: E402
28+
29+DERIVED = dict(DERIVED, **{"kappa/area": lambda d: d["h"] + 0.5 * (d["E"] - d["rho"]) - d["Lx"] - d["Ly"]})
30+
31+
32+def 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)
41+
42+
43+def 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)])
52+
53+
54+def 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}")
87+
88+
89+if __name__ == "__main__":
90+ main()
ml/ddsp.pyadded+285−0View file
@@ -0,0 +1,285 @@
1+"""A small DDSP-style generative model of the struck-sheet clips.
2+
3+ python3 ml/ddsp.py train --data DATASET_DIR --out RUN_DIR [--minutes 55] [--latent 8]
4+ python3 ml/ddsp.py sample --model RUN_DIR/model.pt --out DIR [--n 8] [--seed 1]
5+
6+Encoder: a convolutional network on the cached log-frequency spectrogram
7+(384 x 133) gives the mean and log-variance of an 8-dimensional latent z
8+(a VAE, so that draws from the prior and interpolations are valid sounds).
9+
10+Decoder: z -> the controls of a parametric synthesizer:
11+ - K partials with frequencies f_k (sorted, 30 Hz to 8 kHz), initial
12+ amplitudes a_k, decay rates d_k, and a pitch-glide amount g_k, so that
13+ partial k is a_k exp(-d_k t) sin(phase of f_k (1 + g_k exp(-t / tau)));
14+ - a noise burst with a 48-band spectral shape, an amplitude, and a decay,
15+ for the attack transient and the broadband part of the ring.
16+Synthesis is plain arithmetic at 24 kHz (the simulations carry little energy
17+above 6 kHz), so the same model can be run in JavaScript.
18+
19+Loss: multi-scale log-spectral distance between the synthesized and the
20+original clip (FFT sizes 256 to 4096), plus a small KL term.
21+"""
22+
23+import argparse
24+import json
25+import math
26+import os
27+import sys
28+import time
29+import wave
30+
31+import numpy as np
32+import torch
33+import torch.nn as nn
34+import torch.nn.functional as Fn
35+from scipy.signal import resample_poly
36+
37+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
38+from features import build_cache # noqa: E402
39+
40+torch.backends.cudnn.enabled = False
41+SR = 24000
42+DUR = 3.0
43+N = int(SR * DUR)
44+K = 128
45+NOISE_BANDS = 48
46+F_LO, F_HI = 30.0, 8000.0
47+
48+
49+def read_wav24(path):
50+ with wave.open(path) as w:
51+ x = np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(np.float32) / 32768
52+ y = resample_poly(x, 1, 2).astype(np.float32)
53+ out = np.zeros(N, np.float32)
54+ out[: min(N, len(y))] = y[:N]
55+ return out
56+
57+
58+def write_wav(path, y, sr=SR):
59+ y = np.asarray(y, np.float32)
60+ y = y / (np.abs(y).max() + 1e-9) * 0.95
61+ with wave.open(path, "wb") as w:
62+ w.setnchannels(1); w.setsampwidth(2); w.setframerate(sr)
63+ w.writeframes((y * 32767).astype("<i2").tobytes())
64+
65+
66+class Encoder(nn.Module):
67+ def __init__(self, latent):
68+ super().__init__()
69+ ch, layers, c_in = (16, 32, 64, 96, 128), [], 1
70+ for c in ch:
71+ layers += [nn.Conv2d(c_in, c, 3, padding=1), nn.BatchNorm2d(c), nn.ReLU(inplace=True), nn.MaxPool2d(2)]
72+ c_in = c
73+ self.conv = nn.Sequential(*layers)
74+ self.fc = nn.Sequential(nn.Linear(c_in * 12 * 4, 256), nn.ReLU(inplace=True), nn.Linear(256, 2 * latent))
75+
76+ def forward(self, x): # (B, 1, 384, 133)
77+ h = self.conv(x) # (B, 128, 12, 4)
78+ return self.fc(h.flatten(1)).chunk(2, -1)
79+
80+
81+class Decoder(nn.Module):
82+ def __init__(self, latent):
83+ super().__init__()
84+ self.mlp = nn.Sequential(nn.Linear(latent, 256), nn.LayerNorm(256), nn.ReLU(inplace=True), nn.Linear(256, 512), nn.LayerNorm(512), nn.ReLU(inplace=True))
85+ self.freq = nn.Linear(512, K)
86+ self.amp = nn.Linear(512, K)
87+ self.decay = nn.Linear(512, K)
88+ self.glide = nn.Linear(512, K)
89+ self.noise_shape = nn.Linear(512, NOISE_BANDS)
90+ self.noise_amp = nn.Linear(512, 1)
91+ self.noise_decay = nn.Linear(512, 1)
92+ self.glide_tau = nn.Linear(512, 1)
93+
94+ def forward(self, z):
95+ h = self.mlp(z)
96+ # Sorted frequencies: cumulative positive increments spanning the log range.
97+ inc = Fn.softplus(self.freq(h)) + 1e-3
98+ pos = torch.cumsum(inc, -1) / inc.sum(-1, keepdim=True) # (B, K) in (0, 1]
99+ f = F_LO * (F_HI / F_LO) ** (pos * 0.98 + 0.005)
100+ a = torch.sigmoid(self.amp(h) - 2.0) # amplitudes in (0, 1)
101+ d = 0.3 + 40 * torch.sigmoid(self.decay(h) - 1.0) # 1/s
102+ g = 0.15 * torch.tanh(self.glide(h)) # relative initial detuning
103+ tau = 0.02 + 0.5 * torch.sigmoid(self.glide_tau(h)) # s
104+ ns = torch.sigmoid(self.noise_shape(h) - 2.0)
105+ na = torch.sigmoid(self.noise_amp(h) - 2.0)
106+ nd = 2 + 60 * torch.sigmoid(self.noise_decay(h))
107+ return dict(f=f, a=a, d=d, g=g, tau=tau, noise_shape=ns, noise_amp=na, noise_decay=nd)
108+
109+
110+def synthesize(c, n=N, noise_seed=None):
111+ """Renders a batch of control dicts to audio (B, n) at SR."""
112+ B = c["f"].shape[0]
113+ dev = c["f"].device
114+ t = torch.arange(n, device=dev, dtype=torch.float32) / SR
115+ # Partials: instantaneous frequency with a decaying glide; phase by cumulative sum.
116+ f_inst = c["f"][:, :, None] * (1 + c["g"][:, :, None] * torch.exp(-t[None, None, :] / c["tau"][:, :, None])) # (B, K, n)
117+ phase = 2 * math.pi * torch.cumsum(f_inst, -1) / SR
118+ env = c["a"][:, :, None] * torch.exp(-c["d"][:, :, None] * t[None, None, :])
119+ nyq = (c["f"] < SR / 2 - 200).float()[:, :, None]
120+ y = (env * nyq * torch.sin(phase)).sum(1)
121+ # Noise burst: white noise shaped in the STFT domain by a fixed spectral shape and a decaying envelope.
122+ gen = torch.Generator(device=dev)
123+ if noise_seed is not None:
124+ gen.manual_seed(noise_seed)
125+ noise = torch.randn(B, n, device=dev, generator=gen)
126+ n_fft, hop = 512, 128
127+ win = torch.hann_window(n_fft, device=dev)
128+ S = torch.view_as_real(torch.stft(noise, n_fft, hop, window=win, return_complex=True)) # (B, 257, frames, 2)
129+ shape = Fn.interpolate(c["noise_shape"][:, None, :], size=S.shape[1], mode="linear", align_corners=True)[:, 0, :] # (B, 257)
130+ tf = torch.arange(S.shape[2], device=dev, dtype=torch.float32) * hop / SR
131+ envn = c["noise_amp"] * torch.exp(-c["noise_decay"] * tf[None, :]) # (B, frames)
132+ S = S * (shape[:, :, None] * envn[:, None, :])[..., None]
133+ yn = torch.istft(torch.view_as_complex(S.contiguous()), n_fft, hop, window=win, length=n)
134+ return y + 0.05 * yn
135+
136+
137+def magnitude(x, n_fft):
138+ """STFT magnitude without complex-tensor ops (the office PyTorch cannot JIT-compile those)."""
139+ win = torch.hann_window(n_fft, device=x.device)
140+ S = torch.view_as_real(torch.stft(x, n_fft, n_fft // 4, window=win, return_complex=True))
141+ return torch.sqrt(S[..., 0] ** 2 + S[..., 1] ** 2 + 1e-12)
142+
143+
144+def spectral_loss(y, x):
145+ loss = 0.0
146+ for n_fft in (256, 512, 1024, 2048, 4096):
147+ Y = magnitude(y, n_fft)
148+ X = magnitude(x, n_fft)
149+ loss = loss + (Y - X).abs().mean() + 0.5 * (torch.log(Y + 1e-4) - torch.log(X + 1e-4)).abs().mean()
150+ return loss / 5
151+
152+
153+class Model(nn.Module):
154+ def __init__(self, latent=8):
155+ super().__init__()
156+ self.latent = latent
157+ self.enc = Encoder(latent)
158+ self.dec = Decoder(latent)
159+
160+ def encode(self, spec):
161+ mu, logvar = self.enc(spec)
162+ return mu, logvar.clamp(-8, 4)
163+
164+
165+def load_audio(data, ids):
166+ cache = os.path.join(data, "ml", "audio24.npy")
167+ if os.path.exists(cache):
168+ a = np.load(cache, mmap_mode="r")
169+ if a.shape[0] == len(ids):
170+ return np.array(a)
171+ a = np.stack([read_wav24(os.path.join(data, "clips", i + ".wav")) for i in ids])
172+ np.save(cache, a)
173+ return a
174+
175+
176+def train(args):
177+ os.makedirs(os.path.join(args.out, "samples"), exist_ok=True)
178+ rows, feats = build_cache(args.data, os.path.join(args.data, "ml", "cache-logf"), kind="logf")
179+ ids = [r["id"] for r in rows]
180+ audio = load_audio(args.data, ids)
181+ print(f"{len(ids)} clips, spectrograms {feats.shape}, audio {audio.shape}", flush=True)
182+ dev = torch.device("cuda" if torch.cuda.is_available() else "cpu")
183+ X = torch.from_numpy(feats)
184+ A = torch.from_numpy(audio)
185+ model = Model(args.latent).to(dev)
186+ opt = torch.optim.Adam(model.parameters(), lr=args.lr)
187+ rng = np.random.default_rng(0)
188+ fixed = [0, 1234, 2500, 4321]
189+ t0, step, last_sample = time.time(), 0, 0
190+ log = open(os.path.join(args.out, "train.log"), "a")
191+ while (time.time() - t0) / 60 < args.minutes:
192+ perm = rng.permutation(len(ids))
193+ for k in range(0, len(perm) - args.batch + 1, args.batch):
194+ b = perm[k:k + args.batch]
195+ spec = X[b].to(dev).float().unsqueeze(1)
196+ x = A[b].to(dev)
197+ mu, logvar = model.encode(spec)
198+ z = mu + torch.randn_like(mu) * torch.exp(0.5 * logvar)
199+ y = synthesize(model.dec(z))
200+ rec = spectral_loss(y, x)
201+ kl = -0.5 * torch.mean(1 + logvar - mu ** 2 - logvar.exp())
202+ beta = args.beta * min(1.0, step / 2000)
203+ loss = rec + beta * kl
204+ opt.zero_grad(); loss.backward()
205+ nn.utils.clip_grad_norm_(model.parameters(), 5.0)
206+ opt.step()
207+ step += 1
208+ if step % 50 == 0:
209+ msg = f"step {step} ({(time.time() - t0) / 60:.1f} min): rec {rec.item():.4f} kl {kl.item():.3f}"
210+ print(msg, flush=True); log.write(msg + "\n"); log.flush()
211+ if (time.time() - t0) / 60 - last_sample >= args.sample_every or (time.time() - t0) / 60 >= args.minutes:
212+ last_sample = (time.time() - t0) / 60
213+ torch.save({"state": model.state_dict(), "latent": args.latent, "step": step}, os.path.join(args.out, "model.pt"))
214+ write_samples(model, X, ids, fixed, os.path.join(args.out, "samples", f"m{int(last_sample):03d}"), dev, step)
215+ write_index(args.out)
216+ if (time.time() - t0) / 60 >= args.minutes:
217+ break
218+ print("done", step, "steps")
219+
220+
221+@torch.no_grad()
222+def write_samples(model, X, ids, fixed, out, dev, step, n_prior=6):
223+ model.eval()
224+ os.makedirs(out, exist_ok=True)
225+ names = []
226+ spec = X[fixed].to(dev).float().unsqueeze(1)
227+ mu, _ = model.encode(spec)
228+ y = synthesize(model.dec(mu), noise_seed=1).cpu().numpy()
229+ for i, k in enumerate(fixed):
230+ write_wav(os.path.join(out, f"recon-{ids[k]}.wav"), y[i]); names.append((f"reconstruction of {ids[k]}", f"recon-{ids[k]}.wav"))
231+ g = torch.Generator(device=dev); g.manual_seed(7)
232+ zp = torch.randn(n_prior, model.latent, device=dev, generator=g)
233+ y = synthesize(model.dec(zp), noise_seed=2).cpu().numpy()
234+ for i in range(n_prior):
235+ write_wav(os.path.join(out, f"prior-{i}.wav"), y[i]); names.append((f"prior sample {i}", f"prior-{i}.wav"))
236+ steps = 5
237+ zi = torch.stack([mu[0] + (mu[1] - mu[0]) * s / (steps - 1) for s in range(steps)])
238+ y = synthesize(model.dec(zi), noise_seed=3).cpu().numpy()
239+ for s in range(steps):
240+ write_wav(os.path.join(out, f"interp-{s}.wav"), y[s]); names.append((f"interpolation {ids[fixed[0]]} -> {ids[fixed[1]]}, step {s}", f"interp-{s}.wav"))
241+ with open(os.path.join(out, "list.json"), "w") as f:
242+ json.dump({"step": step, "items": names}, f)
243+ model.train()
244+
245+
246+def write_index(run):
247+ sdir = os.path.join(run, "samples")
248+ snaps = sorted(d for d in os.listdir(sdir) if os.path.isdir(os.path.join(sdir, d)))
249+ blocks = []
250+ for d in reversed(snaps):
251+ try:
252+ meta = json.load(open(os.path.join(sdir, d, "list.json")))
253+ except Exception:
254+ continue
255+ rows = "".join(f"<tr><td class='l'>{label}</td><td class='l'><audio controls preload='none' src='samples/{d}/{fn}'></audio></td></tr>" for label, fn in meta["items"])
256+ blocks.append(f"<h2>after {int(d[1:])} min, step {meta['step']}</h2><table>{rows}</table>")
257+ html = ("<!doctype html><meta charset='utf-8'><title>DDSP sheet model samples</title><style>body{font:14px system-ui,sans-serif;margin:16px} table{border-collapse:collapse} td{padding:3px 8px;border-bottom:1px solid #ddd} td.l{text-align:left} audio{height:28px;width:260px} h2{font-size:15px;margin-top:18px}</style>"
258+ "<h1>DDSP-style sheet model: samples during training</h1><p>Reconstructions of four fixed dataset clips (compare with clips/ID.wav in the gallery), six draws from the prior, and a five-step interpolation between two clips. Newest first.</p>" + "".join(blocks))
259+ with open(os.path.join(run, "samples.html"), "w") as f:
260+ f.write(html)
261+
262+
263+@torch.no_grad()
264+def sample(args):
265+ ck = torch.load(args.model, map_location="cpu", weights_only=False)
266+ model = Model(ck["latent"]); model.load_state_dict(ck["state"]); model.eval()
267+ os.makedirs(args.out, exist_ok=True)
268+ g = torch.Generator(); g.manual_seed(args.seed)
269+ z = torch.randn(args.n, ck["latent"], generator=g) if args.z is None else torch.tensor([[float(v) for v in args.z.split(",")]])
270+ y = synthesize(model.dec(z), noise_seed=args.seed).numpy()
271+ for i in range(y.shape[0]):
272+ write_wav(os.path.join(args.out, f"sample-{i}.wav"), y[i])
273+ print(f"sample-{i}.wav z = " + ", ".join(f"{v:.2f}" for v in z[i].tolist()))
274+
275+
276+if __name__ == "__main__":
277+ ap = argparse.ArgumentParser()
278+ sub = ap.add_subparsers(dest="cmd", required=True)
279+ t = sub.add_parser("train"); t.add_argument("--data", required=True); t.add_argument("--out", required=True)
280+ t.add_argument("--minutes", type=float, default=55); t.add_argument("--latent", type=int, default=8); t.add_argument("--batch", type=int, default=16)
281+ t.add_argument("--lr", type=float, default=1e-3); t.add_argument("--beta", type=float, default=1e-3); t.add_argument("--sample-every", type=float, default=10)
282+ s = sub.add_parser("sample"); s.add_argument("--model", required=True); s.add_argument("--out", required=True)
283+ s.add_argument("--n", type=int, default=8); s.add_argument("--seed", type=int, default=1); s.add_argument("--z", default=None, help="comma-separated latent vector for one clip")
284+ a = ap.parse_args()
285+ (train if a.cmd == "train" else sample)(a)
ml/ddsp2.pyadded+245−0View file
@@ -0,0 +1,245 @@
1+"""Second generative model: an oscillator bank on a fixed log-frequency grid.
2+
3+ python3 ml/ddsp2.py train --data DATASET_DIR --out RUN_DIR [--minutes 20] [--latent 8]
4+ python3 ml/ddsp2.py sample --model RUN_DIR/model.pt --out DIR [--n 8] [--seed 1] [--z v1,...]
5+
6+Instead of letting the decoder choose partial frequencies (which the spectral
7+loss cannot steer, see ddsp.py), the synthesizer has one oscillator at each of
8+M = 96 bins per octave from 30 Hz to 8 kHz (about 770), and the decoder emits
9+for each an amplitude, a decay rate, and a frequency offset within its bin.
10+A partial is then an active bin, its exact frequency the bin plus the offset,
11+and the amplitudes have well-behaved gradients. The decoder is a 1D
12+transposed-convolutional network over the frequency axis, so that patterns of
13+lines can shift along it easily (a change of the sheet's frequency scale is a
14+translation on this axis). A shaped noise burst covers the attack. The loss
15+adds the smoothed log-frequency "line picture" of ml/linesim.py, computed in
16+torch, to the multi-scale spectral loss.
17+"""
18+
19+import argparse
20+import json
21+import math
22+import os
23+import sys
24+import time
25+
26+import numpy as np
27+import torch
28+import torch.nn as nn
29+import torch.nn.functional as Fn
30+
31+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
32+from features import build_cache # noqa: E402
33+from ddsp import SR, DUR, N, NOISE_BANDS, Encoder, magnitude, spectral_loss, load_audio, write_wav, write_index # noqa: E402
34+
35+torch.backends.cudnn.enabled = False
36+F_LO, F_HI = 30.0, 8000.0
37+BPO = 96
38+M = int(math.log2(F_HI / F_LO) * BPO) # oscillators
39+GRID = F_LO * 2 ** (np.arange(M) / BPO)
40+HALF_BIN = 2 ** (0.5 / BPO) - 1 # relative half width of a bin
41+
42+
43+class Decoder(nn.Module):
44+ def __init__(self, latent, base=25, ch=64):
45+ super().__init__()
46+ self.base = base
47+ self.fc = nn.Sequential(nn.Linear(latent, 256), nn.ReLU(inplace=True), nn.Linear(256, ch * base), nn.ReLU(inplace=True))
48+ self.ch = ch
49+ ups = []
50+ c = ch
51+ for k in range(5): # base * 32 = 800 >= M
52+ ups += [nn.ConvTranspose1d(c + 1, ch, 4, stride=2, padding=1), nn.GroupNorm(8, ch), nn.ReLU(inplace=True)]
53+ c = ch
54+ self.ups = nn.ModuleList(ups)
55+ self.head = nn.Conv1d(ch + 1, 3, 5, padding=2)
56+ self.noise = nn.Sequential(nn.Linear(latent, 128), nn.ReLU(inplace=True), nn.Linear(128, NOISE_BANDS + 2))
57+
58+ def forward(self, z):
59+ h = self.fc(z).view(-1, self.ch, self.base)
60+ i = 0
61+ while i < len(self.ups):
62+ coord = torch.linspace(0, 1, h.shape[-1], device=h.device).expand(h.shape[0], 1, -1)
63+ h = self.ups[i](torch.cat([h, coord], 1)); h = self.ups[i + 1](h); h = self.ups[i + 2](h)
64+ i += 3
65+ coord = torch.linspace(0, 1, h.shape[-1], device=h.device).expand(h.shape[0], 1, -1)
66+ out = self.head(torch.cat([h, coord], 1))[:, :, :M] # (B, 3, M)
67+ amp = torch.sigmoid(out[:, 0] - 3.0)
68+ decay = 0.3 + 40 * torch.sigmoid(out[:, 1] - 1.0)
69+ offset = torch.tanh(out[:, 2]) * HALF_BIN
70+ nz = self.noise(z)
71+ return dict(amp=amp, decay=decay, offset=offset, noise_shape=torch.sigmoid(nz[:, :NOISE_BANDS] - 2.0),
72+ noise_amp=torch.sigmoid(nz[:, NOISE_BANDS] - 2.0), noise_decay=2 + 60 * torch.sigmoid(nz[:, NOISE_BANDS + 1]))
73+
74+
75+_GRID_T = None
76+
77+
78+def synthesize(c, n=N, noise_seed=None):
79+ global _GRID_T
80+ B = c["amp"].shape[0]
81+ dev = c["amp"].device
82+ if _GRID_T is None or _GRID_T.device != dev:
83+ _GRID_T = torch.tensor(GRID, dtype=torch.float32, device=dev)
84+ t = torch.arange(n, device=dev, dtype=torch.float32) / SR
85+ f = _GRID_T[None, :] * (1 + c["offset"]) # (B, M)
86+ usable = (f < SR / 2 - 200).float()
87+ phase = 2 * math.pi * f[:, :, None] * t[None, None, :]
88+ env = (c["amp"] * usable)[:, :, None] * torch.exp(-c["decay"][:, :, None] * t[None, None, :])
89+ y = (env * torch.sin(phase)).sum(1)
90+ gen = torch.Generator(device=dev)
91+ if noise_seed is not None:
92+ gen.manual_seed(noise_seed)
93+ noise = torch.randn(B, n, device=dev, generator=gen)
94+ n_fft, hop = 512, 128
95+ win = torch.hann_window(n_fft, device=dev)
96+ S = torch.view_as_real(torch.stft(noise, n_fft, hop, window=win, return_complex=True))
97+ shape = Fn.interpolate(c["noise_shape"][:, None, :], size=S.shape[1], mode="linear", align_corners=True)[:, 0, :]
98+ tf = torch.arange(S.shape[2], device=dev, dtype=torch.float32) * hop / SR
99+ envn = c["noise_amp"][:, None] * torch.exp(-c["noise_decay"][:, None] * tf[None, :])
100+ S = S * (shape[:, :, None] * envn[:, None, :])[..., None]
101+ yn = torch.istft(torch.view_as_complex(S.contiguous()), n_fft, hop, window=win, length=n)
102+ return y + 0.05 * yn
103+
104+
105+# ---------- line-picture loss (linesim.py in torch) ----------
106+
107+_PIC = {}
108+
109+
110+def picture_torch(x, tol=0.02, n_fft=4096, hop=1024, slices=((0, 0.3), (0.3, 1.0), (1.0, 2.0), (2.0, 3.0))):
111+ dev = x.device
112+ key = (tol, str(dev))
113+ if key not in _PIC:
114+ f_lin = np.arange(n_fft // 2 + 1) * SR / n_fft
115+ centers = F_LO * 2 ** (np.arange(int(math.log2(5000 / F_LO) * 48)) / 48)
116+ W = np.exp(-0.5 * ((np.log(np.maximum(f_lin[None, :], 1e-3) / centers[:, None])) / tol) ** 2)
117+ W[:, 0] = 0
118+ W /= np.maximum(W.sum(1, keepdims=True), 1e-12)
119+ _PIC[key] = torch.tensor(W, dtype=torch.float32, device=dev)
120+ W = _PIC[key]
121+ P = magnitude(x, n_fft) ** 2 # (B, bins, frames), hop n_fft // 4
122+ fps = SR / (n_fft // 4)
123+ outs = []
124+ for a, b in slices:
125+ sl = P[:, :, int(a * fps): max(int(a * fps) + 1, int(b * fps))].mean(-1)
126+ outs.append(W @ sl.T) # (n_log, B)
127+ db = 10 * torch.log10(torch.stack(outs, 0) + 1e-12) # (K, n_log, B)
128+ db = db - db.amax((0, 1), keepdim=True)
129+ return torch.clamp(db, min=-70)
130+
131+
132+def picture_loss(y, x):
133+ return (picture_torch(y) - picture_torch(x)).abs().mean() / 10
134+
135+
136+class Model(nn.Module):
137+ def __init__(self, latent=8):
138+ super().__init__()
139+ self.latent = latent
140+ self.enc = Encoder(latent)
141+ self.dec = Decoder(latent)
142+
143+ def encode(self, spec):
144+ mu, logvar = self.enc(spec)
145+ return mu, logvar.clamp(-8, 4)
146+
147+
148+def train(args):
149+ os.makedirs(os.path.join(args.out, "samples"), exist_ok=True)
150+ rows, feats = build_cache(args.data, os.path.join(args.data, "ml", "cache-logf"), kind="logf")
151+ ids = [r["id"] for r in rows]
152+ audio = load_audio(args.data, ids)
153+ print(f"{len(ids)} clips, {M} oscillators, audio {audio.shape}", flush=True)
154+ dev = torch.device("cuda" if torch.cuda.is_available() else "cpu")
155+ X = torch.from_numpy(feats); A = torch.from_numpy(audio)
156+ model = Model(args.latent).to(dev)
157+ print(f"{sum(p.numel() for p in model.parameters()) / 1e6:.2f} M parameters", flush=True)
158+ opt = torch.optim.Adam(model.parameters(), lr=args.lr)
159+ rng = np.random.default_rng(0)
160+ fixed = [0, 1234, 2500, 4321]
161+ t0, step, last_sample = time.time(), 0, 0
162+ log = open(os.path.join(args.out, "train.log"), "a")
163+ done = False
164+ while not done:
165+ perm = rng.permutation(len(ids))
166+ for k in range(0, len(perm) - args.batch + 1, args.batch):
167+ b = perm[k:k + args.batch]
168+ spec = X[b].to(dev).float().unsqueeze(1)
169+ x = A[b].to(dev)
170+ mu, logvar = model.encode(spec)
171+ z = mu + torch.randn_like(mu) * torch.exp(0.5 * logvar)
172+ y = synthesize(model.dec(z))
173+ rec = spectral_loss(y, x)
174+ pic = picture_loss(y, x)
175+ kl = -0.5 * torch.mean(1 + logvar - mu ** 2 - logvar.exp())
176+ beta = args.beta * min(1.0, step / 3000)
177+ loss = rec + pic + beta * torch.clamp(kl, min=args.free_bits)
178+ opt.zero_grad(); loss.backward()
179+ nn.utils.clip_grad_norm_(model.parameters(), 5.0)
180+ opt.step()
181+ step += 1
182+ minutes = (time.time() - t0) / 60
183+ if step % 50 == 0:
184+ msg = f"step {step} ({minutes:.1f} min): rec {rec.item():.4f} pic {pic.item():.4f} ({10 * pic.item():.2f} dB) kl {kl.item():.3f}"
185+ print(msg, flush=True); log.write(msg + "\n"); log.flush()
186+ if minutes - last_sample >= args.sample_every or minutes >= args.minutes:
187+ last_sample = minutes
188+ torch.save({"state": model.state_dict(), "latent": args.latent, "step": step, "kind": "bank"}, os.path.join(args.out, "model.pt"))
189+ write_samples(model, X, ids, fixed, os.path.join(args.out, "samples", f"m{int(minutes):03d}"), dev, step)
190+ write_index(args.out)
191+ if minutes >= args.minutes:
192+ done = True
193+ break
194+ print("done", step, "steps")
195+
196+
197+@torch.no_grad()
198+def write_samples(model, X, ids, fixed, out, dev, step, n_prior=6):
199+ model.eval()
200+ os.makedirs(out, exist_ok=True)
201+ names = []
202+ spec = X[fixed].to(dev).float().unsqueeze(1)
203+ mu, _ = model.encode(spec)
204+ y = synthesize(model.dec(mu), noise_seed=1).cpu().numpy()
205+ for i, k in enumerate(fixed):
206+ write_wav(os.path.join(out, f"recon-{ids[k]}.wav"), y[i]); names.append((f"reconstruction of {ids[k]}", f"recon-{ids[k]}.wav"))
207+ g = torch.Generator(device=dev); g.manual_seed(7)
208+ zp = torch.randn(n_prior, model.latent, device=dev, generator=g)
209+ y = synthesize(model.dec(zp), noise_seed=2).cpu().numpy()
210+ for i in range(n_prior):
211+ write_wav(os.path.join(out, f"prior-{i}.wav"), y[i]); names.append((f"prior sample {i}", f"prior-{i}.wav"))
212+ steps = 5
213+ zi = torch.stack([mu[0] + (mu[1] - mu[0]) * s / (steps - 1) for s in range(steps)])
214+ y = synthesize(model.dec(zi), noise_seed=3).cpu().numpy()
215+ for s in range(steps):
216+ write_wav(os.path.join(out, f"interp-{s}.wav"), y[s]); names.append((f"interpolation {ids[fixed[0]]} -> {ids[fixed[1]]}, step {s}", f"interp-{s}.wav"))
217+ with open(os.path.join(out, "list.json"), "w") as f:
218+ json.dump({"step": step, "items": names}, f)
219+ model.train()
220+
221+
222+@torch.no_grad()
223+def sample(args):
224+ ck = torch.load(args.model, map_location="cpu", weights_only=False)
225+ model = Model(ck["latent"]); model.load_state_dict(ck["state"]); model.eval()
226+ os.makedirs(args.out, exist_ok=True)
227+ g = torch.Generator(); g.manual_seed(args.seed)
228+ z = torch.randn(args.n, ck["latent"], generator=g) if args.z is None else torch.tensor([[float(v) for v in args.z.split(",")]])
229+ y = synthesize(model.dec(z), noise_seed=args.seed).numpy()
230+ for i in range(y.shape[0]):
231+ write_wav(os.path.join(args.out, f"sample-{i}.wav"), y[i])
232+ print(f"sample-{i}.wav z = " + ", ".join(f"{v:.2f}" for v in z[i].tolist()))
233+
234+
235+if __name__ == "__main__":
236+ ap = argparse.ArgumentParser()
237+ sub = ap.add_subparsers(dest="cmd", required=True)
238+ t = sub.add_parser("train"); t.add_argument("--data", required=True); t.add_argument("--out", required=True)
239+ t.add_argument("--minutes", type=float, default=20); t.add_argument("--latent", type=int, default=8); t.add_argument("--batch", type=int, default=16)
240+ t.add_argument("--lr", type=float, default=5e-4); t.add_argument("--beta", type=float, default=2e-3); t.add_argument("--free-bits", type=float, default=0.5)
241+ t.add_argument("--sample-every", type=float, default=5)
242+ s = sub.add_parser("sample"); s.add_argument("--model", required=True); s.add_argument("--out", required=True)
243+ s.add_argument("--n", type=int, default=8); s.add_argument("--seed", type=int, default=1); s.add_argument("--z", default=None)
244+ a = ap.parse_args()
245+ (train if a.cmd == "train" else sample)(a)
ml/features.pyadded+116−0View file
@@ -0,0 +1,116 @@
1+"""Spectrogram features of the dataset clips, cached on disk.
2+
3+Each clip becomes a log-magnitude spectrogram in dB relative to its own
4+maximum, floored at -80 dB and scaled to [0, 1], over 0 to FMAX Hz. The level
5+of the clip therefore does not enter (the WAVs are peak-normalized anyway).
6+The frequency resolution matters more than the time resolution here: the
7+sheet's low partials are tens of hertz apart and ring for seconds, so the
8+window is long (N_FFT / 48000 = 171 ms) and the hop coarse (21 ms).
9+"""
10+
11+import json
12+import os
13+import sys
14+import wave
15+from concurrent.futures import ProcessPoolExecutor
16+
17+import numpy as np
18+
19+FS = 48000
20+N_FFT = 8192
21+HOP = 1024
22+FMAX = 6000.0
23+DB_FLOOR = 80.0
24+N_BINS = int(FMAX / FS * N_FFT) # 1024 bins of 5.86 Hz
25+# Log-frequency variant: 48 bins per octave over 8 octaves from F_LO. On this
26+# axis a change of the sheet's frequency scale (kappa / L^2) translates the
27+# whole pattern of partials, which a convolutional network detects naturally.
28+F_LO = 20.0
29+BINS_PER_OCTAVE = 48
30+N_LOGF = 8 * BINS_PER_OCTAVE
31+
32+
33+def read_wav(path):
34+ with wave.open(path) as w:
35+ assert w.getnchannels() == 1 and w.getsampwidth() == 2 and w.getframerate() == FS, path
36+ x = np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(np.float32) / 32768.0
37+ return x
38+
39+
40+def spectrogram(x):
41+ """(N_BINS, n_frames) float16 in [0, 1]; 1 at the clip's loudest bin."""
42+ win = np.hanning(N_FFT + 1)[:-1].astype(np.float32)
43+ n_frames = 1 + (len(x) - N_FFT) // HOP
44+ frames = np.lib.stride_tricks.as_strided(x, shape=(n_frames, N_FFT), strides=(x.strides[0] * HOP, x.strides[0]))
45+ spec = np.abs(np.fft.rfft(frames * win, axis=1))[:, :N_BINS].T
46+ db = 20 * np.log10(spec + 1e-9)
47+ db -= db.max()
48+ return np.clip(db / DB_FLOOR + 1, 0, 1).astype(np.float16)
49+
50+
51+_FB = None
52+
53+
54+def logf_filterbank():
55+ """(N_LOGF, N_BINS) triangular filters centered on the log-spaced frequencies, at least one linear bin wide."""
56+ global _FB
57+ if _FB is None:
58+ f_lin = np.arange(N_BINS) * FS / N_FFT
59+ centers = F_LO * 2 ** (np.arange(N_LOGF) / BINS_PER_OCTAVE)
60+ fb = np.zeros((N_LOGF, N_BINS), np.float32)
61+ for k, fc in enumerate(centers):
62+ half = max(fc * (2 ** (1 / BINS_PER_OCTAVE) - 1), FS / N_FFT)
63+ w = np.clip(1 - np.abs(f_lin - fc) / half, 0, None)
64+ fb[k] = w / max(w.sum(), 1e-9)
65+ _FB = fb
66+ return _FB
67+
68+
69+def spectrogram_logf(x):
70+ """(N_LOGF, n_frames) float16 in [0, 1] on a log-frequency axis."""
71+ win = np.hanning(N_FFT + 1)[:-1].astype(np.float32)
72+ n_frames = 1 + (len(x) - N_FFT) // HOP
73+ frames = np.lib.stride_tricks.as_strided(x, shape=(n_frames, N_FFT), strides=(x.strides[0] * HOP, x.strides[0]))
74+ power = np.abs(np.fft.rfft(frames * win, axis=1))[:, :N_BINS] ** 2
75+ db = 10 * np.log10(logf_filterbank() @ power.T + 1e-18)
76+ db -= db.max()
77+ return np.clip(db / DB_FLOOR + 1, 0, 1).astype(np.float16)
78+
79+
80+def _one(args):
81+ clip_path, out_path, kind = args
82+ np.save(out_path, (spectrogram_logf if kind == "logf" else spectrogram)(read_wav(clip_path)))
83+ return out_path
84+
85+
86+def load_index(data_dir):
87+ rows = []
88+ with open(os.path.join(data_dir, "index.jsonl")) as f:
89+ for line in f:
90+ if line.strip():
91+ rows.append(json.loads(line))
92+ # Keep only clips whose WAV is present (the index may run ahead of a copy).
93+ return [r for r in rows if os.path.exists(os.path.join(data_dir, "clips", r["id"] + ".wav"))]
94+
95+
96+def build_cache(data_dir, cache_dir, workers=8, kind="lin"):
97+ """Computes the spectrograms (kind "lin" or "logf") of every indexed clip not yet cached; returns the index rows and the stacked features."""
98+ os.makedirs(cache_dir, exist_ok=True)
99+ rows = load_index(data_dir)
100+ todo = []
101+ for r in rows:
102+ out = os.path.join(cache_dir, r["id"] + ".npy")
103+ if not os.path.exists(out):
104+ todo.append((os.path.join(data_dir, "clips", r["id"] + ".wav"), out, kind))
105+ if todo:
106+ with ProcessPoolExecutor(workers) as ex:
107+ for _ in ex.map(_one, todo, chunksize=4):
108+ pass
109+ feats = np.stack([np.load(os.path.join(cache_dir, r["id"] + ".npy")) for r in rows])
110+ return rows, feats
111+
112+
113+if __name__ == "__main__":
114+ data_dir, cache_dir = sys.argv[1], sys.argv[2]
115+ rows, feats = build_cache(data_dir, cache_dir)
116+ print(len(rows), "clips, features", feats.shape, feats.dtype)
ml/fit.pyadded+249−0View file
@@ -0,0 +1,249 @@
1+"""Fits simulation parameters to a target clip by CMA-ES over rendered batches.
2+
3+ python3 ml/fit.py TARGET.wav --out FIT_DIR [--init E=1.7e11,h=8e-4,...] [--generations 12] [--popsize 16]
4+
5+Each generation renders one batch of parameter sets with scripts/render-batch.ts
6+(one simulation per worker in parallel on the GPU), scores every clip against
7+the target, and updates the CMA-ES distribution. The free parameters, in a
8+normalized log (or linear) space, are
9+
10+ E, h, Lx, Ly, sigma0, sigma1, strikeVelocity, strikeFracX, strikeFracY, bow sag (across the width)
11+
12+with rho fixed at the default (E and rho enter the sound only through E / rho
13+at a fixed level, so one of them is redundant) and a single parabolic bow as
14+the shape. The objective combines the log-mel distance of match.py (timbre and
15+decay) with the fine-structure distances of the early and whole-clip spectra
16+(the positions of the partials). Every generation's best clip is kept in
17+FIT_DIR/best-GEN.wav and a listening page FIT_DIR/fit.html is written.
18+"""
19+
20+import argparse
21+import json
22+import math
23+import os
24+import shutil
25+import subprocess
26+import sys
27+import time
28+
29+import numpy as np
30+
31+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
32+from match import read_wav_any, onset, segment, features # noqa: E402
33+import linesim # noqa: E402
34+
35+ROOT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..")
36+NODE = ["node", "--experimental-strip-types", os.path.join(ROOT, "scripts", "render-batch.ts")]
37+
38+# name, kind (log or lin), lower, upper, default
39+SPACE = [
40+ ("E", "log", 60e9, 300e9, 170e9),
41+ ("h", "log", 0.3e-3, 2.0e-3, 0.79e-3),
42+ ("Lx", "log", 6 * 0.0254, 30 * 0.0254, 12 * 0.0254),
43+ ("Ly", "log", 6 * 0.0254, 30 * 0.0254, 18 * 0.0254),
44+ ("sigma0", "log", 0.15, 5.0, 1.15),
45+ ("sigma1", "log", 0.5e-4, 3e-3, 5e-4),
46+ ("strikeVelocity", "log", 0.8, 6.0, 3.0),
47+ ("strikeFracX", "lin", 0.5, 0.92, 0.7),
48+ ("strikeFracY", "lin", 0.12, 0.55, 0.35),
49+ ("bow", "lin", -6e-3, 6e-3, 3e-3),
50+ ("bowY", "lin", -6e-3, 6e-3, 0.0),
51+]
52+FREE = [n for n, *_ in SPACE] # set from --free
53+
54+
55+class SimpleCMA:
56+ """Minimal (mu/mu_w, lambda)-CMA-ES (Hansen's tutorial, without the
57+ rank-mu weights beyond the elite) for when the cma package is absent.
58+ Candidates are clipped to the unit cube."""
59+
60+ def __init__(self, x0, sigma, popsize, seed=1):
61+ self.n = len(x0)
62+ self.mean = np.array(x0, float)
63+ self.sigma = sigma
64+ self.lam = popsize
65+ self.mu = popsize // 2
66+ w = np.log(self.mu + 0.5) - np.log(np.arange(1, self.mu + 1))
67+ self.w = w / w.sum()
68+ self.mueff = 1 / np.sum(self.w ** 2)
69+ n = self.n
70+ self.cc = (4 + self.mueff / n) / (n + 4 + 2 * self.mueff / n)
71+ self.cs = (self.mueff + 2) / (n + self.mueff + 5)
72+ self.c1 = 2 / ((n + 1.3) ** 2 + self.mueff)
73+ self.cmu = min(1 - self.c1, 2 * (self.mueff - 2 + 1 / self.mueff) / ((n + 2) ** 2 + self.mueff))
74+ self.damps = 1 + 2 * max(0, math.sqrt((self.mueff - 1) / (n + 1)) - 1) + self.cs
75+ self.chiN = math.sqrt(n) * (1 - 1 / (4 * n) + 1 / (21 * n * n))
76+ self.pc = np.zeros(n); self.ps = np.zeros(n); self.C = np.eye(n)
77+ self.rng = np.random.default_rng(seed)
78+ self.gen = 0
79+
80+ def ask(self):
81+ self.B, self.D = None, None
82+ vals, vecs = np.linalg.eigh(self.C)
83+ self.D = np.sqrt(np.maximum(vals, 1e-20)); self.B = vecs
84+ self.z = self.rng.standard_normal((self.lam, self.n))
85+ self.y = (self.B @ (self.D[:, None] * self.z.T)).T
86+ X = np.clip(self.mean + self.sigma * self.y, 0, 1)
87+ return [x for x in X]
88+
89+ def tell(self, X, scores):
90+ order = np.argsort(scores)
91+ X = np.array(X)
92+ old = self.mean.copy()
93+ self.mean = self.w @ X[order[: self.mu]]
94+ y_w = (self.mean - old) / self.sigma
95+ C_inv_sqrt = self.B @ np.diag(1 / self.D) @ self.B.T
96+ self.ps = (1 - self.cs) * self.ps + math.sqrt(self.cs * (2 - self.cs) * self.mueff) * (C_inv_sqrt @ y_w)
97+ self.gen += 1
98+ hsig = np.linalg.norm(self.ps) / math.sqrt(1 - (1 - self.cs) ** (2 * self.gen)) / self.chiN < 1.4 + 2 / (self.n + 1)
99+ self.pc = (1 - self.cc) * self.pc + hsig * math.sqrt(self.cc * (2 - self.cc) * self.mueff) * y_w
100+ ys = (X[order[: self.mu]] - old) / self.sigma
101+ rank_mu = sum(wi * np.outer(yi, yi) for wi, yi in zip(self.w, ys))
102+ self.C = (1 - self.c1 - self.cmu) * self.C + self.c1 * (np.outer(self.pc, self.pc) + (1 - hsig) * self.cc * (2 - self.cc) * self.C) + self.cmu * rank_mu
103+ self.sigma *= math.exp((self.cs / self.damps) * (np.linalg.norm(self.ps) / self.chiN - 1))
104+
105+
106+FIXED = {} # name -> value for the parameters that are not free
107+
108+
109+def free_space():
110+ return [e for e in SPACE if e[0] in FREE]
111+
112+
113+def to_unit(vals):
114+ """Physical values of the free parameters -> [0, 1] coordinates."""
115+ u = []
116+ for (name, kind, lo, hi, _), v in zip(free_space(), vals):
117+ u.append((math.log(v / lo) / math.log(hi / lo)) if kind == "log" else (v - lo) / (hi - lo))
118+ return np.array(u)
119+
120+
121+def from_unit(u):
122+ vals = []
123+ for (name, kind, lo, hi, _), x in zip(free_space(), np.clip(u, 0, 1)):
124+ vals.append(lo * (hi / lo) ** x if kind == "log" else lo + (hi - lo) * x)
125+ return vals
126+
127+
128+def params_of(vals):
129+ """Full parameter dict (free values plus the fixed ones) with the shape as a feature string."""
130+ d = {name: d0 for name, _, _, _, d0 in SPACE}
131+ d.update(FIXED)
132+ d.update({name: v for (name, *_), v in zip(free_space(), vals)})
133+ bow, bow_y = d.pop("bow"), d.pop("bowY")
134+ feats = [f"bow:x,{bow:.6g},0"] if abs(bow) > 1e-5 else []
135+ if abs(bow_y) > 1e-5:
136+ feats.append(f"bow:y,{bow_y:.6g},0")
137+ d["features"] = ";".join(feats) if feats else "none"
138+ return d
139+
140+
141+def objective(ref, seg, tol=None):
142+ """Default: the mixed log-mel and fine-structure objective. With tol set: the tolerant line measure of linesim.py."""
143+ if tol is not None:
144+ d, parts = linesim.distance(linesim.picture(ref, tol), linesim.picture(seg, tol))
145+ return d, parts
146+ fe = features(seg)
147+ d_mel = float(np.mean(np.abs(ref[3] - fe[3])))
148+ d_early = 1 - float(ref[0] @ fe[0])
149+ d_whole = 1 - float(ref[1] @ fe[1])
150+ return d_mel + 4 * d_early + 4 * d_whole, (d_mel, d_early, d_whole)
151+
152+
153+def render(jobs, out, duration, workers):
154+ os.makedirs(out, exist_ok=True)
155+ jf = os.path.join(out, "jobs.json")
156+ with open(jf, "w") as f:
157+ json.dump({"jobs": jobs}, f)
158+ env = dict(os.environ)
159+ env.pop("DISPLAY", None)
160+ subprocess.run(NODE + ["--batch", jf, "--out", out, "--workers", str(workers), "--duration", str(duration)], check=True, env=env, cwd=ROOT,
161+ stdout=open(os.path.join(out, "render.log"), "a"), stderr=subprocess.STDOUT)
162+
163+
164+def main():
165+ ap = argparse.ArgumentParser()
166+ ap.add_argument("target")
167+ ap.add_argument("--out", required=True)
168+ ap.add_argument("--init", default="", help="comma-separated name=value in SI for the starting point")
169+ ap.add_argument("--sigma", type=float, default=0.15, help="initial CMA-ES step in the unit cube")
170+ ap.add_argument("--generations", type=int, default=12)
171+ ap.add_argument("--popsize", type=int, default=16)
172+ ap.add_argument("--workers", type=int, default=16)
173+ ap.add_argument("--duration", type=float, default=3.0)
174+ ap.add_argument("--free", default=",".join(FREE), help="comma-separated names of the free parameters")
175+ ap.add_argument("--set", default="", help="comma-separated name=value (SI) for fixed parameters that differ from the defaults")
176+ ap.add_argument("--linesim", action="store_true", help="use the tolerant line measure as the objective")
177+ ap.add_argument("--tol-start", type=float, default=0.05)
178+ ap.add_argument("--tol-end", type=float, default=0.02)
179+ args = ap.parse_args()
180+ FREE[:] = [n.strip() for n in args.free.split(",") if n.strip()]
181+ for kv in args.set.split(","):
182+ if "=" in kv:
183+ k, v = kv.split("="); FIXED[k.strip()] = float(v)
184+ os.makedirs(args.out, exist_ok=True)
185+ x = read_wav_any(args.target)
186+ st = onset(x)
187+ ref_seg = segment(x, st)
188+ ref = ref_seg if args.linesim else features(ref_seg)
189+ init = {name: d for name, _, _, _, d in SPACE}
190+ init.update(FIXED)
191+ for kv in args.init.split(","):
192+ if "=" in kv:
193+ k, v = kv.split("=")
194+ init[k.strip()] = float(v)
195+ x0 = to_unit([init[name] for name, *_ in free_space()])
196+ print(f"target {args.target}: onset {st / 48000:.3f} s; start {params_of(from_unit(x0))}", flush=True)
197+
198+ try:
199+ import cma
200+ es = cma.CMAEvolutionStrategy(x0, args.sigma, {"popsize": args.popsize, "bounds": [0, 1], "seed": 1, "verbose": -9})
201+ except ImportError:
202+ es = SimpleCMA(x0, args.sigma, args.popsize)
203+ log = open(os.path.join(args.out, "fit.jsonl"), "a")
204+ best_overall = (float("inf"), None, None)
205+ history = []
206+ t0 = time.time()
207+ for gen in range(args.generations):
208+ tol = args.tol_start + (args.tol_end - args.tol_start) * gen / max(1, args.generations - 1) if args.linesim else None
209+ X = es.ask()
210+ jobs = [{"id": f"g{gen:02d}-{k:02d}", "params": params_of(from_unit(u))} for k, u in enumerate(X)]
211+ gdir = os.path.join(args.out, f"gen{gen:02d}")
212+ render(jobs, gdir, args.duration, args.workers)
213+ scores = []
214+ for j, u in zip(jobs, X):
215+ p = os.path.join(gdir, "clips", j["id"] + ".wav")
216+ if not os.path.exists(p):
217+ scores.append(1e3); continue
218+ s, parts = objective(ref, segment(read_wav_any(p), 0), tol)
219+ scores.append(s)
220+ log.write(json.dumps({"gen": gen, "id": j["id"], "score": s, "tol": tol, "parts": parts, "params": j["params"]}) + "\n")
221+ log.flush()
222+ es.tell(X, scores)
223+ k = int(np.argmin(scores))
224+ best = (scores[k], jobs[k]["params"], os.path.join(gdir, "clips", jobs[k]["id"] + ".wav"))
225+ shutil.copy(best[2], os.path.join(args.out, f"best-{gen:02d}.wav"))
226+ if best[0] < best_overall[0]:
227+ best_overall = best
228+ history.append({"gen": gen, "best": best[0], "mean": float(np.mean(scores)), "params": best[1]})
229+ print(f"gen {gen}: best {best[0]:.3f} (mean {np.mean(scores):.3f}, overall {best_overall[0]:.3f}{f', tol {tol:.3f}' if tol else ''}) after {(time.time() - t0) / 60:.1f} min; "
230+ + " ".join(f"{n}={v:.4g}" for n, v in best[1].items() if n != "features") + f" {best[1]['features']}", flush=True)
231+ write_page(args, history, best_overall)
232+ print("best:", json.dumps(best_overall[1]))
233+
234+
235+def write_page(args, history, best_overall):
236+ rows = "".join(
237+ f"<tr><td>{h['gen']}</td><td>{h['best']:.3f}</td><td>{h['mean']:.3f}</td><td><audio controls preload='none' src='best-{h['gen']:02d}.wav'></audio></td>"
238+ + "".join(f"<td>{h['params'][n]:.4g}</td>" for n, *_ in SPACE if n not in ("bow", "bowY")) + f"<td class='l'>{h['params']['features']}</td></tr>" for h in history)
239+ html = f"""<!doctype html><meta charset="utf-8"><title>fit to {os.path.basename(args.target)}</title>
240+<style>body{{font:14px system-ui,sans-serif;margin:16px}} table{{border-collapse:collapse}} td,th{{padding:3px 8px;border-bottom:1px solid #ddd;text-align:right;white-space:nowrap}} td.l{{text-align:left}} audio{{height:28px;width:200px}}</style>
241+<h1>Simulation fitted to {os.path.basename(args.target)}</h1>
242+<p>Target: <audio controls src="target.wav"></audio> CMA-ES, population {args.popsize}, objective = log-mel distance + 4 (early fine-structure distance + whole-clip fine-structure distance). Best so far: {best_overall[0]:.3f}.</p>
243+<table><tr><th>gen</th><th>best</th><th>mean</th><th>best clip</th>{"".join(f"<th>{n}</th>" for n, *_ in SPACE if n not in ("bow", "bowY"))}<th class="l">shape</th></tr>{rows}</table>"""
244+ with open(os.path.join(args.out, "fit.html"), "w") as f:
245+ f.write(html)
246+
247+
248+if __name__ == "__main__":
249+ main()
ml/history.pyadded+12−0View file
@@ -0,0 +1,12 @@
1+"""Prints the error of each retraining run against the number of clips.
2+ python3 ml/history.py DATASET_DIR/ml/history.jsonl
3+"""
4+import json
5+import sys
6+
7+COLS = ["E", "rho", "h", "Lx", "Ly", "E/rho", "kappa", "area", "aspect"]
8+rows = [json.loads(l) for l in open(sys.argv[1]) if l.strip()]
9+print(f"{'time':>5} {'clips':>5} " + " ".join(f"{c:>11}" for c in COLS) + " (rel RMS % / R^2)")
10+for h in rows:
11+ m = h["metrics"]
12+ print(f"{h['finished'][11:16]:>5} {h['n_clips']:>5} " + " ".join(f"{m[c]['rel_rms_pct']:5.1f}/{m[c]['r2']:5.2f}" for c in COLS))
ml/linesim.pyadded+129−0View file
@@ -0,0 +1,129 @@
1+"""A similarity score for struck-sheet sounds that compares the partials with a
2+tolerance in frequency, and their levels and decays, without requiring the
3+lines to coincide.
4+
5+Representation. Each clip is cut into K time slices (0-0.3, 0.3-1, 1-2, 2-3 s
6+by default). For each slice the mean power spectrum (8192-point windows,
7+5.9 Hz bins) is mapped onto a log-frequency axis (96 bins per octave, 30 Hz to
8+5 kHz) and smoothed there with a Gaussian whose width is the frequency
9+tolerance (2 percent by default), then converted to dB relative to the clip's
10+loudest value and floored. A line thus becomes a bump of width `tol`: two
11+lines a fraction of `tol` apart overlap almost completely, lines several `tol`
12+apart do not overlap at all, and in between the overlap falls off smoothly.
13+The level of a line is the height of its bump, and its decay is how the bump
14+shrinks from slice to slice, so the level, the position, and the decay of
15+every partial are all in the K x n_bins picture.
16+
17+Score. The mean absolute difference in dB between the two pictures, over the
18+cells where either clip is above the floor (so silent regions do not dilute
19+it), with the later slices weighted equally to the first. Zero for identical
20+clips; smooth in the parameters of a simulation; insensitive to the level and
21+(through the log axis) treating a 2 percent shift at 100 Hz like one at 1 kHz.
22+
23+ python3 ml/linesim.py TARGET.wav CLIP.wav [CLIP.wav ...] [--tol 0.02]
24+ python3 ml/linesim.py TARGET.wav --data DATASET_DIR --top 10 # rank a dataset
25+"""
26+
27+import argparse
28+import json
29+import os
30+import sys
31+
32+import numpy as np
33+
34+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
35+from match import read_wav_any, onset, segment, FS # noqa: E402
36+
37+N_FFT = 8192
38+HOP = 1024
39+F_LO, F_HI = 30.0, 5000.0
40+BINS_PER_OCTAVE = 96
41+SLICES = ((0, 0.3), (0.3, 1.0), (1.0, 2.0), (2.0, 3.0))
42+FLOOR = -70.0
43+
44+_fb_cache = {}
45+
46+
47+def logf_axis():
48+ n = int(np.log2(F_HI / F_LO) * BINS_PER_OCTAVE)
49+ return F_LO * 2 ** (np.arange(n) / BINS_PER_OCTAVE)
50+
51+
52+def smoothing_matrix(tol):
53+ """(n_log, n_lin): maps a linear power spectrum onto the log axis and smooths it with a Gaussian of relative width tol."""
54+ if tol in _fb_cache:
55+ return _fb_cache[tol]
56+ f_lin = np.arange(N_FFT // 2 + 1) * FS / N_FFT
57+ centers = logf_axis()
58+ W = np.exp(-0.5 * ((np.log(f_lin[None, 1:] / centers[:, None])) / tol) ** 2)
59+ W[:, :] /= np.maximum(W.sum(1, keepdims=True), 1e-12)
60+ W = np.concatenate([np.zeros((len(centers), 1)), W], 1).astype(np.float32)
61+ _fb_cache[tol] = W
62+ return W
63+
64+
65+def picture(seg, tol=0.02):
66+ """(K, n_log) dB picture of the clip; see the module docstring."""
67+ win = np.hanning(N_FFT + 1)[:-1].astype(np.float32)
68+ n_frames = 1 + (len(seg) - N_FFT) // HOP
69+ frames = np.lib.stride_tricks.as_strided(seg, shape=(n_frames, N_FFT), strides=(seg.strides[0] * HOP, seg.strides[0]))
70+ P = np.abs(np.fft.rfft(frames * win, axis=1)) ** 2
71+ fps = FS / HOP
72+ W = smoothing_matrix(tol)
73+ out = []
74+ for a, b in SLICES:
75+ sl = P[int(a * fps): max(int(a * fps) + 1, int(b * fps))].mean(0)
76+ out.append(W @ sl)
77+ db = 10 * np.log10(np.array(out) + 1e-20)
78+ db -= db.max()
79+ return np.maximum(db, FLOOR)
80+
81+
82+def distance(A, B):
83+ """Mean |dB difference| over the cells where either picture is above the floor, averaged over slices."""
84+ per_slice = []
85+ for a, b in zip(A, B):
86+ m = (a > FLOOR + 1) | (b > FLOOR + 1)
87+ per_slice.append(float(np.abs(a - b)[m].mean()) if m.any() else 0.0)
88+ return float(np.mean(per_slice)), per_slice
89+
90+
91+def load(path, align):
92+ x = read_wav_any(path)
93+ return segment(x, onset(x) if align else 0)
94+
95+
96+def main():
97+ ap = argparse.ArgumentParser()
98+ ap.add_argument("target")
99+ ap.add_argument("clips", nargs="*")
100+ ap.add_argument("--data", default=None)
101+ ap.add_argument("--top", type=int, default=10)
102+ ap.add_argument("--tol", type=float, default=0.02)
103+ ap.add_argument("--no-align", action="store_true", help="do not detect the onset of the target (simulated clips start at the impact)")
104+ args = ap.parse_args()
105+ T = picture(load(args.target, not args.no_align), args.tol)
106+ if args.clips:
107+ for c in args.clips:
108+ d, parts = distance(T, picture(load(c, False), args.tol))
109+ print(f"{os.path.basename(c):>16}: {d:6.2f} dB slices " + " ".join(f"{p:.2f}" for p in parts))
110+ if args.data:
111+ rows = {json.loads(l)["id"]: json.loads(l) for l in open(os.path.join(args.data, "index.jsonl")) if l.strip()}
112+ ids = sorted(i for i in rows if os.path.exists(os.path.join(args.data, "clips", i + ".wav")))
113+ cache = os.path.join(args.data, "ml", f"linesim-tol{args.tol:g}.npz")
114+ if os.path.exists(cache) and list(np.load(cache)["ids"]) == ids:
115+ pics = np.load(cache)["pics"]
116+ else:
117+ pics = np.stack([picture(load(os.path.join(args.data, "clips", i + ".wav"), False), args.tol) for i in ids]).astype(np.float32)
118+ np.savez(cache, ids=np.array(ids), pics=pics)
119+ d = np.array([distance(T, p)[0] for p in pics])
120+ order = np.argsort(d)
121+ print(f"{len(ids)} clips; median distance {np.median(d):.2f} dB")
122+ for k in order[: args.top]:
123+ r = rows[ids[k]]
124+ print(f"{ids[k]} {d[k]:6.2f} dB | E {r['E'] / 1e9:.0f} rho {r['rho']:.0f} h {r['h'] * 1e3:.2f} {r['Lx'] / 0.0254:.1f}x{r['Ly'] / 0.0254:.1f} in v {r['strikeVelocity']:.1f} sig0 {r['sigma0']:.2f} sig1 {r['sigma1'] * 1e4:.1f} p2p {r['peakToPeak'] * 1e3:.1f}")
125+ return [ids[k] for k in order[: args.top]]
126+
127+
128+if __name__ == "__main__":
129+ main()
ml/match.pyadded+170−0View file
@@ -0,0 +1,170 @@
1+"""Finds the simulated clips closest to a real recording.
2+
3+ python3 ml/match.py REAL.wav [--data DATASET_DIR] [--top 10]
4+
5+The recording is converted to mono at 48 kHz, its onset is found (the first
6+sample of the impact transient), and the 3 s from the onset are compared with
7+every simulated clip, which starts at the impact. Three distances are used,
8+each insensitive to the recording level and to the broad coloring of the
9+microphone and the room:
10+
11+1. the fine structure of the mean log spectrum over the first second (the
12+ log spectrum minus its smooth envelope, over 30 Hz to 5 kHz), which places
13+ the partials;
14+2. the same over the whole clip;
15+3. the decay: the log energy in 24 bands (100 Hz to 5 kHz) at 12 times over
16+ the clip, each band normalized to its own start, which captures how fast
17+ each region of the spectrum dies away.
18+
19+The ranking uses a fourth, perceptual distance: the mean absolute difference
20+in dB between log-mel-like spectrograms (48 log-spaced bands from 60 Hz to
21+6 kHz on 100 ms frames, relative to each clip's maximum, floored at -60 dB),
22+which weighs timbre and decay as they are heard; the three distances above
23+are reported alongside.
24+"""
25+
26+import argparse
27+import json
28+import os
29+import sys
30+import wave
31+
32+import numpy as np
33+from scipy.ndimage import median_filter
34+from scipy.signal import resample_poly
35+
36+FS = 48000
37+N_FFT = 8192
38+HOP = 1024
39+DUR = 3.0
40+
41+
42+def read_wav_any(path):
43+ with wave.open(path) as w:
44+ ch, sw, fs, n = w.getnchannels(), w.getsampwidth(), w.getframerate(), w.getnframes()
45+ raw = w.readframes(n)
46+ if sw == 2:
47+ x = np.frombuffer(raw, dtype="<i2").astype(np.float32) / 32768
48+ elif sw == 3:
49+ b = np.frombuffer(raw, dtype=np.uint8).reshape(-1, 3)
50+ x = ((b[:, 0].astype(np.int32) | (b[:, 1].astype(np.int32) << 8) | (b[:, 2].astype(np.int32) << 16)) << 8 >> 8).astype(np.float32) / 2 ** 23
51+ elif sw == 4:
52+ x = np.frombuffer(raw, dtype="<i4").astype(np.float32) / 2 ** 31
53+ else:
54+ raise ValueError(f"unsupported sample width {sw}")
55+ if ch > 1:
56+ x = x.reshape(-1, ch).mean(1)
57+ if fs != FS:
58+ from math import gcd
59+ g = gcd(fs, FS)
60+ x = resample_poly(x, FS // g, fs // g).astype(np.float32)
61+ return x
62+
63+
64+def onset(x, frac=0.1, pre=0.002):
65+ """First sample where the short-time envelope reaches `frac` of its peak, less `pre` seconds."""
66+ env = np.abs(x)
67+ k = int(0.001 * FS)
68+ env = np.convolve(env, np.ones(k) / k, mode="same")
69+ i = int(np.argmax(env >= frac * env.max()))
70+ return max(0, i - int(pre * FS))
71+
72+
73+def segment(x, start):
74+ n = int(DUR * FS)
75+ seg = np.zeros(n, np.float32)
76+ part = x[start:start + n]
77+ seg[: len(part)] = part
78+ return seg / (np.abs(seg).max() + 1e-9)
79+
80+
81+def stft_db(seg):
82+ win = np.hanning(N_FFT + 1)[:-1].astype(np.float32)
83+ n_frames = 1 + (len(seg) - N_FFT) // HOP
84+ frames = np.lib.stride_tricks.as_strided(seg, shape=(n_frames, N_FFT), strides=(seg.strides[0] * HOP, seg.strides[0]))
85+ return 20 * np.log10(np.abs(np.fft.rfft(frames * win, axis=1)) + 1e-9) # (frames, bins)
86+
87+
88+def melspec(seg, n_bands=48, f_lo=60.0, f_hi=6000.0, frame=0.1, floor=-60.0):
89+ """Log-mel-like spectrogram (log-spaced bands) in dB relative to the clip's
90+ maximum, floored, on 100 ms frames: what the ear hears of timbre and decay."""
91+ n = int(frame * FS)
92+ m = len(seg) // n
93+ frames = seg[: m * n].reshape(m, n) * np.hanning(n)
94+ P = np.abs(np.fft.rfft(frames, axis=1)) ** 2
95+ f = np.fft.rfftfreq(n, 1 / FS)
96+ edges = np.geomspace(f_lo, f_hi, n_bands + 1)
97+ B = np.stack([P[:, (f >= a) & (f < b)].sum(1) for a, b in zip(edges[:-1], edges[1:])], 1)
98+ db = 10 * np.log10(B + 1e-18)
99+ db -= db.max()
100+ return np.maximum(db, floor)
101+
102+
103+def features(seg):
104+ S = stft_db(seg)
105+ f = np.fft.rfftfreq(N_FFT, 1 / FS)
106+ band = (f >= 30) & (f <= 5000)
107+ fps = FS / HOP
108+ def fine(mean_db):
109+ env = median_filter(mean_db, size=int(200 / (FS / N_FFT)) | 1, mode="nearest")
110+ v = (mean_db - env)[band]
111+ return v / (np.linalg.norm(v) + 1e-9)
112+ early = fine(S[: int(1.0 * fps)].mean(0))
113+ whole = fine(S.mean(0))
114+ # Decay profile: 24 log-spaced bands, 12 time steps, each band relative to its first step.
115+ edges = np.geomspace(100, 5000, 25)
116+ P = 10 ** (S / 10)
117+ bands = np.stack([P[:, (f >= a) & (f < b)].sum(1) for a, b in zip(edges[:-1], edges[1:])], 1) # (frames, 24)
118+ steps = np.array_split(np.arange(S.shape[0]), 12)
119+ prof = np.stack([bands[s].mean(0) for s in steps]) # (12, 24)
120+ prof = 10 * np.log10(prof + 1e-12)
121+ prof = prof - prof[0]
122+ return early, whole, prof.ravel(), melspec(seg)
123+
124+
125+def main():
126+ ap = argparse.ArgumentParser()
127+ ap.add_argument("real")
128+ ap.add_argument("--data", default="dataset")
129+ ap.add_argument("--top", type=int, default=10)
130+ ap.add_argument("--onset-frac", type=float, default=0.1)
131+ args = ap.parse_args()
132+ x = read_wav_any(args.real)
133+ st = onset(x, args.onset_frac)
134+ print(f"{args.real}: {len(x) / FS:.2f} s, onset at {st / FS:.3f} s (peak {np.abs(x).max():.3f})")
135+ ref = features(segment(x, st))
136+ rows = {}
137+ with open(os.path.join(args.data, "index.jsonl")) as fh:
138+ for line in fh:
139+ if line.strip():
140+ r = json.loads(line)
141+ rows[r["id"]] = r
142+ ids, D = [], []
143+ for id_ in sorted(rows):
144+ p = os.path.join(args.data, "clips", id_ + ".wav")
145+ if not os.path.exists(p):
146+ continue
147+ seg = segment(read_wav_any(p), 0)
148+ fe = features(seg)
149+ d1 = 1 - float(ref[0] @ fe[0])
150+ d2 = 1 - float(ref[1] @ fe[1])
151+ d3 = float(np.sqrt(np.mean((ref[2] - fe[2]) ** 2)))
152+ d4 = float(np.mean(np.abs(ref[3] - fe[3])))
153+ ids.append(id_); D.append((d1, d2, d3, d4))
154+ D = np.array(D)
155+ # Main score: the mean absolute dB difference of the log-mel spectrograms,
156+ # which weighs timbre and decay the way they are heard; the other three
157+ # distances are reported for reference.
158+ score = D[:, 3]
159+ order = np.argsort(score, kind="stable")
160+ print(f"{len(ids)} simulated clips compared\n")
161+ print(f"{'id':>6} {'mel dB':>6} {'d_early':>7} {'d_whole':>7} {'d_decay':>7} | {'E GPa':>5} {'rho':>5} {'h mm':>5} {'W x H in':>11} {'v':>4} {'strike':>10} {'sig0':>5} {'sig1e4':>6} {'p2p mm':>6} features")
162+ for i in order[: args.top]:
163+ r = rows[ids[i]]
164+ print(f"{ids[i]:>6} {score[i]:6.2f} {D[i, 0]:7.3f} {D[i, 1]:7.3f} {D[i, 2]:7.2f} | {r['E'] / 1e9:5.0f} {r['rho']:5.0f} {r['h'] * 1e3:5.2f} {r['Lx'] / 0.0254:5.1f} x {r['Ly'] / 0.0254:4.1f} {r['strikeVelocity']:4.1f} {r['strikeFracX']:4.2f},{r['strikeFracY']:4.2f} {r['sigma0']:5.2f} {r['sigma1'] * 1e4:6.1f} {r['peakToPeak'] * 1e3:6.1f} {r['features'][:70]}")
165+ print(f"\nmedian over all clips: mel {np.median(D[:, 3]):.2f} dB, d_early {np.median(D[:, 0]):.3f}, d_whole {np.median(D[:, 1]):.3f}, d_decay {np.median(D[:, 2]):.2f}")
166+ return [ids[i] for i in order[: args.top]]
167+
168+
169+if __name__ == "__main__":
170+ main()
ml/predict.pyadded+66−0View file
@@ -0,0 +1,66 @@
1+"""Applies a trained model to recordings.
2+
3+ python3 ml/predict.py --model RUN_DIR/model.pt [--logf] REAL.wav [more.wav ...]
4+
5+Each recording is converted to mono at 48 kHz, cut to 3 s from its onset (see
6+match.py), peak-normalized, turned into the spectrogram the model was trained
7+on, and passed through the network. To show the sensitivity to the alignment,
8+the prediction is repeated with the onset shifted by a few milliseconds and
9+the spread is reported alongside the value at the detected onset.
10+"""
11+
12+import argparse
13+import math
14+import os
15+import sys
16+
17+import numpy as np
18+import torch
19+
20+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
21+from features import spectrogram, spectrogram_logf # noqa: E402
22+from match import read_wav_any, onset, segment, FS # noqa: E402
23+from train import Net, SpecNet, TARGETS, UNITS # noqa: E402
24+
25+torch.backends.cudnn.enabled = False
26+
27+
28+def main():
29+ ap = argparse.ArgumentParser()
30+ ap.add_argument("--model", required=True)
31+ ap.add_argument("--logf", action="store_true")
32+ ap.add_argument("--shifts", default="0,-5,5,10,20", help="onset shifts in ms for the sensitivity check")
33+ ap.add_argument("wavs", nargs="+")
34+ args = ap.parse_args()
35+ ck = torch.load(args.model, map_location="cpu", weights_only=False)
36+ net = SpecNet(len(ck["targets"]), n_freq=ck["n_freq"]) if ck.get("arch") == "spec1d" else Net(len(ck["targets"]))
37+ net.load_state_dict(ck["state"])
38+ net.eval()
39+ mu, sd = ck["mu"], ck["sd"]
40+ Tinv = np.linalg.inv(ck["T"]) if "T" in ck else np.eye(len(ck["targets"]))
41+ logf = ck["input"] == "logf" if "input" in ck else args.logf
42+ spec = spectrogram_logf if logf else spectrogram
43+ shifts = [float(s) for s in args.shifts.split(",")]
44+ for path in args.wavs:
45+ x = read_wav_any(path)
46+ st = onset(x)
47+ preds = []
48+ for sh in shifts:
49+ s0 = max(0, st + int(sh * FS / 1000))
50+ X = torch.from_numpy(spec(segment(x, s0)).astype(np.float32))[None, None]
51+ with torch.no_grad():
52+ z = net(X)[0].numpy()
53+ preds.append((z * sd + mu) @ Tinv.T)
54+ preds = np.array(preds)
55+ p0, lo, hi = preds[0], preds.min(0), preds.max(0)
56+ print(f"{path}: {len(x) / FS:.2f} s, onset at {st / FS:.3f} s")
57+ for i, t in enumerate(TARGETS):
58+ unit, scale = UNITS[t]
59+ print(f" {t:>4}: {math.exp(p0[i]) * scale:8.2f} {unit:<7} (onset shifts give {math.exp(lo[i]) * scale:.2f} to {math.exp(hi[i]) * scale:.2f})")
60+ d = {t: p0[i] for i, t in enumerate(TARGETS)}
61+ kappa = math.exp(d["h"] + 0.5 * (d["E"] - d["rho"])) / math.sqrt(12 * (1 - 0.29 ** 2))
62+ print(f" E/rho {math.exp(d['E'] - d['rho']) / 1e6:.1f} MPa m^3/kg, kappa {kappa:.3f} m^2/s, area {math.exp(d['Lx'] + d['Ly']) / 0.0254 ** 2:.0f} in^2, aspect {math.exp(d['Lx'] - d['Ly']):.2f}")
63+
64+
65+if __name__ == "__main__":
66+ main()
ml/retrain-loop.shadded+17−0View file
@@ -0,0 +1,17 @@
1+#!/bin/bash
2+# Retrains on the current dataset every INTERVAL seconds while clips accumulate.
3+# ml/retrain-loop.sh DATASET_DIR [interval s] [epochs] [tag [extra train.py arguments, e.g. --logf --arch spec1d]]
4+# The tag names the run directories and the report and history files (report-TAG.html).
5+set -u
6+DATA=$1
7+INTERVAL=${2:-3600}
8+EPOCHS=${3:-60}
9+TAG=${4:-}
10+ML="$DATA/ml"
11+mkdir -p "$ML"
12+while true; do
13+ RUN="$ML/run-$(date +%Y%m%d-%H%M%S)${TAG:+-$TAG}"
14+ echo "=== $(date) training into $RUN"
15+ python3 "$(dirname "$0")/train.py" --data "$DATA" --out "$RUN" --epochs "$EPOCHS" ${TAG:+--tag "$TAG"} "${@:5}" 2>&1 | tail -12
16+ sleep "$INTERVAL"
17+done
ml/similar.pyadded+84−0View file
@@ -0,0 +1,84 @@
1+"""Lists the simulated clips most similar to a chosen simulated clip under the
2+perceptual log-mel distance of match.py, and writes a listening page.
3+
4+ python3 ml/similar.py QUERY_ID [--data dataset] [--top 5] [--controls 50,500,2500]
5+
6+The log-mel features of every clip are cached in DATA/ml/mel.npz. The page
7+DATA/similar-QUERY_ID.html has the query, the top matches, and a few control
8+clips at given ranks, so the measure can be judged by ear against the ranks.
9+"""
10+
11+import argparse
12+import json
13+import os
14+import sys
15+
16+import numpy as np
17+
18+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
19+from match import read_wav_any, segment, melspec # noqa: E402
20+
21+
22+def mel_features(data):
23+ cache = os.path.join(data, "ml", "mel.npz")
24+ rows = {json.loads(l)["id"]: json.loads(l) for l in open(os.path.join(data, "index.jsonl")) if l.strip()}
25+ ids = sorted(i for i in rows if os.path.exists(os.path.join(data, "clips", i + ".wav")))
26+ if os.path.exists(cache):
27+ z = np.load(cache, allow_pickle=True)
28+ if list(z["ids"]) == ids:
29+ return rows, ids, z["mel"]
30+ os.makedirs(os.path.join(data, "ml"), exist_ok=True)
31+ mel = np.stack([melspec(segment(read_wav_any(os.path.join(data, "clips", i + ".wav")), 0)) for i in ids]).astype(np.float32)
32+ np.savez(cache, ids=np.array(ids), mel=mel)
33+ return rows, ids, mel
34+
35+
36+def cells(r):
37+ return (f"<td>{r['E'] / 1e9:.0f}</td><td>{r['rho']:.0f}</td><td>{r['h'] * 1e3:.2f}</td><td>{r['Lx'] / 0.0254:.1f}×{r['Ly'] / 0.0254:.1f}</td>"
38+ f"<td>{r['strikeVelocity']:.1f}</td><td>{r['strikeFracX']:.2f}, {r['strikeFracY']:.2f}</td><td>{r['sigma0']:.2f}</td><td>{r['sigma1'] * 1e4:.1f}</td>"
39+ f"<td>{r['peakToPeak'] * 1e3:.1f}</td><td class='l'><div class='feat' title=\"{r['features']}\">{r['features']}</div></td>")
40+
41+
42+def main():
43+ ap = argparse.ArgumentParser()
44+ ap.add_argument("query")
45+ ap.add_argument("--data", default="dataset")
46+ ap.add_argument("--top", type=int, default=5)
47+ ap.add_argument("--controls", default="50,500,2500")
48+ args = ap.parse_args()
49+ rows, ids, mel = mel_features(args.data)
50+ q = ids.index(args.query)
51+ d = np.abs(mel - mel[q]).mean((1, 2))
52+ order = np.argsort(d)
53+ assert order[0] == q
54+ ranked = [(k, ids[i], d[i]) for k, i in enumerate(order)]
55+ sel = [("query", *ranked[0][1:])] + [(f"rank {k}", i, dd) for k, i, dd in ranked[1:args.top + 1]]
56+ for c in [int(x) for x in args.controls.split(",") if x]:
57+ if c < len(ranked):
58+ sel.append((f"control, rank {c}", *ranked[c][1:]))
59+ print(f"query {args.query}; median distance over all clips {np.median(d):.2f} dB")
60+ for label, i, dd in sel:
61+ r = rows[i]
62+ print(f"{label:>18} {i} {dd:5.2f} dB | E {r['E'] / 1e9:.0f} rho {r['rho']:.0f} h {r['h'] * 1e3:.2f} {r['Lx'] / 0.0254:.1f}x{r['Ly'] / 0.0254:.1f} in v {r['strikeVelocity']:.1f} sig0 {r['sigma0']:.2f} sig1 {r['sigma1'] * 1e4:.1f} p2p {r['peakToPeak'] * 1e3:.1f}")
63+ gal = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "scripts", "gallery.html")).read()
64+ style = gal[gal.index("<style>") + 7:gal.index("</style>")]
65+ js = gal[gal.index("// Spectrogram:"):gal.index("$(\"reload\").onclick")]
66+ body = "".join(
67+ f"<tr><td class='l'>{label}</td><td class='l'>{i}</td><td>{dd:.2f}</td><td class='l'><audio controls preload='none' src='clips/{i}.wav'></audio></td>{cells(rows[i])}</tr>"
68+ f"<tr><td colspan='14' class='l'><canvas width='900' height='120' data-src='clips/{i}.wav'></canvas></td></tr>"
69+ for label, i, dd in sel)
70+ html = f"""<!doctype html><meta charset="utf-8"><title>clips similar to {args.query}</title><style>{style}</style>
71+<h1>Simulated clips most similar to {args.query}</h1>
72+<p>Distance: mean absolute dB difference between level-normalized log-mel spectrograms (48 bands, 60 Hz to 6 kHz, 100 ms frames, floored at −60 dB). Median distance over all {len(ids)} clips: {np.median(d):.2f} dB. The control rows are clips at the given ranks, for comparison by ear.</p>
73+<table><thead><tr><th class="l">role</th><th class="l">id</th><th>dist dB</th><th class="l">audio</th><th>E GPa</th><th>ρ</th><th>h mm</th><th>W×H in</th><th>v</th><th>strike</th><th>σ0</th><th>σ1·10⁴</th><th>p-p mm</th><th class="l">features</th></tr></thead><tbody>{body}</tbody></table>
74+<script>const $ = (id) => document.getElementById(id);
75+{js}
76+for (const c of document.querySelectorAll("canvas")) drawSpec(c, c.dataset.src);
77+</script>"""
78+ out = os.path.join(args.data, f"similar-{args.query}.html")
79+ open(out, "w").write(html)
80+ print("wrote", out)
81+
82+
83+if __name__ == "__main__":
84+ main()
ml/train.pyadded+287−0View file
@@ -0,0 +1,287 @@
1+"""Trains a convolutional network that predicts the sheet parameters from the
2+spectrogram of a strike, evaluates it on held-out clips, and writes a report.
3+
4+ python3 ml/train.py --data DATASET_DIR --out RUN_DIR [--cache CACHE_DIR] [--epochs N] [--val-every K] [--device cuda|cpu]
5+
6+Targets are the logarithms of E, rho, h, Lx, Ly, standardized by the training
7+set; the loss is the mean squared error over the five. Every K-th clip by
8+index is held out for validation. The report gives, per target and for a few
9+derived quantities (E/rho, the bending wave constant kappa, the area, the
10+aspect ratio), the relative RMS error of the prediction against that of
11+predicting the training mean, and scatter plots. A line is appended to
12+RUN_DIR/../history.jsonl so that the error can be tracked as clips accumulate.
13+"""
14+
15+import argparse
16+import base64
17+import io
18+import json
19+import math
20+import os
21+import sys
22+import time
23+
24+import numpy as np
25+import torch
26+import torch.nn as nn
27+
28+# The office machine's PyTorch build cannot find a cuDNN engine for the
29+# convolution backward pass; the native kernels work, more slowly.
30+torch.backends.cudnn.enabled = os.environ.get("CUDNN", "0") == "1"
31+
32+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
33+from features import build_cache # noqa: E402
34+
35+TARGETS = ["E", "rho", "h", "Lx", "Ly"]
36+UNITS = {"E": ("GPa", 1e-9), "rho": ("kg/m^3", 1), "h": ("mm", 1e3), "Lx": ("in", 1 / 0.0254), "Ly": ("in", 1 / 0.0254)}
37+DERIVED = {
38+ "E/rho": lambda d: d["E"] - d["rho"],
39+ "kappa": lambda d: d["h"] + 0.5 * (d["E"] - d["rho"]),
40+ "area": lambda d: d["Lx"] + d["Ly"],
41+ "aspect": lambda d: d["Lx"] - d["Ly"],
42+ "kappa/area": lambda d: d["h"] + 0.5 * (d["E"] - d["rho"]) - d["Lx"] - d["Ly"],
43+}
44+
45+
46+class Net(nn.Module):
47+ """Six strided conv blocks, global average pooling, and a small head. A
48+ second input channel carries the frequency coordinate of each row, so the
49+ network can read absolute frequency despite the pooling."""
50+
51+ def __init__(self, n_out=5, channels=(16, 32, 64, 96, 128, 128)):
52+ super().__init__()
53+ layers, c_in = [], 2
54+ for c in channels:
55+ layers += [nn.Conv2d(c_in, c, 3, padding=1), nn.BatchNorm2d(c), nn.ReLU(inplace=True),
56+ nn.Conv2d(c, c, 3, padding=1), nn.BatchNorm2d(c), nn.ReLU(inplace=True), nn.MaxPool2d(2)]
57+ c_in = c
58+ self.features = nn.Sequential(*layers)
59+ self.head = nn.Sequential(nn.Linear(c_in, 128), nn.ReLU(inplace=True), nn.Dropout(0.2), nn.Linear(128, n_out))
60+
61+ def forward(self, x):
62+ b, _, f, t = x.shape
63+ coord = torch.linspace(0, 1, f, device=x.device).view(1, 1, f, 1).expand(b, 1, f, t)
64+ h = self.features(torch.cat([x, coord], 1)).mean((2, 3))
65+ return self.head(h)
66+
67+
68+class SpecNet(nn.Module):
69+ """Full-frequency-resolution model. The spectrogram is collapsed into a few
70+ time slices (mean log spectrum over 0-0.25, 0.25-1, 1-2, 2-3 s), which keep
71+ the partials at the resolution of the bins and the decay of each band as
72+ the ratio between slices. A 1D convolutional stack along frequency (two
73+ stride-2 steps only, dilations for context) is followed by multi-head
74+ attention pooling over frequency positions with a sinusoidal position code,
75+ so that the pooled vector can express where the lines are."""
76+
77+ def __init__(self, n_out=5, n_freq=1024, fps=48000 / 1024, slices=((0, 0.25), (0.25, 1), (1, 2), (2, 3)), heads=4, pos_dim=16):
78+ super().__init__()
79+ self.bounds = [(int(a * fps), int(b * fps)) for a, b in slices]
80+ c_in = len(slices) + 1
81+ def block(ci, co, k, stride=1, dil=1):
82+ return [nn.Conv1d(ci, co, k, stride=stride, padding=dil * (k // 2), dilation=dil), nn.BatchNorm1d(co), nn.ReLU(inplace=True)]
83+ self.conv = nn.Sequential(*block(c_in, 32, 7), *block(32, 32, 7), *block(32, 64, 5, stride=2), *block(64, 64, 5, dil=2),
84+ *block(64, 96, 5, stride=2), *block(96, 96, 5, dil=4), *block(96, 96, 5, dil=8), *block(96, 128, 3))
85+ n_pos = (n_freq + 3) // 4
86+ pos = torch.arange(n_pos, dtype=torch.float32)[:, None] / n_pos
87+ freqs = 2 ** torch.arange(pos_dim // 2, dtype=torch.float32)
88+ self.register_buffer("pos", torch.cat([torch.sin(math.pi * pos * freqs), torch.cos(math.pi * pos * freqs)], 1))
89+ d = 128 + pos_dim
90+ self.att = nn.Linear(d, heads)
91+ self.head = nn.Sequential(nn.Linear(heads * d, 256), nn.ReLU(inplace=True), nn.Dropout(0.2), nn.Linear(256, n_out))
92+
93+ def forward(self, x):
94+ b, _, f, t = x.shape
95+ s = torch.stack([x[:, 0, :, a:min(b_, t)].mean(-1) for a, b_ in self.bounds], 1)
96+ coord = torch.linspace(0, 1, f, device=x.device).view(1, 1, f).expand(b, 1, f)
97+ h = self.conv(torch.cat([s, coord], 1)).transpose(1, 2) # (b, positions, 128)
98+ hp = torch.cat([h, self.pos[: h.shape[1]].expand(b, -1, -1)], -1)
99+ a = torch.softmax(self.att(hp), dim=1) # (b, positions, heads)
100+ pooled = torch.einsum("bpk,bpd->bkd", a, hp).flatten(1)
101+ return self.head(pooled)
102+
103+
104+# Alternative target basis: the combinations the linear physics identifies
105+# (kappa / area up to a constant, aspect, E / rho, area) plus rho, so that the
106+# well-determined directions get their own outputs and clean gradients.
107+COMBOS = np.array([
108+ [0.5, -0.5, 1, -1, -1], # log kappa/area (up to a constant)
109+ [0, 0, 0, 1, -1], # log aspect
110+ [1, -1, 0, 0, 0], # log E/rho
111+ [0, 0, 0, 1, 1], # log area
112+ [0, 1, 0, 0, 0], # log rho
113+], dtype=np.float32)
114+
115+
116+def metrics(pred, true, mean):
117+ """Relative RMS error (percent) of pred and of the constant mean, and R^2, all in log space."""
118+ err = np.sqrt(np.mean((pred - true) ** 2))
119+ base = np.sqrt(np.mean((mean - true) ** 2))
120+ return {"rel_rms_pct": float(100 * err), "baseline_rel_rms_pct": float(100 * base), "r2": float(1 - err ** 2 / max(base ** 2, 1e-12)), "mae_pct": float(100 * np.mean(np.abs(pred - true)))}
121+
122+
123+def scatter_png(true, pred, name, unit, scale):
124+ import matplotlib
125+ matplotlib.use("Agg")
126+ import matplotlib.pyplot as plt
127+ fig, ax = plt.subplots(figsize=(3.2, 3.2), dpi=100)
128+ t, p = np.exp(true) * scale, np.exp(pred) * scale
129+ ax.scatter(t, p, s=6, alpha=0.5)
130+ lo, hi = min(t.min(), p.min()), max(t.max(), p.max())
131+ ax.plot([lo, hi], [lo, hi], "k-", lw=0.8)
132+ ax.set_xlabel(f"true {name} [{unit}]")
133+ ax.set_ylabel(f"predicted {name} [{unit}]")
134+ ax.set_xscale("log"); ax.set_yscale("log")
135+ fig.tight_layout()
136+ buf = io.BytesIO(); fig.savefig(buf, format="png"); plt.close(fig)
137+ return base64.b64encode(buf.getvalue()).decode()
138+
139+
140+def main():
141+ ap = argparse.ArgumentParser()
142+ ap.add_argument("--data", required=True)
143+ ap.add_argument("--out", required=True)
144+ ap.add_argument("--cache", default=None)
145+ ap.add_argument("--epochs", type=int, default=60)
146+ ap.add_argument("--batch", type=int, default=32)
147+ ap.add_argument("--lr", type=float, default=1e-3)
148+ ap.add_argument("--val-every", type=int, default=8, help="every k-th clip by index is held out")
149+ ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu")
150+ ap.add_argument("--workers", type=int, default=8)
151+ ap.add_argument("--logf", action="store_true", help="log-frequency spectrogram input")
152+ ap.add_argument("--no-augment", action="store_true", help="disable the onset-shift and masking augmentation")
153+ ap.add_argument("--arch", default="cnn2d", choices=["cnn2d", "spec1d"], help="pooling 2D CNN (default) or the full-resolution SpecNet")
154+ ap.add_argument("--targets", default="params", choices=["params", "combos"], help="train on log E, rho, h, Lx, Ly or on the identifiable combinations")
155+ ap.add_argument("--tag", default=None, help="suffix of the report and history files (default: 'logf' with --logf, else none)")
156+ args = ap.parse_args()
157+ os.makedirs(args.out, exist_ok=True)
158+ kind = "logf" if args.logf else "lin"
159+ tag = args.tag if args.tag is not None else "-".join(t for t in [("logf" if args.logf else ""), ("" if args.arch == "cnn2d" else args.arch), ("" if args.targets == "params" else args.targets)] if t)
160+ suffix = f"-{tag}" if tag else ""
161+ cache = args.cache or os.path.join(args.data, "ml", "cache" + ("-logf" if args.logf else ""))
162+ t_start = time.time()
163+ rows, feats = build_cache(args.data, cache, args.workers, kind)
164+ n = len(rows)
165+ y = np.stack([[math.log(r[t]) for t in TARGETS] for r in rows]).astype(np.float32)
166+ idx = np.array([r["index"] for r in rows])
167+ val = idx % args.val_every == 0
168+ tr = ~val
169+ if tr.sum() < 20 or val.sum() < 5:
170+ print(f"only {n} clips; need more"); return
171+ T = COMBOS if args.targets == "combos" else np.eye(len(TARGETS), dtype=np.float32)
172+ Tinv = np.linalg.inv(T)
173+ yt = y @ T.T
174+ mu, sd = yt[tr].mean(0), yt[tr].std(0) + 1e-6
175+ z = (yt - mu) / sd
176+ print(f"{n} clips: {tr.sum()} train, {val.sum()} val; features {feats.shape}; device {args.device}")
177+
178+ dev = torch.device(args.device)
179+ X = torch.from_numpy(feats) # float16 (n, F, T), kept on the CPU
180+ Z = torch.from_numpy(z)
181+ net = (SpecNet(len(TARGETS), n_freq=feats.shape[1]) if args.arch == "spec1d" else Net(len(TARGETS))).to(dev)
182+ print(f"arch {args.arch}, targets {args.targets}, {sum(p.numel() for p in net.parameters()) / 1e6:.2f} M parameters")
183+ opt = torch.optim.AdamW(net.parameters(), lr=args.lr, weight_decay=1e-4)
184+ sched = torch.optim.lr_scheduler.OneCycleLR(opt, max_lr=args.lr, total_steps=args.epochs * math.ceil(tr.sum() / args.batch))
185+ tr_idx, val_idx = np.where(tr)[0], np.where(val)[0]
186+ rng = np.random.default_rng(0)
187+
188+ def predict(ids):
189+ net.eval()
190+ out = []
191+ with torch.no_grad():
192+ for k in range(0, len(ids), 64):
193+ xb = X[ids[k:k + 64]].to(dev).float().unsqueeze(1)
194+ out.append(net(xb).cpu().numpy())
195+ return np.concatenate(out)
196+
197+ def augment(xb):
198+ """Label-preserving perturbations against memorization: a random onset
199+ shift of up to 8 frames (170 ms), and one masked frequency band and one
200+ masked time span per clip (the masked cells take the floor value 0)."""
201+ bsz, _, nf, nt = xb.shape
202+ out = torch.zeros_like(xb)
203+ for i in range(bsz):
204+ sh = int(rng.integers(0, 9))
205+ out[i, :, :, sh:] = xb[i, :, :, : nt - sh]
206+ f0, fw = int(rng.integers(0, nf)), int(rng.integers(0, nf // 16))
207+ out[i, :, f0:f0 + fw, :] = 0
208+ t0, tw = int(rng.integers(0, nt)), int(rng.integers(0, nt // 8))
209+ out[i, :, :, t0:t0 + tw] = 0
210+ return out
211+
212+ history, best, best_state = [], float("inf"), None
213+ for epoch in range(args.epochs):
214+ net.train()
215+ perm = rng.permutation(tr_idx)
216+ tot = 0.0
217+ for k in range(0, len(perm), args.batch):
218+ b = perm[k:k + args.batch]
219+ xb = X[b].to(dev).float().unsqueeze(1)
220+ if not args.no_augment:
221+ xb = augment(xb)
222+ zb = Z[b].to(dev)
223+ loss = nn.functional.mse_loss(net(xb), zb)
224+ opt.zero_grad(); loss.backward(); opt.step(); sched.step()
225+ tot += loss.item() * len(b)
226+ pv = predict(val_idx)
227+ vloss = float(np.mean((pv - z[val_idx]) ** 2))
228+ history.append({"epoch": epoch, "train_loss": tot / len(perm), "val_loss": vloss})
229+ if vloss < best:
230+ best, best_state = vloss, {k: v.detach().cpu().clone() for k, v in net.state_dict().items()}
231+ print(f"epoch {epoch}: train {tot / len(perm):.4f} val {vloss:.4f}", flush=True)
232+ net.load_state_dict(best_state)
233+ torch.save({"state": best_state, "mu": mu, "sd": sd, "targets": TARGETS, "T": T, "arch": args.arch, "input": kind, "n_freq": feats.shape[1]}, os.path.join(args.out, "model.pt"))
234+
235+ # Evaluation in log space, back in the natural scale of each target.
236+ pv = (predict(val_idx) * sd + mu) @ Tinv.T
237+ yv = y[val_idx]
238+ ym = np.broadcast_to(mu @ Tinv.T, yv.shape)
239+ res = {t: metrics(pv[:, i], yv[:, i], ym[:, i]) for i, t in enumerate(TARGETS)}
240+ d_pred = {t: pv[:, i] for i, t in enumerate(TARGETS)}
241+ d_true = {t: yv[:, i] for i, t in enumerate(TARGETS)}
242+ d_mean = {t: ym[:, i] for i, t in enumerate(TARGETS)}
243+ for name, f in DERIVED.items():
244+ res[name] = metrics(f(d_pred), f(d_true), f(d_mean))
245+ summary = {"input": kind, "arch": args.arch, "targets": args.targets, "augment": not args.no_augment, "n_clips": n, "n_train": int(tr.sum()), "n_val": int(val.sum()), "best_val_loss": float(best), "epochs": args.epochs,
246+ "minutes": (time.time() - t_start) / 60, "finished": time.strftime("%Y-%m-%d %H:%M:%S"), "metrics": res, "history": history}
247+ with open(os.path.join(args.out, "metrics.json"), "w") as f:
248+ json.dump(summary, f, indent=1)
249+ with open(os.path.join(args.out, "..", f"history{suffix}.jsonl"), "a") as f:
250+ f.write(json.dumps({k: v for k, v in summary.items() if k != "history"}) + "\n")
251+ np.save(os.path.join(args.out, "val_pred.npy"), np.stack([idx[val_idx], *yv.T, *pv.T], 1))
252+ write_report(args.out, summary, yv, pv, suffix)
253+ print("\n".join(f"{t:>8}: rel RMS {m['rel_rms_pct']:.1f}% (mean-baseline {m['baseline_rel_rms_pct']:.1f}%), R^2 {m['r2']:.2f}" for t, m in res.items()))
254+
255+
256+def write_report(out, summary, yv, pv, suffix=""):
257+ hist_lines = []
258+ hp = os.path.join(out, "..", f"history{suffix}.jsonl")
259+ if os.path.exists(hp):
260+ with open(hp) as f:
261+ hist_lines = [json.loads(l) for l in f if l.strip()]
262+ rows = "".join(
263+ f"<tr><td>{t}</td><td>{m['rel_rms_pct']:.1f}</td><td>{m['mae_pct']:.1f}</td><td>{m['baseline_rel_rms_pct']:.1f}</td><td>{m['r2']:.2f}</td></tr>"
264+ for t, m in summary["metrics"].items())
265+ imgs = "".join(f'<img src="data:image/png;base64,{scatter_png(yv[:, i], pv[:, i], t, *UNITS[t])}">' for i, t in enumerate(TARGETS))
266+ hist_rows = "".join(f"<tr><td>{h['finished']}</td><td>{h['n_clips']}</td>" + "".join(f"<td>{h['metrics'][t]['rel_rms_pct']:.1f}</td>" for t in TARGETS) + "</tr>" for h in hist_lines)
267+ curve = "".join(f"<tr><td>{h['epoch']}</td><td>{h['train_loss']:.4f}</td><td>{h['val_loss']:.4f}</td></tr>" for h in summary["history"][::max(1, len(summary["history"]) // 15)])
268+ html = f"""<!doctype html><meta charset="utf-8"><title>steelsim31 prediction report</title>
269+<style>body{{font:14px system-ui,sans-serif;margin:16px;color:#222}} table{{border-collapse:collapse;margin:8px 0}} td,th{{padding:3px 10px;border-bottom:1px solid #ddd;text-align:right}} td:first-child,th:first-child{{text-align:left}} img{{margin:4px}} h2{{font-size:16px;margin-top:20px}}</style>
270+<h1 style="font-size:18px">Prediction of the sheet parameters from the strike sound</h1>
271+<p>Input: {"log-frequency" if summary["input"] == "logf" else "linear-frequency"} spectrogram; architecture {summary.get("arch", "cnn2d")}; targets {summary.get("targets", "params")}. {summary['n_clips']} clips ({summary['n_train']} train, {summary['n_val']} held out), {summary['epochs']} epochs, {summary['minutes']:.1f} min, finished {summary['finished']}. Errors are relative (percent, from the logarithm) on the held-out clips; the baseline predicts the training mean.</p>
272+<table><tr><th>target</th><th>rel RMS %</th><th>rel MAE %</th><th>baseline rel RMS %</th><th>R²</th></tr>{rows}</table>
273+<p>Derived: E/rho (sets the wave speeds), kappa = h sqrt(E / (12 rho (1 - nu²))) (sets the bending frequencies), area = Lx Ly, aspect = Lx / Ly, kappa/area (the scale of the bending spectrum).</p>
274+<h2>Held-out predictions</h2>{imgs}
275+<h2>Error against the number of clips</h2>
276+<table><tr><th>run</th><th>clips</th>{"".join(f"<th>{t} %</th>" for t in TARGETS)}</tr>{hist_rows}</table>
277+<h2>Training curve (standardized MSE)</h2>
278+<table><tr><th>epoch</th><th>train</th><th>val</th></tr>{curve}</table>
279+"""
280+ with open(os.path.join(out, "report.html"), "w") as f:
281+ f.write(html)
282+ with open(os.path.join(out, "..", f"report{suffix}.html"), "w") as f:
283+ f.write(html)
284+
285+
286+if __name__ == "__main__":
287+ main()
reference-light-tap-lower-gain.wavadded+0−0View file
Binary file not shown.
reference-light-tap.wavadded+0−0View file
Binary file not shown.
reference.wavadded+0−0View file
Binary file not shown.
scripts/README-dataset.mdadded+22−0View file
@@ -0,0 +1,22 @@
1+# Dataset generation
2+
3+`scripts/generate.ts` renders many drumstick strikes headlessly and writes one WAV clip and one line of metadata per clip, for training and evaluating models that predict the sheet's material and dimensions from its sound. It serves `dist/` (build first with `npm run build`) to a headless Chromium with WebGPU through Playwright and drives several pages in parallel: one simulation is bound by the latency of its GPU dispatches, so a large GPU has room for many.
4+
5+ node --experimental-strip-types scripts/generate.ts --out DIR --count N [--start I] [--workers W] [--seed S] [--duration D] [--spacing S] [--fps F]
6+
7+Clip `i` is drawn from a generator seeded by `(seed, i)`, so the dataset is reproducible and can be regenerated or extended in pieces; existing clips in `DIR` are skipped. `scripts/show-samples.ts` prints the drawn parameters of a range of clips without rendering them, and `scripts/probe.ts` runs one configuration (a URL query, or `--seed S --index I` for a dataset clip with optional overrides) and prints the status as it runs.
8+
9+## Output
10+
11+- `DIR/clips/NNNNNN.wav`: 16-bit mono at 48 kHz, 3 s from the impact, peak-normalized to 0.99. The pressure in pascals at 1 m is `sample * wavScalePa` with `wavScalePa` from the metadata.
12+- `DIR/index.jsonl`: one JSON object per clip, appended as clips finish. Targets: `E` (Pa), `rho` (kg/m^3), `h` (m), `Lx`, `Ly` (m; the width and the height); `nu` is fixed at 0.29. Derived: `kappa` (bending wave constant sqrt(D / rho h), m^2/s), `mass` (kg). Nuisance variables: `features` (the initial shape in the text form of `shape.ts`, or `none`), `shapeScale`, `peakToPeak`, `maxSlope`, `sigma0`, `sigma1`, `sigmaIn`, `strikeFracX`, `strikeFracY`, `strikeVelocity`, `listenerDist`, `listenerAzimuth`, `listenerElevation`. Numerics: `spacing`, `nx`, `ny`, `dt`, `totalSteps`. Audio: `peakPa`, `rmsPa`, `wavScalePa`, `maxForce` (peak contact force, N), `wmax` (largest displacement, m). Bookkeeping: `wallSeconds`, `link` (the app's query string for the clip, which reproduces it in the browser).
13+- `DIR/failures.log`: clips that did not render (a runaway displacement, a lost GPU).
14+- `DIR/gallery.html`: a page (copy of `scripts/gallery.html`) that lists the newest clips with players, parameters, and optional spectrograms; serve `DIR` with any static server.
15+
16+`scripts/pull-dataset.sh` copies the clips and the index from a remote machine into a local directory every couple of minutes.
17+
18+## Sampling (version 1)
19+
20+All draws are independent. The targets are log-uniform over about plus or minus 30 percent of the defaults: E in 120 to 220 GPa, rho in 6000 to 9500 kg/m^3, h in 0.6 to 1.1 mm, width in 9 to 16 in, height in 13 to 24 in. The nuisance variables are meant to make the model rely on what a real recording would also show: the initial shape is flat with probability 0.08 and otherwise a random feature list (bows, twist, dents, a crease, a smooth random field; see `randomFeatures` in `shape.ts`) with a height scale log-uniform in 0.5 to 5 mm, redrawn at half the scale while the largest slope exceeds 0.1; the strike point is uniform in 0.6 to 0.8 of the width and 0.25 to 0.45 of the height, the stick speed log-uniform in 2 to 4 m/s; sigma0 and sigma1 are the defaults times a log-uniform factor in 0.7 to 1.4; the listener is at 1 m with azimuth and elevation uniform in plus or minus 15 degrees. The grid spacing is fixed at 8 mm, the sheet is a rectangle, and the excitation is a single strike.
21+
22+Two things to keep in mind when interpreting predictions. In the linear regime the sheet's motion depends on E and rho only through E / rho, and its frequencies scale as h sqrt(E / rho) / L^2; separating E from rho, or h from the size, must rely on the amplitude (through the stick's momentum exchange with the sheet's mass) and on the nonlinearity, both of which a real recording shows less reliably than the simulation. With a fixed grid spacing the grid dispersion at a given frequency depends on the bending wave constant, which is a simulator artifact a model could in principle exploit.
scripts/gallery.htmladded+132−0View file
@@ -0,0 +1,132 @@
1+<!doctype html>
2+<meta charset="utf-8">
3+<title>steelsim31 dataset</title>
4+<style>
5+ body { font: 14px system-ui, sans-serif; margin: 16px; color: #222; }
6+ h1 { font-size: 18px; margin: 0 0 8px; }
7+ #summary { margin-bottom: 12px; color: #444; }
8+ table { border-collapse: collapse; width: 100%; }
9+ th, td { padding: 4px 8px; border-bottom: 1px solid #ddd; text-align: right; white-space: nowrap; vertical-align: middle; }
10+ th { position: sticky; top: 0; background: #f4f4f4; cursor: pointer; }
11+ td.l, th.l { text-align: left; }
12+ audio { height: 28px; width: 220px; }
13+ canvas { display: block; }
14+ .feat { max-width: 420px; overflow: hidden; text-overflow: ellipsis; font-family: monospace; font-size: 11px; color: #555; white-space: nowrap; }
15+ .ctl { margin-bottom: 8px; display: flex; gap: 16px; align-items: center; flex-wrap: wrap; }
16+</style>
17+<h1>steelsim31 dataset</h1>
18+<div id="summary">loading…</div>
19+<div class="ctl">
20+ <label><input type="checkbox" id="auto" checked> reload every 30 s</label>
21+ <label>show <input type="number" id="limit" value="60" min="10" step="10" style="width:5em"> newest</label>
22+ <label><input type="checkbox" id="specs"> spectrograms</label>
23+ <button id="reload">reload</button>
24+ <a href="ml/report.html">training report (linear f)</a> <a href="ml/report-logf.html">training report (log f)</a> <a href="ml/report-specnet.html">training report (SpecNet)</a>
25+</div>
26+<table>
27+ <thead><tr>
28+ <th class="l">id</th><th class="l">audio</th><th>E GPa</th><th>ρ kg/m³</th><th>h mm</th><th>W×H in</th>
29+ <th>v m/s</th><th>strike</th><th>σ0</th><th>σ1·10⁴</th><th>listener</th><th>shape p-p mm</th><th class="l">features</th>
30+ <th>peak Pa</th><th>grid</th><th>wall s</th>
31+ </tr></thead>
32+ <tbody id="rows"></tbody>
33+</table>
34+<script>
35+const $ = (id) => document.getElementById(id);
36+let recs = [];
37+async function load() {
38+ const r = await fetch("index.jsonl", { cache: "no-store" });
39+ const text = await r.text();
40+ recs = text.split("\n").filter((l) => l.trim()).map((l) => JSON.parse(l));
41+ render();
42+}
43+function fmt(x, d) { return x == null ? "" : Number(x).toFixed(d); }
44+function render() {
45+ const n = recs.length;
46+ if (!n) { $("summary").textContent = "no clips yet"; return; }
47+ const first = new Date(recs[0].finished).getTime(), last = new Date(recs[n - 1].finished).getTime();
48+ const hours = Math.max(1e-6, (last - first) / 3.6e6);
49+ const rate = n > 1 ? ((n - 1) / hours).toFixed(0) : "?";
50+ const wall = recs.reduce((s, r) => s + r.wallSeconds, 0) / n;
51+ $("summary").textContent = `${n} clips, ${rate} clips/hour, mean ${wall.toFixed(0)} s wall per clip, last ${new Date(recs[n - 1].finished).toLocaleTimeString()}`;
52+ const limit = parseInt($("limit").value) || 60;
53+ const rows = recs.slice(-limit).reverse();
54+ const tb = $("rows");
55+ tb.innerHTML = "";
56+ for (const r of rows) {
57+ const tr = document.createElement("tr");
58+ const feats = r.features === "none" ? "flat" : r.features.split(";").map((f) => f.split(":")[0]).join(", ") + " (" + r.features + ")";
59+ tr.innerHTML = `<td class="l">${r.id}</td><td class="l"><audio controls preload="none" src="clips/${r.id}.wav"></audio></td>` +
60+ `<td>${fmt(r.E / 1e9, 0)}</td><td>${fmt(r.rho, 0)}</td><td>${fmt(r.h * 1e3, 2)}</td><td>${fmt(r.Lx / 0.0254, 1)}×${fmt(r.Ly / 0.0254, 1)}</td>` +
61+ `<td>${fmt(r.strikeVelocity, 2)}</td><td>${fmt(r.strikeFracX, 2)}, ${fmt(r.strikeFracY, 2)}</td><td>${fmt(r.sigma0, 2)}</td><td>${fmt(r.sigma1 * 1e4, 1)}</td>` +
62+ `<td>${fmt(r.listenerAzimuth, 0)}°, ${fmt(r.listenerElevation, 0)}°</td><td>${fmt(r.peakToPeak * 1e3, 1)}</td><td class="l"><div class="feat" title="${r.features}">${feats}</div></td>` +
63+ `<td>${fmt(r.peakPa, 2)}</td><td>${r.nx}×${r.ny}</td><td>${fmt(r.wallSeconds, 0)}</td>`;
64+ tb.appendChild(tr);
65+ if ($("specs").checked) {
66+ const tr2 = document.createElement("tr");
67+ const td = document.createElement("td"); td.colSpan = 16; td.className = "l";
68+ const c = document.createElement("canvas"); c.width = 900; c.height = 120; td.appendChild(c); tr2.appendChild(td); tb.appendChild(tr2);
69+ drawSpec(c, `clips/${r.id}.wav`);
70+ }
71+ }
72+}
73+// Spectrogram: 2048-point Hann windows, hop 512, 0 to 5 kHz, 80 dB range.
74+let actx = null;
75+async function drawSpec(canvas, url) {
76+ actx = actx || new (window.OfflineAudioContext || window.webkitOfflineAudioContext)(1, 48000, 48000);
77+ const buf = await (await fetch(url)).arrayBuffer();
78+ const audio = await actx.decodeAudioData(buf);
79+ const x = audio.getChannelData(0), fs = audio.sampleRate;
80+ const N = 2048, hop = 512, fmax = 5000, kmax = Math.floor((fmax / fs) * N);
81+ const frames = Math.floor((x.length - N) / hop);
82+ const g = canvas.getContext("2d"), W = canvas.width, H = canvas.height;
83+ const img = g.createImageData(W, H);
84+ const win = new Float32Array(N); for (let i = 0; i < N; i++) win[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / N);
85+ const re = new Float32Array(N), im = new Float32Array(N);
86+ const cols = new Array(frames);
87+ let peak = -Infinity;
88+ for (let f = 0; f < frames; f++) {
89+ for (let i = 0; i < N; i++) { re[i] = x[f * hop + i] * win[i]; im[i] = 0; }
90+ fft(re, im);
91+ const col = new Float32Array(kmax);
92+ for (let k = 0; k < kmax; k++) { col[k] = 10 * Math.log10(re[k] * re[k] + im[k] * im[k] + 1e-20); peak = Math.max(peak, col[k]); }
93+ cols[f] = col;
94+ }
95+ for (let px = 0; px < W; px++) {
96+ const f = Math.min(frames - 1, Math.floor((px / W) * frames));
97+ for (let py = 0; py < H; py++) {
98+ const k = Math.floor(((H - 1 - py) / H) * kmax);
99+ const v = Math.max(0, Math.min(1, (cols[f][k] - peak + 80) / 80));
100+ const o = 4 * (py * W + px);
101+ img.data[o] = 255 * v; img.data[o + 1] = 255 * v * v; img.data[o + 2] = 255 * (0.2 + 0.8 * v * v * v) * (1 - v) + 255 * v * v * v; img.data[o + 3] = 255;
102+ }
103+ }
104+ g.putImageData(img, 0, 0);
105+}
106+function fft(re, im) {
107+ const n = re.length;
108+ for (let i = 1, j = 0; i < n; i++) {
109+ let bit = n >> 1;
110+ for (; j & bit; bit >>= 1) j ^= bit;
111+ j ^= bit;
112+ if (i < j) { [re[i], re[j]] = [re[j], re[i]]; [im[i], im[j]] = [im[j], im[i]]; }
113+ }
114+ for (let len = 2; len <= n; len <<= 1) {
115+ const ang = (-2 * Math.PI) / len, wr = Math.cos(ang), wi = Math.sin(ang);
116+ for (let i = 0; i < n; i += len) {
117+ let cr = 1, ci = 0;
118+ for (let k = 0; k < len / 2; k++) {
119+ const a = i + k, b = a + len / 2;
120+ const tr = re[b] * cr - im[b] * ci, ti = re[b] * ci + im[b] * cr;
121+ re[b] = re[a] - tr; im[b] = im[a] - ti; re[a] += tr; im[a] += ti;
122+ const ncr = cr * wr - ci * wi; ci = cr * wi + ci * wr; cr = ncr;
123+ }
124+ }
125+ }
126+}
127+$("reload").onclick = load;
128+$("limit").onchange = render;
129+$("specs").onchange = render;
130+setInterval(() => { if ($("auto").checked) load(); }, 30000);
131+load();
132+</script>
scripts/generate.tsadded+315−0View file
@@ -0,0 +1,315 @@
1+// Dataset generator: renders many drumstick strikes headlessly and writes one
2+// WAV clip plus one line of metadata per clip. Serves dist/ (run `npm run
3+// build` first) to a headless Chromium with WebGPU through Playwright, and
4+// drives several pages in parallel, since one simulation is latency-bound and
5+// leaves a large GPU mostly idle.
6+//
7+// node --experimental-strip-types scripts/generate.ts --out DIR --count N [--start I] [--workers W] [--seed S] [--fps F]
8+//
9+// Clip i (global index) is drawn from a generator seeded by (seed, i), so the
10+// dataset is reproducible and can be extended or regenerated in pieces.
11+// Existing clips in DIR are skipped. Output: DIR/clips/NNNNNN.wav (16-bit,
12+// peak-normalized; the peak pressure in Pa is in the metadata) and DIR/index.jsonl.
13+
14+import { createServer } from "node:http";
15+import { readFileSync, existsSync, mkdirSync, appendFileSync, writeFileSync, readdirSync } from "node:fs";
16+import { fileURLToPath } from "node:url";
17+import { dirname, join, extname } from "node:path";
18+import { chromium, type Page, type Browser } from "playwright";
19+import { defaultParams, makeGrid, makeMaterial, type PhysParams } from "../src/params.ts";
20+import { encodeFeatures, initialShape, randomFeatures, shapeStats, type ShapeFeature } from "../src/shape.ts";
21+
22+const root = join(dirname(fileURLToPath(import.meta.url)), "..");
23+
24+// ---------- Arguments ----------
25+
26+const argv = process.argv.slice(2);
27+function arg(name: string, def: string): string {
28+ const k = argv.indexOf(`--${name}`);
29+ return k >= 0 && k + 1 < argv.length ? argv[k + 1] : def;
30+}
31+const OUT = arg("out", "dataset");
32+const COUNT = parseInt(arg("count", "10"));
33+const START = parseInt(arg("start", "0"));
34+const WORKERS = parseInt(arg("workers", "4"));
35+const SEED = parseInt(arg("seed", "1"));
36+const FPS = parseFloat(arg("fps", "120"));
37+const DURATION = parseFloat(arg("duration", "3"));
38+const JOB_TIMEOUT_MS = parseFloat(arg("timeout", "1800")) * 1000;
39+
40+// ---------- Random numbers ----------
41+
42+/** mulberry32, seeded by a hash of (seed, index). */
43+function rng(seed: number, index: number): () => number {
44+ let a = (Math.imul(seed, 0x9e3779b1) ^ Math.imul(index + 1, 0x85ebca6b)) >>> 0;
45+ for (let k = 0; k < 4; k++) { a = Math.imul(a ^ (a >>> 15), 0x2c1b3c6d); a = Math.imul(a ^ (a >>> 12), 0x297a2d39); }
46+ return () => {
47+ a = (a + 0x6d2b79f5) >>> 0;
48+ let t = a;
49+ t = Math.imul(t ^ (t >>> 15), t | 1);
50+ t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
51+ return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
52+ };
53+}
54+
55+// ---------- Parameter sampling ----------
56+
57+export interface Sample {
58+ id: string;
59+ index: number;
60+ params: PhysParams;
61+ /** Shape statistics on the simulation grid. */
62+ shapeScale: number;
63+ peakToPeak: number;
64+ maxSlope: number;
65+ nx: number;
66+ ny: number;
67+}
68+
69+const NU = 0.29;
70+const SLOPE_LIMIT = 0.1;
71+const SPACING = parseFloat(arg("spacing", "0.008"));
72+
73+/**
74+ * The targets (E, rho, h, Lx, Ly) are drawn log-uniformly and independently
75+ * over ranges of about +-30 percent around the defaults. The nuisance variables
76+ * (initial shape, strike point and speed, losses, listener direction) are
77+ * drawn independently of them. The grid spacing is fixed (8 mm by default),
78+ * so the grid dispersion at a given frequency varies with the bending wave
79+ * constant of the sheet; see the note in the dataset README.
80+ */
81+export function drawSample(seed: number, index: number): Sample {
82+ const rnd = rng(seed, index);
83+ const U = (a: number, b: number) => a + (b - a) * rnd();
84+ const logU = (a: number, b: number) => Math.exp(U(Math.log(a), Math.log(b)));
85+ const IN = 0.0254;
86+ const E = logU(120e9, 220e9);
87+ const rho = logU(6000, 9500);
88+ const h = logU(0.6e-3, 1.1e-3);
89+ const Lx = logU(9 * IN, 16 * IN);
90+ const Ly = logU(13 * IN, 24 * IN);
91+ const p: PhysParams = {
92+ ...defaultParams,
93+ shape: "rectangle",
94+ Lx, Ly, h, E, nu: NU, rho,
95+ features: [],
96+ sigma0: defaultParams.sigma0 * logU(0.7, 1.4),
97+ sigma1: defaultParams.sigma1 * logU(0.7, 1.4),
98+ sigmaIn: 20,
99+ excitation: "strike",
100+ strikeFracX: U(0.6, 0.8),
101+ strikeFracY: U(0.25, 0.45),
102+ strikeVelocity: logU(2, 4),
103+ hitSequence: "single",
104+ listenerDist: 1,
105+ listenerAzimuth: U(-15, 15),
106+ listenerElevation: U(-15, 15),
107+ spacing: SPACING,
108+ fps: FPS,
109+ duration: DURATION,
110+ };
111+ // Initial shape: flat with a small probability, otherwise a random feature
112+ // list with heights of the order of 0.5 to 5 mm (the default is a 3 mm bow);
113+ // redrawn at half the scale while the largest slope leaves the shallow-shell range.
114+ const grid = makeGrid(p);
115+ let shapeScale = rnd() < 0.08 ? 0 : logU(0.5e-3, 5e-3);
116+ let features: ShapeFeature[] = [];
117+ let stats = { peakToPeak: 0, maxSlope: 0 };
118+ for (let attempt = 0; shapeScale > 0 && attempt < 12; attempt++) {
119+ features = randomFeatures(p, shapeScale, rnd);
120+ stats = shapeStats(initialShape({ ...p, features }, grid), grid);
121+ if (stats.maxSlope <= SLOPE_LIMIT) break;
122+ shapeScale *= 0.5;
123+ }
124+ p.features = features;
125+ return { id: String(index).padStart(6, "0"), index, params: p, shapeScale, peakToPeak: stats.peakToPeak, maxSlope: stats.maxSlope, nx: grid.nx, ny: grid.ny };
126+}
127+
128+/** Query string for the app, in the form form.ts applyUrl reads. */
129+export function toQuery(p: PhysParams): string {
130+ const parts: string[] = ["autostart=1"];
131+ for (const k of Object.keys(defaultParams) as (keyof PhysParams)[]) {
132+ if (k === "features") {
133+ parts.push(`features=${encodeURIComponent(p.features.length ? encodeFeatures(p.features) : "none")}`);
134+ } else if (p[k] !== defaultParams[k]) {
135+ const v = p[k];
136+ parts.push(`${k}=${encodeURIComponent(typeof v === "number" ? String(v) : typeof v === "boolean" ? (v ? "1" : "0") : String(v))}`);
137+ }
138+ }
139+ return parts.join("&");
140+}
141+
142+/** Checks that the page parsed the targets as intended (the link is the app's own encoding of its parameters). */
143+function checkLink(link: string, p: PhysParams): void {
144+ const q = new URLSearchParams(link.includes("?") ? link.slice(link.indexOf("?") + 1) : link);
145+ for (const k of ["E", "rho", "h", "Lx", "Ly", "strikeVelocity", "sigma0", "sigma1"] as const) {
146+ const got = parseFloat(q.get(k) ?? String(defaultParams[k]));
147+ if (!(Math.abs(got - p[k]) <= 1e-5 * Math.abs(p[k]))) throw new Error(`parameter ${k} mismatch: page has ${got}, wanted ${p[k]} (link ${link})`);
148+ }
149+}
150+
151+// ---------- WAV ----------
152+
153+export function wavBytes(samples: Float32Array, fs: number, scale: number): Uint8Array {
154+ const buf = new ArrayBuffer(44 + samples.length * 2);
155+ const v = new DataView(buf);
156+ const str = (o: number, s: string) => { for (let k = 0; k < s.length; k++) v.setUint8(o + k, s.charCodeAt(k)); };
157+ str(0, "RIFF"); v.setUint32(4, 36 + samples.length * 2, true); str(8, "WAVE");
158+ str(12, "fmt "); v.setUint32(16, 16, true); v.setUint16(20, 1, true); v.setUint16(22, 1, true);
159+ v.setUint32(24, fs, true); v.setUint32(28, fs * 2, true); v.setUint16(32, 2, true); v.setUint16(34, 16, true);
160+ str(36, "data"); v.setUint32(40, samples.length * 2, true);
161+ for (let k = 0; k < samples.length; k++) {
162+ const s = Math.max(-1, Math.min(1, samples[k] * scale));
163+ v.setInt16(44 + 2 * k, Math.round(s * 32767), true);
164+ }
165+ return new Uint8Array(buf);
166+}
167+
168+// ---------- Static server for dist/ ----------
169+
170+const MIME: Record<string, string> = { ".html": "text/html", ".js": "text/javascript", ".css": "text/css", ".wgsl": "text/plain" };
171+export function serveDist(): Promise<string> {
172+ const dist = join(root, "dist");
173+ if (!existsSync(join(dist, "index.html"))) throw new Error("dist/index.html not found; run `npm run build` first");
174+ const server = createServer((req, res) => {
175+ const path = decodeURIComponent((req.url ?? "/").split("?")[0]);
176+ const file = join(dist, path === "/" ? "index.html" : path);
177+ try {
178+ const data = readFileSync(file);
179+ res.writeHead(200, { "content-type": MIME[extname(file)] ?? "application/octet-stream" });
180+ res.end(data);
181+ } catch {
182+ res.writeHead(404); res.end();
183+ }
184+ });
185+ return new Promise((resolve) => server.listen(0, "127.0.0.1", () => {
186+ const addr = server.address();
187+ resolve(`http://127.0.0.1:${typeof addr === "object" && addr ? addr.port : 0}`);
188+ }));
189+}
190+
191+// ---------- Rendering one clip ----------
192+
193+export const ARGS = ["--enable-unsafe-webgpu", "--enable-features=Vulkan", "--use-angle=vulkan", "--ignore-gpu-blocklist", "--no-sandbox"];
194+
195+interface Status { phase: string; t: number; totalSteps: number; wmax: number; peakPressure: number; samples: number; stepsPerSec: number; dt: number; nx: number; ny: number; stick: { maxForce: number } | null }
196+
197+export async function renderClip(page: Page, url: string, s: Sample): Promise<Record<string, unknown>> {
198+ const errors: string[] = [];
199+ const onErr = (e: Error) => errors.push(e.message);
200+ page.on("pageerror", onErr);
201+ const t0 = Date.now();
202+ await page.goto(`${url}/?${toQuery(s.params)}`);
203+ let st: Status | null = null;
204+ for (;;) {
205+ await page.waitForTimeout(250);
206+ st = await page.evaluate<Status | null>("window.steelsim ? window.steelsim.status() : null");
207+ const noGpu = await page.evaluate<string | null>("window.steelsim ? window.steelsim.noGpu() : null");
208+ if (noGpu) throw new Error(`no WebGPU: ${noGpu}`);
209+ if (errors.length) throw new Error(`page errors: ${errors.join("; ")}`);
210+ if (st && st.phase === "done") break;
211+ if (Date.now() - t0 > JOB_TIMEOUT_MS) throw new Error("timeout");
212+ }
213+ page.off("pageerror", onErr);
214+ const audio = new Float32Array(await page.evaluate<number[]>("window.steelsim.audio()"));
215+ const link = await page.evaluate<string>("window.steelsim.link()");
216+ checkLink(link, s.params);
217+ const wall = (Date.now() - t0) / 1000;
218+ if (!st) throw new Error("no status");
219+ if (!isFinite(st.wmax) || st.wmax > 0.1 || !isFinite(st.peakPressure)) throw new Error(`unstable: wmax ${st.wmax}, peak ${st.peakPressure}`);
220+ let peak = 0, sumsq = 0;
221+ for (let k = 0; k < audio.length; k++) { peak = Math.max(peak, Math.abs(audio[k])); sumsq += audio[k] * audio[k]; }
222+ const scale = peak > 0 ? 0.99 / peak : 0;
223+ const n = Math.min(audio.length, Math.round(DURATION * s.params.fs));
224+ writeFileSync(join(OUT, "clips", `${s.id}.wav`), wavBytes(audio.subarray(0, n), s.params.fs, scale));
225+ const p = s.params;
226+ const mat = makeMaterial(p);
227+ return {
228+ id: s.id,
229+ index: s.index,
230+ seed: SEED,
231+ // Targets
232+ E: p.E, rho: p.rho, h: p.h, Lx: p.Lx, Ly: p.Ly, nu: p.nu,
233+ // Derived
234+ kappa: Math.sqrt(mat.D / mat.rhoh),
235+ mass: mat.rhoh * p.Lx * p.Ly,
236+ // Nuisance
237+ features: p.features.length ? encodeFeatures(p.features) : "none",
238+ shapeScale: s.shapeScale, peakToPeak: s.peakToPeak, maxSlope: s.maxSlope,
239+ sigma0: p.sigma0, sigma1: p.sigma1, sigmaIn: p.sigmaIn,
240+ strikeFracX: p.strikeFracX, strikeFracY: p.strikeFracY, strikeVelocity: p.strikeVelocity,
241+ listenerDist: p.listenerDist, listenerAzimuth: p.listenerAzimuth, listenerElevation: p.listenerElevation,
242+ // Numerics
243+ spacing: p.spacing, nx: st.nx, ny: st.ny, dt: st.dt, fps: p.fps, fs: p.fs, duration: DURATION, totalSteps: st.totalSteps,
244+ // Audio
245+ samples: n, peakPa: peak, rmsPa: Math.sqrt(sumsq / Math.max(1, audio.length)), wavScalePa: peak > 0 ? 1 / scale : 0,
246+ maxForce: st.stick?.maxForce ?? null, wmax: st.wmax,
247+ // Bookkeeping
248+ wallSeconds: wall, stepsPerSec: st.stepsPerSec, link, finished: new Date().toISOString(),
249+ };
250+}
251+
252+// ---------- Main ----------
253+
254+async function main(): Promise<void> {
255+ mkdirSync(join(OUT, "clips"), { recursive: true });
256+ const indexFile = join(OUT, "index.jsonl");
257+ const done = new Set<string>();
258+ if (existsSync(indexFile)) for (const line of readFileSync(indexFile, "utf8").split("\n")) if (line.trim()) done.add(JSON.parse(line).id);
259+ const existing = new Set(readdirSync(join(OUT, "clips")).map((f) => f.replace(/\.wav$/, "")));
260+ const jobs: Sample[] = [];
261+ for (let i = START; i < START + COUNT; i++) {
262+ const s = drawSample(SEED, i);
263+ if (done.has(s.id) && existing.has(s.id)) continue;
264+ jobs.push(s);
265+ }
266+ console.log(`${jobs.length} clips to render (${COUNT - jobs.length} already done), ${WORKERS} workers, out ${OUT}`);
267+ if (!jobs.length) return;
268+ const url = await serveDist();
269+ let next = 0, finished = 0, failed = 0;
270+ const t0 = Date.now();
271+ const sleep = (ms: number) => new Promise((r) => setTimeout(r, ms));
272+ // Each worker owns a browser, so that a lost GPU process costs one worker a
273+ // relaunch rather than every worker its WebGPU; the starts are staggered.
274+ const worker = async (w: number) => {
275+ await sleep(1500 * w);
276+ let browser: Browser = await chromium.launch({ headless: true, args: ARGS });
277+ let page: Page = await (await browser.newContext()).newPage();
278+ const relaunch = async () => {
279+ try { await browser.close(); } catch { /* already gone */ }
280+ await sleep(3000);
281+ browser = await chromium.launch({ headless: true, args: ARGS });
282+ page = await (await browser.newContext()).newPage();
283+ };
284+ while (next < jobs.length) {
285+ const s = jobs[next++];
286+ for (let attempt = 0; ; attempt++) {
287+ try {
288+ const rec = await renderClip(page, url, s);
289+ appendFileSync(indexFile, JSON.stringify(rec) + "\n");
290+ finished++;
291+ const rate = (finished * 3600) / ((Date.now() - t0) / 1000);
292+ console.log(`[w${w}] ${s.id} done in ${(rec.wallSeconds as number).toFixed(0)} s: ${s.nx}x${s.ny}, E ${(s.params.E / 1e9).toFixed(0)} GPa, rho ${s.params.rho.toFixed(0)}, h ${(s.params.h * 1e3).toFixed(2)} mm, ${(s.params.Lx * 100).toFixed(0)}x${(s.params.Ly * 100).toFixed(0)} cm, peak ${(rec.peakPa as number).toFixed(2)} Pa | ${finished} done, ${failed} failed, ${rate.toFixed(0)} clips/h`);
293+ break;
294+ } catch (e) {
295+ const msg = String(e).replace(/\n/g, " ");
296+ if (/no WebGPU|closed|crashed|Protocol error|timeout/i.test(msg) && attempt < 3) {
297+ console.log(`[w${w}] ${s.id} retrying after: ${msg}`);
298+ await relaunch();
299+ continue;
300+ }
301+ failed++;
302+ console.log(`[w${w}] ${s.id} FAILED: ${msg}`);
303+ appendFileSync(join(OUT, "failures.log"), `${s.id}\t${new Date().toISOString()}\t${msg}\n`);
304+ break;
305+ }
306+ }
307+ }
308+ await browser.close();
309+ };
310+ await Promise.all(Array.from({ length: WORKERS }, (_, w) => worker(w)));
311+ console.log(`finished: ${finished} clips, ${failed} failures, ${((Date.now() - t0) / 60000).toFixed(1)} min`);
312+ process.exit(0);
313+}
314+
315+if (process.argv[1] && fileURLToPath(import.meta.url) === process.argv[1]) void main().catch((e) => { console.error(e); process.exit(1); });
scripts/probe.tsadded+51−0View file
@@ -0,0 +1,51 @@
1+// Runs one configuration headlessly and prints the status every half second.
2+// node --experimental-strip-types scripts/probe.ts "duration=1&E=1.6e11&..." [--seed S --index I (parameters of a dataset clip, plus overrides after &)] [--wav out.wav]
3+import { chromium } from "playwright";
4+import { createServer } from "node:http";
5+import { readFileSync } from "node:fs";
6+import { join, extname, dirname } from "node:path";
7+import { fileURLToPath } from "node:url";
8+import { writeFileSync } from "node:fs";
9+import { drawSample, toQuery, wavBytes } from "./generate.ts";
10+
11+const root = join(dirname(fileURLToPath(import.meta.url)), "..");
12+const argv = process.argv.slice(2);
13+const arg = (n: string, d: string) => { const k = argv.indexOf(`--${n}`); return k >= 0 ? argv[k + 1] : d; };
14+let query = argv.find((a) => !a.startsWith("--") && !/^\d+$/.test(a)) ?? "";
15+if (argv.includes("--index")) {
16+ const s = drawSample(parseInt(arg("seed", "1")), parseInt(arg("index", "0")));
17+ query = toQuery(s.params) + (query ? "&" + query : "");
18+} else query = "autostart=1&" + query;
19+const MIME: Record<string, string> = { ".html": "text/html", ".js": "text/javascript", ".css": "text/css" };
20+const server = createServer((req, res) => {
21+ const path = decodeURIComponent((req.url ?? "/").split("?")[0]);
22+ try { const f = join(root, "dist", path === "/" ? "index.html" : path); res.writeHead(200, { "content-type": MIME[extname(f)] ?? "application/octet-stream" }); res.end(readFileSync(f)); } catch { res.writeHead(404); res.end(); }
23+});
24+await new Promise<void>((r) => server.listen(0, "127.0.0.1", () => r()));
25+const port = (server.address() as { port: number }).port;
26+const browser = await chromium.launch({ headless: true, args: ["--enable-unsafe-webgpu", "--enable-features=Vulkan", "--use-angle=vulkan", "--ignore-gpu-blocklist", "--no-sandbox"] });
27+const page = await browser.newPage();
28+page.on("pageerror", (e) => console.log("pageerror", e.message));
29+page.on("console", (m) => { if (m.type() === "error") console.log("console.error", m.text()); });
30+console.log(query);
31+await page.goto(`http://127.0.0.1:${port}/?${query}`);
32+const t0 = Date.now();
33+for (;;) {
34+ await page.waitForTimeout(500);
35+ const st = await page.evaluate<Record<string, unknown> | null>("window.steelsim ? window.steelsim.status() : null");
36+ const noGpu = await page.evaluate<string | null>("window.steelsim ? window.steelsim.noGpu() : null");
37+ if (noGpu) { console.log("no gpu:", noGpu); break; }
38+ if (!st) continue;
39+ const stick = st.stick as { x: number; v: number; force: number; maxForce: number } | null;
40+ console.log(`${((Date.now() - t0) / 1000).toFixed(1)}s phase ${st.phase} t ${(st.t as number).toFixed(3)} ke ${(st.ke as number).toExponential(2)} wmax ${((st.wmax as number) * 1e3).toFixed(2)} mm peak ${(st.peakPressure as number).toFixed(3)} Pa ${(st.stepsPerSec as number / 1e3).toFixed(0)} k/s ${st.nx}x${st.ny}` + (stick ? ` stick x ${(stick.x * 1e3).toFixed(1)} mm v ${stick.v.toFixed(2)} F ${stick.force.toFixed(1)} Fmax ${stick.maxForce.toFixed(1)}` : ""));
41+ if (st.phase === "done") break;
42+}
43+if (argv.includes("--wav")) {
44+ const audio = new Float32Array(await page.evaluate<number[]>("window.steelsim.audio()"));
45+ let peak = 0;
46+ for (const v of audio) peak = Math.max(peak, Math.abs(v));
47+ writeFileSync(arg("wav", "probe.wav"), wavBytes(audio, 48000, peak > 0 ? 0.99 / peak : 0));
48+ console.log(`wrote ${arg("wav", "probe.wav")}, peak ${peak.toFixed(3)} Pa`);
49+}
50+await browser.close();
51+server.close();
scripts/pull-dataset.shadded+28−0View file
@@ -0,0 +1,28 @@
1+#!/bin/bash
2+# Pulls the clips and the index of a dataset rendered on another machine into a
3+# local directory, repeatedly, so that they can be listened to while the run
4+# progresses. The index is snapshotted first and every clip missing locally is
5+# fetched, so the local index never names a clip that is not there yet.
6+# scripts/pull-dataset.sh [remote host] [remote dir] [local dir] [interval s]
7+set -u
8+REMOTE=${1:-office}
9+RDIR=${2:-steelsim31-dataset/v1}
10+LDIR=${3:-dataset}
11+INTERVAL=${4:-120}
12+mkdir -p "$LDIR/clips"
13+cp "$(dirname "$0")/gallery.html" "$LDIR/gallery.html"
14+while true; do
15+ if ssh -x -o BatchMode=yes "$REMOTE" "cat $RDIR/index.jsonl 2>/dev/null" > "$LDIR/index.jsonl.tmp"; then
16+ ssh -x -o BatchMode=yes "$REMOTE" "ls $RDIR/clips" | sort > "$LDIR/.remote-clips"
17+ ls "$LDIR/clips" | sort > "$LDIR/.local-clips"
18+ comm -13 "$LDIR/.local-clips" "$LDIR/.remote-clips" > "$LDIR/.missing"
19+ if [ -s "$LDIR/.missing" ]; then
20+ ssh -x -o BatchMode=yes "$REMOTE" "cd $RDIR/clips && tar cf - $(tr '\n' ' ' < "$LDIR/.missing")" | tar xf - -C "$LDIR/clips"
21+ fi
22+ mv "$LDIR/index.jsonl.tmp" "$LDIR/index.jsonl"
23+ # The training report and history, if a model is being trained there.
24+ ssh -x -o BatchMode=yes "$REMOTE" "cd $RDIR && tar cf - ml/report*.html ml/history*.jsonl 2>/dev/null" | tar xf - -C "$LDIR" 2>/dev/null
25+ echo "$(date +%H:%M:%S) $(wc -l < "$LDIR/index.jsonl") clips in index, $(ls "$LDIR/clips" | wc -l) wav files local"
26+ fi
27+ sleep "$INTERVAL"
28+done
scripts/render-batch.tsadded+75−0View file
@@ -0,0 +1,75 @@
1+// Renders a batch of explicit parameter sets (for the optimizer in ml/fit.py):
2+// node --experimental-strip-types scripts/render-batch.ts --batch JOBS.json --out DIR [--workers W] [--duration D]
3+// JOBS.json is {"jobs": [{"id": "...", "params": {PhysParams overrides in SI, "features" as text}}]}.
4+// Writes DIR/clips/ID.wav and DIR/index.jsonl like generate.ts; existing clips are skipped.
5+
6+import { readFileSync, existsSync, mkdirSync, appendFileSync } from "node:fs";
7+import { join } from "node:path";
8+import { chromium, type Browser, type Page } from "playwright";
9+import { defaultParams, makeGrid, type PhysParams } from "../src/params.ts";
10+import { initialShape, parseFeatures, shapeStats } from "../src/shape.ts";
11+import { renderClip, serveDist, ARGS, type Sample } from "./generate.ts";
12+
13+const argv = process.argv.slice(2);
14+const arg = (n: string, d: string) => { const k = argv.indexOf(`--${n}`); return k >= 0 && k + 1 < argv.length ? argv[k + 1] : d; };
15+const OUT = arg("out", "batch");
16+const WORKERS = parseInt(arg("workers", "16"));
17+
18+interface Job { id: string; params: Record<string, unknown> }
19+
20+function toSample(job: Job): Sample {
21+ const p: PhysParams = { ...defaultParams, excitation: "strike", hitSequence: "single", duration: parseFloat(arg("duration", "3")) };
22+ for (const [k, v] of Object.entries(job.params)) {
23+ if (k === "features") (p as unknown as Record<string, unknown>).features = typeof v === "string" ? (v === "none" ? [] : parseFeatures(v) ?? []) : v;
24+ else if (k in p) (p as unknown as Record<string, unknown>)[k] = v;
25+ }
26+ const grid = makeGrid(p);
27+ const stats = shapeStats(initialShape(p, grid), grid);
28+ return { id: job.id, index: -1, params: p, shapeScale: 0, peakToPeak: stats.peakToPeak, maxSlope: stats.maxSlope, nx: grid.nx, ny: grid.ny };
29+}
30+
31+async function main(): Promise<void> {
32+ mkdirSync(join(OUT, "clips"), { recursive: true });
33+ const jobs = (JSON.parse(readFileSync(arg("batch", "jobs.json"), "utf8")).jobs as Job[]).filter((j) => !existsSync(join(OUT, "clips", `${j.id}.wav`)));
34+ if (!jobs.length) { console.log("nothing to render"); return; }
35+ const url = await serveDist();
36+ let next = 0, done = 0, failed = 0;
37+ const sleep = (ms: number) => new Promise((r) => setTimeout(r, ms));
38+ const worker = async (w: number) => {
39+ await sleep(500 * w);
40+ let browser: Browser = await chromium.launch({ headless: true, args: ARGS });
41+ let page: Page = await (await browser.newContext()).newPage();
42+ while (next < jobs.length) {
43+ const job = jobs[next++];
44+ const s = toSample(job);
45+ for (let attempt = 0; ; attempt++) {
46+ try {
47+ const rec = await renderClip(page, url, s);
48+ appendFileSync(join(OUT, "index.jsonl"), JSON.stringify(rec) + "\n");
49+ done++;
50+ console.log(`[w${w}] ${job.id} done in ${(rec.wallSeconds as number).toFixed(0)} s (${done} done, ${failed} failed)`);
51+ break;
52+ } catch (e) {
53+ const msg = String(e).replace(/\n/g, " ");
54+ if (/no WebGPU|closed|crashed|Protocol error|timeout/i.test(msg) && attempt < 3) {
55+ try { await browser.close(); } catch { /* gone */ }
56+ await sleep(3000);
57+ browser = await chromium.launch({ headless: true, args: ARGS });
58+ page = await (await browser.newContext()).newPage();
59+ continue;
60+ }
61+ failed++;
62+ console.log(`[w${w}] ${job.id} FAILED: ${msg}`);
63+ appendFileSync(join(OUT, "index.jsonl"), JSON.stringify({ id: job.id, failed: msg }) + "\n");
64+ break;
65+ }
66+ }
67+ }
68+ await browser.close();
69+ };
70+ await Promise.all(Array.from({ length: Math.min(WORKERS, jobs.length) }, (_, w) => worker(w)));
71+ console.log(`batch finished: ${done} done, ${failed} failed`);
72+ process.exit(0);
73+}
74+
75+void main().catch((e) => { console.error(e); process.exit(1); });
scripts/show-samples.tsadded+11−0View file
@@ -0,0 +1,11 @@
1+// Prints the sampled parameters of a range of clips without rendering them.
2+// node --experimental-strip-types scripts/show-samples.ts [--seed S] [--start I] [--count N]
3+import { drawSample } from "./generate.ts";
4+const argv = process.argv.slice(2);
5+const arg = (n: string, d: string) => { const k = argv.indexOf(`--${n}`); return k >= 0 ? argv[k + 1] : d; };
6+const seed = parseInt(arg("seed", "1")), start = parseInt(arg("start", "0")), count = parseInt(arg("count", "8"));
7+for (let i = start; i < start + count; i++) {
8+ const s = drawSample(seed, i);
9+ const p = s.params;
10+ console.log(s.id, `E ${(p.E / 1e9).toFixed(0)} rho ${p.rho.toFixed(0)} h ${(p.h * 1e3).toFixed(2)} L ${(p.Lx * 100).toFixed(0)}x${(p.Ly * 100).toFixed(0)} v ${p.strikeVelocity.toFixed(2)} grid ${s.nx}x${s.ny} scale ${(s.shapeScale * 1e3).toFixed(2)}mm p2p ${(s.peakToPeak * 1e3).toFixed(1)} slope ${s.maxSlope.toFixed(3)} sig0 ${p.sigma0.toFixed(2)} sig1 ${p.sigma1.toExponential(1)} mass ${(p.rho * p.h * p.Lx * p.Ly).toFixed(2)} kg`);
11+}