// Web Worker that renders the cymbal stroke 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   -> cymbal_basis (one shell eigensolve per nodal-diameter
//                   order; the expensive stage, cached on geometry, material,
//                   damping and discretization parameters)
//   radiate tier -> cymbal_strikes + cymbal_cascade + cymbal_radiate (cached
//                   on the stick, the gesture, the cascade, the listener and
//                   the synthesis rate)
//   output tier  -> cymbal_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 CymbalParams } from '../params';

export interface RenderRequest {
  type: 'render';
  id: number;
  params: CymbalParams;
  quality: 'preview' | 'full';
  mfiles: MFile[];
}

export interface RenderResponse {
  type: 'result';
  id: number;
  quality: 'preview' | 'full';
  y: Float32Array;               // final signal, [-1, 1]
  sampleRate: number;
  strikeT: Float32Array;         // stroke times, s
  strikeV: Float32Array;         // stroke velocities, m/s
  strikeS: Float32Array;         // stroke arclength fractions
  modeF: Float32Array;           // mode frequencies, Hz
  modeZeta: Float32Array;        // damping ratio of each mode
  modeZetaRad: Float32Array;     // radiation part of the damping
  modeSigma: Float32Array;       // radiation efficiency of each mode
  modeN: Float32Array;           // nodal diameters
  nodeR: Float32Array;           // meridian nodes, m
  nodeZ: Float32Array;
  shapeW: Float32Array;          // normal displacement per mode, nn x n col-major
  cascadeDb: Float32Array;       // energy grid, nBands x nFrames col-major, dB
  excDb: Float32Array;           // first stroke's modal excitation per band, dB
  cascadeBands: number;
  cascadeFrames: number;
  cascadeFLo: number;
  cascadeFHi: number;
  cascadeTEnd: number;
  diag: {
    modesKept: number;
    maxN: number;
    strikeCount: number;
    meanV: number;
    maxV: number;
    meanTauUs: number;
    peakForceN: number;
    durationS: number;
    peakPa: number;
    peakSplDb: number;
    directPeakPa: number;
    listenerDistM: number;
    strikeRadiusM: number;
    strikeRimM: number;
    roll2RadiusM: number;
    peakWrmsUm: number;
    peakStrength: number;
    cascadeLostFrac: number;
    plateMassKg: number;
    warnings: string[];
    solveMs: number;             // 0 when the cached eigensolve 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;
}

interface SolveEntry {
  vars: Record<string, RuntimeValue>;
  modeF: Float32Array;
  modeZeta: Float32Array;
  modeZetaRad: Float32Array;
  modeSigma: Float32Array;
  modeN: Float32Array;
  nodeR: Float32Array;
  nodeZ: Float32Array;
  shapeW: Float32Array;
  maxN: number;
  plateMassKg: number;
}

interface RadiateEntry {
  p: RuntimeValue;
  strikeT: Float32Array;
  strikeV: Float32Array;
  strikeS: Float32Array;
  cascadeDb: Float32Array;
  excDb: Float32Array;
  cascadeBands: number;
  cascadeFrames: number;
  cascadeTEnd: number;
  strikeRadiusM: number;
  strikeRimM: number;
  roll2RadiusM: number;
  strikeCount: number;
  meanV: number;
  maxV: number;
  meanTauUs: number;
  peakForceN: number;
  durationS: number;
  directPeakPa: number;
  listenerDistM: number;
  listenerInside: boolean;
  peakWrmsUm: number;
  peakStrength: number;
  cascadeLostFrac: number;
  oversample: number;
}

// Small LRU caches. The solve cache does not depend on the oversample rate,
// so a preview render right after a full one reuses the same eigensolves.
const solveCache = new Map<string, SolveEntry>();
const radiateCache = new Map<string, RadiateEntry>();
const SOLVE_CACHE_MAX = 2;
const RADIATE_CACHE_MAX = 3;

function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
  const v = cache.get(key);
  if (v !== undefined) {
    cache.delete(key);
    cache.set(key, v);
  }
  return v;
}

function lruSet<T>(cache: Map<string, T>, 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: CymbalParams, keys: readonly (keyof CymbalParams)[]): string {
  return keys.map((k) => String(params[k])).join('|');
}

// MATLAB source lines assigning every parameter onto a struct P.
function pLiteral(params: CymbalParams): 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<string, RuntimeValue>) {
  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;
}

// The basis fields handed from the solve stage into the radiate stage, in the
// order the reassembly script expects.
const BASIS_FIELDS = [
  'f', 'omega', 'zeta', 'zetaRad', 'zetaMount', 'sigma', 'ndia', 'ldip', 'mix',
  'NU', 'NW', 'NT',
  'rnode', 'znode', 'snode', 'drdsnode', 'dzdsnode', 'thknode', 'csnode', 'ctnode',
  'surfR', 'surfZ', 'surfDr', 'surfDz', 'surfLen', 'SW',
  'S', 'area', 'mass', 'thkMean',
];

