Comparing changes
main is 2 commits ahead of compare-interface.
Create pull request
export standalone matlab
Jeremy Magland committed
9b71d5bMerge pull request #10 from concept-collection/compare-interface
Jeremy Magland committed
71b64929 changed files+951−2
.gitignoremodified+1−0View file
@@ -1,3 +1,4 @@
1+tmp/
12 node_modules/
23 dist/
34 *.log
README.mdmodified+6−0View file
@@ -497,6 +497,12 @@ package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
497497 `--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
498498 (`;`-separated) passes Dawn options through.
499499
500+### The same run in MATLAB
501+
502+A run in the page needs a browser and a GPU; further analysis usually wants neither. The app therefore exports the run on screen as one self-contained MATLAB function file: **The same run as a standalone MATLAB script**, under the benchmark command, shows the script for copying and downloads it as `turing_surface_run.m`. The current model and geometry `.m` go in verbatim, edits in the page included, with the parameter values baked in; around them the file carries double-precision ports of everything the host provides: the transforms and their derivative shuffles, the metric weights, the seeded random field, and the run loop ([`src/export/`](src/export/)). It needs base MATLAB only, R2020b or newer, no toolboxes.
503+
504+Two deliberate differences from the page are stated in the script's own header: it runs in f64 where the GPU path is f32, and random draws use MATLAB's own `rng`, so a seed value picks a different member of the same random ensemble than the same value in the app. The script plots the pattern live and writes its initial and final spectral state to HDF5 in the reference-run layout of [docs/ellipsoid-reference-spec.md](docs/ellipsoid-reference-spec.md), so a MATLAB run can be loaded back into the page (**Compare against uploaded data**) or checked with `npm run ref -- --in turing_surface_run.h5`. Exported at the defaults, a 60-step Schnakenberg run on the ellipsoid replayed that way agrees with the app to relative L2 of about 1e-7, which is fp32 accumulation; the exported flux-form and Algorithm-4 models track each other to about 3e-10 in f64.
505+
500506 ## Tests
501507
502508 There is no second implementation of the solver to diff against, so the `.m`
index.htmlmodified+18−0View file
@@ -117,6 +117,14 @@
117117 font: 12px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
118118 color: var(--ink); user-select: all;
119119 }
120+ details.cli > summary { cursor: pointer; }
121+ details.cli > summary::before { content: '▸'; font-size: 10px; color: var(--ink-2); }
122+ details.cli[open] > summary::before { content: '▾'; }
123+ #matlabscript {
124+ margin: 0; padding: 8px 10px; max-height: 30em; overflow: auto;
125+ font: 12px/1.5 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
126+ color: var(--ink); white-space: pre; user-select: text;
127+ }
120128 #blurb { margin-top: 4px; font-size: 13px; color: var(--ink-2); }
121129 #err { color: #b35900; white-space: pre-wrap; font-size: 13px; }
122130 .editor {
@@ -414,6 +422,16 @@
414422 </div>
415423 <code id="cmd"></code>
416424 </div>
425+ <details class="cli" id="matlab">
426+ <summary class="cli-head">
427+ <span>The same run as a standalone MATLAB script</span>
428+ <span>
429+ <button id="copymatlab" type="button">Copy</button>
430+ <button id="downloadmatlab" type="button">Download .m</button>
431+ </span>
432+ </summary>
433+ <pre id="matlabscript"></pre>
434+ </details>
417435 <p id="blurb"></p>
418436 <p id="err"></p>
419437 </main>
scripts/test-node.tsmodified+2−0View file
@@ -20,6 +20,7 @@ import { geometryChecks } from '../test/geometryChecks.ts';
2020 import { fluxChecks } from '../test/fluxChecks.ts';
2121 import { compareChecks } from '../test/compareChecks.ts';
2222 import { referenceChecks, type H5Rt } from '../test/referenceChecks.ts';
23+import { matlabExportChecks } from '../test/matlabExportChecks.ts';
2324
2425 let failures = 0;
2526 const check = (name: string, ok: boolean, detail: string): void => {
@@ -59,6 +60,7 @@ await modelChecks(device, check, log);
5960 await geometryChecks(device, check, log);
6061 await fluxChecks(device, check, log);
6162 await compareChecks(device, check, log);
63+matlabExportChecks(check, log);
6264 await referenceChecks(
6365 h5wasm as unknown as H5Rt,
6466 (name) => join(tmpdir(), `turing-surface-${process.pid}-${name}`),
src/export/matlabScript.tsadded+313−0View file
@@ -0,0 +1,313 @@
1+/**
2+ * The run on screen, as one standalone MATLAB script.
3+ *
4+ * The models and geometries are already MATLAB; what the app supplies around
5+ * them — the transforms, the geometry weights, the seeded field, the driver —
6+ * exists only as TypeScript and WGSL. This module assembles a single function
7+ * file carrying all of it: the current model and geometry sources verbatim as
8+ * local functions, double-precision MATLAB ports of the host-provided
9+ * operations (support.m), and a generated driver with the run's settings
10+ * baked in.
11+ *
12+ * Fidelity is method-for-method, not bit-for-bit: the ports run in f64 where
13+ * the GPU path is f32, and random draws use MATLAB's own rng, so a seed value
14+ * selects a different member of the same random ensemble than the same value
15+ * in the app. The script's results file uses the app's reference-run layout
16+ * (docs/ellipsoid-reference-spec.md), so a MATLAB run can be loaded back into
17+ * the page or checked with `npm run ref`.
18+ */
19+import supportSource from './support.m?raw';
20+import randnfun3Source from '../../tools/randnfun3.m?raw';
21+import randnfunsphereSource from '../../tools/randnfunsphere.m?raw';
22+import type { MModel, Params } from '../mgpu/registry.ts';
23+import type { MGeometry } from '../geom/registry.ts';
24+
25+/** The generated function's name — and therefore the file name to save as. */
26+export const MATLAB_SCRIPT_NAME = 'turing_surface_run';
27+
28+export interface MatlabExportSpec {
29+ model: MModel;
30+ /** Model source as running — the editor's working copy when edited. */
31+ modelSource: string;
32+ params: Params;
33+ geometry: MGeometry;
34+ geometrySource: string;
35+ geometryParams: Params;
36+ lmax: number;
37+ niter: number;
38+ /** Wavelength of the seeded random field. */
39+ lam3: number;
40+ seed: number;
41+ /** Preset key, recorded in the results file's /spec. */
42+ preset: string;
43+ /** The equivalent `npm run bench` command, recorded for provenance. */
44+ command: string;
45+ /** The generated script's own run controls; app defaults when omitted. */
46+ controls?: { nsteps?: number; plotEvery?: number; outFile?: string };
47+}
48+
49+interface Signature {
50+ outputs: string[];
51+ params: string[];
52+}
53+
54+/** First `function [outs] = name(args)` line in a .m — the same contract the
55+ * compiler applies, minus everything it checks later. */
56+function parseSignature(source: string, name: string, file: string): Signature {
57+ const re = new RegExp(
58+ String.raw`^[ \t]*function\s+(?:\[([^\]]*)\]|([A-Za-z]\w*))\s*=\s*${name}\s*\(([^)]*)\)`,
59+ 'm',
60+ );
61+ const m = re.exec(source);
62+ if (!m) {
63+ throw new Error(`cannot export: ${file} defines no function named '${name}'`);
64+ }
65+ const split = (s: string): string[] =>
66+ s.split(',').map((t) => t.trim()).filter((t) => t.length > 0);
67+ return {
68+ outputs: m[1] !== undefined ? split(m[1]) : [m[2]],
69+ params: split(m[3]),
70+ };
71+}
72+
73+/** A number as MATLAB source. JS stringification round-trips doubles exactly
74+ * and every form it produces (0.0004, 1e-21, -3) is a MATLAB literal. */
75+const num = (v: number): string => (Number.isFinite(v) ? String(v) : '0');
76+
77+/** A string as a MATLAB char literal. */
78+const str = (s: string): string => `'${s.replace(/'/g, "''")}'`;
79+
80+const banner = (title: string): string => {
81+ const line = `% ${'='.repeat(72)}`;
82+ return `${line}\n% ${title}\n${line}`;
83+};
84+
85+export function generateMatlabScript(spec: MatlabExportSpec): string {
86+ const { model, geometry } = spec;
87+ const controls = {
88+ nsteps: spec.controls?.nsteps ?? 2000,
89+ plotEvery: spec.controls?.plotEvery ?? 10,
90+ outFile: spec.controls?.outFile ?? `${MATLAB_SCRIPT_NAME}.h5`,
91+ };
92+
93+ const init = parseSignature(spec.modelSource, 'init', `models/${model.key}.m`);
94+ const step = parseSignature(spec.modelSource, 'step', `models/${model.key}.m`);
95+ const shape = parseSignature(spec.geometrySource, 'shape', `geometries/${geometry.key}.m`);
96+
97+ // The driver defines every host-provided name the .m may ask for (lam,
98+ // filt, the geometry fields, jhat, niter, ...) under its canonical name,
99+ // so a model call is its own signature read back. Only the tunable
100+ // parameters live elsewhere — in the mp/gp structs, where the person
101+ // running the script edits them — so those names are mapped.
102+ const modelParamKeys = new Set(model.params.map((p) => p.key));
103+ const modelArg = (a: string): string => (modelParamKeys.has(a) ? `mp.${a}` : a);
104+ const shapeArg = (a: string): string =>
105+ a === 'theta' || a === 'phi' ? a : `gp.${a}`;
106+
107+ const stateOuts = [...model.state, ...model.species];
108+ const outs = `[${stateOuts.join(', ')}]`;
109+ const initCall = `${outs} = init(${init.params.map(modelArg).join(', ')});`;
110+ const stepCall = `${outs} = step(${step.params.map(modelArg).join(', ')});`;
111+ const shapeCall = `[gxr, gyr, gzr] = shape(${shape.params.map(shapeArg).join(', ')});`;
112+
113+ const speciesCell = `{${model.species.join(', ')}}`;
114+ const namesCell = `{${model.species.map((s) => str(s)).join(', ')}}`;
115+
116+ // `noise` is the plain seeded grid perturbation, for a .m that takes it
117+ // instead of calling randnfun3 (none of the shipped models do).
118+ const takesNoise = init.params.includes('noise') || step.params.includes('noise');
119+
120+ const mpBlock = model.params
121+ .map((p) => `mp.${p.key} = ${num(spec.params[p.key] ?? p.value)};`)
122+ .join('\n');
123+ const gpBlock = geometry.params
124+ .map((p) => `gp.${p.key} = ${num(spec.geometryParams[p.key] ?? p.value)};`)
125+ .join('\n');
126+
127+ const driver = `function ${MATLAB_SCRIPT_NAME}()
128+% ${model.label} on ${geometry.label} -- a run captured from the
129+% turing-surface app as one standalone MATLAB script.
130+%
131+% The model and geometry .m below are the app's own, verbatim; around them
132+% this file carries double-precision MATLAB ports of everything the app
133+% provides from the host side: the spherical-harmonic transforms and their
134+% derivative shuffles, the metric weights of the surface Laplace-Beltrami
135+% operator, the seeded random field, and the run loop (src/sht and src/geom
136+% in the repository). The scheme is the app's: IMEX Euler, implicit
137+% diffusion preconditioned on the round sphere, the geometric correction
138+% iterated niter times per step.
139+%
140+% Two deliberate differences from the page. Everything here runs in double
141+% precision, where the app's GPU path is single. And random draws use
142+% MATLAB's own rng, so a seed value selects a different member of the same
143+% random ensemble than the same value in the app.
144+%
145+% Save as ${MATLAB_SCRIPT_NAME}.m and run it. The run plots live, and the
146+% final state is written to an HDF5 file in the app's reference-run layout
147+% (docs/ellipsoid-reference-spec.md in the repository), so it can be loaded
148+% back into the page ("Compare against uploaded data") or checked on the
149+% desktop with \`npm run ref -- --in ${controls.outFile}\`.
150+% Needs base MATLAB, R2020b or newer; no toolboxes.
151+
152+% ---- run controls --------------------------------------------------------
153+nsteps = ${controls.nsteps}; % timesteps to run
154+plot_every = ${controls.plotEvery}; % live-plot interval, in steps; 0 disables plotting
155+out_file = ${str(controls.outFile)}; % results file; '' disables
156+seed = ${num(spec.seed)}; % rng seed for the initial condition
157+
158+% ---- captured from the app -----------------------------------------------
159+lmax = ${spec.lmax}; % spherical-harmonic truncation degree
160+niter = ${spec.niter}; % iterations of the implicit solve's geometric correction
161+lam3 = ${num(spec.lam3)}; % wavelength of the seeded random field
162+${model.params.length ? `% ${model.label} parameters\n${mpBlock}` : `% ${model.label} has no parameters`}
163+${geometry.params.length ? `% ${geometry.label} parameters\n${gpBlock}` : `% ${geometry.label} has no parameters`}
164+
165+% ---- grid and transforms -------------------------------------------------
166+% nlat/nphi follow lmax by the app's dealiasing rule (src/sht/layout.ts) for
167+% a reaction of polynomial degree pdeg.
168+pdeg = ${model.pdeg};
169+mmax = lmax;
170+nlat = 2 * ceil(max(lmax + 1, ((pdeg + 1) * lmax + 1) / 2) / 2);
171+nphi = 2 ^ nextpow2((pdeg + 1) * lmax + 1);
172+npts = nlat * nphi;
173+sht_tables(sht_setup(lmax, mmax, nlat, nphi));
174+S = sht_tables();
175+nlm = S.nlm;
176+lam = S.lam;
177+filt = S.filt;
178+theta = S.theta;
179+phi = S.phi;
180+
181+% ---- the surface ---------------------------------------------------------
182+${shapeCall}
183+% A constant coordinate comes back scalar; spread it over the grid.
184+gxr = gxr + zeros(npts, 1);
185+gyr = gyr + zeros(npts, 1);
186+gzr = gzr + zeros(npts, 1);
187+G = surface_tables(gxr, gyr, gzr);
188+gx = G.gx; gy = G.gy; gz = G.gz;
189+Gx = G.Gx; Gy = G.Gy; Gz = G.Gz;
190+p1 = G.p1; p2 = G.p2; q2 = G.q2; r = G.r;
191+dp1 = G.dp1; dq2 = G.dq2; jinv = G.jinv;
192+Vtx = G.Vtx; Vty = G.Vty; Vtz = G.Vtz;
193+Vpx = G.Vpx; Vpy = G.Vpy; Vpz = G.Vpz;
194+jhat = G.Jhat;
195+radius = sqrt(gx.^2 + gy.^2 + gz.^2);
196+fprintf('grid %d x %d, nlm %d, radius %.3f-%.3f, Jhat %.3f\\n', ...
197+ nlat, nphi, nlm, min(radius), max(radius), jhat);
198+
199+% ---- initial condition ---------------------------------------------------
200+rng(seed);
201+${takesNoise ? `noise = ${num(model.seedAmp)} * randn(npts, 1);\n` : ''}${initCall}
202+${model.state.map((s) => `${s}0 = ${s};`).join('\n')}
203+
204+% ---- time loop -----------------------------------------------------------
205+if plot_every > 0
206+ ph = plot_setup(gx, gy, gz, ${speciesCell}, ${namesCell});
207+ plot_update(ph, ${speciesCell}, 0, 0, nsteps);
208+end
209+report_every = max(1, round(nsteps / 10));
210+t = 0;
211+tstart = tic;
212+for k = 1:nsteps
213+ ${stepCall}
214+ t = t + mp.dt;
215+ if plot_every > 0 && (mod(k, plot_every) == 0 || k == nsteps)
216+ plot_update(ph, ${speciesCell}, t, k, nsteps);
217+ end
218+ if mod(k, report_every) == 0 || k == nsteps
219+ fprintf('step %d/%d t = %.3f (%.1f s)\\n', k, nsteps, t, toc(tstart));
220+ end
221+end
222+
223+% ---- results file --------------------------------------------------------
224+% The app's reference-run layout, plus a /fields group with the final grid
225+% fields, the surface and the grid angles (each field stored nphi x nlat,
226+% ring by ring from the north pole).
227+if ~isempty(out_file)
228+ if exist(out_file, 'file') == 2
229+ delete(out_file);
230+ end
231+${['Gx', 'Gy', 'Gz']
232+ .map((c) => ` write_coeffs(out_file, '/geometry/${c}', ${c});`)
233+ .join('\n')}
234+${model.state
235+ .map((s) => ` write_coeffs(out_file, '/initial/${s}', ${s}0);`)
236+ .join('\n')}
237+${model.state
238+ .map((s) => ` write_coeffs(out_file, '/final/${s}', ${s});`)
239+ .join('\n')}
240+${[...model.species.map((s) => [s, s] as const), (['x', 'gx'] as const), (['y', 'gy'] as const), (['z', 'gz'] as const)]
241+ .map(
242+ ([name, v]) =>
243+ ` h5create(out_file, '/fields/${name}', [nphi nlat]);\n` +
244+ ` h5write(out_file, '/fields/${name}', reshape(${v}, nphi, nlat));`,
245+ )
246+ .join('\n')}
247+ h5create(out_file, '/fields/theta', nlat);
248+ h5write(out_file, '/fields/theta', acos(min(1, max(-1, S.ct))));
249+ h5create(out_file, '/fields/phi', nphi);
250+ h5write(out_file, '/fields/phi', 2*pi*(0:nphi-1)'/nphi);
251+ make_group(out_file, '/backend');
252+ make_group(out_file, '/spec');
253+ make_group(out_file, '/spec/params');
254+ make_group(out_file, '/spec/geometry_params');
255+ make_group(out_file, '/grid');
256+ h5writeatt(out_file, '/', 'model', ${str(model.key)});
257+ h5writeatt(out_file, '/', 'species', [${model.state.map((s) => `"${s}"`).join(' ')}]);
258+ h5writeatt(out_file, '/', 'command', ${str(spec.command)});
259+ h5writeatt(out_file, '/backend', 'runtime', 'matlab');
260+ h5writeatt(out_file, '/backend', 'adapter', ['MATLAB ' version]);
261+ h5writeatt(out_file, '/backend', 'precision', 'double');
262+ h5writeatt(out_file, '/spec', 'preset', ${str(spec.preset)});
263+ h5writeatt(out_file, '/spec', 'geometry', ${str(geometry.key)});
264+ h5writeatt(out_file, '/spec', 'lmax', lmax);
265+ h5writeatt(out_file, '/spec', 'seed', seed);
266+ h5writeatt(out_file, '/spec', 'steps', nsteps);
267+ h5writeatt(out_file, '/spec', 'warmup', 0);
268+ h5writeatt(out_file, '/spec', 'niter', niter);
269+ h5writeatt(out_file, '/spec', 'lam3', lam3);
270+${model.params
271+ .map((p) => ` h5writeatt(out_file, '/spec/params', ${str(p.key)}, mp.${p.key});`)
272+ .join('\n')}
273+${geometry.params
274+ .map((p) => ` h5writeatt(out_file, '/spec/geometry_params', ${str(p.key)}, gp.${p.key});`)
275+ .join('\n')}
276+ h5writeatt(out_file, '/grid', 'lmax', lmax);
277+ h5writeatt(out_file, '/grid', 'mmax', mmax);
278+ h5writeatt(out_file, '/grid', 'nlat', nlat);
279+ h5writeatt(out_file, '/grid', 'nphi', nphi);
280+ h5writeatt(out_file, '/grid', 'nlm', nlm);
281+ fprintf('wrote %s\\n', out_file);
282+end
283+end`;
284+
285+ // tools/randnfun3.m verbatim, renamed: the models call the app's builtin
286+ // `randnfun3(lam3, gx, gy, gz)`, which support.m provides as a dispatcher
287+ // over this mode draw.
288+ const modesSource = randnfun3Source.replace(
289+ /function\s*\[\s*k\s*,\s*c\s*\]\s*=\s*randnfun3\s*\(/,
290+ 'function [k, c] = randnfun3_modes(',
291+ );
292+ if (modesSource === randnfun3Source) {
293+ throw new Error('cannot export: tools/randnfun3.m no longer matches the expected signature');
294+ }
295+
296+ const usesSphere = /\brandnfunsphere\b/.test(spec.geometrySource + spec.modelSource);
297+
298+ const parts = [
299+ driver,
300+ banner(`models/${model.key}.m -- the model, verbatim`),
301+ spec.modelSource.trim(),
302+ banner(`geometries/${geometry.key}.m -- the surface, verbatim`),
303+ spec.geometrySource.trim(),
304+ banner('tools/randnfun3.m -- the random-field mode draw, verbatim'),
305+ modesSource.trim(),
306+ ...(usesSphere
307+ ? [banner('tools/randnfunsphere.m -- verbatim'), randnfunsphereSource.trim()]
308+ : []),
309+ banner('host-provided operations, ported from src/sht and src/geom'),
310+ supportSource.trim(),
311+ ];
312+ return parts.join('\n\n') + '\n';
313+}
src/export/support.madded+388−0View file
@@ -0,0 +1,388 @@
1+% ---------------------------------------------------------------- transforms
2+%
3+% Double-precision MATLAB ports of the operations the app provides to a .m
4+% around its compiled GPU pipeline. Conventions follow src/sht/layout.ts:
5+% orthonormal spherical harmonics with the Condon-Shortley phase, coefficients
6+% stored for m >= 0 only in m-major order (m = 0..mmax, l = m..lmax within
7+% each m) -- here as complex nlm x 1 column vectors where the GPU carries
8+% interleaved [re, im] pairs. Grid fields are npts x 1 columns, phi-fastest:
9+% point (itheta, iphi) sits at row (itheta-1)*nphi + iphi, north row first.
10+
11+% Holds the precomputed tables between calls: set once from the top of the
12+% run, read back by every transform below.
13+function S = sht_tables(S)
14+ persistent stored
15+ if nargin > 0
16+ stored = S;
17+ end
18+ S = stored;
19+end
20+
21+% Everything the transforms need for one grid: Gauss nodes and weights,
22+% per-m Legendre tables, the coefficient layout, the derivative shuffles,
23+% and the eigenvalue/filter vectors the models take as `lam` and `filt`.
24+function S = sht_setup(lmax, mmax, nlat, nphi)
25+ S.lmax = lmax;
26+ S.mmax = mmax;
27+ S.nlat = nlat;
28+ S.nphi = nphi;
29+ S.npts = nlat * nphi;
30+ [ct, wg] = gauss_legendre(nlat);
31+ S.ct = ct;
32+ S.st = sqrt(1 - ct.^2);
33+ S.wg = wg;
34+ S.nlm = (mmax + 1) * (lmax + 1) - mmax * (mmax + 1) / 2;
35+
36+ % The grid angles as npts x 1 fields, phi-fastest like everything else.
37+ S.theta = repelem(acos(min(1, max(-1, ct))), nphi);
38+ S.phi = repmat(2*pi*(0:nphi-1)'/nphi, nlat, 1);
39+ S.stpt = repelem(S.st, nphi);
40+
41+ % Degree and order of each coefficient, and each m block's start.
42+ off = zeros(mmax + 1, 1);
43+ lv = zeros(S.nlm, 1);
44+ mv = zeros(S.nlm, 1);
45+ pos = 1;
46+ for m = 0:mmax
47+ n = lmax - m + 1;
48+ off(m + 1) = pos;
49+ lv(pos:pos + n - 1) = (m:lmax)';
50+ mv(pos:pos + n - 1) = m;
51+ pos = pos + n;
52+ end
53+ S.off = off;
54+ S.lv = lv;
55+ S.mv = mv;
56+ % Laplace-Beltrami eigenvalues l(l+1) and the top-mode filter: 1 below
57+ % lmax-2, 0 at the top two degrees, where the derivative recurrences cannot
58+ % exactly represent a derivative (src/mgpu/model.ts).
59+ S.lam = lv .* (lv + 1);
60+ S.filt = double(lv < lmax - 2);
61+
62+ % Orthonormal Legendre tables ytilde_l^m(theta_i), one nlat x (lmax-m+1)
63+ % block per m, by the standard three-term recurrence (src/sht/coeffs.ts;
64+ % SHTNS normalization, Condon-Shortley phase carried in the seed's sign).
65+ S.Y = cell(mmax + 1, 1);
66+ t = 1 / (4*pi);
67+ amm = sqrt(t);
68+ for m = 0:mmax
69+ if m > 0
70+ t = t * (2*m + 1) / (2*m);
71+ amm = (-1)^m * sqrt(t);
72+ end
73+ n = lmax - m + 1;
74+ Y = zeros(nlat, n);
75+ y0 = amm * S.st.^m;
76+ Y(:, 1) = y0;
77+ if n > 1
78+ y1 = sqrt(2*m + 3) * ct .* y0;
79+ Y(:, 2) = y1;
80+ for l = m + 2:lmax
81+ t1 = (l + m) * (l - m);
82+ a = sqrt((2*l + 1) * (2*l - 1) / t1);
83+ b = -sqrt(((2*l + 1) / (2*l - 3)) * ((l - 1 + m) * (l - 1 - m) / t1));
84+ y2 = a * ct .* y1 + b * y0;
85+ Y(:, l - m + 1) = y2;
86+ y0 = y1;
87+ y1 = y2;
88+ end
89+ end
90+ S.Y{m + 1} = Y;
91+ end
92+
93+ % sin(theta)*dtheta in coefficient space: v_l^m = ap(lm) u_{l-1}^m +
94+ % am(lm) u_{l+1}^m (src/sht/derivCoeffs.ts). Neighbors sit at +-1 within
95+ % each m block; ap/am are zero at the block edges, so the clamped index
96+ % vectors never read across a boundary.
97+ l = lv;
98+ m = mv;
99+ ap = (l - 1) .* sqrt(max(0, (l - m) .* (l + m)) ./ ((2*l - 1) .* (2*l + 1)));
100+ ap(l <= m) = 0;
101+ am = -(l + 2) .* sqrt((l + 1 - m) .* (l + 1 + m) ./ ((2*l + 1) .* (2*l + 3)));
102+ am(l >= lmax) = 0;
103+ S.ap = ap;
104+ S.am = am;
105+ S.iprev = max((1:S.nlm)' - 1, 1);
106+ S.inext = min((1:S.nlm)' + 1, S.nlm);
107+
108+ % dphig's Fourier multiplier: i*m on fft's frequency layout, masked past
109+ % the filter's reach (mcut = lmax-3), mirroring src/sht/wgsl/deriv.ts.
110+ freq = [(0:nphi/2)'; (-nphi/2 + 1:-1)'];
111+ S.dmul = 1i * freq .* (abs(freq) <= max(0, lmax - 3));
112+end
113+
114+% Gauss-Legendre nodes cos(theta), in decreasing order (north pole first),
115+% and weights for integration over cos(theta) -- Newton iteration on P_n,
116+% as src/sht/gauss.ts.
117+function [x, w] = gauss_legendre(n)
118+ x = zeros(n, 1);
119+ w = zeros(n, 1);
120+ half = floor((n + 1) / 2);
121+ for i = 1:half
122+ z = cos(pi * (i - 0.25) / (n + 0.5));
123+ pp = 0;
124+ for it = 1:100
125+ p1 = 1;
126+ p2 = 0;
127+ for j = 1:n
128+ p3 = p2;
129+ p2 = p1;
130+ p1 = ((2*j - 1) * z * p2 - (j - 1) * p3) / j;
131+ end
132+ pp = n * (z * p1 - p2) / (z^2 - 1);
133+ dz = p1 / pp;
134+ z = z - dz;
135+ if abs(dz) < 1e-15 * abs(z) + 1e-300
136+ p1 = 1;
137+ p2 = 0;
138+ for j = 1:n
139+ p3 = p2;
140+ p2 = p1;
141+ p1 = ((2*j - 1) * z * p2 - (j - 1) * p3) / j;
142+ end
143+ pp = n * (z * p1 - p2) / (z^2 - 1);
144+ z = z - p1 / pp;
145+ break;
146+ end
147+ end
148+ x(i) = z;
149+ x(n + 1 - i) = -z;
150+ wi = 2 / ((1 - z^2) * pp^2);
151+ w(i) = wi;
152+ w(n + 1 - i) = wi;
153+ end
154+ if mod(n, 2) == 1
155+ x(half) = 0;
156+ end
157+end
158+
159+% Synthesis, spectral -> grid. Grouped calls -- [a, b] = synth(x, y) -- are
160+% the app's batching hint; here each member simply runs in turn.
161+function varargout = synth(varargin)
162+ S = sht_tables();
163+ varargout = cell(1, nargin);
164+ for k = 1:nargin
165+ varargout{k} = synth_one(S, varargin{k});
166+ end
167+end
168+
169+function f = synth_one(S, Q)
170+ % Legendre stage per m, then one inverse FFT per latitude ring with the
171+ % m < 0 modes filled in by conjugate symmetry (the field is real).
172+ G = zeros(S.nphi, S.nlat);
173+ for m = 0:S.mmax
174+ Fm = (S.Y{m + 1} * Q(S.off(m + 1):S.off(m + 1) + S.lmax - m)).';
175+ G(m + 1, :) = Fm;
176+ if m > 0
177+ G(S.nphi + 1 - m, :) = conj(Fm);
178+ end
179+ end
180+ f = S.nphi * real(ifft(G, [], 1));
181+ f = f(:);
182+end
183+
184+% Analysis, grid -> spectral: forward FFT per ring, then Gauss quadrature
185+% against the same Legendre tables.
186+function varargout = analys(varargin)
187+ S = sht_tables();
188+ varargout = cell(1, nargin);
189+ for k = 1:nargin
190+ varargout{k} = analys_one(S, varargin{k});
191+ end
192+end
193+
194+function Q = analys_one(S, f)
195+ F = fft(reshape(f, S.nphi, S.nlat), [], 1) * (2*pi/S.nphi);
196+ Q = complex(zeros(S.nlm, 1));
197+ for m = 0:S.mmax
198+ Q(S.off(m + 1):S.off(m + 1) + S.lmax - m) = S.Y{m + 1}.' * (S.wg .* F(m + 1, :).');
199+ end
200+end
201+
202+% The coefficients of sin(theta)*dtheta(u): the alpha^+/alpha^- shift by one
203+% degree within each m block.
204+function V = dthetac(Q)
205+ S = sht_tables();
206+ V = S.ap .* Q(S.iprev) + S.am .* Q(S.inext);
207+end
208+
209+% The coefficients of dphi(u): i*m, diagonal.
210+function V = dphic(Q)
211+ S = sht_tables();
212+ V = 1i * (S.mv .* Q);
213+end
214+
215+% Grid-space derivatives, coefficients in: compositions of the shuffles and
216+% the synthesis (src/sht/deriv.ts). dtheta divides by sin(theta) afterwards.
217+function f = dtheta(Q)
218+ S = sht_tables();
219+ f = synth(dthetac(Q)) ./ S.stpt;
220+end
221+
222+function f = dphi(Q)
223+ f = synth(dphic(Q));
224+end
225+
226+% Grid-space phi derivative, grid in: two FFT stages and a pointwise i*m,
227+% no Legendre work -- d/dphi is diagonal in the Fourier index.
228+function g = dphig(f)
229+ S = sht_tables();
230+ F = fft(reshape(f, S.nphi, S.nlat), [], 1);
231+ g = real(ifft(S.dmul .* F, [], 1));
232+ g = g(:);
233+end
234+
235+% ---------------------------------------------------------------- the surface
236+%
237+% What the app precomputes from a shape's raw grid values: the band-limited
238+% embedding and both metric formulations built on it (src/geom/geometry.ts,
239+% src/geom/metric.ts). The solver runs on the synthesis of the coefficients,
240+% not on the raw values -- for a shape with sharp features the two differ.
241+function G = surface_tables(gxr, gyr, gzr)
242+ S = sht_tables();
243+ [G.Gx, G.Gy, G.Gz] = analys(gxr, gyr, gzr);
244+ [G.gx, G.gy, G.gz] = synth(G.Gx, G.Gy, G.Gz);
245+
246+ % Flux-form metric weights, from the sin-weighted theta tangent
247+ % sin(theta)*X_theta and X_phi, both smooth on the sphere:
248+ % gtt~ = sin^2 g_tt, gtp~ = sin g_tp, D = J sin^2(theta).
249+ [sXtx, sXty, sXtz] = synth(dthetac(G.Gx), dthetac(G.Gy), dthetac(G.Gz));
250+ [Xpx, Xpy, Xpz] = synth(dphic(G.Gx), dphic(G.Gy), dphic(G.Gz));
251+ gtt = sXtx.^2 + sXty.^2 + sXtz.^2;
252+ gtp = sXtx.*Xpx + sXty.*Xpy + sXtz.*Xpz;
253+ gpp = Xpx.^2 + Xpy.^2 + Xpz.^2;
254+ D = sqrt(gtt .* gpp - gtp.^2);
255+ G.p1 = gpp ./ D;
256+ G.p2 = -gtp ./ D;
257+ G.q2 = gtt ./ D;
258+ G.r = 1 ./ D;
259+
260+ % The sphere-subtracted weights and the bounded 1/J = r sin^2(theta) --
261+ % what keeps the concentrated division off the round sphere's share of the
262+ % flux divergence. Formed here in f64, as the app forms them.
263+ G.jinv = G.r .* S.stpt.^2;
264+ G.dp1 = G.p1 - 1;
265+ G.dq2 = G.q2 - 1;
266+
267+ % Preconditioner scale Jhat = 2/(muMin + muMax) over the eigenvalues of
268+ % the operator's symbol S = (1/J)[[p1, p2], [p2, q2]].
269+ s11 = G.p1 .* G.jinv;
270+ s12 = G.p2 .* G.jinv;
271+ s22 = G.q2 .* G.jinv;
272+ mn = (s11 + s22) / 2;
273+ disc = sqrt(((s11 - s22) / 2).^2 + s12.^2);
274+ G.Jhat = 2 / (min(mn - disc) + max(mn + disc));
275+
276+ % Inverse metric quantities V_theta/V_phi, for the Algorithm-4 models.
277+ Xtx = sXtx ./ S.stpt;
278+ Xty = sXty ./ S.stpt;
279+ Xtz = sXtz ./ S.stpt;
280+ g11 = Xtx.^2 + Xty.^2 + Xtz.^2;
281+ g12 = Xtx.*Xpx + Xty.*Xpy + Xtz.*Xpz;
282+ g22 = gpp;
283+ det = g11 .* g22 - g12.^2;
284+ G.Vtx = (g22 .* Xtx - g12 .* Xpx) ./ det;
285+ G.Vty = (g22 .* Xty - g12 .* Xpy) ./ det;
286+ G.Vtz = (g22 .* Xtz - g12 .* Xpz) ./ det;
287+ G.Vpx = (g11 .* Xpx - g12 .* Xtx) ./ det;
288+ G.Vpy = (g11 .* Xpy - g12 .* Xty) ./ det;
289+ G.Vpz = (g11 .* Xpz - g12 .* Xtz) ./ det;
290+end
291+
292+% ---------------------------------------------------------------- random field
293+%
294+% chebfun-style smooth random field in 3D, restricted to the surface by
295+% evaluating it at the grid points -- the way surfacefun seeds a run. Two
296+% signatures, as in the app:
297+% [k, c] = randnfun3(lambda, dom) the Fourier-mode draw (tools/randnfun3.m)
298+% f = randnfun3(lambda, gx, gy, gz) that draw, summed at the surface points
299+% Seed with rng(...) before calling.
300+function varargout = randnfun3(lambda, varargin)
301+ if nargin == 2
302+ [k, c] = randnfun3_modes(lambda, varargin{1});
303+ varargout = {k, c};
304+ return;
305+ end
306+ [gx, gy, gz] = deal(varargin{1:3});
307+ dom = [min(gx) max(gx) min(gy) max(gy) min(gz) max(gz)];
308+ [k, c] = randnfun3_modes(lambda, dom);
309+ % Summed in blocks of modes: the full npts x nmodes phase matrix can reach
310+ % hundreds of MB at a fine wavelength.
311+ f = zeros(numel(gx), 1);
312+ blk = 2048;
313+ for j0 = 1:blk:size(k, 1)
314+ j1 = min(j0 + blk - 1, size(k, 1));
315+ t = gx * k(j0:j1, 1)' + gy * k(j0:j1, 2)' + gz * k(j0:j1, 3)';
316+ f = f + cos(t) * c(j0:j1, 1) - sin(t) * c(j0:j1, 2);
317+ end
318+ varargout = {f};
319+end
320+
321+% ---------------------------------------------------------------- display
322+%
323+% The pattern on the surface, one panel per species. The solver grid has no
324+% pole rows and an open phi seam; wrap_grid closes both for display, capping
325+% each pole with the mean of its nearest ring.
326+function h = plot_setup(gx, gy, gz, fields, names)
327+ fig = figure('Name', 'turing-surface', 'Color', 'w');
328+ Xs = wrap_grid(gx);
329+ Ys = wrap_grid(gy);
330+ Zs = wrap_grid(gz);
331+ n = numel(fields);
332+ h.surf = gobjects(1, n);
333+ h.ax = gobjects(1, n);
334+ for k = 1:n
335+ ax = subplot(1, n, k, 'Parent', fig);
336+ h.surf(k) = surf(ax, Xs, Ys, Zs, wrap_grid(fields{k}), 'EdgeColor', 'none');
337+ shading(ax, 'interp');
338+ axis(ax, 'equal');
339+ axis(ax, 'off');
340+ colormap(ax, 'jet');
341+ colorbar(ax);
342+ h.ax(k) = ax;
343+ end
344+ h.names = names;
345+end
346+
347+function plot_update(h, fields, t, k, nsteps)
348+ for i = 1:numel(fields)
349+ C = wrap_grid(fields{i});
350+ set(h.surf(i), 'CData', C);
351+ lo = min(C(:));
352+ hi = max(C(:));
353+ if ~(hi > lo)
354+ hi = lo + 1;
355+ end
356+ caxis(h.ax(i), [lo hi]);
357+ title(h.ax(i), sprintf('%s t = %.3f (step %d/%d)', h.names{i}, t, k, nsteps));
358+ end
359+ drawnow;
360+end
361+
362+function M = wrap_grid(f)
363+ S = sht_tables();
364+ M = reshape(f, S.nphi, S.nlat).';
365+ M = [M, M(:, 1)];
366+ M = [mean(M(1, :)) * ones(1, S.nphi + 1); M; mean(M(end, :)) * ones(1, S.nphi + 1)];
367+end
368+
369+% ---------------------------------------------------------------- results file
370+%
371+% Complex coefficients -> flat float32 [re, im] per (l, m), the layout the
372+% app's reference-file reader expects (docs/ellipsoid-reference-spec.md).
373+function write_coeffs(fname, path, Q)
374+ flat = zeros(2 * numel(Q), 1);
375+ flat(1:2:end) = real(Q);
376+ flat(2:2:end) = imag(Q);
377+ h5create(fname, path, numel(flat), 'Datatype', 'single');
378+ h5write(fname, path, single(flat));
379+end
380+
381+% h5writeatt cannot create a bare group, so the attribute-only groups of the
382+% reference layout are made through the low-level API.
383+function make_group(fname, path)
384+ fid = H5F.open(fname, 'H5F_ACC_RDWR', 'H5P_DEFAULT');
385+ gid = H5G.create(fid, path, 'H5P_DEFAULT', 'H5P_DEFAULT', 'H5P_DEFAULT');
386+ H5G.close(gid);
387+ H5F.close(fid);
388+end
src/main.tsmodified+87−2View file
@@ -41,6 +41,7 @@ import {
4141 } from './compare/variants.ts';
4242 import { loadReferenceFile } from './compare/referenceFile.ts';
4343 import type { ReferenceCase } from './compare/referenceCase.ts';
44+import { generateMatlabScript, MATLAB_SCRIPT_NAME } from './export/matlabScript.ts';
4445
4546 const $ = <T extends HTMLElement>(id: string): T =>
4647 document.getElementById(id) as T;
@@ -86,6 +87,10 @@ const elStats = $('stats');
8687 const elBenchResult = $('benchresult');
8788 const elCmd = $('cmd');
8889 const elCopyCmd = $<HTMLButtonElement>('copycmd');
90+const elMatlab = $<HTMLDetailsElement>('matlab');
91+const elMatlabScript = $('matlabscript');
92+const elCopyMatlab = $<HTMLButtonElement>('copymatlab');
93+const elDownloadMatlab = $<HTMLButtonElement>('downloadmatlab');
8994 const elBlurb = $('blurb');
9095 const elErr = $('err');
9196 const elSource = $<HTMLTextAreaElement>('source');
@@ -349,6 +354,7 @@ function buildGeomParamInputs(): void {
349354 next = spec.min + Math.floor(Math.random() * (span + 1));
350355 }
351356 geomParams[spec.key] = next;
357+ updateCommand();
352358 viewChange = viewChange.then(() => applyGeometry());
353359 });
354360 elGeomParams.append(button);
@@ -365,6 +371,7 @@ function buildGeomParamInputs(): void {
365371 input.addEventListener('change', () => {
366372 const v = Number(input.value);
367373 if (Number.isFinite(v)) geomParams[spec.key] = v;
374+ updateCommand();
368375 viewChange = viewChange.then(() => applyGeometry());
369376 });
370377 label.append(input);
@@ -399,6 +406,8 @@ function applyGeometryChoice(key: string): void {
399406 editedGeomSource = null;
400407 buildGeomParamInputs();
401408 showEditorFile();
409+ // The command line and the MATLAB export both bake the surface in.
410+ updateCommand();
402411 }
403412
404413 /** Load the chosen file into the editor, keeping any unsaved edit to it. */
@@ -440,9 +449,43 @@ function updateCommand(): void {
440449 // equivalent is the ref checker, not the benchmark.
441450 if (compareRun?.refFile) {
442451 elCmd.textContent = `npm run ref -- --in ${compareRun.refFile.label}`;
443- return;
452+ } else {
453+ elCmd.textContent = formatCommand(currentSpec());
454+ }
455+ if (elMatlab.open) refreshMatlabScript();
456+}
457+
458+/** The run on screen as one standalone .m — in a study, its reference
459+ * variant, the same choice `currentSpec` makes. Throws on a working copy
460+ * the export cannot parse (no init/step/shape function). */
461+function matlabScriptText(): string {
462+ const spec = currentSpec();
463+ return generateMatlabScript({
464+ model,
465+ modelSource: source(),
466+ params: spec.params,
467+ geometry,
468+ geometrySource: geomSource(),
469+ geometryParams: geomParams,
470+ lmax: spec.lmax,
471+ niter: spec.niter,
472+ lam3: Number(elLam3.value),
473+ seed,
474+ preset: spec.preset,
475+ command: formatCommand(spec),
476+ });
477+}
478+
479+/** Regenerate the visible script text; on failure, show why in its place. */
480+function refreshMatlabScript(): string | null {
481+ try {
482+ const text = matlabScriptText();
483+ elMatlabScript.textContent = text;
484+ return text;
485+ } catch (e) {
486+ elMatlabScript.textContent = e instanceof Error ? e.message : String(e);
487+ return null;
444488 }
445- elCmd.textContent = formatCommand(currentSpec());
446489 }
447490
448491 elModel.addEventListener('change', () => {
@@ -475,6 +518,8 @@ elGeometry.addEventListener('change', () => {
475518 elLam3.addEventListener('change', () => {
476519 const v = Number(elLam3.value);
477520 if (!Number.isFinite(v) || v <= 0) return;
521+ // Not in the bench command, but the MATLAB export bakes lam3 in.
522+ updateCommand();
478523 // Changing the wavelength redraws the field, which restarts the run — so
479524 // pause first, exactly as the Re-seed button does. Without it the reseed's
480525 // readback races the pump's own, and the two collide on the staging buffer.
@@ -578,6 +623,46 @@ elCopyCmd.addEventListener('click', () => {
578623 navigator.clipboard.writeText(text).then(() => flash('Copied'), selectCommand);
579624 });
580625
626+// The MATLAB export: the same run as one self-contained .m. Regenerated from
627+// the current UI state whenever it is shown, copied or downloaded, so the
628+// text always matches the run on screen.
629+elMatlab.addEventListener('toggle', () => {
630+ if (elMatlab.open) refreshMatlabScript();
631+});
632+elCopyMatlab.addEventListener('click', (e) => {
633+ // The buttons live inside the <summary>; without this a click also toggles.
634+ e.preventDefault();
635+ e.stopPropagation();
636+ const text = refreshMatlabScript();
637+ if (text === null || !navigator.clipboard) {
638+ // Open to show the error, or the text to select by hand.
639+ elMatlab.open = true;
640+ return;
641+ }
642+ navigator.clipboard.writeText(text).then(
643+ () => {
644+ elCopyMatlab.textContent = 'Copied';
645+ setTimeout(() => (elCopyMatlab.textContent = 'Copy'), 1200);
646+ },
647+ () => (elMatlab.open = true),
648+ );
649+});
650+elDownloadMatlab.addEventListener('click', (e) => {
651+ e.preventDefault();
652+ e.stopPropagation();
653+ const text = refreshMatlabScript();
654+ if (text === null) {
655+ elMatlab.open = true;
656+ return;
657+ }
658+ const url = URL.createObjectURL(new Blob([text], { type: 'text/x-matlab' }));
659+ const a = document.createElement('a');
660+ a.href = url;
661+ a.download = `${MATLAB_SCRIPT_NAME}.m`;
662+ a.click();
663+ URL.revokeObjectURL(url);
664+});
665+
581666 // ---------------------------------------------------------------- setup
582667 function disposeView(): void {
583668 for (const s of scenes) s.dispose();
test/matlabExportChecks.tsadded+134−0View file
@@ -0,0 +1,134 @@
1+/**
2+ * The MATLAB export (src/export/matlabScript.ts) is string assembly, so these
3+ * checks are cheap and need no GPU: every preset x geometry combination must
4+ * generate, the assembled file must keep its local-function namespace free of
5+ * collisions, and the driver's calls must match the signatures the .m files
6+ * declare. Whether the generated MATLAB actually reproduces a run is checked
7+ * against MATLAB itself, not here: a run exported at defaults and executed in
8+ * MATLAB R2026b lands within fp32 accumulation error of the app's own replay
9+ * (relL2 ~1e-7 over 60 steps via `npm run ref`), and the flux and Algorithm-4
10+ * exports track each other to ~3e-10 in f64.
11+ */
12+import { generateMatlabScript, MATLAB_SCRIPT_NAME } from '../src/export/matlabScript.ts';
13+import { presets, mModelByKey } from '../src/mgpu/registry.ts';
14+import { mGeometries } from '../src/geom/registry.ts';
15+import { formatCommand, resolvePreset, DEFAULT_WARMUP } from '../src/bench/runSpec.ts';
16+
17+type Check = (name: string, ok: boolean, detail: string) => void;
18+type Log = (s: string) => void;
19+
20+export function matlabExportChecks(check: Check, log: Log): void {
21+ log('--- MATLAB export ---');
22+ for (const preset of presets) {
23+ const { model, params } = resolvePreset(preset.key);
24+ for (const geometry of mGeometries) {
25+ const geometryParams = Object.fromEntries(
26+ geometry.params.map((p) => [p.key, p.value]),
27+ );
28+ const spec = {
29+ preset: preset.key,
30+ lmax: 63,
31+ seed: 1,
32+ steps: 2000,
33+ warmup: DEFAULT_WARMUP,
34+ params,
35+ geometry: geometry.key,
36+ geometryParams,
37+ niter: 8,
38+ };
39+ const name = `matlab-export ${preset.key} on ${geometry.key}`;
40+ let text: string;
41+ try {
42+ text = generateMatlabScript({
43+ model,
44+ modelSource: model.source,
45+ params,
46+ geometry,
47+ geometrySource: geometry.source,
48+ geometryParams,
49+ lmax: 63,
50+ niter: 8,
51+ lam3: 0.5,
52+ seed: 1,
53+ preset: preset.key,
54+ command: formatCommand(spec),
55+ });
56+ } catch (e) {
57+ check(name, false, e instanceof Error ? e.message : String(e));
58+ continue;
59+ }
60+
61+ // One file, one namespace: every local function name must be unique,
62+ // or MATLAB silently shadows one definition with another.
63+ const fnNames = [...text.matchAll(/^[ \t]*function\s+(?:\[[^\]]*\]|\w+)\s*=\s*(\w+)\s*\(/gm)]
64+ .map((m) => m[1]);
65+ const dupes = fnNames.filter((n, i) => fnNames.indexOf(n) !== i);
66+
67+ // The driver must define what it calls: the state it steps, the model
68+ // call mapped through the mp struct, and the transform setup.
69+ const wants = [
70+ `function ${MATLAB_SCRIPT_NAME}()`,
71+ 'sht_tables(sht_setup(lmax, mmax, nlat, nphi));',
72+ `= init(`,
73+ `= step(${model.state.join(', ')}, `,
74+ 'surface_tables(gxr, gyr, gzr)',
75+ `'/final/${model.state[0]}'`,
76+ ];
77+ const missing = wants.filter((w) => !text.includes(w));
78+
79+ // randnfunsphere rides along exactly when the geometry draws on it.
80+ const wantsSphereTool = /\brandnfunsphere\b/.test(geometry.source);
81+ const carriesSphereTool = /function f = randnfunsphere\(/.test(text);
82+
83+ const problems = [
84+ ...(dupes.length ? [`duplicate local functions: ${[...new Set(dupes)].join(', ')}`] : []),
85+ ...(missing.length ? [`missing: ${missing.join(' | ')}`] : []),
86+ ...(wantsSphereTool !== carriesSphereTool
87+ ? [`randnfunsphere ${wantsSphereTool ? 'missing' : 'included needlessly'}`]
88+ : []),
89+ ];
90+ check(name, problems.length === 0, problems.join('; ') || `${fnNames.length} local functions`);
91+ }
92+ }
93+
94+ // An edited working copy that dropped a required function is refused with a
95+ // message naming the file, not exported broken.
96+ const { model, params } = resolvePreset(presets[0].key);
97+ const geometry = mGeometries[0];
98+ try {
99+ generateMatlabScript({
100+ model,
101+ modelSource: '% nothing here',
102+ params,
103+ geometry,
104+ geometrySource: geometry.source,
105+ geometryParams: {},
106+ lmax: 63,
107+ niter: 8,
108+ lam3: 0.5,
109+ seed: 1,
110+ preset: presets[0].key,
111+ command: '',
112+ });
113+ check('matlab-export refuses a source without init', false, 'no error thrown');
114+ } catch (e) {
115+ const msg = e instanceof Error ? e.message : String(e);
116+ check(
117+ 'matlab-export refuses a source without init',
118+ msg.includes("'init'") && msg.includes(model.key),
119+ msg,
120+ );
121+ }
122+ // Guard the assumption the model registry makes for stateFor: state names
123+ // are used as `<name>0` initial-capture variables, which must not collide
124+ // with the species names.
125+ for (const p of presets) {
126+ const m = mModelByKey(p.modelKey)!;
127+ const all = new Set([...m.state, ...m.species]);
128+ check(
129+ `matlab-export names disjoint for ${m.key}`,
130+ all.size === m.state.length + m.species.length,
131+ [...all].join(', '),
132+ );
133+ }
134+}
test/test-page.tsmodified+2−0View file
@@ -31,6 +31,7 @@ import { geometryChecks } from './geometryChecks.ts';
3131 import { fluxChecks } from './fluxChecks.ts';
3232 import { compareChecks } from './compareChecks.ts';
3333 import { referenceChecks, type H5Rt } from './referenceChecks.ts';
34+import { matlabExportChecks } from './matlabExportChecks.ts';
3435
3536 declare global {
3637 interface Window {
@@ -226,6 +227,7 @@ async function main(): Promise<void> {
226227 await compareChecks(device, check, log);
227228 // '/' is the wasm module's in-memory filesystem — nothing touches disk.
228229 await referenceChecks(h5wasm as unknown as H5Rt, (name) => `/${name}`, check, log);
230+ matlabExportChecks(check, log);
229231
230232 window.__RESULTS__ = { ok: failures === 0, lines };
231233 log(failures === 0 ? 'ALL PASS' : `${failures} FAILURE(S)`);