/**
 * GPU validation, run in a browser (scripts/test-gpu.mjs drives it headlessly):
 *
 *  1. near weights against the brute-force double precision reference, for the
 *     self pair and the near pairs of a few nodes;
 *  2. Green's identity S_k[du/dn] - D_k[u] = u/2 at all wall nodes through the
 *     full operator (far field plus near weights);
 *  3. the self-test problem: the field of a source outside the room, solved
 *     with the impedance equation and compared with the exact field inside.
 *
 * The same three checks then run on a coarse cymbal (178 patches), which uses
 * the interpolated geometry and the exterior equation: the near weights of a face
 * patch and a rim patch, Green's identity for a source outside the cymbal, and
 * the exterior self-test with sources spread inside the cymbal (?cymbal=0 skips it,
 * ?room=0 skips the rounded cube).
 *
 * ?bench=small|medium|large instead runs a named cymbal case of the CUDA proof of
 * concept (same surface, same unknowns) and reports its self-test error, the
 * scattered field at the listener, the reciprocity check and the timings, for a
 * direct comparison with run.py --problem cymbal.
 *
 * The surfaces are kept small so that the software WebGPU adapter can run them.
 */
import { roundedCube, type Surface } from '../src/surface.ts';
import { cymbal, CymbalShape, CYMBAL_DEFAULTS } from '../src/cymbal.ts';
import { selftestSources, selftestPoints as cymbalPoints, sourceSum, LISTENER, C_SOUND } from '../src/cymbalProblem.ts';
import { green, greenDn, selftestRhs, selftestPoints, relError, ALPHA, type CVec } from '../src/problem.ts';
import { refSelf, refLeaf, momentsToWeights } from '../src/reference.ts';
import { adapterDescription, requestSolverDevice } from '../src/gpu/device.ts';
import { Operator } from '../src/gpu/operator.ts';

declare global {
  interface Window {
    __RESULTS__?: { passed: number; failed: number; fatal?: string; lines: string[] };
  }
}

const logEl = document.getElementById('log')!;
const lines: string[] = [];
let passed = 0;
let failed = 0;
const log = (s: string, cls = '') => {
  lines.push(s);
  const span = document.createElement('span');
  span.textContent = s + '\n';
  if (cls) span.className = cls;
  logEl.append(span);
  console.log(s);
};
const check = (name: string, val: number, tol: number) => {
  const ok = val < tol;
  ok ? passed++ : failed++;
  log(`${ok ? 'ok  ' : 'FAIL'} ${name}: ${val.toExponential(2)} (tol ${tol.toExponential(0)})`, ok ? 'ok' : 'fail');
};

/** Near weights against the brute-force reference for three nodes of each patch in tps. */
async function checkWeights(op: Operator, surf: Surface, tps: number[], k: number, label: string) {
  const npts = surf.npts;
  const npols = surf.npols;
  const xs = surf.srcvals.subarray(0, 3 * npts);
  const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
  const { pairs, W } = await op.readWallWeights();
  const worst = new Float64Array(6);
  for (const tp of tps) {
    for (const loc of [0, Math.floor(npols / 2), npols - 2]) {
      const i = tp * npols + loc;
      const xt = [xs[i], xs[npts + i], xs[2 * npts + i]];
      const nt = [nrm[i], nrm[npts + i], nrm[2 * npts + i]];
      for (let pp = pairs.rowPtr[i]; pp < pairs.rowPtr[i + 1]; pp++) {
        const j = pairs.pairP[pp];
        const self = j === tp;
        const I = self ? refSelf(surf, j, surf.uvNodes[i], surf.uvNodes[npts + i], xt, nt, k) : refLeaf(surf, j, xt, nt, k);
        const Wref = momentsToWeights(surf, I);
        const errs: string[] = [];
        for (let mm = 0; mm < 6; mm++) {
          let num = 0;
          let den = 0;
          for (let l = 0; l < 2 * npols; l++) {
            const g = W[2 * ((pp * 6 + mm) * npols) + l];
            num += (g - Wref[mm][l]) ** 2;
            den += Wref[mm][l] ** 2;
          }
          const e = Math.sqrt(num / den);
          worst[mm] = Math.max(worst[mm], e);
          errs.push(e.toExponential(1));
        }
        const kind = ['plain', 'preimage', 'self'][pairs.pairF[8 * pp + 2]];
        const d = Math.hypot(xt[0] - surf.centroids[j], xt[1] - surf.centroids[surf.npatches + j], xt[2] - surf.centroids[2 * surf.npatches + j]) / surf.rads[j];
        log(`     node ${i} patch ${j} (${kind}, dist/rad ${d.toFixed(2)}): rel err per kernel ${errs.join(' ')}`);
      }
    }
  }
  for (let mm = 0; mm < 6; mm++) check(`${label}near weights, kernel ${mm}, worst relative error`, worst[mm], 1e-4);
}

