1/**
2 * The fp32 transform round-trip floor, swept in lmax — the measurement
3 * docs/reduced-transforms.md Sec 5b prescribes before and
4 * after any change to summation order in the Legendre kernels.
5 *
6 * npx vite-node scripts/sht-accuracy.ts [--lmax 63,127,255] [--seed 42]
7 *
8 * For a band-limited spectrum q, analys(synth(q)) = q exactly (Gauss
9 * quadrature is exact for the band), so the relative round-trip error is the
10 * transforms' own fp32 round-off with no reference implementation in the
11 * loop. Sequential accumulation over l in the synthesis shows this floor
12 * growing roughly linearly in lmax; pairwise/compensated accumulation shows
13 * it near-flat. The per-degree profile says *where* the error lives (the
14 * high-l coefficients are the ones the alpha shifts and the l(l+1)
15 * eigenvalues amplify).
16 */
17import { ShtPlan, requestShtDevice, describeAdapter } from '../src/sht/sht.ts';
18import { gridForLmax, lmIndex, nlmCalc } from '../src/sht/layout.ts';
19import { randomSpectrum } from '../src/sht/reference.ts';
20import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
22const arg = (name: string): string | undefined => {
23 const i = process.argv.indexOf(`--${name}`);
24 return i >= 0 ? process.argv[i + 1] : undefined;
25};
26const LMAXES = (arg('lmax') ?? '63,127,255').split(',').map(Number);
27const SEED = Number(arg('seed') ?? 42);
29let runtime: string;
30try {
31 runtime = await installWebGpu();
32} catch (e) {
33 console.error(`sht-accuracy: ${errMsg(e)}\n${NO_ADAPTER_HINT}`);
34 process.exit(1);
35}
36const device = await requestShtDevice();
37console.log(`sht-accuracy — ${runtime}, ${await describeAdapter(device)}\n`);
38console.log(' lmax grid rel L2 roundtrip worst degree (rel)');
40for (const lmax of LMAXES) {
41 const { nlat, nphi } = gridForLmax(lmax, 1);
42 const cfg = { lmax, mmax: lmax, nlat, nphi };
43 const plan = await ShtPlan.create(device, cfg);
44 const q = randomSpectrum(cfg, SEED);
46 const grid = await plan.synth(q);
47 const back = await plan.analys(grid);
49 // Overall relative L2, and the same per degree — errors concentrate in l.
50 let num = 0;
51 let den = 0;
52 const nlm = nlmCalc(lmax, lmax);
53 const errL = new Float64Array(lmax + 1);
54 const magL = new Float64Array(lmax + 1);
55 for (let m = 0; m <= lmax; m++) {
56 for (let l = m; l <= lmax; l++) {
57 const i = lmIndex(lmax, l, m);
58 const dr = back[2 * i] - q[2 * i];
59 const di = back[2 * i + 1] - q[2 * i + 1];
60 const d2 = dr * dr + di * di;
61 const m2 = q[2 * i] ** 2 + q[2 * i + 1] ** 2;
62 num += d2;
63 den += m2;
64 errL[l] += d2;
65 magL[l] += m2;
66 }
67 }
68 let worstL = 0;
69 let worstRel = 0;
70 for (let l = 0; l <= lmax; l++) {
71 const rel = Math.sqrt(errL[l] / Math.max(magL[l], 1e-300));
72 if (rel > worstRel) {
73 worstRel = rel;
74 worstL = l;
75 }
76 }
77 console.log(
78 ` ${String(lmax).padEnd(6)} ${`${nlat}x${nphi}`.padEnd(11)} ` +
79 `${Math.sqrt(num / den).toExponential(3).padEnd(18)} ` +
80 `l=${worstL}: ${worstRel.toExponential(3)} (${nlm} coefficients)`,
81 );
82 plan.destroy();
83}