const CASCADE_BANDS = 48;
const CASCADE_FRAMES = 240;

const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;

self.onmessage = (ev: MessageEvent<RenderRequest>) => {
  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) {
    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 solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);

  const fullOversample = params.oversample;
  const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
  if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {
    quality = 'full'; // full-rate radiation 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: one shell eigensolve per nodal-diameter order ----
  let solveMs = 0;
  let solveCached = true;
  let solveEntry = lruGet(solveCache, solveKey);
  if (!solveEntry) {
    solveCached = false;
    post({ type: 'progress', id, stage: 'solve', frac: 0.12 });
    const t0 = performance.now();
    const r = runStage(mfiles, `
${pLiteral(params)}
basis = cymbal_basis(P);
${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
maxN = max(basis.ndia);
`);
    solveMs = performance.now() - t0;

    const vars: Record<string, RuntimeValue> = {};
    for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
    solveEntry = {
      vars,
      modeF: toF32(r.variableValues.b_f),
      modeZeta: toF32(r.variableValues.b_zeta),
      modeZetaRad: toF32(r.variableValues.b_zetaRad),
      modeSigma: toF32(r.variableValues.b_sigma),
      modeN: toF32(r.variableValues.b_ndia),
      nodeR: toF32(r.variableValues.b_rnode),
      nodeZ: toF32(r.variableValues.b_znode),
      shapeW: toF32(r.variableValues.b_NW),
      maxN: asScalar(r.variableValues.maxN),
      plateMassKg: asScalar(r.variableValues.b_mass),
    };
    lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
  }

  // ---- radiate stage: strikes + cascade + modal synthesis ----
  const radiateKey = radiateParamsKey + `|os=${oversample}`;
  let radiateMs = 0;
  let radiateCached = true;
  let radEntry = lruGet(radiateCache, radiateKey);
  if (!radEntry) {
    radiateCached = false;
    post({ type: 'progress', id, stage: 'radiate', frac: 0.5 });
    const t0 = performance.now();
    const r = runStage(mfiles, `
${pLiteral(params)}
${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
basis.n = length(basis.f);
[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis);
[env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, sinfo.durationS);
[p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, ${fsHi});
% collapse the modal energies onto a log-frequency by time grid for the plot
nb = ${CASCADE_BANDS};
nt = ${CASCADE_FRAMES};
fLo = max(1, min(basis.f));
fHi = max(fLo * 1.2, max(basis.f));
bandOf = min(nb, max(1, 1 + floor(nb * (log(basis.f) - log(fLo)) / (log(fHi) - log(fLo) + 1e-12))));
nf = size(env, 2);
grid = zeros(nb, nt);
for q = 1:nt
  c0 = 1 + floor((q - 1) * nf / nt);
  c1 = max(c0, floor(q * nf / nt));
  e = sum(env(:, c0:c1).^2, 2) / (c1 - c0 + 1);
  grid(:, q) = accumarray(bandOf, e, [nb, 1]);
end
cascadeDb = 10 * log10(grid + 1e-30);
excDb = 10 * log10(accumarray(bandOf, 0.5 * DV(:, 1).^2, [nb, 1]) + 1e-30);
cascadeTEnd = tframe(length(tframe));
[strikeRadiusM, ~] = cymbal_point(basis, P.strike_r);
strikeRimM = P.radius_m - strikeRadiusM;
[roll2RadiusM, ~] = cymbal_point(basis, P.roll_r2);
strikeCount = sinfo.count;
meanV = sinfo.meanV;
maxV = sinfo.maxV;
meanTauUs = sinfo.meanTauUs;
peakForceN = sinfo.peakForceN;
durationS = sinfo.durationS;
directPeakPa = radinfo.directPeakPa;
listenerDistM = radinfo.listenerDistM;
listenerInside = radinfo.listenerInside;
peakWrmsUm = cinfo.peakWrmsUm;
peakStrength = cinfo.peakStrength;
cascadeLostFrac = cinfo.cascadeLostJ / max(1e-30, cinfo.injectedJ);
`, solveEntry.vars);
    radiateMs = performance.now() - t0;
    radEntry = {
      p: r.variableValues.p,
      strikeT: toF32(r.variableValues.tstr),
      strikeV: toF32(r.variableValues.vstr),
      strikeS: toF32(r.variableValues.sstr),
      cascadeDb: toF32(r.variableValues.cascadeDb),
      excDb: toF32(r.variableValues.excDb),
      cascadeBands: CASCADE_BANDS,
      cascadeFrames: CASCADE_FRAMES,
      cascadeTEnd: asScalar(r.variableValues.cascadeTEnd),
      strikeRadiusM: asScalar(r.variableValues.strikeRadiusM),
      strikeRimM: asScalar(r.variableValues.strikeRimM),
      roll2RadiusM: asScalar(r.variableValues.roll2RadiusM),
      strikeCount: asScalar(r.variableValues.strikeCount),
      meanV: asScalar(r.variableValues.meanV),
      maxV: asScalar(r.variableValues.maxV),
      meanTauUs: asScalar(r.variableValues.meanTauUs),
      peakForceN: asScalar(r.variableValues.peakForceN),
      durationS: asScalar(r.variableValues.durationS),
      directPeakPa: asScalar(r.variableValues.directPeakPa),
      listenerDistM: asScalar(r.variableValues.listenerDistM),
      listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
      peakWrmsUm: asScalar(r.variableValues.peakWrmsUm),
      peakStrength: asScalar(r.variableValues.peakStrength),
      cascadeLostFrac: asScalar(r.variableValues.cascadeLostFrac),
      oversample,
    };
    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] = cymbal_output(P, p, ${radEntry.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 within 2 cm of the plate; the dipole field is evaluated at a clamped minimum distance',
    );
  }
  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)}`,
      );
    }
  }
  if (params.n_azimuth < 2.5 * solveEntry.maxN) {
    warnings.push(
      `only ${params.n_azimuth} azimuth samples for ${solveEntry.maxN} nodal diameters; the listener transfer of the high-order modes is aliased`,
    );
  }
  if (radEntry.cascadeLostFrac > 0.25) {
    warnings.push(
      `${(100 * radEntry.cascadeLostFrac).toFixed(0)}% of the struck energy left the top of the retained band; raise max_mode_hz or lower cascade_leak`,
    );
  }
  const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
  if (topMode < 0.6 * params.max_mode_hz) {
    warnings.push(
      `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max nodal diameters to fill the band`,
    );
  }

  const resp: RenderResponse = {
    type: 'result',
    id,
    quality,
    y: y32,
    sampleRate: params.sample_rate,
    strikeT: radEntry.strikeT.slice(),
    strikeV: radEntry.strikeV.slice(),
    strikeS: radEntry.strikeS.slice(),
    modeF: solveEntry.modeF.slice(),
    modeZeta: solveEntry.modeZeta.slice(),
    modeZetaRad: solveEntry.modeZetaRad.slice(),
    modeSigma: solveEntry.modeSigma.slice(),
    modeN: solveEntry.modeN.slice(),
    nodeR: solveEntry.nodeR.slice(),
    nodeZ: solveEntry.nodeZ.slice(),
    shapeW: solveEntry.shapeW.slice(),
    cascadeDb: radEntry.cascadeDb.slice(),
    excDb: radEntry.excDb.slice(),
    cascadeBands: radEntry.cascadeBands,
    cascadeFrames: radEntry.cascadeFrames,
    cascadeFLo: solveEntry.modeF[0],
    cascadeFHi: solveEntry.modeF[solveEntry.modeF.length - 1],
    cascadeTEnd: radEntry.cascadeTEnd,
    diag: {
      modesKept: solveEntry.modeF.length,
      maxN: solveEntry.maxN,
      strikeCount: radEntry.strikeCount,
      meanV: radEntry.meanV,
      maxV: radEntry.maxV,
      meanTauUs: radEntry.meanTauUs,
      peakForceN: radEntry.peakForceN,
      durationS: radEntry.durationS,
      peakPa,
      peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
      directPeakPa: radEntry.directPeakPa,
      listenerDistM: radEntry.listenerDistM,
      strikeRadiusM: radEntry.strikeRadiusM,
      strikeRimM: radEntry.strikeRimM,
      roll2RadiusM: radEntry.roll2RadiusM,
      peakWrmsUm: radEntry.peakWrmsUm,
      peakStrength: radEntry.peakStrength,
      cascadeLostFrac: radEntry.cascadeLostFrac,
      plateMassKg: solveEntry.plateMassKg,
      warnings,
      solveMs,
      radiateMs,
      solveCached,
      radiateCached,
      oversample: radEntry.oversample,
    },
  };
  (self as unknown as Worker).postMessage(resp, [
    resp.y.buffer,
    resp.strikeT.buffer,
    resp.strikeV.buffer,
    resp.strikeS.buffer,
    resp.modeF.buffer,
    resp.modeZeta.buffer,
    resp.modeZetaRad.buffer,
    resp.modeSigma.buffer,
    resp.modeN.buffer,
    resp.nodeR.buffer,
    resp.nodeZ.buffer,
    resp.shapeW.buffer,
    resp.cascadeDb.buffer,
    resp.excDb.buffer,
  ]);
}