async function main() {
  const t0 = performance.now();
  if (logEl.textContent === 'starting…') logEl.textContent = '';
  const gpu = await requestSolverDevice();
  const device = gpu.device;
  const lim = device.limits;
  log(`device: ${adapterDescription(gpu.adapter ?? null)}; maxStorageBufferBindingSize ${(lim.maxStorageBufferBindingSize / 2 ** 20).toFixed(0)} MB, workgroup storage ${lim.maxComputeWorkgroupStorageSize} B`);

  const params = new URLSearchParams(location.search);
  const bench = params.get('bench');
  if (bench) {
    await cymbalBench(gpu, bench);
    log(`done, ${((performance.now() - t0) / 1000).toFixed(1)} s`);
    window.__RESULTS__ = { passed, failed, lines };
    return;
  }
  if (params.get('room') !== '0') await roomTests(gpu, params);
  if (params.get('cymbal') !== '0') await cymbalTests(gpu);
  log(`${passed} passed, ${failed} failed, ${((performance.now() - t0) / 1000).toFixed(1)} s`, failed ? 'fail' : 'ok');
  window.__RESULTS__ = { passed, failed, lines };
}

/** The rounded cube: exact geometry, interior equation. */
async function roomTests(gpu: Awaited<ReturnType<typeof requestSolverDevice>>, params: URLSearchParams) {
  const m = Number(params.get('m') ?? 3);
  const norder = Number(params.get('norder') ?? 6);
  const k = Number(params.get('k') ?? 4);
  const a = 1;
  const p = Number(params.get('p') ?? 8);
  const surf = roundedCube(a, p, m, norder);
  const npts = surf.npts;
  const npols = surf.npols;
  log(`surface: p = ${p}, ${surf.npatches} patches of order ${norder}, ${npts} nodes; k = ${k}`);
  const xs = surf.srcvals.subarray(0, 3 * npts);
  const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);

  const lam = -0.2;
  const op = new Operator(gpu, surf, k, ALPHA, lam, 0, { progress: () => {} });
  const info = await op.setup();
  log(`operator: ${info.npairs} near pairs, ${info.nleaves} leaves, ${info.noversampled} oversampled nodes, weights ${(info.weightsBytes / 2 ** 20).toFixed(1)} MB, setup ${info.setupSeconds.toFixed(2)} s`);

  // 1. near weights against the reference, for three nodes of one patch
  await checkWeights(op, surf, [4], k, '');

  // 2. Green's identity through the operator
  {
    const xe = [1.7, 0.6, -0.5];
    const u = green(xs, xe, k);
    const un = greenDn(xs, nrm, xe, k);
    const S = await op.applySingle(4, un);
    const D = await op.applySingle(5, u);
    const res = new Float64Array(2 * npts);
    for (let i = 0; i < 2 * npts; i++) res[i] = S[i] - D[i] - u[i] / 2;
    let num = 0;
    let den = 0;
    for (let i = 0; i < npts; i++) {
      num += (res[2 * i] ** 2 + res[2 * i + 1] ** 2) * surf.wts[i];
      den += (u[2 * i] ** 2 + u[2 * i + 1] ** 2) * surf.wts[i];
    }
    check("Green's identity S[du/dn] - D[u] - u/2 on the walls (relative, weighted)", Math.sqrt(num / den), 1e-4);
  }

  // 3. the self-test problem
  {
    const xe = [1.7, 0.6, -0.5];
    const rhs = selftestRhs(xs, nrm, xe, k, lam);
    const iters: string[] = [];
    const sol = await op.solve(rhs, 1e-7, 100, (it, r) => iters.push(`${it}:${r.toExponential(1)}`));
    log(`     GMRES: ${sol.niter} iterations, residual ${sol.rres.toExponential(2)}, ${sol.seconds.toFixed(2)} s`);
    const xt = selftestPoints(30, a);
    const ut = await op.eval(xt, sol.x, sol.phi1);
    const uex = green(xt, xe, k);
    check('self-test: relative error at 30 interior points', relError(ut, uex), 2e-4);
    const ub = await op.evalOnSurface(sol.x, sol.phi1);
    const ubex = green(xs, xe, k);
    check('self-test: relative error on the walls', relError(ub, ubex), 2e-4);
  }

  op.destroy();
}

