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