1// Convert the `sys` struct produced by a numbl script into a BeadSystem.
2// The script uses 1-based bead indices; they become 0-based here.
4import type { RuntimeValue } from "numbl";
5import type { BeadSystem, Hit, Method } from "./types";
7type Mat = { rows: number; cols: number; data: Float64Array };
9class SysError extends Error {}
11function asMat(v: RuntimeValue | undefined, name: string): Mat | null {
12 if (v === undefined) return null;
13 if (typeof v === "number") return { rows: 1, cols: 1, data: Float64Array.of(v) };
14 if (typeof v === "boolean") return { rows: 1, cols: 1, data: Float64Array.of(v ? 1 : 0) };
15 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor") {
16 const t = v as unknown as { data: Float64Array; shape: number[]; imag?: Float64Array };
17 if (t.imag) throw new SysError(`${name} must be real`);
18 const rows = t.shape[0] ?? 1;
19 const cols = t.shape.slice(1).reduce((a, b) => a * b, 1);
20 return { rows, cols, data: t.data };
21 }
22 throw new SysError(`${name} must be numeric`);
23}
25function scalar(v: RuntimeValue | undefined, name: string, def: number): number {
26 const m = asMat(v, name);
27 if (!m || m.data.length === 0) return def;
28 if (m.data.length !== 1) throw new SysError(`${name} must be a scalar`);
29 return m.data[0];
30}
32function vec3(v: RuntimeValue | undefined, name: string, def: [number, number, number]): [number, number, number] {
33 const m = asMat(v, name);
34 if (!m || m.data.length === 0) return def;
35 if (m.data.length !== 3) throw new SysError(`${name} must have 3 elements`);
36 return [m.data[0], m.data[1], m.data[2]];
37}
39function str(v: RuntimeValue | undefined, name: string, def: string): string {
40 if (v === undefined) return def;
41 if (typeof v === "string") return v;
42 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "char") {
43 return (v as unknown as { value: string }).value;
44 }
45 throw new SysError(`${name} must be a string`);
46}
48function perBead(v: RuntimeValue | undefined, name: string, n: number, def: number): Float64Array {
49 const out = new Float64Array(n);
50 const m = asMat(v, name);
51 if (!m || m.data.length === 0) return out.fill(def);
52 if (m.data.length === 1) return out.fill(m.data[0]);
53 if (m.data.length !== n) throw new SysError(`${name} must be a scalar or have one entry per bead (${n})`);
54 out.set(m.data);
55 return out;
56}
58function indices(v: RuntimeValue | undefined, name: string, n: number): Int32Array {
59 const m = asMat(v, name);
60 if (!m) return new Int32Array(0);
61 const isLogical = (v as { _isLogical?: boolean })._isLogical || typeof v === "boolean";
62 if (isLogical) {
63 if (m.data.length !== n) throw new SysError(`${name}: a logical mask must have one entry per bead (${n})`);
64 const out: number[] = [];
65 m.data.forEach((b, i) => b && out.push(i));
66 return Int32Array.from(out);
67 }
68 const out = new Int32Array(m.data.length);
69 m.data.forEach((b, q) => {
70 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`${name}: bead index ${b} is not in 1..${n}`);
71 out[q] = b - 1;
72 });
73 return out;
74}
76function fields(v: RuntimeValue | undefined, name: string): Map<string, RuntimeValue>[] {
77 if (v === undefined) return [];
78 const k = (v as { kind?: string }).kind;
79 if (k === "struct") return [(v as unknown as { fields: Map<string, RuntimeValue> }).fields];
80 if (k === "struct_array") {
81 return (v as unknown as { elements: { fields: Map<string, RuntimeValue> }[] }).elements.map(e => e.fields);
82 }
83 throw new SysError(`${name} must be a struct or struct array`);
84}
86const KNOWN = new Set([
87 "pos", "mass", "pinned", "springs", "dashpots", "friction", "hits", "listen",
88 "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes",
89]);
91export function parseSystem(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
92 try {
93 return parseInner(sysValue);
94 } catch (e) {
95 if (e instanceof SysError) throw new Error(e.message);
96 throw e;
97 }
98}
100function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
101 if (sysValue === undefined) throw new SysError("The script must define a struct named sys.");
102 const [f] = fields(sysValue, "sys");
103 if (!f) throw new SysError("sys must be a single struct");
104 const warnings: string[] = [];
105 for (const key of f.keys()) if (!KNOWN.has(key)) warnings.push(`sys.${key} is not used by the solver.`);
107 const P = asMat(f.get("pos"), "sys.pos");
108 if (!P || P.cols !== 3 || P.rows < 1) throw new SysError("sys.pos must be an n-by-3 matrix of bead positions");
109 const n = P.rows;
110 const pos0 = new Float64Array(3 * n);
111 for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) pos0[3 * i + c] = P.data[i + c * n];
113 if (!f.has("mass")) throw new SysError("sys.mass is required (scalar or one per bead)");
114 const mass = perBead(f.get("mass"), "sys.mass", n, 0);
115 for (let i = 0; i < n; i++) if (!(mass[i] > 0)) throw new SysError(`sys.mass must be positive (bead ${i + 1})`);
117 const pinned = new Uint8Array(n);
118 for (const b of indices(f.get("pinned"), "sys.pinned", n)) pinned[b] = 1;
120 const S = asMat(f.get("springs"), "sys.springs");
121 const ns = S && S.data.length ? S.rows : 0;
122 if (S && ns && S.cols !== 4 && S.cols !== 5) {
123 throw new SysError("sys.springs must have 4 or 5 columns: [i j k L0] or [i j k L0 k3]");
124 }
125 const si = new Int32Array(ns);
126 const sj = new Int32Array(ns);
127 const sk = new Float64Array(ns);
128 const sL0 = new Float64Array(ns);
129 const sk3 = new Float64Array(ns);
130 for (let s = 0; s < ns; s++) {
131 const col = (c: number) => S!.data[s + c * ns];
132 const i = col(0);
133 const j = col(1);
134 for (const b of [i, j]) {
135 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.springs row ${s + 1}: bead index ${b} is not in 1..${n}`);
136 }
137 if (i === j) throw new SysError(`sys.springs row ${s + 1} joins bead ${i} to itself`);
138 si[s] = i - 1;
139 sj[s] = j - 1;
140 sk[s] = col(2);
141 sL0[s] = col(3);
142 sk3[s] = S!.cols === 5 ? col(4) : 0;
143 if (sk[s] < 0 || sk3[s] < 0 || sL0[s] < 0) throw new SysError(`sys.springs row ${s + 1}: k, L0 and k3 must be nonnegative`);
144 }
146 const D = asMat(f.get("dashpots"), "sys.dashpots");
147 const nd = D && D.data.length ? D.rows : 0;
148 if (D && nd && D.cols !== 3) throw new SysError("sys.dashpots must have 3 columns: [i j c]");
149 const di = new Int32Array(nd);
150 const dj = new Int32Array(nd);
151 const dc = new Float64Array(nd);
152 for (let d = 0; d < nd; d++) {
153 const col = (c: number) => D!.data[d + c * nd];
154 for (const b of [col(0), col(1)]) {
155 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.dashpots row ${d + 1}: bead index ${b} is not in 1..${n}`);
156 }
157 di[d] = col(0) - 1;
158 dj[d] = col(1) - 1;
159 dc[d] = col(2);
160 if (dc[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: c must be nonnegative`);
161 }
162 const friction = perBead(f.get("friction"), "sys.friction", n, 0);
164 const hits: Hit[] = fields(f.get("hits"), "sys.hits").map((h, q) => {
165 const nm = `sys.hits(${q + 1})`;
166 const beads = indices(h.get("beads"), `${nm}.beads`, n);
167 if (beads.length === 0) throw new SysError(`${nm}.beads is required`);
168 const duration = scalar(h.get("duration"), `${nm}.duration`, 0);
169 const hit: Hit = {
170 t: scalar(h.get("t"), `${nm}.t`, 0),
171 beads,
172 velocity: vec3(h.get("velocity"), `${nm}.velocity`, [0, 0, 0]),
173 force: vec3(h.get("force"), `${nm}.force`, [0, 0, 0]),
174 duration,
175 };
176 // In a struct array every hit has every field, so treat empty [] as absent.
177 const given = (k: string) => (asMat(h.get(k), `${nm}.${k}`)?.data.length ?? 0) > 0;
178 if (duration > 0 && given("velocity")) throw new SysError(`${nm}: give velocity (impulse) or force with duration (pulse), not both`);
179 if (duration === 0 && given("force")) throw new SysError(`${nm}.force needs a positive duration`);
180 return hit;
181 });
182 if (hits.length === 0) warnings.push("sys.hits is empty: nothing excites the system.");
184 const [L] = fields(f.get("listen"), "sys.listen");
185 if (!L) throw new SysError("sys.listen is required, e.g. struct('beads', 10, 'dir', [0 0 1])");
186 const listenBeads = indices(L.get("beads") ?? L.get("bead"), "sys.listen.beads", n);
187 if (listenBeads.length === 0) throw new SysError("sys.listen.beads is required");
188 const dir = vec3(L.get("dir"), "sys.listen.dir", [0, 0, 1]);
189 const dn = Math.hypot(...dir);
190 if (dn === 0) throw new SysError("sys.listen.dir must be nonzero");
192 const method = str(f.get("method"), "sys.method", "verlet") as Method;
193 if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
195 const sys: BeadSystem = {
196 n, pos0, mass, pinned, si, sj, sk, sL0, sk3, di, dj, dc, friction, hits, listenBeads,
197 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
198 duration: scalar(f.get("duration"), "sys.duration", 2),
199 sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
200 dt: scalar(f.get("dt"), "sys.dt", 0),
201 safety: scalar(f.get("safety"), "sys.safety", 0.5),
202 method,
203 frameRate: scalar(f.get("frameRate"), "sys.frameRate", 500),
204 modes: scalar(f.get("modes"), "sys.modes", 1) !== 0,
205 };
206 if (!(sys.duration > 0)) throw new SysError("sys.duration must be positive");
207 if (!(sys.safety > 0 && sys.safety <= 1)) throw new SysError("sys.safety must be in (0, 1]");
208 if (ns === 0) warnings.push("sys.springs is empty: the beads do not interact.");
209 return { sys, warnings };
210}