main().catch((e) => {
  const msg = e instanceof Error ? `${e.message}\n${e.stack ?? ''}` : String(e);
  log(`FATAL ${msg}`, 'fail');
  window.__RESULTS__ = { passed, failed: failed + 1, fatal: msg, lines };
});

/** The same checks on a coarse cymbal: interpolated geometry, exterior equation, lambda = 0. */
async function cymbalTests(gpu: Awaited<ReturnType<typeof requestSolverDevice>>) {
  const shape = new CymbalShape({ ...CYMBAL_DEFAULTS, tRim: 0.02 });
  const surf = cymbal(shape, 0.1, 5);
  const npts = surf.npts;
  const k = (2 * Math.PI * 1000) / 343;
  log(`cymbal surface: rim ${shape.par.tRim * 1000} mm, ${surf.npatches} patches of order ${surf.norder}, ${npts} nodes; k = ${k.toFixed(2)}`);
  const xs = surf.srcvals.subarray(0, 3 * npts);
  const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
  const op = new Operator(gpu, surf, k, ALPHA, 0, 1, { progress: () => {} });
  const info = await op.setup();
  log(`operator: ${info.npairs} near pairs, ${info.nleaves} leaves, weights ${(info.weightsBytes / 2 ** 20).toFixed(1)} MB, setup ${info.setupSeconds.toFixed(2)} s`);

  // a patch on the top face near the bell, and the first rim patch
  const nface = (surf.npatches - surf.layout.nrim) / 2;
  await checkWeights(op, surf, [surf.layout.mC ** 2 + 2, 2 * nface], k, 'cymbal: ');

  // Green's identity for a source outside the body: S[du/dn] - D[u] = u/2 inside it
  {
    const xe = [0.05, -0.1, 0.3];
    const u = green(xs, xe, k);
    const un = greenDn(xs, nrm, xe, k);
    const S = await op.applySingle(4, un);
    const D = await op.applySingle(5, u);
    let num = 0;
    let den = 0;
    for (let i = 0; i < npts; i++) {
      const rr = S[2 * i] - D[2 * i] - u[2 * i] / 2;
      const ri = S[2 * i + 1] - D[2 * i + 1] - u[2 * i + 1] / 2;
      num += (rr * rr + ri * ri) * surf.wts[i];
      den += (u[2 * i] ** 2 + u[2 * i + 1] ** 2) * surf.wts[i];
    }
    check("cymbal: Green's identity on the surface (relative, weighted)", Math.sqrt(num / den), 1e-3);
  }

  // the exterior self-test: sources spread inside the cymbal
  {
    const { x0s, ws } = selftestSources(shape);
    const rhs = sourceSum(xs, x0s, ws, k, nrm);
    const sol = await op.solve(rhs, 1e-7, 100);
    log(`     GMRES: ${sol.niter} iterations, residual ${sol.rres.toExponential(2)}, ${sol.seconds.toFixed(2)} s`);
    const xt = cymbalPoints(shape, 30);
    const ut = await op.eval(xt, sol.x, sol.phi1);
    check('cymbal: self-test, relative error at 30 points around the cymbal', relError(ut, sourceSum(xt, x0s, ws, k)), 2e-3);
    const ub = await op.evalOnSurface(sol.x, sol.phi1);
    check('cymbal: self-test, relative error on the surface', relError(ub, sourceSum(xs, x0s, ws, k)), 5e-3);
  }
  op.destroy();
}

