1// Web Worker that renders the woodblock strike 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 -> woodblock_solve (cached on solve-tier params + oversample)
6// radiate tier -> woodblock_radiate (cached on radiate-tier params + M)
7// output tier -> woodblock_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 WoodblockParams } from '../params';
14export interface RenderRequest {
15 type: 'render';
16 id: number;
17 params: WoodblockParams;
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 force: Float32Array; // contact force, first forceMs at fsHi
29 forceFsHz: number;
30 modeShapes: Float32Array; // n x shapePts, concatenated
31 shapePts: number;
32 diag: {
33 modesKept: number;
34 partialsHz: number[];
35 ratios: number[];
36 helmholtzHz: number;
37 contactMs: number;
38 peakForceN: number;
39 peakPa: number;
40 peakSplDb: number;
41 listenerDistM: number;
42 warnings: string[];
43 solveMs: number; // 0 when the cached solve was reused
44 radiateMs: number; // 0 when the cached radiation was reused
45 solveCached: boolean;
46 radiateCached: boolean;
47 oversample: number;
48 elements: number;
49 };
50}
52export interface ProgressMessage {
53 type: 'progress';
54 id: number;
55 stage: 'solve' | 'radiate' | 'output';
56 frac: number;
57}
59// Raised by a fault in the .m code (the user can edit it), so it is reported
60// rather than left to stall the request.
61export interface ErrorMessage {
62 type: 'error';
63 id: number;
64 message: string;
65 file?: string;
66 line?: number;
67}
69const SHAPE_PTS = 128;
70const FORCE_MS = 8;
72interface SolveEntry {
73 Qdd: RuntimeValue;
74 Xdd: RuntimeValue;
75 fsHi: number;
76 f: number[];
77 modeShapes: Float32Array;
78 force: Float32Array;
79 helmholtzHz: number;
80 contactMs: number;
81 peakForceN: number;
82}
84interface RadiateEntry {
85 p: RuntimeValue;
86 warnings: string[];
87}
89// Small LRU caches so preview- and full-quality intermediates coexist. A
90// radiate-tier change right after a full render then reuses the full solve,
91// and the "preview" is silently upgraded to full quality.
92const solveCache = new Map<string, SolveEntry>();
93const radiateCache = new Map<string, RadiateEntry>();
94const SOLVE_CACHE_MAX = 2;
95const RADIATE_CACHE_MAX = 3;
97function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
98 const v = cache.get(key);
99 if (v !== undefined) {
100 cache.delete(key);
101 cache.set(key, v);
102 }
103 return v;
104}
106function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
107 cache.delete(key);
108 cache.set(key, value);
109 while (cache.size > max) {
110 const oldest = cache.keys().next().value as string;
111 cache.delete(oldest);
112 }
113}
115function keyFor(params: WoodblockParams, keys: readonly (keyof WoodblockParams)[]): string {
116 return keys.map((k) => String(params[k])).join('|');
117}
119// MATLAB source lines assigning every parameter onto a struct P.
120function pLiteral(params: WoodblockParams): string {
121 const lines: string[] = [];
122 for (const [k, v] of Object.entries(params)) {
123 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
124 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
125 }
126 return lines.join('\n');
127}
129function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
130 return executeCode(
131 source,
132 { initialVariableValues: inputs, displayResults: false },
133 mfiles,
134 'stage.m',
135 );
136}
138function asTensor(v: RuntimeValue): RuntimeTensor {
139 const t = v as RuntimeTensor;
140 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
141 return t;
142}
144function asScalar(v: RuntimeValue): number {
145 if (typeof v === 'number') return v;
146 const t = v as RuntimeTensor;
147 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
148 if (typeof v === 'boolean') return v ? 1 : 0;
149 throw new Error('expected scalar result from numbl stage');
150}
152const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
154self.onmessage = (ev: MessageEvent<RenderRequest>) => {
155 const msg = ev.data;
156 if (msg.type !== 'render') return;
158 // numbl reports a file it could not parse through console.warn and then
159 // skips it, so the failure surfaces at the call site as "Undefined
160 // function". Capture the warning to report the actual syntax error instead.
161 const warnings: string[] = [];
162 const origWarn = console.warn;
163 console.warn = (...args: unknown[]) => {
164 warnings.push(args.map(String).join(' '));
165 };
167 try {
168 render(msg);
169 } catch (err) {
170 // Most likely a fault in the (user-editable) .m code. Report it so the
171 // request completes and the message reaches the code panel.
172 const e = err as { message?: string; file?: string; line?: number };
173 let m: ErrorMessage = {
174 type: 'error',
175 id: msg.id,
176 message: e?.message ? String(e.message) : String(err),
177 file: typeof e?.file === 'string' ? e.file : undefined,
178 line: typeof e?.line === 'number' ? e.line : undefined,
179 };
180 for (const w of warnings) {
181 const hit = SKIP_WARNING.exec(w);
182 if (hit) {
183 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
184 break;
185 }
186 }
187 (self as unknown as Worker).postMessage(m);
188 } finally {
189 console.warn = origWarn;
190 }
191};
193function render(msg: RenderRequest): void {
194 const { id, params, mfiles } = msg;
195 let quality = msg.quality;
197 // The .m sources are part of the cache identity: editing the code must
198 // invalidate everything computed from the previous version.
199 const codeKey = hashMFiles(mfiles);
200 const solveParamsKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
201 const fullSolveKey = solveParamsKey + `|os=${params.oversample}`;
202 if (quality === 'preview' && solveCache.has(fullSolveKey)) {
203 quality = 'full'; // full solve already cached; only cheap stages remain
204 }
206 const oversample = quality === 'preview' ? Math.min(8, params.oversample) : params.oversample;
207 const elements = quality === 'preview'
208 ? Math.min(24, params.radiating_elements)
209 : params.radiating_elements;
210 const fsHi = params.sample_rate * oversample;
212 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
214 // ---- solve stage ----
215 const solveKey = solveParamsKey + `|os=${oversample}`;
216 let solveMs = 0;
217 let solveCached = true;
218 let solveEntry = lruGet(solveCache, solveKey);
219 if (!solveEntry) {
220 solveCached = false;
221 post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
222 const t0 = performance.now();
223 const r = runStage(mfiles, `
224${pLiteral(params)}
225basis = woodblock_basis(P);
226[Qdd, Xdd, force, helmholtzHz, contactMs, peakForceN, maxRealLam] = woodblock_solve(P, basis, ${oversample});
227f = basis.f;
228shapes = basis.shapes(1:32:end, :);
229`);
230 solveMs = performance.now() - t0;
232 const fT = asTensor(r.variableValues.f);
233 const f = Array.from(fT.data);
234 const n = f.length;
236 // decimated mode shapes, normalized per mode for the viewer
237 const shapesT = asTensor(r.variableValues.shapes);
238 const rows = shapesT.shape[0];
239 const modeShapes = new Float32Array(n * SHAPE_PTS);
240 for (let m = 0; m < n; m++) {
241 let peak = 0;
242 for (let i = 0; i < rows; i++) peak = Math.max(peak, Math.abs(shapesT.data[m * rows + i]));
243 for (let i = 0; i < SHAPE_PTS; i++) {
244 const j = Math.round((i / (SHAPE_PTS - 1)) * (rows - 1));
245 modeShapes[m * SHAPE_PTS + i] = peak > 0 ? shapesT.data[m * rows + j] / peak : 0;
246 }
247 }
249 const forceT = asTensor(r.variableValues.force);
250 const forceLen = Math.min(forceT.data.length, Math.round((FORCE_MS / 1000) * fsHi));
251 const force = new Float32Array(forceLen);
252 for (let i = 0; i < forceLen; i++) force[i] = forceT.data[i];
254 solveEntry = {
255 Qdd: r.variableValues.Qdd,
256 Xdd: r.variableValues.Xdd,
257 fsHi,
258 f,
259 modeShapes,
260 force,
261 helmholtzHz: asScalar(r.variableValues.helmholtzHz),
262 contactMs: asScalar(r.variableValues.contactMs),
263 peakForceN: asScalar(r.variableValues.peakForceN),
264 };
265 lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
266 }
268 // ---- radiate stage ----
269 const radiateKey = solveKey + '||' + keyFor(params, RADIATE_KEYS) + `|M=${elements}`;
270 let radiateMs = 0;
271 let radiateCached = true;
272 let radEntry = lruGet(radiateCache, radiateKey);
273 if (!radEntry) {
274 radiateCached = false;
275 post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
276 const t0 = performance.now();
277 const r = runStage(mfiles, `
278${pLiteral(params)}
279basis = woodblock_basis(P);
280[p, spacingNeed, listenerInside] = woodblock_radiate(P, basis, Qdd, Xdd, ${solveEntry.fsHi}, ${elements});
281`, { Qdd: solveEntry.Qdd, Xdd: solveEntry.Xdd });
282 radiateMs = performance.now() - t0;
284 const warnings: string[] = [];
285 const spacingNeed = asScalar(r.variableValues.spacingNeed);
286 if (spacingNeed > 0) {
287 const fMax = solveEntry.f[solveEntry.f.length - 1];
288 warnings.push(
289 `element spacing is coarse at ${fMax.toFixed(0)} Hz; consider ${spacingNeed} radiating elements`,
290 );
291 }
292 if (asScalar(r.variableValues.listenerInside) !== 0) {
293 const d = Math.hypot(params.listener_x_m, params.listener_y_m, params.listener_z_m);
294 warnings.push(
295 `listener is ${(d * 100).toFixed(0)} cm from the block centre, inside the source region; point-source elements are crude here`,
296 );
297 }
298 radEntry = { p: r.variableValues.p, warnings };
299 lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
300 }
302 // ---- output stage ----
303 post({ type: 'progress', id, stage: 'output', frac: 0.9 });
304 const rOut = runStage(mfiles, `
305${pLiteral(params)}
306[y, peakPa] = woodblock_output(P, p, ${oversample});
307`, { p: radEntry.p });
309 const yT = asTensor(rOut.variableValues.y);
310 const y32 = new Float32Array(yT.data.length);
311 for (let i = 0; i < yT.data.length; i++) y32[i] = yT.data[i];
312 const peakPa = asScalar(rOut.variableValues.peakPa);
314 const warnings = [...radEntry.warnings];
315 if (!params.normalize) {
316 let peak = 0;
317 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
318 if (peak > 1.0) {
319 warnings.push(
320 `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
321 );
322 }
323 }
324 const lx = params.listener_x_m, ly = params.listener_y_m, lz = params.listener_z_m;
326 const force = solveEntry.force.slice();
327 const modeShapes = solveEntry.modeShapes.slice();
328 const resp: RenderResponse = {
329 type: 'result',
330 id,
331 quality,
332 y: y32,
333 sampleRate: params.sample_rate,
334 force,
335 forceFsHz: solveEntry.fsHi,
336 modeShapes,
337 shapePts: SHAPE_PTS,
338 diag: {
339 modesKept: solveEntry.f.length,
340 partialsHz: solveEntry.f,
341 ratios: solveEntry.f.map((f) => f / solveEntry.f[0]),
342 helmholtzHz: solveEntry.helmholtzHz,
343 contactMs: solveEntry.contactMs,
344 peakForceN: solveEntry.peakForceN,
345 peakPa,
346 peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
347 listenerDistM: Math.sqrt(lx * lx + ly * ly + lz * lz),
348 warnings,
349 solveMs,
350 radiateMs,
351 solveCached,
352 radiateCached,
353 oversample,
354 elements,
355 },
356 };
357 (self as unknown as Worker).postMessage(resp, [y32.buffer, force.buffer, modeShapes.buffer]);
358}