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 */
19import supportSource from './support.m?raw';
20import randnfun3Source from '../../tools/randnfun3.m?raw';
21import randnfunsphereSource from '../../tools/randnfunsphere.m?raw';
22import type { MModel, Params } from '../mgpu/registry.ts';
23import type { MGeometry } from '../geom/registry.ts';
25/** The generated function's name — and therefore the file name to save as. */
26export const MATLAB_SCRIPT_NAME = 'turing_surface_run';
28export 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}
49interface Signature {
50 outputs: string[];
51 params: string[];
52}
54/** First `function [outs] = name(args)` line in a .m — the same contract the
55 * compiler applies, minus everything it checks later. */
56function 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}
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. */
75const num = (v: number): string => (Number.isFinite(v) ? String(v) : '0');
77/** A string as a MATLAB char literal. */
78const str = (s: string): string => `'${s.replace(/'/g, "''")}'`;
80const banner = (title: string): string => {
81 const line = `% ${'='.repeat(72)}`;
82 return `${line}\n% ${title}\n${line}`;
83};
85export 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 };
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`);
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}`;
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(', ')});`;
113 const speciesCell = `{${model.species.join(', ')}}`;
114 const namesCell = `{${model.species.map((s) => str(s)).join(', ')}}`;
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');
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');
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.
152% ---- run controls --------------------------------------------------------
153nsteps = ${controls.nsteps}; % timesteps to run
154plot_every = ${controls.plotEvery}; % live-plot interval, in steps; 0 disables plotting
155out_file = ${str(controls.outFile)}; % results file; '' disables
156seed = ${num(spec.seed)}; % rng seed for the initial condition
158% ---- captured from the app -----------------------------------------------
159lmax = ${spec.lmax}; % spherical-harmonic truncation degree
160niter = ${spec.niter}; % iterations of the implicit solve's geometric correction
161lam3 = ${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`}
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.
168pdeg = ${model.pdeg};
169mmax = lmax;
170nlat = 2 * ceil(max(lmax + 1, ((pdeg + 1) * lmax + 1) / 2) / 2);
171nphi = 2 ^ nextpow2((pdeg + 1) * lmax + 1);
172npts = nlat * nphi;
173sht_tables(sht_setup(lmax, mmax, nlat, nphi));
174S = sht_tables();
175nlm = S.nlm;
176lam = S.lam;
177filt = S.filt;
178theta = S.theta;
179phi = S.phi;
181% ---- the surface ---------------------------------------------------------
182${shapeCall}
183% A constant coordinate comes back scalar; spread it over the grid.
184gxr = gxr + zeros(npts, 1);
185gyr = gyr + zeros(npts, 1);
186gzr = gzr + zeros(npts, 1);
187G = surface_tables(gxr, gyr, gzr);
188gx = G.gx; gy = G.gy; gz = G.gz;
189Gx = G.Gx; Gy = G.Gy; Gz = G.Gz;
190p1 = G.p1; p2 = G.p2; q2 = G.q2; r = G.r;
191dp1 = G.dp1; dq2 = G.dq2; jinv = G.jinv;
192Vtx = G.Vtx; Vty = G.Vty; Vtz = G.Vtz;
193Vpx = G.Vpx; Vpy = G.Vpy; Vpz = G.Vpz;
194jhat = G.Jhat;
195radius = sqrt(gx.^2 + gy.^2 + gz.^2);
196fprintf('grid %d x %d, nlm %d, radius %.3f-%.3f, Jhat %.3f\\n', ...
197 nlat, nphi, nlm, min(radius), max(radius), jhat);
199% ---- initial condition ---------------------------------------------------
200rng(seed);
201${takesNoise ? `noise = ${num(model.seedAmp)} * randn(npts, 1);\n` : ''}${initCall}
202${model.state.map((s) => `${s}0 = ${s};`).join('\n')}
204% ---- time loop -----------------------------------------------------------
205if plot_every > 0
206 ph = plot_setup(gx, gy, gz, ${speciesCell}, ${namesCell});
207 plot_update(ph, ${speciesCell}, 0, 0, nsteps);
208end
209report_every = max(1, round(nsteps / 10));
210t = 0;
211tstart = tic;
212for 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
221end
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).
227if ~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);
282end
283end`;
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 }
296 const usesSphere = /\brandnfunsphere\b/.test(spec.geometrySource + spec.modelSource);
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}