/** A named cymbal case of the CUDA proof of concept, for comparison with run.py. */
async function cymbalBench(gpu: Awaited<ReturnType<typeof requestSolverDevice>>, name: string) {
  const cases: Record<string, { freq: number; h: number; norder: number }> = {
    small: { freq: 1000, h: 0.06, norder: 6 },
    medium: { freq: 3000, h: 0.04, norder: 8 },
    large: { freq: 8000, h: 0.03, norder: 8 },
  };
  const c = cases[name];
  if (!c) throw new Error(`no case ${name}`);
  const shape = new CymbalShape();
  const surf = cymbal(shape, c.h, c.norder);
  const npts = surf.npts;
  const k = (2 * Math.PI * c.freq) / C_SOUND;
  log(`cymbal ${name}: ${surf.npatches} patches of order ${c.norder}, ${npts} points, k = ${k.toFixed(2)}`);
  const xs = surf.srcvals.subarray(0, 3 * npts);
  const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
  const op = new Operator(gpu, surf, k, ALPHA, 0, 1, { progress: () => {} });
  const info = await op.setup();
  log(`setup (near quadrature) ${info.setupSeconds.toFixed(2)} s: ${info.npairs} near pairs, ${info.nleaves} leaves, weights ${(info.weightsBytes / 2 ** 20).toFixed(0)} MB`);
  const eps = 1e-6;
  {
    const { x0s, ws } = selftestSources(shape);
    const sol = await op.solve(sourceSum(xs, x0s, ws, k, nrm), eps, 400);
    const xt = cymbalPoints(shape, 50);
    const e1 = relError(await op.eval(xt, sol.x, sol.phi1), sourceSum(xt, x0s, ws, k));
    const e2 = relError(await op.evalOnSurface(sol.x, sol.phi1), sourceSum(xs, x0s, ws, k));
    log(`self-test: ${sol.niter} iterations, ${sol.seconds.toFixed(2)} s; relative error ${e1.toExponential(2)} off the surface, ${e2.toExponential(2)} on it`);
  }
  const B = [-0.15, 0.05, -0.12];
  const rhs = (x0: number[]) => {
    const g = greenDn(xs, nrm, x0, k);
    for (let i = 0; i < g.length; i++) g[i] = -g[i];
    return g;
  };
  const solL = await op.solve(rhs(LISTENER), eps, 400);
  log(`main problem: ${solL.niter} iterations, residual ${solL.rres.toExponential(2)}, solve ${solL.seconds.toFixed(2)} s`);
  const vL = await op.eval(Float64Array.from([LISTENER[0], B[0], LISTENER[1], B[1], LISTENER[2], B[2]]), solL.x, solL.phi1);
  log(`scattered field at the listener ${vL[0].toExponential(6)} ${vL[1] < 0 ? '-' : '+'} ${Math.abs(vL[1]).toExponential(6)}i`);
  const solB = await op.solve(rhs(B), eps, 400);
  const vB = await op.eval(Float64Array.from(LISTENER), solB.x, solB.phi1);
  const rec = Math.hypot(vL[2] - vB[0], vL[3] - vB[1]) / Math.hypot(vL[2], vL[3]);
  log(`reciprocity: at B due to the listener ${vL[2].toExponential(6)} ${vL[3] < 0 ? '-' : '+'} ${Math.abs(vL[3]).toExponential(6)}i, relative difference ${rec.toExponential(2)}`);
  op.destroy();
}