// Web Worker that renders the maraca shake by executing the MATLAB // implementation in ../matlab through numbl. The stages run as separate // executeCode calls so intermediates can be cached and cheap parameter // changes skip the expensive ones: // solve tier -> maraca_basis + maraca_impacts (cached on solve-tier params) // radiate tier -> maraca_radiate (cached on the listener and direct_gain) // output tier -> maraca_output (always cheap, always re-run) // Tensors cross the boundary as numbl RuntimeValues (column-major Float64). import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl'; import { hashMFiles, type MFile } from '../mfiles'; import { SOLVE_KEYS, RADIATE_KEYS, type MaracaParams } from '../params'; export interface RenderRequest { type: 'render'; id: number; params: MaracaParams; quality: 'preview' | 'full'; mfiles: MFile[]; } export interface RenderResponse { type: 'result'; id: number; quality: 'preview' | 'full'; y: Float32Array; // final signal, [-1, 1] sampleRate: number; impactT: Float32Array; // impact times, s impactV: Float32Array; // impact speeds, m/s impactPos: Float32Array; // 3 x I unit positions on the shell rateT: Float32Array; // expected-rate envelope, s rateHz: Float32Array; // expected impacts per second modeF: Float32Array; // mode frequencies, Hz modeL: Float32Array; // spherical-harmonic degree of each mode modeZeta: Float32Array; // damping ratio of each mode modeSigma: Float32Array; // radiation efficiency of each mode diag: { modesKept: number; breathingHz: number; breathingDryHz: number; impactCount: number; meanV: number; maxV: number; meanTauUs: number; durationS: number; peakPa: number; peakSplDb: number; directPeakPa: number; listenerDistM: number; warnings: string[]; solveMs: number; // 0 when the cached solve was reused radiateMs: number; // 0 when the cached radiation was reused solveCached: boolean; radiateCached: boolean; oversample: number; }; } export interface ProgressMessage { type: 'progress'; id: number; stage: 'solve' | 'radiate' | 'output'; frac: number; } // Raised by a fault in the .m code (the user can edit it), so it is reported // rather than left to stall the request. export interface ErrorMessage { type: 'error'; id: number; message: string; file?: string; line?: number; } // Impacts drawn for the raster plot. The full process runs to thousands, and // plotting every one of them is neither visible nor worth the transfer. const MAX_PLOT_IMPACTS = 4000; interface SolveEntry { timp: RuntimeValue; jimp: RuntimeValue; tauimp: RuntimeValue; fimp: RuntimeValue; ximp: RuntimeValue; fsHi: number; impactT: Float32Array; impactV: Float32Array; impactPos: Float32Array; rateT: Float32Array; rateHz: Float32Array; modeF: Float32Array; modeL: Float32Array; modeZeta: Float32Array; modeSigma: Float32Array; breathingHz: number; breathingDryHz: number; impactCount: number; meanV: number; maxV: number; meanTauUs: number; durationS: number; } interface RadiateEntry { p: RuntimeValue; directPeakPa: number; listenerDistM: number; listenerInside: boolean; } // Small LRU caches so preview- and full-quality intermediates coexist. A // radiate-tier change right after a full render then reuses the full solve, // and the "preview" is silently upgraded to full quality. const solveCache = new Map(); const radiateCache = new Map(); const SOLVE_CACHE_MAX = 2; const RADIATE_CACHE_MAX = 3; function lruGet(cache: Map, key: string): T | undefined { const v = cache.get(key); if (v !== undefined) { cache.delete(key); cache.set(key, v); } return v; } function lruSet(cache: Map, key: string, value: T, max: number): void { cache.delete(key); cache.set(key, value); while (cache.size > max) { const oldest = cache.keys().next().value as string; cache.delete(oldest); } } function keyFor(params: MaracaParams, keys: readonly (keyof MaracaParams)[]): string { return keys.map((k) => String(params[k])).join('|'); } // MATLAB source lines assigning every parameter onto a struct P. function pLiteral(params: MaracaParams): string { const lines: string[] = []; for (const [k, v] of Object.entries(params)) { if (typeof v === 'number') lines.push(`P.${k} = ${v};`); else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`); } return lines.join('\n'); } function runStage(mfiles: MFile[], source: string, inputs?: Record) { return executeCode( source, { initialVariableValues: inputs, displayResults: false }, mfiles, 'stage.m', ); } function asTensor(v: RuntimeValue): RuntimeTensor { const t = v as RuntimeTensor; if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage'); return t; } function asScalar(v: RuntimeValue): number { if (typeof v === 'number') return v; const t = v as RuntimeTensor; if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0]; if (typeof v === 'boolean') return v ? 1 : 0; throw new Error('expected scalar result from numbl stage'); } function toF32(v: RuntimeValue): Float32Array { const t = asTensor(v); const out = new Float32Array(t.data.length); for (let i = 0; i < t.data.length; i++) out[i] = t.data[i]; return out; } const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/; self.onmessage = (ev: MessageEvent) => { const msg = ev.data; if (msg.type !== 'render') return; // numbl reports a file it could not parse through console.warn and then // skips it, so the failure surfaces at the call site as "Undefined // function". Capture the warning to report the actual syntax error instead. const warnings: string[] = []; const origWarn = console.warn; console.warn = (...args: unknown[]) => { warnings.push(args.map(String).join(' ')); }; try { render(msg); } catch (err) { // Most likely a fault in the (user-editable) .m code. Report it so the // request completes and the message reaches the code panel. const e = err as { message?: string; file?: string; line?: number }; let m: ErrorMessage = { type: 'error', id: msg.id, message: e?.message ? String(e.message) : String(err), file: typeof e?.file === 'string' ? e.file : undefined, line: typeof e?.line === 'number' ? e.line : undefined, }; for (const w of warnings) { const hit = SKIP_WARNING.exec(w); if (hit) { m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) }; break; } } (self as unknown as Worker).postMessage(m); } finally { console.warn = origWarn; } }; function render(msg: RenderRequest): void { const { id, params, mfiles } = msg; let quality = msg.quality; // The .m sources are part of the cache identity: editing the code must // invalidate everything computed from the previous version. const codeKey = hashMFiles(mfiles); const solveParamsKey = codeKey + '|' + keyFor(params, SOLVE_KEYS); const fullSolveKey = solveParamsKey + `|os=${params.oversample}`; if (quality === 'preview' && solveCache.has(fullSolveKey)) { quality = 'full'; // full-rate impact process already cached } const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample; const fsHi = params.sample_rate * oversample; const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m); // ---- solve stage: mode basis and the stochastic impact process ---- const solveKey = solveParamsKey + `|os=${oversample}`; let solveMs = 0; let solveCached = true; let solveEntry = lruGet(solveCache, solveKey); if (!solveEntry) { solveCached = false; post({ type: 'progress', id, stage: 'solve', frac: 0.15 }); const t0 = performance.now(); const r = runStage(mfiles, ` ${pLiteral(params)} basis = maraca_basis(P); [timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P); modeF = basis.f; modeL = basis.l; modeZeta = basis.zeta; modeSigma = basis.sigma; breathingHz = basis.breathingHz; breathingDryHz = basis.breathingDryHz; impactCount = ginfo.count; meanV = ginfo.meanV; maxV = ginfo.maxV; meanTauUs = ginfo.meanTauUs; durationS = ginfo.durationS; rateT = ginfo.rateT; rateHz = ginfo.rateHz; vimp = ginfo.vimp; `); solveMs = performance.now() - t0; const impactCount = asScalar(r.variableValues.impactCount); const allT = toF32(r.variableValues.timp); const allV = toF32(r.variableValues.vimp); const allX = toF32(r.variableValues.ximp); // thin uniformly for the raster plot, keeping the time ordering const stride = Math.max(1, Math.ceil(impactCount / MAX_PLOT_IMPACTS)); const nPlot = Math.ceil(impactCount / stride); const impactT = new Float32Array(nPlot); const impactV = new Float32Array(nPlot); const impactPos = new Float32Array(nPlot * 3); for (let i = 0, j = 0; j < nPlot; i += stride, j++) { impactT[j] = allT[i]; impactV[j] = allV[i]; impactPos[j * 3] = allX[i * 3]; impactPos[j * 3 + 1] = allX[i * 3 + 1]; impactPos[j * 3 + 2] = allX[i * 3 + 2]; } solveEntry = { timp: r.variableValues.timp, jimp: r.variableValues.jimp, tauimp: r.variableValues.tauimp, fimp: r.variableValues.fimp, ximp: r.variableValues.ximp, fsHi, impactT, impactV, impactPos, rateT: toF32(r.variableValues.rateT), rateHz: toF32(r.variableValues.rateHz), modeF: toF32(r.variableValues.modeF), modeL: toF32(r.variableValues.modeL), modeZeta: toF32(r.variableValues.modeZeta), modeSigma: toF32(r.variableValues.modeSigma), breathingHz: asScalar(r.variableValues.breathingHz), breathingDryHz: asScalar(r.variableValues.breathingDryHz), impactCount, meanV: asScalar(r.variableValues.meanV), maxV: asScalar(r.variableValues.maxV), meanTauUs: asScalar(r.variableValues.meanTauUs), durationS: asScalar(r.variableValues.durationS), }; lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX); } // ---- radiate stage ---- const radiateKey = solveKey + '||' + keyFor(params, RADIATE_KEYS); let radiateMs = 0; let radiateCached = true; let radEntry = lruGet(radiateCache, radiateKey); if (!radEntry) { radiateCached = false; post({ type: 'progress', id, stage: 'radiate', frac: 0.55 }); const t0 = performance.now(); const r = runStage(mfiles, ` ${pLiteral(params)} basis = maraca_basis(P); [p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, ${solveEntry.fsHi}); directPeakPa = radinfo.directPeakPa; listenerDistM = radinfo.listenerDistM; listenerInside = radinfo.listenerInside; `, { timp: solveEntry.timp, jimp: solveEntry.jimp, tauimp: solveEntry.tauimp, fimp: solveEntry.fimp, ximp: solveEntry.ximp, }); radiateMs = performance.now() - t0; radEntry = { p: r.variableValues.p, directPeakPa: asScalar(r.variableValues.directPeakPa), listenerDistM: asScalar(r.variableValues.listenerDistM), listenerInside: asScalar(r.variableValues.listenerInside) !== 0, }; lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX); } // ---- output stage ---- post({ type: 'progress', id, stage: 'output', frac: 0.9 }); const rOut = runStage(mfiles, ` ${pLiteral(params)} [y, peakPa] = maraca_output(P, p, ${oversample}); `, { p: radEntry.p }); const y32 = toF32(rOut.variableValues.y); const peakPa = asScalar(rOut.variableValues.peakPa); const warnings: string[] = []; if (radEntry.listenerInside) { warnings.push( 'the listener is at or inside the shell surface; the multipole field is evaluated at the nearest valid radius instead', ); } if (!params.normalize) { let peak = 0; for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i])); if (peak > 1.0) { warnings.push( `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`, ); } } const topMode = solveEntry.modeF[solveEntry.modeF.length - 1]; if (topMode < 0.4 * params.max_mode_hz) { warnings.push( `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max degree l to keep the shell from sounding truncated`, ); } const resp: RenderResponse = { type: 'result', id, quality, y: y32, sampleRate: params.sample_rate, impactT: solveEntry.impactT.slice(), impactV: solveEntry.impactV.slice(), impactPos: solveEntry.impactPos.slice(), rateT: solveEntry.rateT.slice(), rateHz: solveEntry.rateHz.slice(), modeF: solveEntry.modeF.slice(), modeL: solveEntry.modeL.slice(), modeZeta: solveEntry.modeZeta.slice(), modeSigma: solveEntry.modeSigma.slice(), diag: { modesKept: solveEntry.modeF.length, breathingHz: solveEntry.breathingHz, breathingDryHz: solveEntry.breathingDryHz, impactCount: solveEntry.impactCount, meanV: solveEntry.meanV, maxV: solveEntry.maxV, meanTauUs: solveEntry.meanTauUs, durationS: solveEntry.durationS, peakPa, peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30), directPeakPa: radEntry.directPeakPa, listenerDistM: radEntry.listenerDistM, warnings, solveMs, radiateMs, solveCached, radiateCached, oversample, }, }; (self as unknown as Worker).postMessage(resp, [ resp.y.buffer, resp.impactT.buffer, resp.impactV.buffer, resp.impactPos.buffer, resp.rateT.buffer, resp.rateHz.buffer, resp.modeF.buffer, resp.modeL.buffer, resp.modeZeta.buffer, resp.modeSigma.buffer, ]); }