1/**
2 * The checks, against the real pipeline: SceneSpec -> generated MATLAB ->
3 * numbl lowering -> generated WGSL -> GPU. Physics first — does each kind of
4 * object sound the pitch it is tuned to, does the sealed box actually seal,
5 * do two objects in one scene both reach the microphone — and the planner's
6 * mechanics after.
7 *
8 * Everything runs on deliberately small grids: these are correctness checks,
9 * and a 32-point air grid already carries every code path the 128-point one
10 * does.
11 */
12import { StudioSession, specParams } from '../src/studio/session.ts';
13import { cloneSpec, presetByKey } from '../src/studio/presets.ts';
14import { makeSceneFile, parseSceneFile, serializeSceneFile } from '../src/studio/scenefile.ts';
15import { generate } from '../src/studio/codegen.ts';
16import { fastRender, prepareRoom, roomContributions } from '../src/studio/fast.ts';
17import { synthesize } from '../src/studio/modal.ts';
18import { makeAirGrid, stableDt } from '../src/grids.ts';
19import { C_AIR, CFL } from '../src/units.ts';
20import type { Params, SceneSpec } from '../src/studio/types.ts';
21import { planPlayback } from '../src/audio/play.ts';
23type Check = (name: string, ok: boolean, detail: string) => void;
24type Log = (s: string) => void;
26const peakOf = (a: Float32Array): number => {
27 let m = 0;
28 for (const v of a) m = Math.max(m, Math.abs(v));
29 return m;
30};
32const spec = (key: string): SceneSpec => cloneSpec(presetByKey(key)!.spec);
34async function makeSession(
35 device: GPUDevice,
36 sceneSpec: SceneSpec,
37 nx: number,
38 params: Params = {},
39 sceneParams: Params = {},
40 operandBudget?: number,
41): Promise<StudioSession> {
42 return StudioSession.create({
43 device, spec: sceneSpec, nx, params, sceneParams, operandBudget,
44 });
45}
47/** Frequency content of a series sampled every `sdt` seconds. */
48function dftAt(series: number[], sdt: number, f: number): number {
49 let re = 0;
50 let im = 0;
51 for (let i = 0; i < series.length; i++) {
52 const w = 2 * Math.PI * f * i * sdt;
53 re += series[i] * Math.cos(w);
54 im -= series[i] * Math.sin(w);
55 }
56 return Math.hypot(re, im) / series.length;
57}
59/** Run a session and sample one state element every `every` steps. */
60async function sampleField(
61 s: StudioSession,
62 field: string,
63 index: number,
64 totalSteps: number,
65 every: number,
66): Promise<{ series: number[]; sdt: number }> {
67 const series: number[] = [];
68 for (let done = 0; done < totalSteps; done += every) {
69 s.step(every);
70 const u = await s.read(field);
71 series.push(u[index]);
72 }
73 return { series, sdt: every * s.dt };
74}
76/** The compiled plan: per-object external ops present, excitation in place,
77 * mic riding along. */
78export async function planChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
79 const s = await makeSession(device, spec('dulcimer'), 32);
80 const { step } = s.describe();
81 for (const l of step) log(` ${l}`);
82 const externals = step.filter((l) => l.startsWith('external'));
83 const wanted = ['dxx_s1(u_s1)', 'dxx_s1(um_s1)', 'dxxxx_s1(u_s1)', 'rad_s1(acc_s1)', 'brg_s1(un_s1)', 'lapw(p, wall)'];
84 for (const w of wanted) {
85 check(
86 `the step uses ${w}`,
87 externals.some((l) => l.includes(w)),
88 externals.length ? externals.join('; ') : 'no external ops planned',
89 );
90 }
92 s.excite();
93 const u = await s.read('u_s1');
94 const amp = s.params.amp_s1;
95 const peak = peakOf(u);
96 check(
97 'the pluck draws the string to its set height',
98 Math.abs(peak - amp) < 0.15 * amp,
99 `peak |u| = ${peak.toExponential(3)}, amp = ${amp}`,
100 );
101 check('the ends are pinned', u[0] === 0 && u[u.length - 1] === 0, `u[0]=${u[0]}, u[end]=${u[u.length - 1]}`);
103 s.step(100);
104 check(
105 'the microphone samples once per timestep',
106 s.recorder.count === 100,
107 `${s.recorder.count} samples after 100 steps`,
108 );
109 const plan = planPlayback(s.recorder.count, s.dt);
110 check(
111 'playback is real time at the solver rate',
112 plan.realTime && Math.abs(plan.rate * s.dt - 1) < 1e-9,
113 `rate ${plan.rate.toFixed(0)} Hz, dt ${s.dt.toExponential(3)}`,
114 );
115 s.destroy();
116}
118/** The generated sources look like what the Lab view promises. */
119export async function codegenChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
120 const s = await makeSession(device, spec('duet'), 32);
121 const m = s.generated.modelSource;
122 const sc = s.generated.sceneSource;
123 check(
124 'the generated model carries both strings',
125 m.includes('un_s1 = pin_s1') && m.includes('un_s2 = pin_s2'),
126 'looked for un_s1 and un_s2 updates',
127 );
128 check(
129 'the generated scene builds a profile per object',
130 sc.includes('prof_s1 =') && sc.includes('prof_s2 =') && sc.includes('bprof_s2 ='),
131 'looked for prof_s1, prof_s2, bprof_s2',
132 );
133 log(` model: ${m.split('\n').length} lines, scene: ${sc.split('\n').length} lines`);
134 s.destroy();
135}
137/** A saved scene file round-trips: same spec, same generated MATLAB, and a
138 * malformed file is refused with a named error rather than a compile
139 * failure three layers down. */
140export async function sceneFileChecks(_device: GPUDevice, check: Check, log: Log): Promise<void> {
141 const original = spec('duet');
142 const defaults = specParams(original);
143 const file = makeSceneFile({
144 spec: original,
145 params: { ...defaults.params, f0_s1: 300 },
146 sceneParams: { ...defaults.sceneParams, oy_s1: -0.05 },
147 mic: { x: 0.1, y: 0.05, z: 0.08 },
148 nx: 96,
149 });
150 const parsed = parseSceneFile(serializeSceneFile(file));
151 check(
152 'a scene file round-trips its values',
153 parsed.params.f0_s1 === 300 && parsed.sceneParams.oy_s1 === -0.05 &&
154 parsed.nx === 96 && parsed.spec.objects.length === 2,
155 `f0_s1 = ${parsed.params.f0_s1}, oy_s1 = ${parsed.sceneParams.oy_s1}`,
156 );
158 const air = makeAirGrid(32, 1.0);
159 const dt = stableDt(air.h, C_AIR, CFL);
160 const a = generate(original, air, dt);
161 const b = generate(parsed.spec, air, dt);
162 check(
163 'a round-tripped spec generates identical MATLAB',
164 a.modelSource === b.modelSource && a.sceneSource === b.sceneSource,
165 `model ${a.modelSource === b.modelSource ? 'matches' : 'DIFFERS'}, ` +
166 `scene ${a.sceneSource === b.sceneSource ? 'matches' : 'DIFFERS'}`,
167 );
169 const bad = [
170 ['{}', 'format marker'],
171 ['{"soundbox":1,"spec":{"objects":[],"body":{"type":"none"}},"mic":{"x":0,"y":0,"z":0},"nx":32}', 'non-empty'],
172 ['{"soundbox":1,"spec":{"objects":[{"id":"s1","type":"gong","size":{"L1":0.3}}],"body":{"type":"none"}},"mic":{"x":0,"y":0,"z":0},"nx":32}', 'unknown type'],
173 ] as const;
174 for (const [text, wants] of bad) {
175 let msg = '';
176 try {
177 parseSceneFile(text);
178 } catch (e) {
179 msg = e instanceof Error ? e.message : String(e);
180 }
181 check(
182 `a malformed scene file is refused (${wants})`,
183 msg.includes(wants),
184 msg || 'no error was thrown',
185 );
186 }
187 log(' scene file round-trip and refusals exercised');
188}
190/** A string sounds the pitch it is tuned to. */
191export async function stringChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
192 // An ideal string (no stiffness, no damping to speak of), watched at a
193 // point for a dozen periods; the zero crossings say the frequency.
194 const f0 = 294;
195 const s = await makeSession(device, spec('dulcimer'), 32, {
196 f0_s1: f0, B_s1: 0, sig1_s1: 0, t60_s1: 8,
197 });
198 s.excite();
199 const g = s.grids.s1;
200 const node = Math.round(g.n1 * 0.4);
201 const periods = 12;
202 const stepsTotal = Math.round(periods / f0 / s.dt);
203 const { series } = await sampleField(s, 'u_s1', node, stepsTotal, 10);
204 let crossings = 0;
205 for (let i = 1; i < series.length; i++) {
206 if ((series[i - 1] < 0 && series[i] >= 0) || (series[i - 1] >= 0 && series[i] < 0)) crossings++;
207 }
208 const measured = crossings / 2 / (stepsTotal * s.dt);
209 log(` ${crossings} zero crossings over ${(stepsTotal * s.dt * 1000).toFixed(1)} ms -> ${measured.toFixed(1)} Hz`);
210 check(
211 `the string sounds its fundamental (${f0} Hz)`,
212 Math.abs(measured - f0) < 0.04 * f0,
213 `measured ${measured.toFixed(1)} Hz`,
214 );
216 // The decay knob means what it says: at half of t60 the string should be
217 // down 30 dB.
218 const t60 = 0.5;
219 const s2 = await makeSession(device, spec('dulcimer'), 32, {
220 t60_s1: t60, B_s1: 0, sig1_s1: 0,
221 });
222 s2.excite();
223 const before = peakOf(await s2.read('u_s1'));
224 s2.step(Math.round(t60 / 2 / s2.dt));
225 const after = peakOf(await s2.read('u_s1'));
226 const db = 20 * Math.log10(after / before);
227 log(` |u| ${before.toExponential(2)} -> ${after.toExponential(2)} in ${t60 / 2} s (${db.toFixed(1)} dB)`);
228 check(
229 't60 decays the string on schedule',
230 Math.abs(db + 30) < 4,
231 `${db.toFixed(1)} dB over t60/2, want -30`,
232 );
233 s.destroy();
234 s2.destroy();
235}
237/**
238 * A bar rings at the frequencies its own discrete dispersion relation
239 * predicts. The bar grid is coarse by construction (an explicit stiff scheme
240 * wants coarse grids), so the continuous n² f1 comb is audibly detuned; what
241 * must hold exactly is the scheme's own relation: mode n of the pinned
242 * fourth-difference operator has eigenvalue lam_n = (4/h²) sin²(n π h / 2L),
243 * and rings at f_n = kap·lam_n / 2π. This is the "you can hear the
244 * dispersion" check: the same numbers the Lab spectrum overlay draws.
245 */
246export async function barChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
247 const f1 = 780;
248 const s = await makeSession(device, spec('chime'), 32, {
249 f1_b1: f1, sig1_b1: 0, t60_b1: 12,
250 });
251 s.excite();
252 const g = s.grids.b1;
253 const kap = (2 * f1 * g.L1 * g.L1) / Math.PI;
254 const lam = (n: number): number =>
255 (4 / (g.h1 * g.h1)) * Math.sin((n * Math.PI * g.h1) / (2 * g.L1)) ** 2;
256 const fd = (n: number): number => (kap * lam(n) * lam(n) === 0 ? 0 : (kap * lam(n)) / (2 * Math.PI));
257 const node = Math.round(g.n1 * 0.37);
258 const total = Math.round(0.06 / s.dt);
259 const { series, sdt } = await sampleField(s, 'u_b1', node, total, 4);
260 const comb = (dftAt(series, sdt, fd(1)) + dftAt(series, sdt, fd(2))) / 2;
261 const off =
262 (dftAt(series, sdt, 1.4 * fd(1)) + dftAt(series, sdt, Math.sqrt(fd(1) * fd(2)))) / 2;
263 log(
264 ` discrete modes ${fd(1).toFixed(0)}, ${fd(2).toFixed(0)} Hz ` +
265 `(continuous would be ${f1}, ${4 * f1}); comb ${comb.toExponential(2)}, off ${off.toExponential(2)}`,
266 );
267 check(
268 'the bar rings where its discrete dispersion relation says',
269 comb > 5 * off,
270 `comb/off = ${(comb / off).toFixed(1)}`,
271 );
272 s.destroy();
273}
275/** A membrane struck at its centre sounds its (1,1) fundamental. */
276export async function membraneChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
277 const f0 = 110;
278 const s = await makeSession(device, spec('drum'), 32, {
279 f0_m1: f0, sig1_m1: 0, t60_m1: 6, strikex_m1: 0.5, strikey_m1: 0.5,
280 });
281 s.excite();
282 const g = s.grids.m1;
283 const centre = Math.round(g.n1 / 2) + g.n1 * Math.round(g.n2 / 2);
284 const total = Math.round(10 / f0 / s.dt);
285 const { series, sdt } = await sampleField(s, 'u_m1', centre, total, 10);
286 const at = dftAt(series, sdt, f0);
287 const off = (dftAt(series, sdt, 1.31 * f0) + dftAt(series, sdt, 0.71 * f0)) / 2;
288 log(` membrane at f0 ${at.toExponential(2)}, off ${off.toExponential(2)}`);
289 check(
290 `the membrane sounds its (1,1) mode (${f0} Hz)`,
291 at > 4 * off,
292 `on/off = ${(at / off).toFixed(1)}`,
293 );
294 s.destroy();
295}
297/** The wall mask is a wall: a sealed box keeps the sound out. */
298export async function wallChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
299 // Bridge drive only (the direct route off), so the source sits above the
300 // plate; a microphone inside a sealed box should hear far less than one
301 // outside beside it.
302 const seconds = 0.06;
303 const run = async (mic: [number, number, number], holer: number): Promise<number> => {
304 const s = await makeSession(
305 device,
306 spec('dulcimer'),
307 64,
308 { gair_s1: 0, gbridge_s1: 1 },
309 { holer, absorb: 0, boxd: 0.12, thick: 0.02 },
310 );
311 s.setMic(...mic);
312 s.excite();
313 s.step(Math.round(seconds / s.dt));
314 const trace = await s.recorder.read();
315 s.destroy();
316 return peakOf(trace);
317 };
318 const inside = await run([0, 0, -0.06], 0);
319 const outside = await run([0, 0, 0.08], 0);
320 log(` |p| inside the sealed box ${inside.toExponential(2)}, outside ${outside.toExponential(2)}`);
321 check(
322 'a sealed box keeps the sound out',
323 inside < 0.05 * outside,
324 `inside/outside = ${(inside / outside).toExponential(2)}`,
325 );
327 const open = await run([0, 0, -0.06], 0.05);
328 log(` |p| in the cavity: sealed ${inside.toExponential(2)}, open ${open.toExponential(2)}`);
329 check(
330 'the sound hole lets the cavity speak',
331 open > 5 * inside,
332 `open/sealed = ${(open / inside).toExponential(2)}`,
333 );
334}
336/** Two objects in one scene both reach the microphone. */
337export async function duetChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
338 const s = await makeSession(device, spec('duet'), 64, {
339 B_s1: 0, B_s2: 0, sig1_s1: 0, sig1_s2: 0,
340 });
341 s.setMic(0.12, 0.08, 0.1);
342 s.excite();
343 const seconds = 0.12;
344 const total = Math.round(seconds / s.dt);
345 for (let done = 0; done < total; done += 2048) {
346 s.step(Math.min(2048, total - done));
347 await s.sync();
348 }
349 const trace = await s.recorder.read();
350 const series = [...trace];
351 const at = (f: number): number => dftAt(series, s.dt, f);
352 const s1 = at(294);
353 const s2 = at(441);
354 const off = (at(360) + at(520)) / 2;
355 log(` duet: 294 Hz ${s1.toExponential(2)}, 441 Hz ${s2.toExponential(2)}, off ${off.toExponential(2)}`);
356 check('the first string reaches the microphone', s1 > 4 * off, `294/off = ${(s1 / off).toFixed(1)}`);
357 check('the second string reaches the microphone', s2 > 4 * off, `441/off = ${(s2 / off).toFixed(1)}`);
358 s.destroy();
359}
361/** Moving an object is a scene edit: no recompile, and the profile follows. */
362export async function placementChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
363 const s = await makeSession(device, spec('bare'), 32);
364 const before = s.scene.fields.prof_s1.slice();
365 s.setScene({ ...s.sceneParams, oy_s1: 0.1 });
366 const after = s.scene.fields.prof_s1;
367 let moved = 0;
368 for (let i = 0; i < before.length; i++) {
369 if (Math.abs(before[i] - after[i]) > 1e-6) moved++;
370 }
371 log(` ${moved} profile cells changed after moving the string`);
372 check('moving an object moves its coupling profile', moved > 0, `${moved} cells changed`);
373 s.excite();
374 s.step(50);
375 const p = await s.read('p');
376 let bad = 0;
377 for (const v of p) if (!Number.isFinite(v)) bad++;
378 check('the moved scene still runs', bad === 0, `${bad} non-finite values`);
379 s.destroy();
380}
382/** Air sound sits on the string's partial comb, not somewhere else. */
383export async function spectrumChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
384 const f0 = 294;
385 const s = await makeSession(device, spec('dulcimer'), 64, { f0_s1: f0 });
386 s.setMic(0.12, 0.08, 0.1);
387 s.excite();
388 const seconds = 0.12;
389 const total = Math.round(seconds / s.dt);
390 for (let done = 0; done < total; done += 2048) {
391 s.step(Math.min(2048, total - done));
392 await s.sync();
393 }
394 const trace = await s.recorder.read();
395 const series = [...trace];
396 const comb = (dftAt(series, s.dt, f0) + dftAt(series, s.dt, 2 * f0) + dftAt(series, s.dt, 3 * f0)) / 3;
397 const off = (dftAt(series, s.dt, 1.41 * f0) + dftAt(series, s.dt, 2.53 * f0)) / 2;
398 log(` comb ${comb.toExponential(2)}, off-comb ${off.toExponential(2)}`);
399 check(
400 'the microphone hears the string’s partials',
401 comb > 8 * off,
402 `comb/off = ${(comb / off).toFixed(1)}`,
403 );
404 s.destroy();
405}
407/** A starved operand budget splits kernels instead of failing. */
408export async function splitChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
409 const s = await makeSession(device, spec('dulcimer'), 32, {}, {}, 2);
410 const { step } = s.describe();
411 const parts = step.filter((l) => l.includes('_part')).length;
412 log(` ${parts} split kernels at budget 2`);
413 check('a tight binding budget splits the update', parts > 0, `${parts} split kernels`);
414 s.excite();
415 s.step(50);
416 const p = await s.read('p');
417 let bad = 0;
418 for (const v of p) if (!Number.isFinite(v)) bad++;
419 check('the split model still runs', bad === 0, `${bad} non-finite values`);
421 // The split result must agree with the unsplit one.
422 const s2 = await makeSession(device, spec('dulcimer'), 32);
423 s2.excite();
424 s2.step(50);
425 const p2 = await s2.read('p');
426 let worst = 0;
427 const scale = Math.max(peakOf(p), 1e-30);
428 for (let i = 0; i < p.length; i++) worst = Math.max(worst, Math.abs(p[i] - p2[i]));
429 check(
430 'split and unsplit agree',
431 worst < 1e-4 * scale,
432 `worst |diff| = ${worst.toExponential(2)} of peak ${scale.toExponential(2)}`,
433 );
434 s.destroy();
435 s2.destroy();
436}
438/** The modal fast path (src/studio/fast.ts): every preset's mic trace,
439 * synthesized from closed-form modes and one reciprocal room recording per
440 * scene, must match the timestepped solver to its own f32 rounding noise —
441 * across all three object types, both string routes, and every body. And
442 * it must keep matching from the *same* room when pitch, decay and strike
443 * change, since the room depends only on geometry. */
444export async function modalChecks(device: GPUDevice, check: Check, log: Log): Promise<void> {
445 const relErr = (a: Float32Array, b: Float32Array, n: number): number => {
446 let d = 0;
447 let r = 0;
448 for (let i = 0; i < n; i++) {
449 d += (a[i] - b[i]) * (a[i] - b[i]);
450 r += b[i] * b[i];
451 }
452 return Math.sqrt(d / Math.max(r, 1e-300));
453 };
454 const steps = 6000;
456 const run = async (sceneSpec: SceneSpec, mic: { x: number; y: number; z: number }, params: Params): Promise<Float32Array> => {
457 const s = await makeSession(device, sceneSpec, 32, params);
458 s.setMic(mic.x, mic.y, mic.z);
459 s.excite();
460 s.step(steps);
461 await s.sync();
462 const trace = await s.recorder.read();
463 s.destroy();
464 return trace;
465 };
467 // One scene per object type: the drum (membrane over a shell), the
468 // dulcimer (string, both routes, over the box), and the kit (membrane and
469 // bar together — two objects in one room).
470 for (const key of ['drum', 'dulcimer', 'kit'] as const) {
471 const preset = presetByKey(key)!;
472 const sceneSpec = cloneSpec(preset.spec);
473 const defaults = specParams(sceneSpec);
474 const fast = await fastRender({
475 device, spec: sceneSpec, params: defaults.params, sceneParams: defaults.sceneParams,
476 nx: 32, mic: preset.mic, steps, irSeconds: 1e9, modeOpts: { eps: 1e-9 },
477 });
478 const gpu = await run(sceneSpec, preset.mic, {});
479 const err = relErr(fast.trace, gpu, steps);
480 log(` ${key}: ${fast.room.rows} kernels, ${fast.used} modes; rel err ${err.toExponential(2)}`);
481 check(`the modal path reproduces the ${key}`, err < 5e-3, `rel err ${err.toExponential(2)}`);
482 }
484 // Retune and restrike from one kept room: pitch, decay, brightness and
485 // excitation move, the room does not. Membranes are recorded a decade
486 // deeper than synthesis keeps, so a moved strike stays covered.
487 {
488 const preset = presetByKey('drum')!;
489 const sceneSpec = cloneSpec(preset.spec);
490 const defaults = specParams(sceneSpec);
491 const fr = await prepareRoom({
492 device, spec: sceneSpec, params: defaults.params, sceneParams: defaults.sceneParams,
493 nx: 32, mic: preset.mic, irSeconds: 1e9, maxKernelSteps: steps, modeOpts: { eps: 1e-12 },
494 });
495 const edits: Params = { f0_m1: 150, t60_m1: 0.7, sig1_m1: 0.002, strikex_m1: 0.62, striker_m1: 0.03 };
496 const paramsB = { ...defaults.params, ...edits };
497 const { contribs, used, wanted, missing } = roomContributions(fr, sceneSpec, paramsB, { eps: 1e-8 });
498 const traceB = synthesize(contribs, fr.room.steps, steps);
499 const gpuB = await run(sceneSpec, preset.mic, edits);
500 const errB = relErr(traceB, gpuB, steps);
501 log(` retuned drum: ${used}/${wanted} modes (${missing} unrecorded); rel err ${errB.toExponential(2)}`);
502 // A stray borderline mode may sit outside the recorded set (a mode near
503 // a node of the first strike); the error bound is the gate, the count
504 // just guards against wholesale mismatch.
505 check(
506 'the same room serves a retuned, restruck drum',
507 errB < 5e-3 && missing <= 0.02 * wanted,
508 `rel err ${errB.toExponential(2)}, ${missing}/${wanted} modes missing from the room`,
509 );
510 }
511}