1/**
2 * The WGSL spherical-harmonic transforms against the f64 CPU reference.
3 *
4 * This is the one place a second implementation is still the right oracle: the
5 * transforms are vendored shtns-webgpu, and `src/sht/reference.ts` is its direct-
6 * summation f64 twin. Everything above them (the .m models) is checked against
7 * closed-form answers instead — see analyticChecks.ts.
8 */
9import { ShtPlan } from '../src/sht/sht.ts';
10import { ShtReference, randomSpectrum } from '../src/sht/reference.ts';
11import { gridForLmax } from '../src/sht/layout.ts';
12import type { Check, Log } from './analyticChecks.ts';
14function relL2(a: ArrayLike<number>, b: ArrayLike<number>): number {
15 let num = 0;
16 let den = 0;
17 for (let i = 0; i < a.length; i++) {
18 const d = a[i] - b[i];
19 num += d * d;
20 den += b[i] * b[i];
21 }
22 return Math.sqrt(num / Math.max(den, 1e-300));
23}
25export async function transformChecks(
26 device: GPUDevice,
27 check: Check,
28 _log: Log,
29): Promise<void> {
30 const lmax = 31;
31 const { nlat, nphi } = gridForLmax(lmax, 1);
32 const cfg = { lmax, mmax: lmax, nlat, nphi };
34 const plan = await ShtPlan.create(device, cfg);
35 const ref = new ShtReference(cfg);
37 const q = randomSpectrum(cfg, 42);
38 const q64 = new Float64Array(q);
40 const spatGpu = await plan.synth(new Float32Array(q64));
41 const spatCpu = ref.synth(q64);
42 const errSynth = relL2(spatGpu, spatCpu);
44 const qGpu = await plan.analys(new Float32Array(spatCpu));
45 const qCpu = ref.analys(new Float64Array(spatCpu));
46 const errAnalys = relL2(qGpu, qCpu);
48 check(
49 'transforms: WGSL fp32 vs f64 CPU reference',
50 errSynth < 1e-4 && errAnalys < 1e-4,
51 `synth ${errSynth.toExponential(2)}, analys ${errAnalys.toExponential(2)}`,
52 );
54 plan.destroy();
55}