2// Runs the MATLAB layer (driver.m + a method script + the helpers) through the
3// numbl CLI outside the browser, and checks the numbers against the paper.
4//
5// node scripts/matlab-test.mjs # everything except the slow cases
6// node scripts/matlab-test.mjs --full # also n = 640 in the convergence study
7//
8// NUMBL_DIR points at a clone of https://github.com/flatironinstitute/numbl
9// (default ~/src/numbl); the CLI is run with npx tsx, no global install.
10import { execFileSync } from 'node:child_process'
11import { mkdtempSync, readFileSync, writeFileSync, rmSync } from 'node:fs'
12import { tmpdir } from 'node:os'
13import { join, dirname } from 'node:path'
14import { fileURLToPath } from 'node:url'
16const here = dirname(fileURLToPath(import.meta.url))
17const root = join(here, '..')
18const NUMBL = process.env.NUMBL_DIR ?? join(process.env.HOME, 'src', 'numbl')
19const FULL = process.argv.includes('--full')
21const HELPERS = ['nodes_of.m', 'testfun.m', 'evalexpr.m', 'cubic_spline.m', 'classical_rational.m']
22const read = (p) => readFileSync(join(root, p), 'utf8')
24function run(methodFile, params) {
25 const dir = mkdtempSync(join(tmpdir(), 'bary-'))
26 try {
27 for (const h of HELPERS) writeFileSync(join(dir, h), read(`src/matlab/lib/${h}`))
28 writeFileSync(join(dir, 'main.m'), read('src/matlab/driver.m') + '\n' + read(`src/methods/${methodFile}`))
29 writeFileSync(join(dir, 'params.json'), JSON.stringify(params))
30 execFileSync('npx', ['tsx', join(NUMBL, 'src', 'cli.ts'), 'run', 'main.m'], {
31 cwd: dir,
32 stdio: ['ignore', 'pipe', 'pipe'],
33 encoding: 'utf8',
34 })
35 return JSON.parse(readFileSync(join(dir, 'out.json'), 'utf8'))
36 } finally {
37 rmSync(dir, { recursive: true, force: true })
38 }
39}
41const explore = (over = {}) => ({
42 mode: 'explore',
43 f: 'runge',
44 fexpr: '',
45 a: -5,
46 b: 5,
47 n: 20,
48 d: 3,
49 nodes: 'uniform',
50 seed: 1,
51 ngrid: 801,
52 ngridwide: 801,
53 rootsMaxN: 40,
54 want: { poly: false, spline: false, blend: false, poles: false, classical: false },
55 ...over,
56})
58let failures = 0
59function check(name, ok, detail = '') {
60 console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
61 if (!ok) failures++
62}
63const close = (a, b, tol) => Math.abs(a - b) <= tol
65// ── the integer weight patterns of Section 4 ───────────────────────────────
66console.log('\nSection 4: integer weights on a uniform mesh')
67const PATTERNS = {
68 0: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
69 1: [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1],
70 2: [1, 3, 4, 4, 4, 4, 4, 4, 4, 3, 1],
71 3: [1, 4, 7, 8, 8, 8, 8, 8, 7, 4, 1],
72 4: [1, 5, 11, 15, 16, 16, 16, 15, 11, 5, 1],
73}
74for (const [d, want] of Object.entries(PATTERNS)) {
75 for (const method of ['fh.m', 'uniform-integer.m']) {
76 const o = run(method, explore({ n: 10, d: Number(d) }))
77 const got = o.wscaled.map((v) => Math.round(v * 1e6) / 1e6)
78 check(
79 `d = ${d} ${method.padEnd(18)} delta_k`,
80 o.wIsInteger && JSON.stringify(got) === JSON.stringify(want) && o.wAlternates,
81 `[${got.join(' ')}]`,
82 )
83 }
84}
86// ── Table 1: Runge with d = 3, sine with d = 4, |x| with d = 3 ────────────
87console.log('\nTable 1: error in the rational interpolant')
88const ns = FULL ? [10, 20, 40, 80, 160, 320, 640] : [10, 20, 40, 80, 160, 320]
89const TABLE1 = {
90 runge: { d: 3, err: [6.9e-2, 2.8e-3, 4.3e-6, 5.1e-8, 3.0e-9, 1.8e-10, 1.1e-11] },
91 // The n = 20 entry is printed as 3.9e-05 in Table 1, but the paper's own
92 // order column next to it (5.5) says 1.7e-2 / 2^5.5 = 3.8e-04, and every
93 // other entry in the row matches to two figures. We take it as a misprint.
94 sine: { d: 4, err: [1.7e-2, 3.9e-4, 7.1e-6, 1.3e-7, 2.7e-9, 6.0e-11, 1.5e-12] },
95 abs: { d: 3, err: [1.9e-1, 9.5e-2, 4.8e-2, 2.4e-2, 1.2e-2, 5.9e-3, 3.0e-3] },
96}
97for (const [f, spec] of Object.entries(TABLE1)) {
98 const o = run('fh.m', {
99 mode: 'converge',
100 f,
101 fexpr: '',
102 a: -5,
103 b: 5,
104 nodes: 'uniform',
105 seed: 1,
106 ngrid: 4001,
107 ns,
108 ds: [spec.d],
109 want: { poly: false, spline: true },
110 })
111 const got = o.E[0]
112 const ok = got.every((e, i) => close(Math.log10(e), Math.log10(spec.err[i]), 0.05))
113 check(
114 `${f.padEnd(6)} d = ${spec.d}`,
115 ok,
116 got.map((e) => e.toExponential(1)).join(' '),
117 )
118 const ord = o.orders[0].slice(1)
119 console.log(` orders ${ord.map((v) => v.toFixed(1)).join(' ')}`)
120}
122// ── Table 3: rational (d = 3) against the clamped cubic spline, Runge ─────
123console.log('\nTable 3: rational d = 3 vs clamped cubic spline (Runge)')
124{
125 const o = run('fh.m', {
126 mode: 'converge',
127 f: 'runge',
128 fexpr: '',
129 a: -5,
130 b: 5,
131 nodes: 'uniform',
132 seed: 1,
133 ngrid: 4001,
134 ns,
135 ds: [3],
136 want: { poly: false, spline: true },
137 })
138 const WANT = [2.2e-2, 3.2e-3, 2.8e-4, 1.6e-5, 9.5e-7, 5.9e-8, 3.7e-9]
139 const ok = o.splineErr.every((e, i) => close(Math.log10(e), Math.log10(WANT[i]), 0.06))
140 check('spline error', ok, o.splineErr.map((e) => e.toExponential(1)).join(' '))
141 const last = o.splineErr.length - 1
142 check(
143 'rational beats spline by >100x at the largest n',
144 o.splineErr[last] / o.E[0][last] > 100,
145 `ratio ${(o.splineErr[last] / o.E[0][last]).toFixed(0)}x`,
146 )
147}
149// ── Table 4: the sine function, where the spline wins ─────────────────────
150console.log('\nTable 4: the sine function, where the spline is the better one')
151{
152 const o = run('fh.m', {
153 mode: 'converge',
154 f: 'sine',
155 fexpr: '',
156 a: -5,
157 b: 5,
158 nodes: 'uniform',
159 seed: 1,
160 ngrid: 4001,
161 ns,
162 ds: [3],
163 want: { poly: false, spline: true },
164 })
165 const RAT = [1.3e-2, 1.2e-3, 8.4e-5, 5.4e-6, 3.4e-7, 2.1e-8, 1.3e-9]
166 const SPL = [3.3e-3, 1.7e-4, 1.0e-5, 6.4e-7, 4.0e-8, 2.5e-9, 1.6e-10]
167 check('rational d = 3', o.E[0].every((e, i) => close(Math.log10(e), Math.log10(RAT[i]), 0.06)),
168 o.E[0].map((e) => e.toExponential(1)).join(' '))
169 check('spline', o.splineErr.every((e, i) => close(Math.log10(e), Math.log10(SPL[i]), 0.06)),
170 o.splineErr.map((e) => e.toExponential(1)).join(' '))
171}
173// ── Theorem 1: no real poles, for any d and any node distribution ─────────
174console.log('\nTheorem 1: no poles in R')
175for (const nodes of ['uniform', 'chebyshev', 'random', 'paired', 'graded']) {
176 for (const d of [0, 1, 3, 6]) {
177 const n = 16
178 const o = run('fh.m', explore({ n, d, nodes, want: { ...explore().want, poles: true } }))
179 const p = o.poles
180 const minIm = Math.min(...p.rootsIm.map(Math.abs))
181 // The denominator s of equation (10) has degree at most n - d. The
182 // leading coefficient of mu_i is (-1)^(n-i-d), so the leading coefficient
183 // of s is +-sum_{i=0}^{n-d} (-1)^i, which is 1 when n - d is even and 0
184 // when it is odd: the same parity that splits Theorem 2 into two cases.
185 // Theorem 1 then puts every one of those roots off the real axis.
186 const deg = (n - d) % 2 === 0 ? n - d : n - d - 1
187 check(
188 `${nodes.padEnd(10)} d = ${d}`,
189 p.realPoles.length === 0 && p.rootsShown && p.rootsRe.length === deg && minIm > 1e-8,
190 `${p.rootsRe.length} roots (want ${deg}), min |Im| = ${minIm.toExponential(1)}`,
191 )
192 }
193}
195// ── the counter-examples: weights that do not alternate ───────────────────
196console.log('\nWeights that do not alternate in sign do have poles')
197for (const method of ['equal.m', 'random.m']) {
198 const o = run(method, explore({ n: 12, d: 3, want: { ...explore().want, poles: true } }))
199 check(`${method.padEnd(10)} real poles found`, o.poles.realPoles.length > 0,
200 `${o.poles.realPoles.length} poles`)
201}
202{
203 // the Lagrange weights are the degenerate case: the denominator is constant
204 const o = run('lagrange.m', explore({ n: 12, d: 3, want: { ...explore().want, poles: true } }))
205 check('lagrange.m no roots at all (denominator is 1)',
206 o.poles.realPoles.length === 0 && o.poles.rootsRe.length === 0)
207}
209// ── the blend of equations (4) and (5) reproduces r ───────────────────────
210console.log('\nEquations (4) and (5): the blend equals the barycentric form')
211for (const d of [0, 1, 3, 5]) {
212 const o = run('fh.m', explore({ n: 14, d, want: { ...explore().want, blend: true } }))
213 check(`d = ${d} hasBlend`, o.hasBlend === true, o.blendError ?? '')
214 if (!o.hasBlend) continue
215 let maxDiff = 0
216 let maxPU = 0
217 for (let j = 0; j < o.t.length; j++) {
218 let s = 0
219 let pu = 0
220 for (let i = 0; i < o.L.length; i++) {
221 s += o.L[i][j] * o.P[i][j]
222 pu += o.L[i][j]
223 }
224 maxDiff = Math.max(maxDiff, Math.abs(s - o.r[j]))
225 maxPU = Math.max(maxPU, Math.abs(pu - 1))
226 }
227 check(`d = ${d} sum_i L_i p_i == r`, maxDiff < 1e-9, `max diff ${maxDiff.toExponential(1)}`)
228 check(`d = ${d} sum_i L_i == 1`, maxPU < 1e-12, `max dev ${maxPU.toExponential(1)}`)
229}
231// ── d = n is the polynomial interpolant ───────────────────────────────────
232console.log('\nd = n is the polynomial interpolant of equation (2)')
233{
234 const o = run('fh.m', explore({ n: 12, d: 12, want: { ...explore().want, poly: true } }))
235 const m = Math.max(...o.r.map((v, i) => Math.abs(v - o.rpoly[i])))
236 check('r (d = n) == barycentric Lagrange', m < 1e-10, `max diff ${m.toExponential(1)}`)
237}
239// ── Berrut on a badly graded mesh: the point of Theorem 3's beta ──────────
240console.log("\nTheorem 3: d = 0 needs a bounded mesh ratio, d >= 1 does not")
241{
242 const o = run('fh.m', {
243 mode: 'converge',
244 f: 'runge',
245 fexpr: '',
246 a: -5,
247 b: 5,
248 nodes: 'paired',
249 seed: 1,
250 ngrid: 2001,
251 ns: [20, 40, 80, 160],
252 ds: [0, 1, 3],
253 want: { poly: false, spline: false },
254 })
255 const ord = (row) => o.orders[row].slice(1)
256 console.log(` d = 0 orders ${ord(0).map((v) => v.toFixed(1)).join(' ')}`)
257 console.log(` d = 1 orders ${ord(1).map((v) => v.toFixed(1)).join(' ')}`)
258 console.log(` d = 3 orders ${ord(2).map((v) => v.toFixed(1)).join(' ')}`)
259 const last = o.ns.length - 1
260 check('d = 0 does much worse than d = 1 on the paired mesh',
261 o.E[0][last] / o.E[1][last] > 10, `ratio ${(o.E[0][last] / o.E[1][last]).toExponential(1)}`)
262 check('d = 3 still converges near h^4 on the paired mesh',
263 ord(2).slice(-1)[0] > 3.0, `last order ${ord(2).slice(-1)[0].toFixed(1)}`)
264}
266// ── the classical alternative does put poles in the interval ──────────────
267console.log('\nThe classical rational interpolant p_M/q_N')
268{
269 const o = run('fh.m', explore({ n: 12, d: 3, want: { ...explore().want, poles: true, classical: true } }))
270 const inside = o.poles.classicalPoles.filter((p) => p >= -5 && p <= 5)
271 check('has real poles', o.poles.classicalPoles.length > 0,
272 `${o.poles.classicalPoles.length} real poles, ${inside.length} inside [-5, 5]`)
273}
275console.log(`\n${failures === 0 ? 'all checks passed' : `${failures} FAILURES`}\n`)
276process.exit(failures === 0 ? 0 : 1)