aa03146Maraca shake simulator: MATLAB model and interactive web appmagland 1// Web Worker that renders the maraca shake by executing the MATLAB
2// implementation in ../matlab through numbl. The stages run as separate
3// executeCode calls so intermediates can be cached and cheap parameter
4// changes skip the expensive ones:
5// solve tier -> maraca_basis + maraca_impacts (cached on solve-tier params)
6// radiate tier -> maraca_radiate (cached on the listener and direct_gain)
7// output tier -> maraca_output (always cheap, always re-run)
8// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
10import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
11import { hashMFiles, type MFile } from '../mfiles';
12import { SOLVE_KEYS, RADIATE_KEYS, type MaracaParams } from '../params';
14export interface RenderRequest {
15 type: 'render';
16 id: number;
17 params: MaracaParams;
18 quality: 'preview' | 'full';
19 mfiles: MFile[];
20}
22export interface RenderResponse {
23 type: 'result';
24 id: number;
25 quality: 'preview' | 'full';
26 y: Float32Array; // final signal, [-1, 1]
27 sampleRate: number;
28 impactT: Float32Array; // impact times, s
29 impactV: Float32Array; // impact speeds, m/s
30 impactPos: Float32Array; // 3 x I unit positions on the shell
31 rateT: Float32Array; // expected-rate envelope, s
32 rateHz: Float32Array; // expected impacts per second
33 modeF: Float32Array; // mode frequencies, Hz
34 modeL: Float32Array; // spherical-harmonic degree of each mode
35 modeZeta: Float32Array; // damping ratio of each mode
36 modeSigma: Float32Array; // radiation efficiency of each mode
37 diag: {
38 modesKept: number;
39 breathingHz: number;
40 breathingDryHz: number;
41 impactCount: number;
42 meanV: number;
43 maxV: number;
44 meanTauUs: number;
45 durationS: number;
46 peakPa: number;
47 peakSplDb: number;
48 directPeakPa: number;
49 listenerDistM: number;
50 warnings: string[];
51 solveMs: number; // 0 when the cached solve was reused
52 radiateMs: number; // 0 when the cached radiation was reused
53 solveCached: boolean;
54 radiateCached: boolean;
55 oversample: number;
56 };
57}
59export interface ProgressMessage {
60 type: 'progress';
61 id: number;
62 stage: 'solve' | 'radiate' | 'output';
63 frac: number;
64}
66// Raised by a fault in the .m code (the user can edit it), so it is reported
67// rather than left to stall the request.
68export interface ErrorMessage {
69 type: 'error';
70 id: number;
71 message: string;
72 file?: string;
73 line?: number;
74}
76// Impacts drawn for the raster plot. The full process runs to thousands, and
77// plotting every one of them is neither visible nor worth the transfer.
78const MAX_PLOT_IMPACTS = 4000;
80interface SolveEntry {
81 timp: RuntimeValue;
82 jimp: RuntimeValue;
83 tauimp: RuntimeValue;
84 fimp: RuntimeValue;
85 ximp: RuntimeValue;
86 fsHi: number;
87 impactT: Float32Array;
88 impactV: Float32Array;
89 impactPos: Float32Array;
90 rateT: Float32Array;
91 rateHz: Float32Array;
92 modeF: Float32Array;
93 modeL: Float32Array;
94 modeZeta: Float32Array;
95 modeSigma: Float32Array;
96 breathingHz: number;
97 breathingDryHz: number;
98 impactCount: number;
99 meanV: number;
100 maxV: number;
101 meanTauUs: number;
102 durationS: number;
103}
105interface RadiateEntry {
106 p: RuntimeValue;
107 directPeakPa: number;
108 listenerDistM: number;
109 listenerInside: boolean;
110}
112// Small LRU caches so preview- and full-quality intermediates coexist. A
113// radiate-tier change right after a full render then reuses the full solve,
114// and the "preview" is silently upgraded to full quality.
115const solveCache = new Map<string, SolveEntry>();
116const radiateCache = new Map<string, RadiateEntry>();
117const SOLVE_CACHE_MAX = 2;
118const RADIATE_CACHE_MAX = 3;
120function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
121 const v = cache.get(key);
122 if (v !== undefined) {
123 cache.delete(key);
124 cache.set(key, v);
125 }
126 return v;
127}
129function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
130 cache.delete(key);
131 cache.set(key, value);
132 while (cache.size > max) {
133 const oldest = cache.keys().next().value as string;
134 cache.delete(oldest);
135 }
136}
138function keyFor(params: MaracaParams, keys: readonly (keyof MaracaParams)[]): string {
139 return keys.map((k) => String(params[k])).join('|');
140}
142// MATLAB source lines assigning every parameter onto a struct P.
143function pLiteral(params: MaracaParams): string {
144 const lines: string[] = [];
145 for (const [k, v] of Object.entries(params)) {
146 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
147 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
148 }
149 return lines.join('\n');
150}
152function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
153 return executeCode(
154 source,
155 { initialVariableValues: inputs, displayResults: false },
156 mfiles,
157 'stage.m',
158 );
159}
161function asTensor(v: RuntimeValue): RuntimeTensor {
162 const t = v as RuntimeTensor;
163 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
164 return t;
165}
167function asScalar(v: RuntimeValue): number {
168 if (typeof v === 'number') return v;
169 const t = v as RuntimeTensor;
170 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
171 if (typeof v === 'boolean') return v ? 1 : 0;
172 throw new Error('expected scalar result from numbl stage');
173}
175function toF32(v: RuntimeValue): Float32Array {
176 const t = asTensor(v);
177 const out = new Float32Array(t.data.length);
178 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
179 return out;
180}
182const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
184self.onmessage = (ev: MessageEvent<RenderRequest>) => {
185 const msg = ev.data;
186 if (msg.type !== 'render') return;
188 // numbl reports a file it could not parse through console.warn and then
189 // skips it, so the failure surfaces at the call site as "Undefined
190 // function". Capture the warning to report the actual syntax error instead.
191 const warnings: string[] = [];
192 const origWarn = console.warn;
193 console.warn = (...args: unknown[]) => {
194 warnings.push(args.map(String).join(' '));
195 };
197 try {
198 render(msg);
199 } catch (err) {
200 // Most likely a fault in the (user-editable) .m code. Report it so the
201 // request completes and the message reaches the code panel.
202 const e = err as { message?: string; file?: string; line?: number };
203 let m: ErrorMessage = {
204 type: 'error',
205 id: msg.id,
206 message: e?.message ? String(e.message) : String(err),
207 file: typeof e?.file === 'string' ? e.file : undefined,
208 line: typeof e?.line === 'number' ? e.line : undefined,
209 };
210 for (const w of warnings) {
211 const hit = SKIP_WARNING.exec(w);
212 if (hit) {
213 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
214 break;
215 }
216 }
217 (self as unknown as Worker).postMessage(m);
218 } finally {
219 console.warn = origWarn;
220 }
221};
223function render(msg: RenderRequest): void {
224 const { id, params, mfiles } = msg;
225 let quality = msg.quality;
227 // The .m sources are part of the cache identity: editing the code must
228 // invalidate everything computed from the previous version.
229 const codeKey = hashMFiles(mfiles);
230 const solveParamsKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
231 const fullSolveKey = solveParamsKey + `|os=${params.oversample}`;
232 if (quality === 'preview' && solveCache.has(fullSolveKey)) {
233 quality = 'full'; // full-rate impact process already cached
234 }
236 const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
237 const fsHi = params.sample_rate * oversample;
239 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
241 // ---- solve stage: mode basis and the stochastic impact process ----
242 const solveKey = solveParamsKey + `|os=${oversample}`;
243 let solveMs = 0;
244 let solveCached = true;
245 let solveEntry = lruGet(solveCache, solveKey);
246 if (!solveEntry) {
247 solveCached = false;
248 post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
249 const t0 = performance.now();
250 const r = runStage(mfiles, `
251${pLiteral(params)}
252basis = maraca_basis(P);
253[timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P);
254modeF = basis.f;
255modeL = basis.l;
256modeZeta = basis.zeta;
257modeSigma = basis.sigma;
258breathingHz = basis.breathingHz;
259breathingDryHz = basis.breathingDryHz;
260impactCount = ginfo.count;
261meanV = ginfo.meanV;
262maxV = ginfo.maxV;
263meanTauUs = ginfo.meanTauUs;
264durationS = ginfo.durationS;
265rateT = ginfo.rateT;
266rateHz = ginfo.rateHz;
267vimp = ginfo.vimp;
268`);
269 solveMs = performance.now() - t0;
271 const impactCount = asScalar(r.variableValues.impactCount);
272 const allT = toF32(r.variableValues.timp);
273 const allV = toF32(r.variableValues.vimp);
274 const allX = toF32(r.variableValues.ximp);
275 // thin uniformly for the raster plot, keeping the time ordering
276 const stride = Math.max(1, Math.ceil(impactCount / MAX_PLOT_IMPACTS));
277 const nPlot = Math.ceil(impactCount / stride);
278 const impactT = new Float32Array(nPlot);
279 const impactV = new Float32Array(nPlot);
280 const impactPos = new Float32Array(nPlot * 3);
281 for (let i = 0, j = 0; j < nPlot; i += stride, j++) {
282 impactT[j] = allT[i];
283 impactV[j] = allV[i];
284 impactPos[j * 3] = allX[i * 3];
285 impactPos[j * 3 + 1] = allX[i * 3 + 1];
286 impactPos[j * 3 + 2] = allX[i * 3 + 2];
287 }
289 solveEntry = {
290 timp: r.variableValues.timp,
291 jimp: r.variableValues.jimp,
292 tauimp: r.variableValues.tauimp,
293 fimp: r.variableValues.fimp,
294 ximp: r.variableValues.ximp,
295 fsHi,
296 impactT,
297 impactV,
298 impactPos,
299 rateT: toF32(r.variableValues.rateT),
300 rateHz: toF32(r.variableValues.rateHz),
301 modeF: toF32(r.variableValues.modeF),
302 modeL: toF32(r.variableValues.modeL),
303 modeZeta: toF32(r.variableValues.modeZeta),
304 modeSigma: toF32(r.variableValues.modeSigma),
305 breathingHz: asScalar(r.variableValues.breathingHz),
306 breathingDryHz: asScalar(r.variableValues.breathingDryHz),
307 impactCount,
308 meanV: asScalar(r.variableValues.meanV),
309 maxV: asScalar(r.variableValues.maxV),
310 meanTauUs: asScalar(r.variableValues.meanTauUs),
311 durationS: asScalar(r.variableValues.durationS),
312 };
313 lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
314 }
316 // ---- radiate stage ----
317 const radiateKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
318 let radiateMs = 0;
319 let radiateCached = true;
320 let radEntry = lruGet(radiateCache, radiateKey);
321 if (!radEntry) {
322 radiateCached = false;
323 post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
324 const t0 = performance.now();
325 const r = runStage(mfiles, `
326${pLiteral(params)}
327basis = maraca_basis(P);
328[p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, ${solveEntry.fsHi});
329directPeakPa = radinfo.directPeakPa;
330listenerDistM = radinfo.listenerDistM;
331listenerInside = radinfo.listenerInside;
332`, {
333 timp: solveEntry.timp,
334 jimp: solveEntry.jimp,
335 tauimp: solveEntry.tauimp,
336 fimp: solveEntry.fimp,
337 ximp: solveEntry.ximp,
338 });
339 radiateMs = performance.now() - t0;
340 radEntry = {
341 p: r.variableValues.p,
342 directPeakPa: asScalar(r.variableValues.directPeakPa),
343 listenerDistM: asScalar(r.variableValues.listenerDistM),
344 listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
345 };
346 lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
347 }
349 // ---- output stage ----
350 post({ type: 'progress', id, stage: 'output', frac: 0.9 });
351 const rOut = runStage(mfiles, `
352${pLiteral(params)}
353[y, peakPa] = maraca_output(P, p, ${oversample});
354`, { p: radEntry.p });
356 const y32 = toF32(rOut.variableValues.y);
357 const peakPa = asScalar(rOut.variableValues.peakPa);
359 const warnings: string[] = [];
360 if (radEntry.listenerInside) {
361 warnings.push(
362 'the listener is at or inside the shell surface; the multipole field is evaluated at the nearest valid radius instead',
363 );
364 }
365 if (!params.normalize) {
366 let peak = 0;
367 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
368 if (peak > 1.0) {
369 warnings.push(
370 `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
371 );
372 }
373 }
374 const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
375 if (topMode < 0.4 * params.max_mode_hz) {
376 warnings.push(
377 `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max degree l to keep the shell from sounding truncated`,
378 );
379 }
381 const resp: RenderResponse = {
382 type: 'result',
383 id,
384 quality,
385 y: y32,
386 sampleRate: params.sample_rate,
387 impactT: solveEntry.impactT.slice(),
388 impactV: solveEntry.impactV.slice(),
389 impactPos: solveEntry.impactPos.slice(),
390 rateT: solveEntry.rateT.slice(),
391 rateHz: solveEntry.rateHz.slice(),
392 modeF: solveEntry.modeF.slice(),
393 modeL: solveEntry.modeL.slice(),
394 modeZeta: solveEntry.modeZeta.slice(),
395 modeSigma: solveEntry.modeSigma.slice(),
396 diag: {
397 modesKept: solveEntry.modeF.length,
398 breathingHz: solveEntry.breathingHz,
399 breathingDryHz: solveEntry.breathingDryHz,
400 impactCount: solveEntry.impactCount,
401 meanV: solveEntry.meanV,
402 maxV: solveEntry.maxV,
403 meanTauUs: solveEntry.meanTauUs,
404 durationS: solveEntry.durationS,
405 peakPa,
406 peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
407 directPeakPa: radEntry.directPeakPa,
408 listenerDistM: radEntry.listenerDistM,
409 warnings,
410 solveMs,
411 radiateMs,
412 solveCached,
413 radiateCached,
414 oversample,
415 },
416 };
417 (self as unknown as Worker).postMessage(resp, [
418 resp.y.buffer,
419 resp.impactT.buffer,
420 resp.impactV.buffer,
421 resp.impactPos.buffer,
422 resp.rateT.buffer,
423 resp.rateHz.buffer,
424 resp.modeF.buffer,
425 resp.modeL.buffer,
426 resp.modeZeta.buffer,
427 resp.modeSigma.buffer,
428 ]);
429}