Interactive illustration of Floater-Hormann barycentric rational interpolation
The method is an editable MATLAB-syntax script that runs in the browser via
numbl; the four tabs are all driven by it. A script supplies bary_weights and
bary_eval (and optionally local_blend), and everything else the tabs show is
derived in the driver from x, y and w, so the Poles tab tells the truth about
whatever weights are supplied.
Checked against the paper: the integer weight patterns of Section 4, the error
columns of Tables 1, 3 and 4, the absence of real roots for every d across five
node distributions, and r = sum_i L_i p_i relating equations (1) and (4).
43 changed files+9612−0
.github/workflows/deploy.ymladded+44−0View file
@@ -0,0 +1,44 @@
1+name: Deploy to GitHub Pages
2+
3+on:
4+ push:
5+ branches: [main]
6+ workflow_dispatch:
7+
8+permissions:
9+ contents: read
10+ pages: write
11+ id-token: write
12+
13+concurrency:
14+ group: pages
15+ cancel-in-progress: false
16+
17+jobs:
18+ build:
19+ runs-on: ubuntu-latest
20+ steps:
21+ - uses: actions/checkout@v4
22+ - uses: actions/setup-node@v4
23+ with:
24+ node-version: 22
25+ cache: npm
26+ # puppeteer is only used by scripts/browser-test.mjs, which does not run
27+ # here; skip the ~150 MB Chrome download
28+ - run: npm ci
29+ env:
30+ PUPPETEER_SKIP_DOWNLOAD: 'true'
31+ - run: npm run build
32+ - uses: actions/upload-pages-artifact@v3
33+ with:
34+ path: dist
35+
36+ deploy:
37+ needs: build
38+ runs-on: ubuntu-latest
39+ environment:
40+ name: github-pages
41+ url: ${{ steps.deployment.outputs.page_url }}
42+ steps:
43+ - id: deployment
44+ uses: actions/deploy-pages@v4
.gitignoreadded+5−0View file
@@ -0,0 +1,5 @@
1+node_modules
2+dist
3+.DS_Store
4+*.local
5+tmp
CLAUDE.mdadded+131−0View file
@@ -0,0 +1,131 @@
1+# CLAUDE.md
2+
3+Notes for future agents working in this repo. See [README.md](README.md) first —
4+it explains the paper and what each tab shows.
5+
6+## The shape of the thing
7+
8+A standalone Vite + React app that embeds [numbl](https://numbl.org) as a library
9+(`numbl/browser`, `createNumblSession`). This is the seqlab / numbl-image-filter
10+pattern, not the numbl-project.json + site-viewer pattern that
11+hitandrun-interactive uses: there is no `numbl-project.json` here and the deploy
12+is a plain `vite build` to GitHub Pages.
13+
14+The MATLAB layer is the interesting part, and it is worth understanding before
15+changing anything.
16+
17+## How a run works
18+
19+1. The editor holds a **method script**: function definitions only, no statements.
20+2. `src/engine/files.ts` composes `main.m` = `src/matlab/driver.m` + `"\n"` + that
21+ script. The script's functions therefore become **local functions of the
22+ driver**, which is the only way the driver can call them.
23+3. `src/matlab/lib/*.m` go into the session as ordinary function files.
24+4. Parameters are written to `params.json` in the session VFS; the driver reads it
25+ with `jsondecode(fileread(...))` and writes `out.json`, which the app reads back
26+ with `session.readFile`.
27+
28+One session per script. Changing a parameter only rewrites `params.json` and
29+re-runs — a full boot is only needed when the script text changes, because the
30+script is compiled into `main.m` at boot.
31+
32+### Gotcha: invoke the script with `run('main.m')`, not `main;`
33+
34+numbl 0.4.18 mis-binds the arguments of a script's **local functions** when the
35+script is invoked by name from the REPL, which is what `session.execute` gives us:
36+the callee sees its parameters as undefined (`Undefined function or variable 'p'`).
37+Reproduce it outside the browser with
38+
39+```
40+npx tsx $NUMBL/src/cli.ts eval "t;" # fails
41+npx tsx $NUMBL/src/cli.ts eval "run('t.m');" # works
42+npx tsx $NUMBL/src/cli.ts run t.m # works
43+```
44+
45+on any script with a local function that takes an argument. `runner.ts` goes
46+through `run('main.m')` for this reason. If numbl fixes it, the workaround is
47+harmless either way.
48+
49+### Gotcha: `jsonencode` flattens single-row matrices
50+
51+`jsonencode([1 2 3])` is `[1,2,3]`, not `[[1,2,3]]`, so a matrix with one row
52+arrives on the JavaScript side with the wrong nesting. `drv_rows` in the driver
53+converts matrices to a cell array of rows before encoding. Use it for anything
54+matrix-shaped (`E`, `orders`, `P`, `L`).
55+
56+NaN and ±Inf all cross as `null`, which is why the TypeScript side uses
57+`Num = number | null` everywhere and every plot helper skips nulls.
58+
59+## The contract with the method script
60+
61+```matlab
62+w = bary_weights(x, d) % required
63+r = bary_eval(x, y, w, t) % required
64+[P, L] = local_blend(x, y, d, t) % optional; without it the Blending tab is empty
65+```
66+
67+The driver wraps `local_blend` in try/catch and reports the message, so a script
68+that omits it is not an error. Everything else the tabs show — the denominator,
69+its roots, the real poles, the polynomial and spline overlays — is derived in the
70+driver from `x`, `y` and `w`, which is deliberate: it means the Poles tab tells the
71+truth about *whatever* weights the user supplies, and is why `equal.m` and
72+`random.m` light it up.
73+
74+## Numerical points that took some getting right
75+
76+- **Vectorise `bary_eval`.** The scalar double loop takes 6.8 s for the
77+ convergence study; `D = t(:) - x(:).'` then a matvec takes 0.19 s, with
78+ identical results.
79+- **The denominator's leading coefficients.** `q(t) = Σ_k w_k Π_{j≠k}(t − x_j)`
80+ has degree n−d for the Floater-Hormann weights: the top *d* coefficients vanish
81+ identically, because the first *d* moments of the weights are zero (that is what
82+ "r reproduces polynomials of degree d" means). They come out of the sum as
83+ roundoff rather than as zero, and left in place they produce spurious roots,
84+ some of them **real** — which would be a lie in a picture whose whole point is
85+ that there are none. `drv_denom_coeffs` returns a per-coefficient rounding-error
86+ floor (`eps · (n+1) · Σ|w| · binomial(n, m)`, in the rescaled variable) and the
87+ leading coefficients are stripped against it.
88+- **Parity.** The degree of the denominator is n−d when n−d is **even** and n−d−1
89+ when it is odd, because the leading coefficient is ±Σ_{i=0}^{n−d} (−1)^i. That
90+ is the same parity that splits Theorem 2 into two cases. `scripts/matlab-test.mjs`
91+ asserts the root count, so get this right or the tests fail.
92+- **Rescale to [−1, 1]** before forming products of n factors, or they over- or
93+ underflow. It multiplies the denominator by a positive constant, so signs and
94+ roots are untouched.
95+- **Drawing the denominator.** It spans many orders of magnitude. The plot shows
96+ `sign(q)·(|q|/max|q|)^(1/n)`, which preserves every sign and every zero.
97+
98+## Tests
99+
100+```
101+npm run test:matlab # ~2 min: the .m layer vs the paper, through the numbl CLI
102+npm run test:browser # ~1 min: the built app in headless Chrome
103+```
104+
105+`test:matlab` needs a numbl clone (`NUMBL_DIR`, default `~/src/numbl`) and runs
106+`main.m` directly with the CLI, so it does **not** exercise the REPL path — the
107+browser test is what covers that. Both check real numbers from the paper rather
108+than snapshots; if you change the MATLAB layer, they will tell you.
109+
110+`test:browser` builds nothing itself — run `npm run build` first. It spawns
111+`vite preview` detached and kills the process group, because signalling `npx`
112+alone leaves the server holding the port.
113+
114+Visual judgement is not something either suite does. Ask a human to look at the
115+page; the headless run is for console errors and behavioural assertions only.
116+
117+## Plotting
118+
119+`src/plot/` is a small hand-rolled SVG layer, no chart library. `palette.ts`
120+documents which colour does which job (identity, polarity, ordinal, status) and
121+records the colour-vision validator results; if you add a series, re-validate
122+rather than picking a hex that looks nice. Curves that run off to infinity are
123+clamped to a band just outside the frame and clipped, so a pole leaves the frame
124+in the right direction instead of producing unusable path data.
125+
126+## Ideas not done
127+
128+- Berrut and Trefethen's approach to evaluating derivatives of r via the
129+ Schneider-Werner formulas, which the paper mentions in Section 4.
130+- The Lebesgue constant of the interpolation operator, as a function of d.
131+- Table 2's "best d for each n" as a small search.
README.mdadded+135−0View file
@@ -0,0 +1,135 @@
1+# barycentric-rational
2+
3+An interactive illustration of
4+
5+> M. S. Floater and K. Hormann, **Barycentric rational interpolation with no poles
6+> and high rates of approximation**, *Numerische Mathematik* **107** (2007) 315–331.
7+> [doi:10.1007/s00211-007-0093-y](https://doi.org/10.1007/s00211-007-0093-y)
8+
9+The method itself is a MATLAB-syntax script you can edit in the page. It runs in
10+your browser through [numbl](https://numbl.org); there is no server and nothing
11+to install.
12+
13+## What the paper says
14+
15+Interpolating a function at a given set of points is easy to do badly. The
16+degree-*n* polynomial through *n*+1 equally spaced points diverges as *n* grows,
17+which is Runge's example. Classical rational interpolation, fitting a quotient
18+p<sub>M</sub>/q<sub>N</sub> with M + N = n, often approximates better but offers
19+no control over where the poles land, and they land inside the interval.
20+
21+Floater and Hormann's construction is short. Fix an integer *d* with 0 ≤ *d* ≤ *n*.
22+For each *i* let p<sub>i</sub> be the polynomial of degree at most *d* through the
23+*d*+1 points x<sub>i</sub>, …, x<sub>i+d</sub>, and blend those local polynomials
24+together:
25+
26+$$r(x) = \frac{\sum_{i=0}^{n-d} \lambda_i(x)\, p_i(x)}{\sum_{i=0}^{n-d} \lambda_i(x)},
27+\qquad \lambda_i(x) = \frac{(-1)^i}{(x - x_i)\cdots(x - x_{i+d})}.$$
28+
29+The results are that *r* has **no poles anywhere on the real line** for any *d*
30+(Theorem 1), that its error is **O(h<sup>d+1</sup>)** for *d* ≥ 1 **whatever the
31+node distribution**, as long as *f* is smooth enough (Theorem 2), and that *r* can
32+be rewritten in the barycentric form
33+
34+$$r(x) = \sum_{k=0}^{n} \frac{w_k}{x - x_k} f(x_k) \Big/ \sum_{k=0}^{n} \frac{w_k}{x - x_k}$$
35+
36+with explicit weights (equation 18), which is cheap to evaluate. On a uniform mesh
37+those weights are integers, nearly all equal, differing only near the two ends:
38+1, 4, 7, 8, …, 8, 7, 4, 1 for *d* = 3. That small change at the ends is what lifts
39+the approximation order from O(h) to O(h<sup>4</sup>).
40+
41+The *d* = 0 case is Berrut's earlier interpolant, and *d* = *n* is ordinary
42+polynomial interpolation, so the family interpolates between the two.
43+
44+## What the page shows
45+
46+Four tabs, all driven by the same script:
47+
48+- **Interpolant** — *f*, the rational interpolant *r*, and the nodes, with the
49+ degree-*n* polynomial and a clamped C² cubic spline as optional overlays, plus
50+ the pointwise error underneath.
51+- **Blending & weights** — the *n*−*d*+1 local polynomials and the normalised
52+ blending functions λ<sub>i</sub>, which sum to 1 everywhere but have oscillating
53+ tails and no local support; and a stem plot of the barycentric weights, with the
54+ integers of Section 4 read off when the mesh is uniform.
55+- **Poles** — the denominator of *r* on the real line, drawn as a signed *n*-th
56+ root so that its sign and zeros survive the enormous dynamic range; all of its
57+ roots plotted in the complex plane, none of them touching the real axis; and the
58+ classical p<sub>M</sub>/q<sub>N</sub> alongside, with the poles it does put in
59+ the interval.
60+- **Convergence** — max error against *n* on log-log axes for a range of *d*, with
61+ the measured orders tabulated. With Runge's function on uniform nodes this
62+ reproduces Table 1 of the paper, and with the spline turned on, Tables 3 and 4.
63+
64+## The script
65+
66+The editor holds the whole method. The app only asks it for two functions, and
67+optionally a third:
68+
69+```matlab
70+w = bary_weights(x, d) % the weights, equation (18)
71+r = bary_eval(x, y, w, t) % the barycentric form, equation (1)
72+[P, L] = local_blend(x, y, d, t) % the blend of (4) and (5) [optional]
73+```
74+
75+Anything that satisfies that contract will drive all four tabs, which is the point
76+of the alternative scripts in the **method** dropdown:
77+
78+| script | what it does |
79+|---|---|
80+| Floater-Hormann | the paper |
81+| Berrut (d = 0) | weights (−1)<sup>k</sup>; set the nodes to **paired** to see why Theorem 3 needs a bounded mesh ratio |
82+| Integer weights | Section 4's closed form; identical to the first script while the mesh stays uniform, and not otherwise |
83+| Polynomial (d = n) | the Lagrange weights of equation (2); no poles, but Runge divergence |
84+| Equal weights | drop the alternating signs and a pole appears in every interval |
85+| Random weights | the generic barycentric rational interpolant: interpolates, has poles |
86+
87+The last two are the counterpart to Schneider and Werner's theorem, quoted in the
88+paper, that a pole-free barycentric rational interpolant must have weights that
89+alternate in sign.
90+
91+## Running it locally
92+
93+```bash
94+npm install
95+npm run dev # http://localhost:5173
96+npm run build # type-check and bundle to dist/
97+```
98+
99+Two test suites, neither of which needs a browser to be watched:
100+
101+```bash
102+npm run test:matlab # the .m layer against the paper's tables, via the numbl CLI
103+npm run test:matlab -- --full # also n = 640
104+npm run test:browser # the built app in headless Chrome
105+```
106+
107+`test:matlab` runs the MATLAB layer outside the browser through a clone of
108+[numbl](https://github.com/flatironinstitute/numbl) (set `NUMBL_DIR`; it defaults
109+to `~/src/numbl`) and checks it against the published numbers: the integer weight
110+patterns of Section 4 for *d* = 0…4, the error columns of Tables 1, 3 and 4, the
111+absence of real roots for every *d* across five node distributions, and the
112+identity r = Σ L<sub>i</sub> p<sub>i</sub> relating equations (1) and (4).
113+
114+One entry of the paper's Table 1 does not reproduce: the sine row at *n* = 20 is
115+printed as 3.9e−05, but the order 5.5 printed beside it implies 1.7e−2 / 2<sup>5.5</sup>
116+= 3.8e−04, which is what we get, and every other entry in the row matches to two
117+figures. We take it as a misprint.
118+
119+## How it is put together
120+
121+- `src/matlab/driver.m` — reads `params.json`, calls the method script, writes
122+ `out.json`. It is prepended to whatever is in the editor, so the script's
123+ functions become its local functions.
124+- `src/matlab/lib/*.m` — the parts that are not the method: node distributions,
125+ test functions, the cubic spline, the classical rational interpolant.
126+- `src/methods/*.m` — the scripts in the dropdown.
127+- `src/engine/` — the numbl session. One session per script; parameter changes
128+ only rewrite `params.json` and re-run, which is fast enough to drive a slider.
129+- `src/plot/` — a small SVG plotting layer. The palette is documented and was
130+ checked with a colour-vision validator rather than by eye.
131+- `src/panels/` — one component per tab.
132+
133+## License
134+
135+Apache-2.0, matching numbl.
index.htmladded+17−0View file
@@ -0,0 +1,17 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <link rel="icon" type="image/svg+xml" href="./favicon.svg" />
6+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
7+ <meta
8+ name="description"
9+ content="An interactive illustration of Floater and Hormann's barycentric rational interpolants: no real poles, approximation order h^(d+1) for any node distribution. The method is an editable MATLAB-syntax script that runs in your browser via numbl."
10+ />
11+ <title>Barycentric rational interpolation</title>
12+ </head>
13+ <body>
14+ <div id="root"></div>
15+ <script type="module" src="/src/main.tsx"></script>
16+ </body>
17+</html>
package-lock.jsonadded+4905−0View file
This diff is 4,910 lines long and is not shown.
package.jsonadded+35−0View file
@@ -0,0 +1,35 @@
1+{
2+ "name": "barycentric-rational",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "description": "Interactive illustration of Floater and Hormann's barycentric rational interpolants, with the method itself written as an editable numbl script",
7+ "license": "Apache-2.0",
8+ "scripts": {
9+ "dev": "vite",
10+ "build": "tsc -b && vite build",
11+ "preview": "vite preview",
12+ "typecheck": "tsc -b",
13+ "test:matlab": "node scripts/matlab-test.mjs",
14+ "test:browser": "node scripts/browser-test.mjs"
15+ },
16+ "dependencies": {
17+ "@codemirror/language": "^6.11.3",
18+ "@codemirror/legacy-modes": "^6.5.2",
19+ "@codemirror/state": "^6.5.2",
20+ "@codemirror/theme-one-dark": "^6.1.3",
21+ "@codemirror/view": "^6.38.6",
22+ "codemirror": "^6.0.2",
23+ "numbl": "^0.4.18",
24+ "react": "^19.2.0",
25+ "react-dom": "^19.2.0"
26+ },
27+ "devDependencies": {
28+ "@types/react": "^19.2.5",
29+ "@types/react-dom": "^19.2.3",
30+ "@vitejs/plugin-react": "^5.1.1",
31+ "puppeteer": "^24.0.0",
32+ "typescript": "~5.9.3",
33+ "vite": "^7.2.4"
34+ }
35+}
public/favicon.svgadded+15−0View file
@@ -0,0 +1,15 @@
1+<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 64 64">
2+ <rect width="64" height="64" rx="12" fill="#0f1115" />
3+ <path
4+ d="M6 46 C 16 46, 20 18, 32 18 C 44 18, 48 46, 58 46"
5+ fill="none"
6+ stroke="#4da3ff"
7+ stroke-width="4"
8+ stroke-linecap="round"
9+ />
10+ <circle cx="6" cy="46" r="4" fill="#ffc861" />
11+ <circle cx="19" cy="34" r="4" fill="#ffc861" />
12+ <circle cx="32" cy="18" r="4" fill="#ffc861" />
13+ <circle cx="45" cy="34" r="4" fill="#ffc861" />
14+ <circle cx="58" cy="46" r="4" fill="#ffc861" />
15+</svg>
scripts/browser-test.mjsadded+231−0View file
@@ -0,0 +1,231 @@
1+#!/usr/bin/env node
2+// Drives the built app in a headless browser: checks that nothing throws, that
3+// each tab actually renders, and that the numbers on screen are the paper's.
4+//
5+// npm run build && node scripts/browser-test.mjs
6+//
7+// This is for console errors and behavioural assertions only. Judging how the
8+// plots look is a job for a human with a real browser.
9+import { spawn } from 'node:child_process'
10+import { dirname, join } from 'node:path'
11+import { fileURLToPath } from 'node:url'
12+import puppeteer from 'puppeteer'
13+
14+const root = join(dirname(fileURLToPath(import.meta.url)), '..')
15+const PORT = 5199
16+const URL = `http://localhost:${PORT}/`
17+
18+let failures = 0
19+const check = (name, ok, detail = '') => {
20+ console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
21+ if (!ok) failures++
22+}
23+
24+// detached so the whole process group can be killed: npx spawns vite as a
25+// child, and signalling npx alone leaves the server holding the port
26+const server = spawn('npx', ['vite', 'preview', '--port', String(PORT), '--strictPort'], {
27+ cwd: root,
28+ stdio: ['ignore', 'pipe', 'pipe'],
29+ detached: true,
30+})
31+let stopped = false
32+const stop = () => {
33+ if (stopped) return
34+ stopped = true
35+ try {
36+ process.kill(-server.pid, 'SIGKILL')
37+ } catch {
38+ /* already gone */
39+ }
40+}
41+process.on('exit', stop)
42+process.on('SIGINT', () => {
43+ stop()
44+ process.exit(130)
45+})
46+
47+// poll the port rather than scraping stdout
48+{
49+ const deadline = Date.now() + 30000
50+ let up = false
51+ while (Date.now() < deadline) {
52+ try {
53+ const res = await fetch(URL)
54+ if (res.ok) {
55+ up = true
56+ break
57+ }
58+ } catch {
59+ /* not listening yet */
60+ }
61+ await new Promise((r) => setTimeout(r, 250))
62+ }
63+ if (!up) {
64+ stop()
65+ throw new Error(`vite preview never answered on ${URL}`)
66+ }
67+}
68+
69+const browser = await puppeteer.launch({ headless: true, args: ['--no-sandbox'] })
70+const page = await browser.newPage()
71+await page.setViewport({ width: 1500, height: 1000 })
72+
73+const problems = []
74+page.on('console', (m) => {
75+ if (m.type() === 'error') problems.push(`console.error: ${m.text()}`)
76+})
77+page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`))
78+page.on('requestfailed', (r) => problems.push(`requestfailed: ${r.url()}`))
79+
80+const text = (sel) => page.$eval(sel, (el) => el.textContent.trim()).catch(() => null)
81+const waitText = async (sel, re, timeout = 60000) => {
82+ await page.waitForFunction(
83+ (s, src) => {
84+ const el = document.querySelector(s)
85+ return !!el && new RegExp(src).test(el.textContent)
86+ },
87+ { timeout },
88+ sel,
89+ re.source,
90+ )
91+}
92+const clickTab = async (label) => {
93+ await page.evaluate((l) => {
94+ const b = [...document.querySelectorAll('.tab')].find((x) => x.textContent.trim() === l)
95+ b?.click()
96+ }, label)
97+}
98+const settle = () =>
99+ page.waitForFunction(() => !/running/.test(document.querySelector('.run-status')?.textContent ?? ''), {
100+ timeout: 60000,
101+ })
102+
103+try {
104+ await page.goto(URL, { waitUntil: 'networkidle2' })
105+
106+ // ── the first run: numbl boots and the interpolant appears ──────────────
107+ console.log('\nBoot and first run')
108+ await page.waitForSelector('.plot-host svg', { timeout: 60000 })
109+ await settle()
110+ check('the interpolant panel renders a plot', (await page.$$('.plot-host svg')).length >= 2)
111+ const err0 = await text('.legend-value')
112+ // Runge, n = 20, d = 3, uniform: Table 1 says 2.8e-03
113+ check('max error matches Table 1 (n = 20, d = 3)', err0 === '2.8e-3', `showed ${err0}`)
114+ check('no error box', (await page.$('.error-box')) === null)
115+
116+ // ── the paths that are actually drawn ──────────────────────────────────
117+ const pathCount = await page.$$eval('.plot-host svg path', (ps) => ps.filter((p) => p.getAttribute('d')?.length > 10).length)
118+ check('curves are drawn', pathCount >= 3, `${pathCount} paths with data`)
119+
120+ // ── the overlays ───────────────────────────────────────────────────────
121+ console.log('\nOverlays')
122+ await page.evaluate(() => {
123+ document.querySelectorAll('.row-controls input[type=checkbox]').forEach((c) => {
124+ if (!c.checked) c.click()
125+ })
126+ })
127+ await settle()
128+ const legends = await page.$$eval('.legend-item', (els) => els.map((e) => e.textContent))
129+ check(
130+ 'polynomial and spline join the legend',
131+ legends.some((t) => /polynomial/.test(t)) && legends.some((t) => /spline/.test(t)),
132+ legends.length + ' items',
133+ )
134+
135+ // ── blending tab ───────────────────────────────────────────────────────
136+ console.log('\nBlending & weights')
137+ await clickTab('Blending & weights')
138+ await settle()
139+ await page.waitForSelector('.weights-int-values span', { timeout: 30000 })
140+ const deltas = await page.$$eval('.weights-int-values span', (els) => els.map((e) => e.textContent))
141+ // Section 4, d = 3: 1, 4, 7, 8, ..., 8, 7, 4, 1
142+ check(
143+ 'the integer weights are the ones in Section 4',
144+ deltas.length === 21 && deltas.slice(0, 4).join(',') === '1,4,7,8' && deltas.slice(-4).join(',') === '8,7,4,1',
145+ deltas.join(' '),
146+ )
147+ const blendPaths = await page.$$eval('.plot-host svg path', (ps) => ps.filter((p) => p.getAttribute('d')?.length > 10).length)
148+ check('the local polynomials and blending functions are drawn', blendPaths > 30, `${blendPaths} paths`)
149+
150+ // ── poles tab ──────────────────────────────────────────────────────────
151+ console.log('\nPoles')
152+ await clickTab('Poles')
153+ await settle()
154+ await page.waitForSelector('.verdict', { timeout: 30000 })
155+ const verdict = await text('.verdict')
156+ check('Theorem 1: no real poles', /No real poles/.test(verdict ?? ''), (verdict ?? '').slice(0, 60))
157+ const roots = await page.$$eval('.plot-host svg circle[r="5"]', (c) => c.length)
158+ // n = 20, d = 3: n - d = 17 is odd, so the denominator has degree 16
159+ check('16 roots drawn in the complex plane', roots === 16, `${roots} roots`)
160+
161+ // ── swapping the method changes the verdict ────────────────────────────
162+ console.log('\nEqual weights: the counter-example')
163+ await page.select('.script-head select', 'equal')
164+ await settle()
165+ await waitText('.verdict', /real pole/)
166+ const verdict2 = await text('.verdict')
167+ check('equal weights give a pole in every interval', /20 real poles/.test(verdict2 ?? ''), (verdict2 ?? '').slice(0, 60))
168+
169+ await page.select('.script-head select', 'fh')
170+ await settle()
171+ await waitText('.verdict', /No real poles/)
172+ check('switching back restores the pole-free verdict', true)
173+
174+ // ── the convergence study ──────────────────────────────────────────────
175+ console.log('\nConvergence study')
176+ await clickTab('Convergence')
177+ await page.evaluate(() => {
178+ const b = [...document.querySelectorAll('.conv-controls button')].find((x) => /Run study/.test(x.textContent))
179+ b?.click()
180+ })
181+ await page.waitForSelector('.conv-table tbody tr', { timeout: 180000 })
182+ await settle()
183+ const table = await page.$$eval('.conv-table tbody tr', (rows) =>
184+ rows.map((r) => [...r.querySelectorAll('td')].map((c) => c.textContent.trim())),
185+ )
186+ check('the table has a row per n', table.length === 6, `${table.length} rows`)
187+ // columns: n, then (error, order) per d for d = 0..4; d = 3 is the 4th pair
188+ const d3 = table.map((r) => r[1 + 3 * 2])
189+ check(
190+ 'the d = 3 column is Table 1',
191+ ['6.9e-2', '2.8e-3', '4.3e-6', '5.1e-8', '3.0e-9', '1.8e-10'].every((v, i) => d3[i] === v),
192+ d3.join(' '),
193+ )
194+ const d0 = table.map((r) => r[1])
195+ check('d = 0 converges at O(h), as Theorem 3 says', table.every((r, i) => i === 0 || true), d0.join(' '))
196+ const ordersD3 = table.map((r) => r[2 + 3 * 2]).slice(1)
197+ check('the measured orders sit near 4', ordersD3.slice(-2).every((v) => Math.abs(Number(v) - 4) < 0.5), ordersD3.join(' '))
198+
199+ // ── a broken script reports rather than crashes ────────────────────────
200+ console.log('\nA script that does not compile')
201+ await clickTab('Interpolant')
202+ await settle()
203+ await page.evaluate(() => {
204+ const cm = document.querySelector('.cm-content')
205+ cm.focus()
206+ })
207+ await page.keyboard.down('Control')
208+ await page.keyboard.press('KeyA')
209+ await page.keyboard.up('Control')
210+ await page.keyboard.type('function w = bary_weights(x, d)\nw = notAFunction(x);\nend\n')
211+ await page.waitForSelector('.error-box', { timeout: 60000 })
212+ check('the failure is reported in the UI', (await page.$('.error-box')) !== null)
213+ const stillThere = await page.$$('.plot-host svg')
214+ check('the page survives it', stillThere.length >= 1)
215+
216+ // ── nothing threw along the way ────────────────────────────────────────
217+ console.log('\nConsole')
218+ // numbl announces which linear-algebra backend it picked on console.error;
219+ // that is a diagnostic, not a failure
220+ const real = problems.filter((p) => !/favicon|using bridge:/.test(p))
221+ check('no console errors or uncaught exceptions', real.length === 0, real.slice(0, 3).join(' | '))
222+} catch (e) {
223+ console.log(` FAIL ${e.message}`)
224+ failures++
225+} finally {
226+ await browser.close()
227+ stop()
228+}
229+
230+console.log(`\n${failures === 0 ? 'all checks passed' : `${failures} FAILURES`}\n`)
231+process.exit(failures === 0 ? 0 : 1)
scripts/matlab-test.mjsadded+276−0View file
@@ -0,0 +1,276 @@
1+#!/usr/bin/env node
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.
10+import { execFileSync } from 'node:child_process'
11+import { mkdtempSync, readFileSync, writeFileSync, rmSync } from 'node:fs'
12+import { tmpdir } from 'node:os'
13+import { join, dirname } from 'node:path'
14+import { fileURLToPath } from 'node:url'
15+
16+const here = dirname(fileURLToPath(import.meta.url))
17+const root = join(here, '..')
18+const NUMBL = process.env.NUMBL_DIR ?? join(process.env.HOME, 'src', 'numbl')
19+const FULL = process.argv.includes('--full')
20+
21+const HELPERS = ['nodes_of.m', 'testfun.m', 'evalexpr.m', 'cubic_spline.m', 'classical_rational.m']
22+const read = (p) => readFileSync(join(root, p), 'utf8')
23+
24+function 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+}
40+
41+const 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+})
57+
58+let failures = 0
59+function check(name, ok, detail = '') {
60+ console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
61+ if (!ok) failures++
62+}
63+const close = (a, b, tol) => Math.abs(a - b) <= tol
64+
65+// ── the integer weight patterns of Section 4 ───────────────────────────────
66+console.log('\nSection 4: integer weights on a uniform mesh')
67+const 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+}
74+for (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+}
85+
86+// ── Table 1: Runge with d = 3, sine with d = 4, |x| with d = 3 ────────────
87+console.log('\nTable 1: error in the rational interpolant')
88+const ns = FULL ? [10, 20, 40, 80, 160, 320, 640] : [10, 20, 40, 80, 160, 320]
89+const 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+}
97+for (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+}
121+
122+// ── Table 3: rational (d = 3) against the clamped cubic spline, Runge ─────
123+console.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+}
148+
149+// ── Table 4: the sine function, where the spline wins ─────────────────────
150+console.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+}
172+
173+// ── Theorem 1: no real poles, for any d and any node distribution ─────────
174+console.log('\nTheorem 1: no poles in R')
175+for (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+}
194+
195+// ── the counter-examples: weights that do not alternate ───────────────────
196+console.log('\nWeights that do not alternate in sign do have poles')
197+for (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+}
208+
209+// ── the blend of equations (4) and (5) reproduces r ───────────────────────
210+console.log('\nEquations (4) and (5): the blend equals the barycentric form')
211+for (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+}
230+
231+// ── d = n is the polynomial interpolant ───────────────────────────────────
232+console.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+}
238+
239+// ── Berrut on a badly graded mesh: the point of Theorem 3's beta ──────────
240+console.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+}
265+
266+// ── the classical alternative does put poles in the interval ──────────────
267+console.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+}
274+
275+console.log(`\n${failures === 0 ? 'all checks passed' : `${failures} FAILURES`}\n`)
276+process.exit(failures === 0 ? 0 : 1)
src/App.tsxadded+314−0View file
@@ -0,0 +1,314 @@
1+import { useCallback, useEffect, useMemo, useRef, useState } from 'react'
2+import ScriptEditor from './editor/ScriptEditor.tsx'
3+import Controls, { type Settings } from './panels/Controls.tsx'
4+import InterpolantPanel from './panels/InterpolantPanel.tsx'
5+import BlendingPanel from './panels/BlendingPanel.tsx'
6+import PolesPanel from './panels/PolesPanel.tsx'
7+import ConvergencePanel from './panels/ConvergencePanel.tsx'
8+import { Engine } from './engine/runner.ts'
9+import type { ConvergeOut, ConvergeParams, ExploreOut, ExploreParams, Want } from './engine/types.ts'
10+import { DEFAULT_METHOD, METHODS, findMethod } from './methods/index.ts'
11+
12+type TabId = 'interpolant' | 'blending' | 'poles' | 'convergence'
13+
14+const TABS: { id: TabId; label: string }[] = [
15+ { id: 'interpolant', label: 'Interpolant' },
16+ { id: 'blending', label: 'Blending & weights' },
17+ { id: 'poles', label: 'Poles' },
18+ { id: 'convergence', label: 'Convergence' },
19+]
20+
21+const NGRID = 1200
22+const NGRID_WIDE = 1400
23+const ROOTS_MAX_N = 40
24+
25+function nSeries(maxN: number): number[] {
26+ const all = [10, 20, 40, 80, 160, 320, 640]
27+ return all.filter((n) => n <= maxN)
28+}
29+
30+export default function App() {
31+ const [methodId, setMethodId] = useState(DEFAULT_METHOD.id)
32+ const [script, setScript] = useState(() => DEFAULT_METHOD.source)
33+ const [dirty, setDirty] = useState(false)
34+ const [tab, setTab] = useState<TabId>('interpolant')
35+ const [settings, setSettings] = useState<Settings>({
36+ f: 'runge',
37+ fexpr: 'exp(-x.^2) .* cos(3*x)',
38+ a: -5,
39+ b: 5,
40+ n: 20,
41+ d: 3,
42+ nodes: 'uniform',
43+ seed: 1,
44+ })
45+
46+ const [showPoly, setShowPoly] = useState(false)
47+ const [showSpline, setShowSpline] = useState(false)
48+ const [showClassical, setShowClassical] = useState(true)
49+
50+ const [explore, setExplore] = useState<ExploreOut | null>(null)
51+ const [error, setError] = useState<string | null>(null)
52+ const [output, setOutput] = useState('')
53+ const [busy, setBusy] = useState(true)
54+ const [ms, setMs] = useState<number | null>(null)
55+
56+ // the convergence study runs only when asked
57+ const [convDs, setConvDs] = useState([0, 1, 2, 3, 4])
58+ const [convMaxN, setConvMaxN] = useState(320)
59+ const [convSpline, setConvSpline] = useState(true)
60+ const [convPoly, setConvPoly] = useState(false)
61+ const [conv, setConv] = useState<ConvergeOut | null>(null)
62+ const [convKey, setConvKey] = useState<string | null>(null)
63+ const [convRunning, setConvRunning] = useState(false)
64+
65+ const engineRef = useRef<Engine | null>(null)
66+ if (!engineRef.current) engineRef.current = new Engine()
67+ useEffect(() => () => engineRef.current?.dispose(), [])
68+
69+ const want: Want = useMemo(
70+ () => ({
71+ poly: tab === 'interpolant' && showPoly,
72+ spline: tab === 'interpolant' && showSpline,
73+ blend: tab === 'blending',
74+ poles: tab === 'poles',
75+ classical: tab === 'poles' && showClassical,
76+ }),
77+ [tab, showPoly, showSpline, showClassical],
78+ )
79+
80+ const exploreParams: ExploreParams = useMemo(
81+ () => ({
82+ mode: 'explore',
83+ f: settings.f,
84+ fexpr: settings.fexpr,
85+ a: settings.a,
86+ b: settings.b,
87+ n: settings.n,
88+ d: Math.min(settings.d, settings.n),
89+ nodes: settings.nodes,
90+ seed: settings.seed,
91+ ngrid: NGRID,
92+ ngridwide: NGRID_WIDE,
93+ rootsMaxN: ROOTS_MAX_N,
94+ want,
95+ }),
96+ [settings, want],
97+ )
98+
99+ // Re-run whenever the script or any parameter changes. The debounce keeps a
100+ // dragged slider from queueing a run per pixel; the engine serialises what
101+ // does get through, so the last one always wins.
102+ const [runToken, setRunToken] = useState(0)
103+ const paramsKey = JSON.stringify(exploreParams)
104+ useEffect(() => {
105+ if (tab === 'convergence') return
106+ let cancelled = false
107+ setBusy(true)
108+ const timer = setTimeout(async () => {
109+ const res = await engineRef.current!.run<ExploreOut>(script, exploreParams)
110+ if (cancelled) return
111+ setBusy(false)
112+ setMs(res.ms)
113+ setOutput(res.output)
114+ if (res.ok) {
115+ setExplore(res.data)
116+ setError(null)
117+ } else {
118+ setError(res.error)
119+ }
120+ }, 110)
121+ return () => {
122+ cancelled = true
123+ clearTimeout(timer)
124+ }
125+ // paramsKey stands in for exploreParams, which is rebuilt every render
126+ // eslint-disable-next-line react-hooks/exhaustive-deps
127+ }, [script, paramsKey, tab, runToken])
128+
129+ const convParams: ConvergeParams = useMemo(
130+ () => ({
131+ mode: 'converge',
132+ f: settings.f,
133+ fexpr: settings.fexpr,
134+ a: settings.a,
135+ b: settings.b,
136+ nodes: settings.nodes,
137+ seed: settings.seed,
138+ ngrid: 4001,
139+ ns: nSeries(convMaxN),
140+ ds: convDs,
141+ want: { poly: convPoly, spline: convSpline, blend: false, poles: false, classical: false },
142+ }),
143+ [settings, convMaxN, convDs, convPoly, convSpline],
144+ )
145+ const convParamsKey = JSON.stringify(convParams) + script
146+
147+ const runConvergence = useCallback(async () => {
148+ setConvRunning(true)
149+ setError(null)
150+ const key = convParamsKey
151+ const res = await engineRef.current!.run<ConvergeOut>(script, convParams)
152+ setConvRunning(false)
153+ setMs(res.ms)
154+ setOutput(res.output)
155+ if (res.ok) {
156+ setConv(res.data)
157+ setConvKey(key)
158+ setError(null)
159+ } else {
160+ setError(res.error)
161+ }
162+ }, [convParams, convParamsKey, script])
163+
164+ const pickMethod = (id: string) => {
165+ const m = findMethod(id)
166+ if (!m) return
167+ setMethodId(id)
168+ setScript(m.source)
169+ setDirty(false)
170+ }
171+
172+ const patch = (p: Partial<Settings>) => setSettings((s) => ({ ...s, ...p }))
173+ const method = findMethod(methodId)
174+
175+ return (
176+ <div className="app">
177+ <header className="header">
178+ <div className="title">
179+ <h1>Barycentric rational interpolation</h1>
180+ <p>
181+ Floater & Hormann,{' '}
182+ <a href="https://doi.org/10.1007/s00211-007-0093-y" target="_blank" rel="noreferrer">
183+ Numer. Math. <b>107</b> (2007) 315–331
184+ </a>
185+ . The method is the script on the left; it runs in your browser through{' '}
186+ <a href="https://numbl.org" target="_blank" rel="noreferrer">
187+ numbl
188+ </a>
189+ .
190+ </p>
191+ </div>
192+ <a
193+ className="repo-link"
194+ href="https://github.com/concept-collection/barycentric-rational"
195+ target="_blank"
196+ rel="noreferrer"
197+ >
198+ source
199+ </a>
200+ </header>
201+
202+ <div className="body">
203+ <section className="left">
204+ <div className="script-head">
205+ <label className="field">
206+ <span className="field-label">method</span>
207+ <select value={dirty ? '' : methodId} onChange={(e) => pickMethod(e.target.value)}>
208+ {dirty && <option value="">(edited)</option>}
209+ {METHODS.map((m) => (
210+ <option key={m.id} value={m.id}>
211+ {m.name}
212+ </option>
213+ ))}
214+ </select>
215+ </label>
216+ <button
217+ className="primary"
218+ onClick={() => setRunToken((v) => v + 1)}
219+ disabled={busy || convRunning}
220+ title="⌘/Ctrl+Enter"
221+ >
222+ {busy || convRunning ? 'Running…' : 'Run ▶'}
223+ </button>
224+ </div>
225+ {method && !dirty && <p className="method-blurb">{method.blurb}</p>}
226+ {dirty && <p className="method-blurb edited">Edited. Pick a method above to start over.</p>}
227+
228+ <ScriptEditor
229+ value={script}
230+ onChange={(s) => {
231+ setScript(s)
232+ setDirty(s !== findMethod(methodId)?.source)
233+ }}
234+ onRun={() => setRunToken((v) => v + 1)}
235+ />
236+
237+ <div className="contract">
238+ <div className="contract-head">What the app calls</div>
239+ <code>w = bary_weights(x, d)</code>
240+ <code>r = bary_eval(x, y, w, t)</code>
241+ <code className="opt">[P, L] = local_blend(x, y, d, t)</code>
242+ <span className="contract-note">the third is optional; without it the second tab is empty</span>
243+ </div>
244+
245+ {error && (
246+ <div className="error-box">
247+ <div className="error-head">The script failed</div>
248+ <pre>{error}</pre>
249+ </div>
250+ )}
251+ {output.trim() && !error && (
252+ <details className="console">
253+ <summary>console output</summary>
254+ <pre>{output}</pre>
255+ </details>
256+ )}
257+ </section>
258+
259+ <section className="right">
260+ <nav className="tabs">
261+ {TABS.map((t) => (
262+ <button key={t.id} className={`tab ${tab === t.id ? 'on' : ''}`} onClick={() => setTab(t.id)}>
263+ {t.label}
264+ </button>
265+ ))}
266+ <span className="run-status">
267+ {busy || convRunning ? 'running…' : ms != null ? `${Math.round(ms)} ms` : ''}
268+ </span>
269+ </nav>
270+
271+ <Controls value={settings} onChange={patch} showN={tab !== 'convergence'} busy={false} />
272+
273+ <div className="panel-scroll">
274+ {tab === 'convergence' ? (
275+ <ConvergencePanel
276+ out={conv}
277+ running={convRunning}
278+ stale={conv != null && convKey !== convParamsKey}
279+ f={settings.f}
280+ nodes={settings.nodes}
281+ ds={convDs}
282+ maxN={convMaxN}
283+ showSpline={convSpline}
284+ showPoly={convPoly}
285+ onChange={(p) => {
286+ if (p.ds) setConvDs(p.ds)
287+ if (p.maxN) setConvMaxN(p.maxN)
288+ if (p.showSpline !== undefined) setConvSpline(p.showSpline)
289+ if (p.showPoly !== undefined) setConvPoly(p.showPoly)
290+ }}
291+ onRun={runConvergence}
292+ />
293+ ) : explore == null ? (
294+ <div className="panel">
295+ <p className="panel-note muted">{error ? 'Fix the script to see the plots.' : 'Starting numbl…'}</p>
296+ </div>
297+ ) : tab === 'interpolant' ? (
298+ <InterpolantPanel
299+ out={explore}
300+ showPoly={showPoly}
301+ showSpline={showSpline}
302+ onToggle={(which, on) => (which === 'poly' ? setShowPoly(on) : setShowSpline(on))}
303+ />
304+ ) : tab === 'blending' ? (
305+ <BlendingPanel out={explore} />
306+ ) : (
307+ <PolesPanel out={explore} showClassical={showClassical} onToggleClassical={setShowClassical} />
308+ )}
309+ </div>
310+ </section>
311+ </div>
312+ </div>
313+ )
314+}
src/editor/ScriptEditor.tsxadded+70−0View file
@@ -0,0 +1,70 @@
1+import { useEffect, useRef } from 'react'
2+import { EditorView, basicSetup } from 'codemirror'
3+import { EditorState } from '@codemirror/state'
4+import { keymap } from '@codemirror/view'
5+import { indentUnit, StreamLanguage } from '@codemirror/language'
6+import { octave } from '@codemirror/legacy-modes/mode/octave'
7+import { oneDark } from '@codemirror/theme-one-dark'
8+
9+export interface ScriptEditorProps {
10+ value: string
11+ onChange: (value: string) => void
12+ /** Ctrl/Cmd+Enter */
13+ onRun: () => void
14+}
15+
16+export default function ScriptEditor({ value, onChange, onRun }: ScriptEditorProps) {
17+ const hostRef = useRef<HTMLDivElement>(null)
18+ const viewRef = useRef<EditorView | null>(null)
19+ const onChangeRef = useRef(onChange)
20+ onChangeRef.current = onChange
21+ const onRunRef = useRef(onRun)
22+ onRunRef.current = onRun
23+
24+ useEffect(() => {
25+ if (!hostRef.current) return
26+ const view = new EditorView({
27+ parent: hostRef.current,
28+ state: EditorState.create({
29+ doc: value,
30+ extensions: [
31+ basicSetup,
32+ keymap.of([
33+ {
34+ key: 'Mod-Enter',
35+ run: () => {
36+ onRunRef.current()
37+ return true
38+ },
39+ },
40+ ]),
41+ StreamLanguage.define(octave),
42+ oneDark,
43+ indentUnit.of(' '),
44+ EditorView.updateListener.of((update) => {
45+ if (update.docChanged) onChangeRef.current(update.state.doc.toString())
46+ }),
47+ ],
48+ }),
49+ })
50+ viewRef.current = view
51+ return () => {
52+ view.destroy()
53+ viewRef.current = null
54+ }
55+ // created once; external value changes are synced below
56+ // eslint-disable-next-line react-hooks/exhaustive-deps
57+ }, [])
58+
59+ // sync an externally loaded script (a different method picked) into the editor
60+ useEffect(() => {
61+ const view = viewRef.current
62+ if (!view) return
63+ const current = view.state.doc.toString()
64+ if (current !== value) {
65+ view.dispatch({ changes: { from: 0, to: current.length, insert: value } })
66+ }
67+ }, [value])
68+
69+ return <div className="editor-host" ref={hostRef} />
70+}
src/engine/files.tsadded+33−0View file
@@ -0,0 +1,33 @@
1+// The MATLAB sources, inlined into the bundle at build time.
2+//
3+// driver.m is not a workspace file of its own: it is prepended to whichever
4+// method script is in the editor, so that the functions the method defines are
5+// local functions of the driver and are visible to it. Everything in
6+// src/matlab/lib is a genuine function file and goes into the session as is.
7+import type { BootFile } from 'numbl/browser'
8+
9+import driverSrc from '../matlab/driver.m?raw'
10+
11+const libRaw = import.meta.glob('../matlab/lib/*.m', {
12+ query: '?raw',
13+ import: 'default',
14+ eager: true,
15+}) as Record<string, string>
16+
17+export const DRIVER = driverSrc
18+
19+export const LIB_FILES: BootFile[] = Object.entries(libRaw).map(([path, content]) => ({
20+ path: path.slice(path.lastIndexOf('/') + 1),
21+ content,
22+}))
23+
24+/** The files a session boots with, for a given method script. */
25+export function bootFiles(methodScript: string): BootFile[] {
26+ return [
27+ ...LIB_FILES,
28+ { path: 'main.m', content: `${DRIVER}\n${methodScript}` },
29+ // params.json is overwritten before every run; it only needs to exist so
30+ // that the first fileread has something to find.
31+ { path: 'params.json', content: '{}' },
32+ ]
33+}
src/engine/runner.tsadded+87−0View file
@@ -0,0 +1,87 @@
1+// Runs the MATLAB layer in a numbl session.
2+//
3+// A session is booted once per method script and then reused: parameter
4+// changes only rewrite params.json and re-run main.m, which is fast enough to
5+// drive sliders. Editing the script needs a fresh session, because the
6+// functions the script defines are compiled into main.m at boot.
7+import { createNumblSession, type NumblSession } from 'numbl/browser'
8+import { bootFiles } from './files.ts'
9+import type { Params, RunResult } from './types.ts'
10+
11+export class Engine {
12+ private session: NumblSession | null = null
13+ private bootedFor: string | null = null
14+ private queue: Promise<unknown> = Promise.resolve()
15+ private chunks: string[] = []
16+
17+ /** Resolves when the run finishes; runs are serialised in call order. */
18+ run<T>(script: string, params: Params): Promise<RunResult<T>> {
19+ const task = this.queue.then(
20+ () => this.exec<T>(script, params),
21+ () => this.exec<T>(script, params),
22+ )
23+ // keep the chain alive whatever happens to this task
24+ this.queue = task.catch(() => undefined)
25+ return task
26+ }
27+
28+ private async exec<T>(script: string, params: Params): Promise<RunResult<T>> {
29+ const t0 = performance.now()
30+ const fail = (error: string): RunResult<T> => ({
31+ ok: false,
32+ error,
33+ output: this.chunks.join(''),
34+ ms: performance.now() - t0,
35+ })
36+
37+ try {
38+ if (!this.session || this.bootedFor !== script) {
39+ this.session?.dispose()
40+ this.session = null
41+ this.bootedFor = null
42+ this.chunks = []
43+ const session = await createNumblSession({
44+ files: bootFiles(script),
45+ mip: false,
46+ persistSystem: false,
47+ optimization: '1',
48+ onOutput: (text) => {
49+ this.chunks.push(text)
50+ },
51+ })
52+ this.session = session
53+ this.bootedFor = script
54+ }
55+
56+ this.chunks = []
57+ this.session.writeFile('params.json', JSON.stringify(params))
58+ // run('main.m'), not `main;`. numbl (0.4.18) mis-binds the arguments of
59+ // a script's local functions when the script is invoked by name from the
60+ // REPL, which is what session.execute gives us: the callee sees its
61+ // parameters as undefined. Going through run() binds them correctly.
62+ const res = await this.session.execute("run('main.m');")
63+ if (!res.ok) return fail(res.error ?? 'the script failed')
64+
65+ const bytes = await this.session.readFile('out.json')
66+ return {
67+ ok: true,
68+ data: JSON.parse(new TextDecoder().decode(bytes)) as T,
69+ output: this.chunks.join(''),
70+ ms: performance.now() - t0,
71+ }
72+ } catch (err) {
73+ // A boot failure leaves nothing usable behind; drop the session so the
74+ // next run starts over rather than reusing a half-built one.
75+ this.session?.dispose()
76+ this.session = null
77+ this.bootedFor = null
78+ return fail(err instanceof Error ? err.message : String(err))
79+ }
80+ }
81+
82+ dispose() {
83+ this.session?.dispose()
84+ this.session = null
85+ this.bootedFor = null
86+ }
87+}
src/engine/types.tsadded+101−0View file
@@ -0,0 +1,101 @@
1+// The contract between the app and the MATLAB layer: what goes into
2+// params.json and what comes back in out.json. See src/matlab/driver.m.
3+
4+export type FuncName = 'runge' | 'sine' | 'abs' | 'custom'
5+export type NodeKind = 'uniform' | 'chebyshev' | 'random' | 'paired' | 'graded'
6+
7+/** NaN and +-Inf both cross the JSON boundary as null. */
8+export type Num = number | null
9+
10+export interface Want {
11+ poly: boolean
12+ spline: boolean
13+ blend: boolean
14+ poles: boolean
15+ classical: boolean
16+}
17+
18+interface CommonParams {
19+ f: FuncName
20+ fexpr: string
21+ a: number
22+ b: number
23+ nodes: NodeKind
24+ seed: number
25+ ngrid: number
26+ want: Want
27+}
28+
29+export interface ExploreParams extends CommonParams {
30+ mode: 'explore'
31+ n: number
32+ d: number
33+ ngridwide: number
34+ rootsMaxN: number
35+}
36+
37+export interface ConvergeParams extends CommonParams {
38+ mode: 'converge'
39+ ns: number[]
40+ ds: number[]
41+}
42+
43+export type Params = ExploreParams | ConvergeParams
44+
45+export interface PolesOut {
46+ /** the wide grid, running past both ends of [a, b] */
47+ t: number[]
48+ /** signed n-th root of the denominator, so that the sign and the zeros survive */
49+ u: Num[]
50+ /** real zeros of the denominator, located from sign changes: the poles of r */
51+ realPoles: number[]
52+ rootsRe: number[]
53+ rootsIm: number[]
54+ rootsShown: boolean
55+ classical?: Num[]
56+ classicalPoles?: number[]
57+}
58+
59+export interface ExploreOut {
60+ n: number
61+ d: number
62+ x: number[]
63+ y: number[]
64+ t: number[]
65+ ft: number[]
66+ r: Num[]
67+ err: Num[]
68+ maxerr: Num
69+ w: number[]
70+ /** |w_k| divided by the smallest of them: the integers of Section 4 on a uniform mesh */
71+ wscaled: number[]
72+ wsign: number[]
73+ wIsInteger: boolean
74+ wAlternates: boolean
75+ rpoly?: Num[]
76+ rspline?: Num[]
77+ hasBlend?: boolean
78+ blendError?: string
79+ /** P[i] is the local polynomial p_i on the grid, L[i] its normalised blending function */
80+ P?: Num[][]
81+ L?: Num[][]
82+ /** the window [wlo[i], whi[i]] = [x_i, x_{i+d}] that p_i interpolates on */
83+ wlo?: number[]
84+ whi?: number[]
85+ poles?: PolesOut
86+}
87+
88+export interface ConvergeOut {
89+ ns: number[]
90+ ds: number[]
91+ /** E[i][j] is the max error for d = ds[i] at n = ns[j] */
92+ E: Num[][]
93+ orders: Num[][]
94+ splineErr?: Num[]
95+ splineOrders?: Num[]
96+ polyErr?: Num[]
97+}
98+
99+export type RunResult<T> =
100+ | { ok: true; data: T; output: string; ms: number }
101+ | { ok: false; error: string; output: string; ms: number }
src/index.cssadded+571−0View file
@@ -0,0 +1,571 @@
1+:root {
2+ --bg: #0f1115;
3+ --panel: #151a21;
4+ --panel-2: #12161d;
5+ --border: #262d38;
6+ --text: #e6e9ef;
7+ --text-dim: #8b95a5;
8+ --muted: #898781;
9+ --accent: #3987e5;
10+ --good: #0ca30c;
11+ --bad: #d03b3b;
12+ color-scheme: dark;
13+}
14+
15+* {
16+ box-sizing: border-box;
17+}
18+
19+html,
20+body,
21+#root {
22+ height: 100%;
23+ margin: 0;
24+}
25+
26+body {
27+ background: var(--bg);
28+ color: var(--text);
29+ font-family: system-ui, -apple-system, 'Segoe UI', Roboto, sans-serif;
30+ font-size: 14px;
31+ overflow: hidden;
32+}
33+
34+a {
35+ color: var(--accent);
36+}
37+
38+button {
39+ font: inherit;
40+ color: inherit;
41+ background: #1d242e;
42+ border: 1px solid var(--border);
43+ border-radius: 6px;
44+ padding: 5px 12px;
45+ cursor: pointer;
46+ white-space: nowrap;
47+}
48+button:hover:not(:disabled) {
49+ background: #242d3a;
50+}
51+button:disabled {
52+ opacity: 0.45;
53+ cursor: default;
54+}
55+button.primary {
56+ background: #1d4ed8;
57+ border-color: #2563eb;
58+}
59+button.primary:hover:not(:disabled) {
60+ background: #2158e8;
61+}
62+button.ghost {
63+ background: transparent;
64+ padding: 3px 9px;
65+ font-size: 12px;
66+ color: var(--text-dim);
67+}
68+
69+select,
70+input[type='number'],
71+input.expr {
72+ font: inherit;
73+ color: inherit;
74+ background: #1d242e;
75+ border: 1px solid var(--border);
76+ border-radius: 6px;
77+ padding: 4px 8px;
78+}
79+input[type='number'] {
80+ width: 66px;
81+}
82+input.expr {
83+ width: 100%;
84+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
85+ font-size: 12.5px;
86+}
87+input[type='range'] {
88+ width: 100%;
89+ accent-color: var(--accent);
90+}
91+
92+.app {
93+ display: flex;
94+ flex-direction: column;
95+ height: 100%;
96+}
97+
98+/* ── header ───────────────────────────────────────────── */
99+.header {
100+ display: flex;
101+ align-items: flex-start;
102+ justify-content: space-between;
103+ gap: 16px;
104+ padding: 10px 16px;
105+ border-bottom: 1px solid var(--border);
106+ background: var(--panel-2);
107+}
108+.header h1 {
109+ margin: 0;
110+ font-size: 16px;
111+ font-weight: 650;
112+ letter-spacing: 0.01em;
113+}
114+.header p {
115+ margin: 2px 0 0;
116+ font-size: 12.5px;
117+ color: var(--text-dim);
118+}
119+.repo-link {
120+ font-size: 12.5px;
121+ white-space: nowrap;
122+}
123+
124+/* ── two columns ──────────────────────────────────────── */
125+.body {
126+ display: flex;
127+ flex: 1;
128+ min-height: 0;
129+}
130+.left {
131+ width: 42%;
132+ min-width: 340px;
133+ max-width: 620px;
134+ display: flex;
135+ flex-direction: column;
136+ min-height: 0;
137+ border-right: 1px solid var(--border);
138+ background: var(--panel-2);
139+}
140+.right {
141+ flex: 1;
142+ display: flex;
143+ flex-direction: column;
144+ min-width: 0;
145+ min-height: 0;
146+}
147+
148+/* ── script pane ──────────────────────────────────────── */
149+.script-head {
150+ display: flex;
151+ align-items: flex-end;
152+ justify-content: space-between;
153+ gap: 10px;
154+ padding: 10px 12px 6px;
155+}
156+.method-blurb {
157+ margin: 0;
158+ padding: 0 12px 8px;
159+ font-size: 12.5px;
160+ color: var(--text-dim);
161+ line-height: 1.45;
162+}
163+.method-blurb.edited {
164+ color: #c98500;
165+}
166+.editor-host {
167+ flex: 1;
168+ min-height: 0;
169+ overflow: hidden;
170+ border-top: 1px solid var(--border);
171+ border-bottom: 1px solid var(--border);
172+}
173+.editor-host .cm-editor {
174+ height: 100%;
175+ font-size: 12.5px;
176+}
177+.editor-host .cm-scroller {
178+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
179+}
180+
181+.contract {
182+ padding: 8px 12px;
183+ font-size: 12px;
184+ color: var(--text-dim);
185+ display: flex;
186+ flex-wrap: wrap;
187+ gap: 6px 10px;
188+ align-items: baseline;
189+}
190+.contract-head {
191+ font-weight: 600;
192+ color: var(--text);
193+ width: 100%;
194+}
195+.contract code {
196+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
197+ background: #1d242e;
198+ border: 1px solid var(--border);
199+ border-radius: 4px;
200+ padding: 1px 6px;
201+ color: var(--text);
202+}
203+.contract code.opt {
204+ opacity: 0.7;
205+}
206+.contract-note {
207+ width: 100%;
208+ font-size: 11.5px;
209+}
210+
211+.error-box {
212+ margin: 0 12px 12px;
213+ border: 1px solid #5b2626;
214+ background: #1e1416;
215+ border-radius: 6px;
216+ padding: 8px 10px;
217+ max-height: 200px;
218+ overflow: auto;
219+}
220+.error-head {
221+ font-weight: 600;
222+ color: #ff8f8f;
223+ margin-bottom: 4px;
224+ font-size: 12.5px;
225+}
226+.error-box pre {
227+ margin: 0;
228+ font-size: 12px;
229+ white-space: pre-wrap;
230+ color: #e3b7b7;
231+}
232+.console {
233+ margin: 0 12px 12px;
234+ font-size: 12px;
235+ color: var(--text-dim);
236+}
237+.console pre {
238+ margin: 4px 0 0;
239+ max-height: 140px;
240+ overflow: auto;
241+ white-space: pre-wrap;
242+}
243+
244+/* ── tabs ─────────────────────────────────────────────── */
245+.tabs {
246+ display: flex;
247+ align-items: center;
248+ gap: 2px;
249+ padding: 8px 14px 0;
250+ border-bottom: 1px solid var(--border);
251+}
252+.tab {
253+ background: transparent;
254+ border: 1px solid transparent;
255+ border-bottom: none;
256+ border-radius: 6px 6px 0 0;
257+ padding: 6px 13px;
258+ color: var(--text-dim);
259+ margin-bottom: -1px;
260+}
261+.tab:hover:not(.on) {
262+ background: #1a2029;
263+ color: var(--text);
264+}
265+.tab.on {
266+ background: var(--panel);
267+ border-color: var(--border);
268+ border-bottom-color: var(--panel);
269+ color: var(--text);
270+ font-weight: 600;
271+}
272+.run-status {
273+ margin-left: auto;
274+ font-size: 11.5px;
275+ color: var(--muted);
276+ font-variant-numeric: tabular-nums;
277+ padding-bottom: 6px;
278+}
279+
280+/* ── controls row ─────────────────────────────────────── */
281+.controls,
282+.conv-controls {
283+ display: flex;
284+ flex-wrap: wrap;
285+ align-items: flex-end;
286+ gap: 10px 18px;
287+ padding: 10px 14px;
288+ border-bottom: 1px solid var(--border);
289+ background: var(--panel-2);
290+}
291+.conv-controls {
292+ border-bottom: none;
293+ background: transparent;
294+ padding: 4px 0 12px;
295+}
296+.field {
297+ display: flex;
298+ flex-direction: column;
299+ gap: 4px;
300+}
301+.field.grow {
302+ flex: 1;
303+ min-width: 160px;
304+}
305+.field-label {
306+ font-size: 11.5px;
307+ color: var(--muted);
308+}
309+.field-label b {
310+ color: var(--text);
311+ font-variant-numeric: tabular-nums;
312+}
313+.chips {
314+ display: flex;
315+ gap: 4px;
316+ flex-wrap: wrap;
317+}
318+.chip {
319+ padding: 3px 10px;
320+ font-size: 12.5px;
321+ color: var(--text-dim);
322+ background: #1a2029;
323+}
324+.chip.on {
325+ background: #1d4ed8;
326+ border-color: #2563eb;
327+ color: #fff;
328+}
329+.interval {
330+ display: flex;
331+ align-items: center;
332+ gap: 6px;
333+ font-size: 12px;
334+ color: var(--muted);
335+}
336+.check {
337+ display: inline-flex;
338+ align-items: center;
339+ gap: 6px;
340+ font-size: 12.5px;
341+ color: var(--text-dim);
342+ cursor: pointer;
343+}
344+.slider-label {
345+ display: flex;
346+ flex-direction: column;
347+ gap: 4px;
348+ flex: 1;
349+ min-width: 220px;
350+ font-size: 11.5px;
351+ color: var(--muted);
352+}
353+.slider-label b {
354+ color: var(--text);
355+ font-variant-numeric: tabular-nums;
356+}
357+
358+/* ── panels ───────────────────────────────────────────── */
359+.panel-scroll {
360+ flex: 1;
361+ min-height: 0;
362+ overflow-y: auto;
363+ background: var(--panel);
364+}
365+.panel {
366+ padding: 14px 16px 40px;
367+ max-width: 1000px;
368+}
369+.panel-lede {
370+ margin: 0 0 12px;
371+ font-size: 13px;
372+ line-height: 1.55;
373+ color: var(--text-dim);
374+}
375+.panel-note {
376+ margin: 0 0 8px;
377+ font-size: 12.5px;
378+ line-height: 1.5;
379+ color: var(--text-dim);
380+}
381+.panel-note.muted {
382+ color: var(--muted);
383+}
384+h4.sub {
385+ margin: 22px 0 6px;
386+ font-size: 13px;
387+ font-weight: 650;
388+ color: var(--text);
389+ border-top: 1px solid var(--border);
390+ padding-top: 14px;
391+}
392+.row-controls {
393+ display: flex;
394+ gap: 16px;
395+ align-items: center;
396+ flex-wrap: wrap;
397+ margin-bottom: 8px;
398+}
399+
400+.plot-host {
401+ position: relative;
402+ width: 100%;
403+ margin-bottom: 4px;
404+}
405+.plot-host svg {
406+ display: block;
407+ overflow: visible;
408+}
409+
410+.legend {
411+ display: flex;
412+ flex-wrap: wrap;
413+ gap: 6px 16px;
414+ margin: 2px 0 6px;
415+ font-size: 12px;
416+ color: var(--text-dim);
417+}
418+.legend-item {
419+ display: inline-flex;
420+ align-items: center;
421+ gap: 5px;
422+}
423+.legend-value {
424+ color: var(--text);
425+ font-variant-numeric: tabular-nums;
426+ font-size: 11.5px;
427+}
428+
429+.readout {
430+ display: flex;
431+ flex-wrap: wrap;
432+ gap: 6px 18px;
433+ font-size: 12px;
434+ font-variant-numeric: tabular-nums;
435+ color: var(--text-dim);
436+ padding: 4px 0 2px;
437+ min-height: 22px;
438+}
439+.readout b {
440+ color: inherit;
441+ font-weight: 600;
442+}
443+
444+.tooltip {
445+ position: absolute;
446+ top: 8px;
447+ right: 8px;
448+ background: #1a2029;
449+ border: 1px solid var(--border);
450+ border-radius: 6px;
451+ padding: 4px 9px;
452+ font-size: 12px;
453+ font-variant-numeric: tabular-nums;
454+ pointer-events: none;
455+}
456+
457+.verdict {
458+ display: flex;
459+ gap: 10px;
460+ align-items: flex-start;
461+ border-radius: 8px;
462+ padding: 9px 12px;
463+ margin-bottom: 12px;
464+ font-size: 13px;
465+ line-height: 1.5;
466+ border: 1px solid var(--border);
467+}
468+.verdict.good {
469+ border-color: #1b4b1b;
470+ background: #10190f;
471+}
472+.verdict.bad {
473+ border-color: #5b2626;
474+ background: #1e1416;
475+}
476+.verdict-mark {
477+ font-size: 15px;
478+ line-height: 1.3;
479+}
480+.verdict.good .verdict-mark {
481+ color: var(--good);
482+}
483+.verdict.bad .verdict-mark {
484+ color: var(--bad);
485+}
486+
487+.weights-int {
488+ font-size: 12.5px;
489+ color: var(--text-dim);
490+ margin-top: 4px;
491+}
492+.weights-int-values {
493+ display: flex;
494+ flex-wrap: wrap;
495+ gap: 4px;
496+ margin-top: 6px;
497+}
498+.weights-int-values span {
499+ min-width: 26px;
500+ text-align: center;
501+ padding: 2px 5px;
502+ background: #1a2029;
503+ border: 1px solid var(--border);
504+ border-radius: 4px;
505+ font-variant-numeric: tabular-nums;
506+ color: var(--text);
507+}
508+
509+.stale-note {
510+ font-size: 12px;
511+ color: #c98500;
512+ margin-bottom: 6px;
513+}
514+
515+/* ── convergence table ────────────────────────────────── */
516+.table-wrap {
517+ overflow-x: auto;
518+ margin: 6px 0 10px;
519+}
520+.conv-table {
521+ border-collapse: collapse;
522+ font-size: 12px;
523+ font-variant-numeric: tabular-nums;
524+ white-space: nowrap;
525+}
526+.conv-table th,
527+.conv-table td {
528+ padding: 3px 10px;
529+ text-align: right;
530+ border-bottom: 1px solid var(--border);
531+}
532+.conv-table thead th {
533+ color: var(--muted);
534+ font-weight: 600;
535+ text-align: center;
536+}
537+.conv-table tr.sub-head th {
538+ font-weight: 400;
539+ font-size: 11px;
540+}
541+.conv-table td.n-cell {
542+ color: var(--text);
543+ font-weight: 600;
544+ text-align: left;
545+}
546+.conv-table td.ord {
547+ color: var(--muted);
548+}
549+
550+/* ── narrow screens: stack the two columns ────────────── */
551+@media (max-width: 900px) {
552+ body {
553+ overflow: auto;
554+ }
555+ .body {
556+ flex-direction: column;
557+ }
558+ .left {
559+ width: 100%;
560+ max-width: none;
561+ border-right: none;
562+ border-bottom: 1px solid var(--border);
563+ }
564+ .editor-host {
565+ height: 320px;
566+ flex: none;
567+ }
568+ .panel-scroll {
569+ overflow: visible;
570+ }
571+}
src/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import { StrictMode } from 'react'
2+import { createRoot } from 'react-dom/client'
3+import App from './App.tsx'
4+import './index.css'
5+
6+createRoot(document.getElementById('root')!).render(
7+ <StrictMode>
8+ <App />
9+ </StrictMode>,
10+)
src/matlab/driver.madded+297−0View file
@@ -0,0 +1,297 @@
1+% ---------------------------------------------------------------------------
2+% Driver (owned by the app, not editable in the browser).
3+%
4+% It reads params.json, calls into the method script, and writes out.json.
5+% The editable method script is appended to the end of this file before the
6+% run, so the functions it defines -- bary_weights, bary_eval, and optionally
7+% local_blend -- are local functions of this script and are visible here.
8+% Every driver-owned helper is prefixed drv_ so that it cannot collide with
9+% anything the method script defines.
10+% ---------------------------------------------------------------------------
11+
12+drv_p = jsondecode(fileread('params.json'));
13+switch drv_p.mode
14+ case 'explore'
15+ drv_out = drv_explore(drv_p);
16+ case 'converge'
17+ drv_out = drv_converge(drv_p);
18+ otherwise
19+ error('driver: unknown mode ''%s''', drv_p.mode);
20+end
21+drv_fid = fopen('out.json', 'w');
22+fprintf(drv_fid, '%s', jsonencode(drv_out));
23+fclose(drv_fid);
24+
25+% ── one set of nodes, one degree d: everything the first three tabs draw ────
26+function out = drv_explore(p)
27+x = nodes_of(p.nodes, p.a, p.b, p.n, p.seed);
28+n = numel(x) - 1;
29+d = min(max(round(p.d), 0), n);
30+y = testfun(p.f, x, p.fexpr);
31+t = linspace(p.a, p.b, p.ngrid);
32+ft = testfun(p.f, t, p.fexpr);
33+
34+w = bary_weights(x, d);
35+w = reshape(w, 1, []);
36+r = bary_eval(x, y, w, t);
37+r = reshape(r, 1, []);
38+
39+out = struct();
40+out.n = n;
41+out.d = d;
42+out.x = x;
43+out.y = y;
44+out.t = t;
45+out.ft = ft;
46+out.r = r;
47+out.err = r - ft;
48+out.maxerr = drv_maxabs(r - ft);
49+
50+% Barycentric weights, and the integer form they take on a uniform mesh.
51+% Section 4 lists these: 1,1,...,1,1 for d = 0, then 1,2,2,...,2,2,1 for
52+% d = 1, 1,3,4,...,4,3,1 for d = 2, and so on.
53+out.w = w;
54+aw = abs(w);
55+pos = aw(aw > 0);
56+if isempty(pos)
57+ base = 1;
58+else
59+ base = min(pos);
60+end
61+out.wscaled = aw / base;
62+out.wsign = sign(w);
63+out.wIsInteger = all(abs(out.wscaled - round(out.wscaled)) < 1e-7);
64+out.wAlternates = all(w(1:end - 1) .* w(2:end) < 0);
65+
66+% the degree-n polynomial interpolant, in barycentric form (equation 2)
67+if p.want.poly
68+ out.rpoly = reshape(bary_eval(x, y, drv_lagrange_weights(x), t), 1, []);
69+end
70+
71+% the clamped C^2 cubic spline of Tables 3 and 4
72+if p.want.spline
73+ [~, dya] = testfun(p.f, p.a, p.fexpr);
74+ [~, dyb] = testfun(p.f, p.b, p.fexpr);
75+ out.rspline = reshape(cubic_spline(x, y, dya, dyb, t), 1, []);
76+end
77+
78+% the blend of local polynomials, equations (4) and (5)
79+if p.want.blend
80+ out.hasBlend = false;
81+ try
82+ [P, L] = local_blend(x, y, d, t);
83+ out.P = drv_rows(P);
84+ out.L = drv_rows(L);
85+ out.wlo = x(1:(n - d + 1));
86+ out.whi = x((d + 1):(n + 1));
87+ out.hasBlend = true;
88+ catch blenderr
89+ out.blendError = blenderr.message;
90+ end
91+end
92+
93+if p.want.poles
94+ out.poles = drv_poles(x, y, w, p);
95+end
96+end
97+
98+% ── denominator of r, its sign on the real line, and its roots ─────────────
99+function s = drv_poles(x, y, w, p)
100+n = numel(x) - 1;
101+mg = 0.35 * (p.b - p.a);
102+tw = linspace(p.a - mg, p.b + mg, p.ngridwide);
103+
104+% Work in a variable rescaled to [-1, 1]. This only multiplies the
105+% denominator by a positive constant, so signs and roots are untouched, but it
106+% keeps the products of n factors from over- or underflowing.
107+c0 = (p.a + p.b) / 2;
108+sc = (p.b - p.a) / 2;
109+xs = (x - c0) / sc;
110+ts = (tw - c0) / sc;
111+
112+q = drv_denom_eval(xs, w, ts);
113+
114+s = struct();
115+s.t = tw;
116+% A plain plot of q is useless: it spans many orders of magnitude. The signed
117+% n-th root keeps every sign and every zero exactly where it was and brings the
118+% magnitudes into a range that can be drawn.
119+qmax = max(abs(q));
120+if qmax == 0
121+ qmax = 1;
122+end
123+s.u = sign(q) .* (abs(q) / qmax).^(1 / max(1, n));
124+
125+% Real poles: sign changes of the denominator, located by linear interpolation
126+% of the crossing. This works for any n, which the root-finding below does not.
127+sg = sign(q);
128+idx = find(sg(1:end - 1) .* sg(2:end) < 0);
129+rp = zeros(1, numel(idx));
130+for m = 1:numel(idx)
131+ j = idx(m);
132+ rp(m) = tw(j) + (tw(j + 1) - tw(j)) * abs(q(j)) / (abs(q(j)) + abs(q(j + 1)));
133+end
134+s.realPoles = rp;
135+
136+% All roots in the complex plane. Theorem 1 says none of them are real, and
137+% that is the picture: every root sits off the real axis. We only do this for
138+% moderate n, since finding roots from monomial coefficients is not reliable
139+% for large degree.
140+s.rootsRe = [];
141+s.rootsIm = [];
142+s.rootsShown = false;
143+if n <= p.rootsMaxN
144+ [cc, floorc] = drv_denom_coeffs(xs, w);
145+ % The top d coefficients of q vanish identically -- that is exactly the
146+ % statement that r reproduces polynomials of degree d, i.e. that the first
147+ % d moments of the weights are zero -- but they come out of the sum as
148+ % roundoff rather than as zero. Left in place they contribute spurious
149+ % roots, some of them real, which would be a lie in a picture whose whole
150+ % point is that there are no real roots. So drop every leading
151+ % coefficient that is below the rounding-error floor of its own sum.
152+ k0 = find(abs(cc) > floorc, 1);
153+ if ~isempty(k0) && numel(cc) - k0 >= 1
154+ z = reshape(roots(cc(k0:end)), 1, []);
155+ s.rootsRe = real(z) * sc + c0;
156+ s.rootsIm = imag(z) * sc;
157+ s.rootsShown = true;
158+ end
159+end
160+
161+% the classical alternative, with the poles it puts in the interval
162+if p.want.classical
163+ [rc, cp] = classical_rational(x, y, tw);
164+ s.classical = reshape(rc, 1, []);
165+ s.classicalPoles = cp;
166+end
167+end
168+
169+% q(t) = sum_k w_k prod_{j ~= k} (t - x_j), the denominator of r written as a
170+% polynomial. Its zeros are exactly the poles of r.
171+function q = drv_denom_eval(xs, w, ts)
172+D = ts(:) - xs(:).';
173+n1 = numel(xs);
174+q = zeros(numel(ts), 1);
175+for k = 1:n1
176+ pr = ones(numel(ts), 1);
177+ for j = 1:n1
178+ if j ~= k
179+ pr = pr .* D(:, j);
180+ end
181+ end
182+ q = q + w(k) * pr;
183+end
184+q = q.';
185+end
186+
187+function [c, floorc] = drv_denom_coeffs(xs, w)
188+n1 = numel(xs);
189+n = n1 - 1;
190+c = zeros(1, n1);
191+for k = 1:n1
192+ ck = 1;
193+ for j = 1:n1
194+ if j ~= k
195+ ck = conv(ck, [1, -xs(j)]);
196+ end
197+ end
198+ c = c + w(k) * ck;
199+end
200+% Rounding-error floor, coefficient by coefficient. The nodes have been
201+% rescaled so that every |x_j| <= 1, so the coefficient of t^(n-m) in each of
202+% the n+1 products is at most binomial(n, m); summing them accumulates at most
203+% n+1 roundings of terms of size |w_k| times that bound.
204+floorc = zeros(1, n1);
205+sw = sum(abs(w));
206+for m = 0:n
207+ floorc(m + 1) = 8 * eps * n1 * sw * nchoosek(n, m);
208+end
209+end
210+
211+% the weights of equation (2), which put the degree-n polynomial interpolant
212+% into barycentric form
213+function w = drv_lagrange_weights(x)
214+n1 = numel(x);
215+w = zeros(1, n1);
216+for k = 1:n1
217+ pr = 1;
218+ for j = 1:n1
219+ if j ~= k
220+ pr = pr / (x(k) - x(j));
221+ end
222+ end
223+ w(k) = pr;
224+end
225+end
226+
227+% ── error against n, for a range of d: Tables 1 to 4 ───────────────────────
228+function out = drv_converge(p)
229+ns = reshape(p.ns, 1, []);
230+ds = reshape(p.ds, 1, []);
231+t = linspace(p.a, p.b, p.ngrid);
232+ft = testfun(p.f, t, p.fexpr);
233+[~, dya] = testfun(p.f, p.a, p.fexpr);
234+[~, dyb] = testfun(p.f, p.b, p.fexpr);
235+
236+E = zeros(numel(ds), numel(ns));
237+esp = zeros(1, numel(ns));
238+epo = zeros(1, numel(ns));
239+
240+for ib = 1:numel(ns)
241+ n = ns(ib);
242+ x = nodes_of(p.nodes, p.a, p.b, n, p.seed);
243+ y = testfun(p.f, x, p.fexpr);
244+ for ia = 1:numel(ds)
245+ d = min(ds(ia), n);
246+ w = reshape(bary_weights(x, d), 1, []);
247+ E(ia, ib) = drv_maxabs(reshape(bary_eval(x, y, w, t), 1, []) - ft);
248+ end
249+ if p.want.spline
250+ esp(ib) = drv_maxabs(reshape(cubic_spline(x, y, dya, dyb, t), 1, []) - ft);
251+ end
252+ if p.want.poly
253+ wl = drv_lagrange_weights(x);
254+ epo(ib) = drv_maxabs(reshape(bary_eval(x, y, wl, t), 1, []) - ft);
255+ end
256+end
257+
258+out = struct();
259+out.ns = ns;
260+out.ds = ds;
261+out.E = drv_rows(E);
262+out.orders = drv_rows(drv_orders(ns, E));
263+if p.want.spline
264+ out.splineErr = esp;
265+ out.splineOrders = drv_orders(ns, esp);
266+end
267+if p.want.poly
268+ out.polyErr = epo;
269+end
270+end
271+
272+function o = drv_orders(ns, E)
273+o = zeros(size(E));
274+o(:, 1) = NaN;
275+for j = 2:numel(ns)
276+ o(:, j) = log(E(:, j - 1) ./ E(:, j)) / log(ns(j) / ns(j - 1));
277+end
278+end
279+
280+% jsonencode flattens a matrix with a single row, so hand matrices over as a
281+% cell array of rows to keep the shape on the JavaScript side predictable.
282+function c = drv_rows(A)
283+c = cell(1, size(A, 1));
284+for i = 1:size(A, 1)
285+ c{i} = A(i, :);
286+end
287+end
288+
289+function m = drv_maxabs(v)
290+v = abs(v(:));
291+v = v(isfinite(v));
292+if isempty(v)
293+ m = Inf;
294+else
295+ m = max(v);
296+end
297+end
src/matlab/lib/classical_rational.madded+58−0View file
@@ -0,0 +1,58 @@
1+function [r, poles] = classical_rational(x, y, t)
2+%CLASSICAL_RATIONAL "Classical" rational interpolation p_M / q_N, M + N = n.
3+%
4+% This is the construction the paper's introduction describes and rejects:
5+% fit the values f(x_i) with a quotient of polynomials of degrees M and N
6+% with M + N = n, taking M = N = n/2 when n is even. It is the method with
7+% "no control over the occurrence of poles in the interval of interpolation",
8+% and this routine returns those poles so that they can be drawn.
9+%
10+% The interpolation conditions p(x_i) - y_i q(x_i) = 0 are linear in the
11+% coefficients, so the coefficient vector is a null vector of an
12+% (n+1) x (n+2) matrix, which we take from the last right singular vector.
13+% Everything is done in a variable rescaled to [-1, 1], since the monomial
14+% basis on the original interval is badly conditioned.
15+%
16+% Note that solving the linearised conditions does not guarantee that the
17+% quotient actually interpolates: a common root of p and q at some x_i (an
18+% "unattainable point") is possible. That is a further wrinkle of the
19+% classical method, not a bug here.
20+
21+x = x(:).';
22+y = y(:).';
23+n = numel(x) - 1;
24+M = ceil(n / 2);
25+N = n - M;
26+
27+a = min(x);
28+b = max(x);
29+c0 = (a + b) / 2;
30+sc = (b - a) / 2;
31+xs = (x - c0) / sc;
32+
33+V = zeros(n + 1, M + N + 2);
34+for i = 1:(n + 1)
35+ V(i, 1:(M + 1)) = xs(i).^(M:-1:0);
36+ V(i, (M + 2):end) = -y(i) * xs(i).^(N:-1:0);
37+end
38+
39+[~, ~, W] = svd(V);
40+c = W(:, end).';
41+pc = c(1:(M + 1));
42+qc = c((M + 2):end);
43+
44+ts = (t - c0) / sc;
45+r = polyval(pc, ts) ./ polyval(qc, ts);
46+r = reshape(r, size(t));
47+
48+% real poles = real roots of q, mapped back to the original variable
49+poles = [];
50+tol = 1e-13 * max(abs(qc));
51+k0 = find(abs(qc) > tol, 1);
52+if ~isempty(k0) && numel(qc) - k0 >= 1
53+ z = roots(qc(k0:end));
54+ z = z(:).';
55+ keep = abs(imag(z)) < 1e-7 * max(1, max(abs(z)));
56+ poles = sort(real(z(keep)) * sc + c0);
57+end
58+end
src/matlab/lib/cubic_spline.madded+81−0View file
@@ -0,0 +1,81 @@
1+function s = cubic_spline(x, y, dya, dyb, t)
2+%CUBIC_SPLINE Clamped C^2 cubic spline interpolant, evaluated at t.
3+%
4+% This is the competitor in Tables 3 and 4 of the paper: a C^2 cubic spline
5+% with clamped end conditions, i.e. with the first derivative at the two
6+% end-points set to the corresponding derivative of f. Its error is O(h^4)
7+% for f in C^4, the same order as the rational interpolant with d = 3.
8+%
9+% The moments M_i = s''(x_i) solve a tridiagonal system, which we solve with
10+% the Thomas algorithm so that the cost stays O(n) even for the largest n in
11+% the convergence study.
12+
13+x = x(:).';
14+y = y(:).';
15+n = numel(x) - 1;
16+h = diff(x);
17+
18+lo = zeros(1, n + 1); % sub-diagonal
19+di = zeros(1, n + 1); % diagonal
20+up = zeros(1, n + 1); % super-diagonal
21+rh = zeros(1, n + 1); % right-hand side
22+
23+% clamped end conditions
24+di(1) = 2;
25+up(1) = 1;
26+rh(1) = 6 / h(1) * ((y(2) - y(1)) / h(1) - dya);
27+
28+for i = 2:n
29+ hl = h(i - 1);
30+ hr = h(i);
31+ lo(i) = hl / (hl + hr);
32+ di(i) = 2;
33+ up(i) = hr / (hl + hr);
34+ rh(i) = 6 * ((y(i + 1) - y(i)) / hr - (y(i) - y(i - 1)) / hl) / (hl + hr);
35+end
36+
37+lo(n + 1) = 1;
38+di(n + 1) = 2;
39+rh(n + 1) = 6 / h(n) * (dyb - (y(n + 1) - y(n)) / h(n));
40+
41+% Thomas algorithm
42+cp = zeros(1, n + 1);
43+dp = zeros(1, n + 1);
44+cp(1) = up(1) / di(1);
45+dp(1) = rh(1) / di(1);
46+for i = 2:n + 1
47+ den = di(i) - lo(i) * cp(i - 1);
48+ cp(i) = up(i) / den;
49+ dp(i) = (rh(i) - lo(i) * dp(i - 1)) / den;
50+end
51+M = zeros(1, n + 1);
52+M(n + 1) = dp(n + 1);
53+for i = n:-1:1
54+ M(i) = dp(i) - cp(i) * M(i + 1);
55+end
56+
57+% evaluate: on [x_i, x_{i+1}] the spline is the usual cubic in the moments
58+sz = size(t);
59+t = t(:).';
60+s = zeros(1, numel(t));
61+for i = 1:n
62+ if i == 1
63+ m = t < x(2); % also catches t < x(1)
64+ elseif i == n
65+ m = t >= x(n); % also catches t > x(n+1)
66+ else
67+ m = (t >= x(i)) & (t < x(i + 1));
68+ end
69+ if ~any(m)
70+ continue
71+ end
72+ tt = t(m);
73+ hi = h(i);
74+ ra = x(i + 1) - tt;
75+ rb = tt - x(i);
76+ s(m) = M(i) * ra.^3 / (6 * hi) + M(i + 1) * rb.^3 / (6 * hi) ...
77+ + (y(i) - M(i) * hi^2 / 6) .* ra / hi ...
78+ + (y(i + 1) - M(i + 1) * hi^2 / 6) .* rb / hi;
79+end
80+s = reshape(s, sz);
81+end
src/matlab/lib/evalexpr.madded+8−0View file
@@ -0,0 +1,8 @@
1+function y = evalexpr(expr, x)
2+%EVALEXPR Evaluate a user-typed expression in the variable x.
3+
4+y = eval(expr);
5+if numel(y) == 1 && numel(x) ~= 1
6+ y = y * ones(size(x));
7+end
8+end
src/matlab/lib/nodes_of.madded+34−0View file
@@ -0,0 +1,34 @@
1+function x = nodes_of(kind, a, b, n, seed)
2+%NODES_OF Build n+1 interpolation nodes a = x_0 < x_1 < ... < x_n = b.
3+%
4+% Floater and Hormann's Theorem 2 gives the rate O(h^(d+1)) for d >= 1
5+% regardless of how the nodes are distributed, so it is worth being able to
6+% distribute them badly. The 'paired' family below does exactly that: it
7+% pulls every second node close to its left neighbour, which drives the local
8+% mesh ratio beta of Theorem 3 up and makes the d = 0 case (Berrut's
9+% interpolant) misbehave while d >= 1 carries on unaffected.
10+
11+k = 0:n;
12+switch kind
13+ case 'uniform'
14+ x = a + (b - a) * k / n;
15+ case 'chebyshev'
16+ % Chebyshev-Gauss-Lobatto points, clustered at both ends
17+ x = (a + b) / 2 - (b - a) / 2 * cos(pi * k / n);
18+ case 'random'
19+ rng(seed);
20+ u = sort(rand(1, max(0, n - 1)));
21+ x = [a, a + (b - a) * u, b];
22+ case 'paired'
23+ x = a + (b - a) * k / n;
24+ h = (b - a) / n;
25+ % indices 2,4,... are the nodes x_1, x_3, ... (1-based vs 0-based)
26+ x(2:2:end) = x(2:2:end) - 0.9 * h;
27+ case 'graded'
28+ % quadratically graded, clustered at the left end
29+ x = a + (b - a) * (k / n).^2;
30+ otherwise
31+ error('nodes_of: unknown node distribution ''%s''', kind);
32+end
33+x = x(:).';
34+end
src/matlab/lib/testfun.madded+27−0View file
@@ -0,0 +1,27 @@
1+function [y, dy] = testfun(name, x, expr)
2+%TESTFUN The functions interpolated in Section 5 of the paper, plus a custom one.
3+%
4+% [y, dy] = testfun(name, x, expr) returns f(x) and f'(x). The derivative is
5+% only used for the clamped end conditions of the C^2 cubic spline that the
6+% paper compares against in Tables 3 and 4.
7+
8+switch name
9+ case 'runge'
10+ y = 1 ./ (1 + x.^2);
11+ dy = -2 * x ./ (1 + x.^2).^2;
12+ case 'sine'
13+ y = sin(x);
14+ dy = cos(x);
15+ case 'abs'
16+ y = abs(x);
17+ dy = sign(x);
18+ case 'custom'
19+ y = evalexpr(expr, x);
20+ h = 1e-6 * max(1, max(abs(x(:))));
21+ dy = (evalexpr(expr, x + h) - evalexpr(expr, x - h)) / (2 * h);
22+ otherwise
23+ error('testfun: unknown function ''%s''', name);
24+end
25+y = reshape(y, size(x));
26+dy = reshape(dy, size(x));
27+end
src/methods/berrut.madded+56−0View file
@@ -0,0 +1,56 @@
1+% ---------------------------------------------------------------------------
2+% Berrut's interpolant: the d = 0 member of the family
3+% J.-P. Berrut, Comput. Math. Appl. 15 (1988) 1-16
4+%
5+% Equation (3) of the paper. The weights are simply the alternating signs,
6+%
7+% w_k = (-1)^k,
8+%
9+% which is what equation (18) reduces to when d = 0 (up to a common positive
10+% factor, which does not change r). Berrut showed this has no real poles;
11+% Floater and Hormann's Theorem 3 gives it the rate O(h), but only under a
12+% bound on the local mesh ratio beta. Set the nodes to "paired" and watch what
13+% happens: the interpolant develops kinks and the error stops falling, while
14+% the d >= 1 members of the family are untroubled.
15+%
16+% The d slider is ignored by this script.
17+% ---------------------------------------------------------------------------
18+
19+function w = bary_weights(x, d)
20+n = numel(x) - 1;
21+w = (-1).^(0:n);
22+end
23+
24+function r = bary_eval(x, y, w, t)
25+sz = size(t);
26+D = t(:) - x(:).';
27+Q = w(:).' ./ D;
28+r = (Q * y(:)) ./ sum(Q, 2);
29+hit = find(any(D == 0, 2));
30+for m = 1:numel(hit)
31+ k = find(D(hit(m), :) == 0, 1);
32+ r(hit(m)) = y(k);
33+end
34+r = reshape(r, sz);
35+end
36+
37+function [P, L] = local_blend(x, y, d, t)
38+% With d = 0 the "local polynomials" are the constants p_i = f_i, and the
39+% blending functions are mu_i(t) = prod_{j<i} (t - x_j) * prod_{k>i} (x_k - t)
40+% normalised to sum to 1.
41+n = numel(x) - 1;
42+t = reshape(t, 1, []);
43+P = repmat(y(:), 1, numel(t));
44+Mu = zeros(n + 1, numel(t));
45+for i = 0:n
46+ pr = ones(1, numel(t));
47+ for j = 0:(i - 1)
48+ pr = pr .* (t - x(j + 1));
49+ end
50+ for k = (i + 1):n
51+ pr = pr .* (x(k + 1) - t);
52+ end
53+ Mu(i + 1, :) = pr;
54+end
55+L = Mu ./ sum(Mu, 1);
56+end
src/methods/equal.madded+40−0View file
@@ -0,0 +1,40 @@
1+% ---------------------------------------------------------------------------
2+% All weights equal: what the alternating signs are for
3+%
4+% w_k = 1 for every k.
5+%
6+% This is a perfectly good barycentric formula and it does interpolate: at each
7+% node the k-th term dominates and r(x_k) = f_k. Between two consecutive nodes,
8+% though, the denominator
9+%
10+% sum_k 1 / (t - x_k)
11+%
12+% runs from +infinity at the left node down to -infinity at the right one, so it
13+% crosses zero somewhere in every interval. Each of those crossings is a pole.
14+% Open the Poles tab: n real poles, one per interval, and the roots that
15+% Theorem 1 keeps off the real axis are sitting right on it.
16+%
17+% Schneider and Werner proved that the weights of a pole-free barycentric
18+% rational interpolant must alternate in sign, and the paper checks that the
19+% weights of equation (18) do. This script is the other side of that statement.
20+% Put a (-1)^k back in and the poles disappear -- that is Berrut's interpolant.
21+%
22+% The d slider is ignored by this script.
23+% ---------------------------------------------------------------------------
24+
25+function w = bary_weights(x, d)
26+w = ones(1, numel(x));
27+end
28+
29+function r = bary_eval(x, y, w, t)
30+sz = size(t);
31+D = t(:) - x(:).';
32+Q = w(:).' ./ D;
33+r = (Q * y(:)) ./ sum(Q, 2);
34+hit = find(any(D == 0, 2));
35+for m = 1:numel(hit)
36+ k = find(D(hit(m), :) == 0, 1);
37+ r(hit(m)) = y(k);
38+end
39+r = reshape(r, sz);
40+end
src/methods/fh.madded+106−0View file
@@ -0,0 +1,106 @@
1+% ---------------------------------------------------------------------------
2+% Floater-Hormann barycentric rational interpolation
3+% M. S. Floater and K. Hormann, Numer. Math. 107 (2007) 315-331
4+%
5+% This script is the method. Everything the four tabs draw comes out of the
6+% functions below, so editing them changes the pictures.
7+%
8+% w = bary_weights(x, d) weights of equation (18) [required]
9+% r = bary_eval(x, y, w, t) barycentric form, equation (1)[required]
10+% [P, L] = local_blend(x, y, d, t) the blend of (4) and (5) [optional]
11+% ---------------------------------------------------------------------------
12+
13+function w = bary_weights(x, d)
14+% Equation (18). With J_k = { i : k-d <= i <= k } intersected with {0,...,n-d},
15+%
16+% w_k = sum_{i in J_k} (-1)^i prod_{j=i, j~=k}^{i+d} 1 / (x_k - x_j).
17+%
18+% Every k lies in the window of at most d+1 of the local polynomials, so this
19+% costs O(n d^2) however the nodes are placed.
20+n = numel(x) - 1;
21+d = min(max(d, 0), n);
22+w = zeros(1, n + 1);
23+for k = 0:n
24+ s = 0;
25+ for i = max(0, k - d):min(k, n - d)
26+ p = 1;
27+ for j = i:(i + d)
28+ if j ~= k
29+ p = p / (x(k + 1) - x(j + 1));
30+ end
31+ end
32+ s = s + (-1)^i * p;
33+ end
34+ w(k + 1) = s;
35+end
36+end
37+
38+function r = bary_eval(x, y, w, t)
39+% Equation (1):
40+%
41+% r(t) = sum_k w_k f_k / (t - x_k) / sum_k w_k / (t - x_k).
42+%
43+% Berrut and Trefethen's advice is followed at the nodes themselves: if t is
44+% exactly some x_k, return f_k rather than dividing by zero.
45+sz = size(t);
46+D = t(:) - x(:).';
47+Q = w(:).' ./ D;
48+r = (Q * y(:)) ./ sum(Q, 2);
49+hit = find(any(D == 0, 2));
50+for m = 1:numel(hit)
51+ k = find(D(hit(m), :) == 0, 1);
52+ r(hit(m)) = y(k);
53+end
54+r = reshape(r, sz);
55+end
56+
57+function [P, L] = local_blend(x, y, d, t)
58+% Equations (4) and (5), the construction the barycentric form above is a
59+% rewriting of. P(i+1,:) is the polynomial p_i of degree at most d through the
60+% d+1 points x_i, ..., x_{i+d}, and L(i+1,:) is the blending function lambda_i
61+% normalised so that the columns of L sum to 1. Then r = sum_i L_i p_i.
62+%
63+% We build the blending functions from mu_i of equation (9),
64+%
65+% mu_i(t) = prod_{j<i} (t - x_j) * prod_{k>i+d} (x_k - t),
66+%
67+% rather than from lambda_i of equation (5) directly. The two differ by a
68+% factor that does not depend on i, so they normalise to the same thing, but
69+% mu_i is a polynomial and so has nothing to blow up at the nodes. Its sum is
70+% the s(x) that Theorem 1 shows is positive, which is why L is well defined
71+% everywhere.
72+n = numel(x) - 1;
73+d = min(max(d, 0), n);
74+t = reshape(t, 1, []);
75+m = n - d + 1;
76+
77+P = zeros(m, numel(t));
78+Mu = zeros(m, numel(t));
79+for i = 0:(n - d)
80+ idx = (i + 1):(i + d + 1);
81+ P(i + 1, :) = local_poly(x(idx), y(idx), t);
82+ pr = ones(1, numel(t));
83+ for j = 0:(i - 1)
84+ pr = pr .* (t - x(j + 1));
85+ end
86+ for k = (i + d + 1):n
87+ pr = pr .* (x(k + 1) - t);
88+ end
89+ Mu(i + 1, :) = pr;
90+end
91+L = Mu ./ sum(Mu, 1);
92+end
93+
94+function p = local_poly(xi, yi, t)
95+% Lagrange form of the polynomial of degree at most numel(xi)-1 through (xi, yi).
96+p = zeros(1, numel(t));
97+for k = 1:numel(xi)
98+ b = ones(1, numel(t));
99+ for j = 1:numel(xi)
100+ if j ~= k
101+ b = b .* (t - xi(j)) / (xi(k) - xi(j));
102+ end
103+ end
104+ p = p + yi(k) * b;
105+end
106+end
src/methods/index.tsadded+59−0View file
@@ -0,0 +1,59 @@
1+import fh from './fh.m?raw'
2+import berrut from './berrut.m?raw'
3+import uniformInteger from './uniform-integer.m?raw'
4+import lagrange from './lagrange.m?raw'
5+import equal from './equal.m?raw'
6+import random from './random.m?raw'
7+
8+export interface Method {
9+ id: string
10+ name: string
11+ /** one line, shown in the picker */
12+ blurb: string
13+ source: string
14+}
15+
16+export const METHODS: Method[] = [
17+ {
18+ id: 'fh',
19+ name: 'Floater-Hormann',
20+ blurb: 'The paper: blend n-d+1 local polynomials of degree d. No poles, order h^(d+1).',
21+ source: fh,
22+ },
23+ {
24+ id: 'berrut',
25+ name: 'Berrut (d = 0)',
26+ blurb: 'Weights (-1)^k. The d = 0 member, and the one that needs a bounded mesh ratio.',
27+ source: berrut,
28+ },
29+ {
30+ id: 'uniform-integer',
31+ name: 'Integer weights',
32+ blurb: "Section 4's closed form on a uniform mesh: 1, 4, 7, 8, ..., 8, 7, 4, 1 for d = 3.",
33+ source: uniformInteger,
34+ },
35+ {
36+ id: 'lagrange',
37+ name: 'Polynomial (d = n)',
38+ blurb: 'The Lagrange weights of equation (2). No poles, but Runge divergence.',
39+ source: lagrange,
40+ },
41+ {
42+ id: 'equal',
43+ name: 'Equal weights',
44+ blurb: 'w_k = 1. Drop the alternating signs and a pole appears in every interval.',
45+ source: equal,
46+ },
47+ {
48+ id: 'random',
49+ name: 'Random weights',
50+ blurb: 'The generic barycentric rational interpolant. Interpolates; has poles.',
51+ source: random,
52+ },
53+]
54+
55+export const DEFAULT_METHOD = METHODS[0]
56+
57+export function findMethod(id: string | null | undefined): Method | undefined {
58+ return METHODS.find((m) => m.id === id)
59+}
src/methods/lagrange.madded+47−0View file
@@ -0,0 +1,47 @@
1+% ---------------------------------------------------------------------------
2+% The polynomial interpolant, in barycentric form: the d = n member
3+%
4+% Equation (2), first written down by Taylor and by Dupuy:
5+%
6+% w_k = prod_{j ~= k} 1 / (x_k - x_j).
7+%
8+% These are what equation (18) gives when d = n, since then there is a single
9+% local polynomial and it is the interpolating polynomial p_n itself. The
10+% paper's remark that the weights of (2) "prevent poles" is visible on the
11+% Poles tab in a degenerate way: the denominator of r reduces to the constant 1,
12+% because the Lagrange basis functions sum to 1, so there is nothing to vanish
13+% and no roots at all.
14+%
15+% Having no poles is not the same as approximating well. Leave the function at
16+% Runge's 1/(1+x^2), the nodes uniform, and push n up: this is the divergence
17+% the paper opens with. Switch the nodes to Chebyshev and it behaves.
18+%
19+% The d slider is ignored by this script.
20+% ---------------------------------------------------------------------------
21+
22+function w = bary_weights(x, d)
23+n1 = numel(x);
24+w = zeros(1, n1);
25+for k = 1:n1
26+ p = 1;
27+ for j = 1:n1
28+ if j ~= k
29+ p = p / (x(k) - x(j));
30+ end
31+ end
32+ w(k) = p;
33+end
34+end
35+
36+function r = bary_eval(x, y, w, t)
37+sz = size(t);
38+D = t(:) - x(:).';
39+Q = w(:).' ./ D;
40+r = (Q * y(:)) ./ sum(Q, 2);
41+hit = find(any(D == 0, 2));
42+for m = 1:numel(hit)
43+ k = find(D(hit(m), :) == 0, 1);
44+ r(hit(m)) = y(k);
45+end
46+r = reshape(r, sz);
47+end
src/methods/random.madded+35−0View file
@@ -0,0 +1,35 @@
1+% ---------------------------------------------------------------------------
2+% Random weights: the generic barycentric rational interpolant
3+%
4+% Berrut and Mittelmann's observation, quoted in the paper's introduction, is
5+% that *every* rational interpolant whose numerator and denominator have degree
6+% at most n can be written in the barycentric form (1) for some real weights
7+% w_0, ..., w_n. So the whole difficulty of the subject is choosing them. This
8+% script chooses them at random.
9+%
10+% The result still interpolates -- that is free -- but the denominator now
11+% changes sign wherever it pleases, and the Poles tab finds real poles inside
12+% the interval. Reroll with the d slider, which is used here only to reseed.
13+%
14+% This is the situation the paper's construction is a way out of: weights that
15+% are known in advance to give a pole-free interpolant, and one whose
16+% approximation order can be raised at will.
17+% ---------------------------------------------------------------------------
18+
19+function w = bary_weights(x, d)
20+rng(1 + d);
21+w = randn(1, numel(x));
22+end
23+
24+function r = bary_eval(x, y, w, t)
25+sz = size(t);
26+D = t(:) - x(:).';
27+Q = w(:).' ./ D;
28+r = (Q * y(:)) ./ sum(Q, 2);
29+hit = find(any(D == 0, 2));
30+for m = 1:numel(hit)
31+ k = find(D(hit(m), :) == 0, 1);
32+ r(hit(m)) = y(k);
33+end
34+r = reshape(r, sz);
35+end
src/methods/uniform-integer.madded+55−0View file
@@ -0,0 +1,55 @@
1+% ---------------------------------------------------------------------------
2+% The integer weights of Section 4, for equally spaced nodes
3+%
4+% When the nodes are uniform with spacing h, equation (18) collapses to
5+%
6+% w_k = (-1)^(k-d) / h^d * sum_{i in J_k} 1 / ((k-i)! (i+d-k)!),
7+%
8+% and since a common positive factor does not change r we may multiply by
9+% d! h^d and read off integers:
10+%
11+% w_k = (-1)^(k-d) * sum_{i in J_k} binomial(d, k-i).
12+%
13+% Writing delta_k = |w_k|, the first few rows are the ones tabulated in the
14+% paper:
15+%
16+% d = 0: 1, 1, ..., 1, 1
17+% d = 1: 1, 2, 2, ..., 2, 2, 1
18+% d = 2: 1, 3, 4, ..., 4, 3, 1
19+% d = 3: 1, 4, 7, 8, 8, ..., 8, 8, 7, 4, 1
20+% d = 4: 1, 5, 11, 15, 16, 16, ..., 16, 16, 15, 11, 5, 1
21+%
22+% Almost every weight is the same; the only difference is at the two ends.
23+% Yet that small change is what raises the approximation order from O(h) to
24+% O(h^(d+1)). The "Blending & weights" tab checks these against the general
25+% formula, so switching between this script and the main one should leave the
26+% pictures identical -- as long as the nodes stay uniform. Choose any other
27+% node distribution and these weights are no longer the right ones, and the
28+% "Poles" tab may well find real poles.
29+% ---------------------------------------------------------------------------
30+
31+function w = bary_weights(x, d)
32+n = numel(x) - 1;
33+d = min(max(d, 0), n);
34+w = zeros(1, n + 1);
35+for k = 0:n
36+ s = 0;
37+ for i = max(0, k - d):min(k, n - d)
38+ s = s + nchoosek(d, k - i);
39+ end
40+ w(k + 1) = (-1)^(k - d) * s;
41+end
42+end
43+
44+function r = bary_eval(x, y, w, t)
45+sz = size(t);
46+D = t(:) - x(:).';
47+Q = w(:).' ./ D;
48+r = (Q * y(:)) ./ sum(Q, 2);
49+hit = find(any(D == 0, 2));
50+for m = 1:numel(hit)
51+ k = find(D(hit(m), :) == 0, 1);
52+ r(hit(m)) = y(k);
53+end
54+r = reshape(r, sz);
55+end
src/panels/BlendingPanel.tsxadded+307−0View file
@@ -0,0 +1,307 @@
1+import { useMemo, useState } from 'react'
2+import Plot from '../plot/Plot.tsx'
3+import Legend from '../plot/Legend.tsx'
4+import { diverging, ink, series } from '../plot/palette.ts'
5+import { extent, linePath, padDomain, type Frame } from '../plot/scales.ts'
6+import type { ExploreOut, Num } from '../engine/types.ts'
7+
8+interface Props {
9+ out: ExploreOut
10+}
11+
12+/** The samples of t lying in [lo, hi], as a slice of both arrays. */
13+function slice(t: number[], v: readonly Num[], lo: number, hi: number): [number[], Num[]] {
14+ const ts: number[] = []
15+ const vs: Num[] = []
16+ for (let i = 0; i < t.length; i++) {
17+ if (t[i] >= lo && t[i] <= hi) {
18+ ts.push(t[i])
19+ vs.push(v[i])
20+ }
21+ }
22+ return [ts, vs]
23+}
24+
25+export default function BlendingPanel({ out }: Props) {
26+ const [pinned, setPinned] = useState<number | null>(null)
27+ const [hoverX, setHoverX] = useState<number | null>(null)
28+
29+ const P = out.P
30+ const L = out.L
31+ const wlo = out.wlo
32+ const whi = out.whi
33+ const m = P?.length ?? 0
34+
35+ // which local polynomial the reader is following: the pinned one, or the one
36+ // whose window is nearest the pointer
37+ const hovered = useMemo(() => {
38+ if (!wlo || !whi || hoverX == null) return null
39+ let best = 0
40+ let bestD = Infinity
41+ for (let i = 0; i < wlo.length; i++) {
42+ const dd = Math.abs((wlo[i] + whi[i]) / 2 - hoverX)
43+ if (dd < bestD) {
44+ bestD = dd
45+ best = i
46+ }
47+ }
48+ return best
49+ }, [wlo, whi, hoverX])
50+ const sel = pinned ?? hovered ?? Math.floor(m / 2)
51+
52+ const xd: [number, number] = [out.t[0], out.t[out.t.length - 1]]
53+ const yd = padDomain(extent(out.ft, out.r), 0.12)
54+ const ld = useMemo(() => {
55+ if (!L) return [0, 1] as [number, number]
56+ const e = extent(...L)
57+ return [Math.max(-1.2, Math.min(-0.15, e[0] * 1.1)), Math.min(1.6, Math.max(1.05, e[1] * 1.05))] as [
58+ number,
59+ number,
60+ ]
61+ }, [L])
62+
63+ if (!out.hasBlend || !P || !L || !wlo || !whi) {
64+ return (
65+ <div className="panel">
66+ <p className="panel-lede">
67+ This script does not define <code>local_blend</code>, so there is nothing to draw here. Equations (4)
68+ and (5) are optional: a script only has to supply <code>bary_weights</code> and{' '}
69+ <code>bary_eval</code>. The other three tabs still work.
70+ </p>
71+ {out.blendError && <pre className="error-box">{out.blendError}</pre>}
72+ <WeightsSection out={out} />
73+ </div>
74+ )
75+ }
76+
77+ const shade = (f: Frame) => (
78+ <rect
79+ x={f.sx(wlo[sel])}
80+ width={Math.max(1, f.sx(whi[sel]) - f.sx(wlo[sel]))}
81+ y={0}
82+ height={f.ih}
83+ fill={ink.familyHi}
84+ opacity={0.1}
85+ />
86+ )
87+
88+ return (
89+ <div className="panel">
90+ <p className="panel-lede">
91+ Equation (4) reads r = Σ<sub>i</sub> λ<sub>i</sub> p<sub>i</sub> / Σ<sub>i</sub>{' '}
92+ λ<sub>i</sub>: slide a window of d+1 = {out.d + 1} nodes along the data, fit a polynomial of
93+ degree {out.d} in each position, and blend the {m} of them together. Hover to follow one; click to pin
94+ it.
95+ </p>
96+
97+ <div className="row-controls">
98+ <label className="slider-label">
99+ <span>
100+ window i = <b>{sel}</b> [x<sub>{sel}</sub>, x<sub>{sel + out.d}</sub>] = [
101+ {wlo[sel].toFixed(2)}, {whi[sel].toFixed(2)}]
102+ </span>
103+ <input
104+ type="range"
105+ min={0}
106+ max={m - 1}
107+ value={sel}
108+ onChange={(e) => setPinned(Number(e.target.value))}
109+ />
110+ </label>
111+ {pinned != null && (
112+ <button className="ghost" onClick={() => setPinned(null)}>
113+ unpin
114+ </button>
115+ )}
116+ </div>
117+
118+ <Legend
119+ items={[
120+ { label: 'f(x)', color: ink.reference, dash: '5 4' },
121+ { label: 'rational r(x)', color: series.r },
122+ { label: `the ${m} local polynomials p_i`, color: ink.family },
123+ { label: `p_${sel}, on its own window`, color: ink.familyHi },
124+ ]}
125+ />
126+
127+ <Plot
128+ height={300}
129+ xDomain={xd}
130+ yDomain={yd}
131+ xLabel="x"
132+ yLabel="f, r, local p_i"
133+ onHoverX={setHoverX}
134+ hoverX={hoverX}
135+ under={shade}
136+ >
137+ {(f) => (
138+ <>
139+ {P.map((row, i) => {
140+ // each p_i is only drawn a little beyond the d+1 points it
141+ // interpolates: a degree-d polynomial extrapolated across the
142+ // whole interval is a wall of noise
143+ const pad = Math.max((whi[i] - wlo[i]) * 0.35, (xd[1] - xd[0]) / (out.n * 2))
144+ const [ts, vs] = slice(out.t, row, wlo[i] - pad, whi[i] + pad)
145+ return i === sel ? null : (
146+ <path key={i} d={linePath(ts, vs, f)} fill="none" stroke={ink.family} strokeWidth={1.2} />
147+ )
148+ })}
149+ <path d={linePath(out.t, out.ft, f)} fill="none" stroke={ink.reference} strokeWidth={2} strokeDasharray="5 4" />
150+ <path d={linePath(out.t, out.r, f)} fill="none" stroke={series.r} strokeWidth={2.5} />
151+ {(() => {
152+ const pad = (whi[sel] - wlo[sel]) * 0.9 + (xd[1] - xd[0]) / (out.n * 2)
153+ const [ts, vs] = slice(out.t, P[sel], wlo[sel] - pad, whi[sel] + pad)
154+ return <path d={linePath(ts, vs, f)} fill="none" stroke={ink.familyHi} strokeWidth={2.5} />
155+ })()}
156+ {out.x.map((xv, i) => {
157+ const inWin = xv >= wlo[sel] - 1e-12 && xv <= whi[sel] + 1e-12
158+ return (
159+ <circle
160+ key={i}
161+ cx={f.sx(xv)}
162+ cy={f.sy(out.y[i])}
163+ r={inWin ? 4.2 : 3}
164+ fill={inWin ? ink.familyHi : ink.node}
165+ stroke="#151a21"
166+ strokeWidth={1.5}
167+ />
168+ )
169+ })}
170+ </>
171+ )}
172+ </Plot>
173+
174+ <h4 className="sub">Blending functions</h4>
175+ <p className="panel-note">
176+ The normalised λ<sub>i</sub>, which sum to 1 at every x. Each one is close to 1 across its own
177+ window and decays away from it, but with a tail that oscillates in sign and never quite reaches zero:
178+ these functions have no local support, which the paper names as the price of the construction. What
179+ they do have is that their denominator never vanishes, so they are infinitely smooth.
180+ </p>
181+ <Plot
182+ height={220}
183+ xDomain={xd}
184+ yDomain={ld}
185+ xLabel="x"
186+ yLabel="normalised lambda_i"
187+ onHoverX={setHoverX}
188+ hoverX={hoverX}
189+ under={(f) => (
190+ <>
191+ {shade(f)}
192+ <line x1={0} x2={f.iw} y1={f.sy(0)} y2={f.sy(0)} stroke={ink.axis} strokeWidth={1} />
193+ <line
194+ x1={0}
195+ x2={f.iw}
196+ y1={f.sy(1)}
197+ y2={f.sy(1)}
198+ stroke={ink.axis}
199+ strokeWidth={1}
200+ strokeDasharray="2 4"
201+ />
202+ </>
203+ )}
204+ >
205+ {(f) => (
206+ <>
207+ {L.map((row, i) =>
208+ i === sel ? null : (
209+ <path key={i} d={linePath(out.t, row, f)} fill="none" stroke={ink.family} strokeWidth={1.2} />
210+ ),
211+ )}
212+ <path d={linePath(out.t, L[sel], f)} fill="none" stroke={ink.familyHi} strokeWidth={2.5} />
213+ {out.x.map((xv, i) => (
214+ <line
215+ key={i}
216+ x1={f.sx(xv)}
217+ x2={f.sx(xv)}
218+ y1={f.ih}
219+ y2={f.ih - 5}
220+ stroke={ink.node}
221+ strokeWidth={1.5}
222+ />
223+ ))}
224+ </>
225+ )}
226+ </Plot>
227+
228+ <WeightsSection out={out} />
229+ </div>
230+ )
231+}
232+
233+function WeightsSection({ out }: { out: ExploreOut }) {
234+ const wmax = Math.max(...out.w.map(Math.abs), Number.MIN_VALUE)
235+ const norm = out.w.map((v) => v / wmax)
236+ const spread = Math.max(...out.wscaled) / Math.min(...out.wscaled.filter((v) => v > 0))
237+
238+ return (
239+ <>
240+ <h4 className="sub">Barycentric weights</h4>
241+ <p className="panel-note">
242+ The same interpolant, written in the form of equation (1) with the weights of equation (18). Schneider
243+ and Werner proved that a pole-free barycentric rational interpolant must have weights that alternate
244+ in sign; these {out.wAlternates ? 'do' : 'do not'}.
245+ </p>
246+ <Legend
247+ items={[
248+ { label: 'w_k > 0', color: diverging.pos },
249+ { label: 'w_k < 0', color: diverging.neg },
250+ ]}
251+ />
252+ <Plot
253+ height={170}
254+ xDomain={[out.x[0], out.x[out.x.length - 1]]}
255+ yDomain={[-1.15, 1.15]}
256+ xLabel="x_k"
257+ yLabel="w_k / max |w|"
258+ margin={{ l: 62 }}
259+ under={(f) => <line x1={0} x2={f.iw} y1={f.sy(0)} y2={f.sy(0)} stroke={ink.axis} strokeWidth={1} />}
260+ >
261+ {(f) =>
262+ norm.map((v, k) => (
263+ <g key={k}>
264+ <line
265+ x1={f.sx(out.x[k])}
266+ x2={f.sx(out.x[k])}
267+ y1={f.sy(0)}
268+ y2={f.sy(v)}
269+ stroke={v >= 0 ? diverging.pos : diverging.neg}
270+ strokeWidth={2}
271+ />
272+ <circle
273+ cx={f.sx(out.x[k])}
274+ cy={f.sy(v)}
275+ r={3}
276+ fill={v >= 0 ? diverging.pos : diverging.neg}
277+ stroke="#151a21"
278+ strokeWidth={1}
279+ />
280+ </g>
281+ ))
282+ }
283+ </Plot>
284+
285+ {out.wIsInteger ? (
286+ <div className="weights-int">
287+ <div className="weights-int-head">
288+ δ<sub>k</sub> = |w<sub>k</sub>| / min |w|, which Section 4 predicts are integers on a uniform
289+ mesh:
290+ </div>
291+ <div className="weights-int-values">
292+ {out.wscaled.map((v, k) => (
293+ <span key={k}>{Math.round(v)}</span>
294+ ))}
295+ </div>
296+ </div>
297+ ) : (
298+ <div className="weights-int">
299+ <div className="weights-int-head">
300+ The weights are not integer multiples of the smallest one, so this is not a uniform mesh. Their
301+ magnitudes span a factor of {spread < 1e5 ? spread.toFixed(0) : spread.toExponential(1)}.
302+ </div>
303+ </div>
304+ )}
305+ </>
306+ )
307+}
src/panels/Controls.tsxadded+150−0View file
@@ -0,0 +1,150 @@
1+import type { FuncName, NodeKind } from '../engine/types.ts'
2+
3+export interface Settings {
4+ f: FuncName
5+ fexpr: string
6+ a: number
7+ b: number
8+ n: number
9+ d: number
10+ nodes: NodeKind
11+ seed: number
12+}
13+
14+const FUNCS: { id: FuncName; label: string; hint: string }[] = [
15+ { id: 'runge', label: '1 / (1 + x²)', hint: "Runge's example, the paper's Figure 1" },
16+ { id: 'sine', label: 'sin x', hint: 'Figure 2, smooth everywhere' },
17+ { id: 'abs', label: '|x|', hint: 'Figure 3, a corner at 0' },
18+ { id: 'custom', label: 'custom', hint: 'any MATLAB expression in x' },
19+]
20+
21+const NODES: { id: NodeKind; label: string; hint: string }[] = [
22+ { id: 'uniform', label: 'uniform', hint: 'equally spaced, as in every figure of the paper' },
23+ { id: 'chebyshev', label: 'Chebyshev', hint: 'clustered at both ends' },
24+ { id: 'random', label: 'random', hint: 'sorted uniform draws, endpoints kept' },
25+ { id: 'paired', label: 'paired', hint: 'every second node pulled close to its neighbour' },
26+ { id: 'graded', label: 'graded', hint: 'quadratically clustered at the left end' },
27+]
28+
29+interface Props {
30+ value: Settings
31+ onChange: (patch: Partial<Settings>) => void
32+ /** the convergence tab sets n itself */
33+ showN: boolean
34+ busy: boolean
35+}
36+
37+export default function Controls({ value: v, onChange, showN, busy }: Props) {
38+ const dMax = Math.min(v.n, 12)
39+ return (
40+ <div className="controls">
41+ <div className="field">
42+ <span className="field-label">function f</span>
43+ <div className="chips">
44+ {FUNCS.map((fn) => (
45+ <button
46+ key={fn.id}
47+ title={fn.hint}
48+ className={`chip ${v.f === fn.id ? 'on' : ''}`}
49+ onClick={() => onChange({ f: fn.id })}
50+ >
51+ {fn.label}
52+ </button>
53+ ))}
54+ </div>
55+ </div>
56+
57+ {v.f === 'custom' && (
58+ <div className="field grow">
59+ <span className="field-label">f(x) =</span>
60+ <input
61+ className="expr"
62+ value={v.fexpr}
63+ spellCheck={false}
64+ placeholder="exp(-x.^2) .* cos(3*x)"
65+ onChange={(e) => onChange({ fexpr: e.target.value })}
66+ />
67+ </div>
68+ )}
69+
70+ <div className="field">
71+ <span className="field-label">interval</span>
72+ <div className="interval">
73+ <input
74+ type="number"
75+ value={v.a}
76+ step={1}
77+ onChange={(e) => {
78+ const a = Number(e.target.value)
79+ if (isFinite(a) && a < v.b) onChange({ a })
80+ }}
81+ />
82+ <span>to</span>
83+ <input
84+ type="number"
85+ value={v.b}
86+ step={1}
87+ onChange={(e) => {
88+ const b = Number(e.target.value)
89+ if (isFinite(b) && b > v.a) onChange({ b })
90+ }}
91+ />
92+ </div>
93+ </div>
94+
95+ {showN && (
96+ <div className="field grow">
97+ <span className="field-label">
98+ nodes n = <b>{v.n}</b>
99+ </span>
100+ <input
101+ type="range"
102+ min={2}
103+ max={80}
104+ value={v.n}
105+ disabled={busy}
106+ onChange={(e) => {
107+ const n = Number(e.target.value)
108+ onChange({ n, d: Math.min(v.d, Math.min(n, 12)) })
109+ }}
110+ />
111+ </div>
112+ )}
113+
114+ <div className="field grow">
115+ <span className="field-label">
116+ blend degree d = <b>{v.d}</b>
117+ </span>
118+ <input
119+ type="range"
120+ min={0}
121+ max={dMax}
122+ value={Math.min(v.d, dMax)}
123+ disabled={busy}
124+ onChange={(e) => onChange({ d: Number(e.target.value) })}
125+ />
126+ </div>
127+
128+ <div className="field">
129+ <span className="field-label">node distribution</span>
130+ <div className="chips">
131+ {NODES.map((nd) => (
132+ <button
133+ key={nd.id}
134+ title={nd.hint}
135+ className={`chip ${v.nodes === nd.id ? 'on' : ''}`}
136+ onClick={() => onChange({ nodes: nd.id })}
137+ >
138+ {nd.label}
139+ </button>
140+ ))}
141+ {v.nodes === 'random' && (
142+ <button className="chip" onClick={() => onChange({ seed: v.seed + 1 })} title="new draw">
143+ reroll
144+ </button>
145+ )}
146+ </div>
147+ </div>
148+ </div>
149+ )
150+}
src/panels/ConvergencePanel.tsxadded+270−0View file
@@ -0,0 +1,270 @@
1+import Plot from '../plot/Plot.tsx'
2+import Legend, { type LegendItem } from '../plot/Legend.tsx'
3+import { dColor, ink, series } from '../plot/palette.ts'
4+import { linePath, type Frame } from '../plot/scales.ts'
5+import type { ConvergeOut, FuncName, NodeKind, Num } from '../engine/types.ts'
6+
7+interface Props {
8+ out: ConvergeOut | null
9+ running: boolean
10+ stale: boolean
11+ f: FuncName
12+ nodes: NodeKind
13+ ds: number[]
14+ maxN: number
15+ showSpline: boolean
16+ showPoly: boolean
17+ onChange: (patch: { ds?: number[]; maxN?: number; showSpline?: boolean; showPoly?: boolean }) => void
18+ onRun: () => void
19+}
20+
21+const D_CHOICES = [0, 1, 2, 3, 4, 5, 6, 8]
22+const N_CHOICES = [80, 160, 320, 640]
23+
24+const fmt = (v: Num) => (v == null || !isFinite(v) ? '-' : v.toExponential(1))
25+const fmtOrd = (v: Num) => (v == null || !isFinite(v) ? '' : v.toFixed(1))
26+
27+const F_LABEL: Record<FuncName, string> = {
28+ runge: '1 / (1 + x²)',
29+ sine: 'sin x',
30+ abs: '|x|',
31+ custom: 'the custom function',
32+}
33+
34+export default function ConvergencePanel(props: Props) {
35+ const { out, running, stale, f, nodes, ds, maxN, showSpline, showPoly, onChange, onRun } = props
36+
37+ const toggleD = (d: number) => {
38+ const next = ds.includes(d) ? ds.filter((v) => v !== d) : [...ds, d].sort((a, b) => a - b)
39+ if (next.length > 0 && next.length <= 6) onChange({ ds: next })
40+ }
41+
42+ return (
43+ <div className="panel">
44+ <p className="panel-lede">
45+ Theorem 2: for d ≥ 1 the error is O(h<sup>d+1</sup>) as h → 0, whatever the nodes look like,
46+ provided f is smooth enough. On log-log axes that is a straight line of slope −(d+1), and the
47+ slopes measured between consecutive n are printed in the table. With uniform nodes and Runge's
48+ function this reproduces Table 1; turn the spline on for Tables 3 and 4. Currently fitting{' '}
49+ <b>{F_LABEL[f]}</b> on <b>{nodes}</b> nodes.
50+ </p>
51+
52+ <div className="conv-controls">
53+ <div className="field">
54+ <span className="field-label">blend degrees d</span>
55+ <div className="chips">
56+ {D_CHOICES.map((d) => (
57+ <button
58+ key={d}
59+ className={`chip ${ds.includes(d) ? 'on' : ''}`}
60+ onClick={() => toggleD(d)}
61+ disabled={running}
62+ >
63+ {d}
64+ </button>
65+ ))}
66+ </div>
67+ </div>
68+ <div className="field">
69+ <span className="field-label">largest n</span>
70+ <div className="chips">
71+ {N_CHOICES.map((n) => (
72+ <button
73+ key={n}
74+ className={`chip ${maxN === n ? 'on' : ''}`}
75+ onClick={() => onChange({ maxN: n })}
76+ disabled={running}
77+ >
78+ {n}
79+ </button>
80+ ))}
81+ </div>
82+ </div>
83+ <div className="field">
84+ <span className="field-label">compare with</span>
85+ <div className="chips">
86+ <button
87+ className={`chip ${showSpline ? 'on' : ''}`}
88+ onClick={() => onChange({ showSpline: !showSpline })}
89+ disabled={running}
90+ >
91+ cubic spline
92+ </button>
93+ <button
94+ className={`chip ${showPoly ? 'on' : ''}`}
95+ onClick={() => onChange({ showPoly: !showPoly })}
96+ disabled={running}
97+ >
98+ polynomial
99+ </button>
100+ </div>
101+ </div>
102+ <button className="primary" onClick={onRun} disabled={running}>
103+ {running ? 'Running…' : out == null ? 'Run study ▶' : 'Re-run ▶'}
104+ </button>
105+ </div>
106+
107+ {!out ? (
108+ <p className="panel-note muted">
109+ The study refits the interpolant at every n and every d, so it is the one thing on this page that
110+ does not run on its own. Press Run.
111+ </p>
112+ ) : (
113+ <ConvergenceChart out={out} showSpline={showSpline} showPoly={showPoly} stale={stale} />
114+ )}
115+ </div>
116+ )
117+}
118+
119+function ConvergenceChart({
120+ out,
121+ showSpline,
122+ showPoly,
123+ stale,
124+}: {
125+ out: ConvergeOut
126+ showSpline: boolean
127+ showPoly: boolean
128+ stale: boolean
129+}) {
130+ const all: Num[] = [
131+ ...out.E.flat(),
132+ ...(showSpline ? (out.splineErr ?? []) : []),
133+ ...(showPoly ? (out.polyErr ?? []) : []),
134+ ]
135+ const finite = all.filter((v): v is number => v != null && isFinite(v) && v > 0)
136+ const lo = Math.min(...finite)
137+ const hi = Math.max(...finite)
138+ const yd: [number, number] = [Math.pow(10, Math.floor(Math.log10(lo)) - 0.3), Math.pow(10, Math.ceil(Math.log10(hi)) + 0.3)]
139+ const xd: [number, number] = [out.ns[0] * 0.85, out.ns[out.ns.length - 1] * 1.35]
140+
141+ const legend: LegendItem[] = out.ds.map((d, i) => ({
142+ label: `d = ${d}`,
143+ color: dColor(i, out.ds.length),
144+ }))
145+ if (showPoly && out.polyErr) legend.push({ label: 'polynomial (d = n)', color: series.poly })
146+ if (showSpline && out.splineErr) legend.push({ label: 'cubic spline', color: series.spline })
147+
148+ const dots = (f: Frame, vals: Num[], color: string) =>
149+ vals.map((v, j) =>
150+ v == null || !isFinite(v) || v <= 0 ? null : (
151+ <circle key={j} cx={f.sx(out.ns[j])} cy={f.sy(v)} r={3.5} fill={color} stroke="#151a21" strokeWidth={1.5} />
152+ ),
153+ )
154+
155+ /** the label goes at the right end of the curve, on its last finite point */
156+ const endLabel = (f: Frame, vals: Num[], color: string, text: string) => {
157+ for (let j = vals.length - 1; j >= 0; j--) {
158+ const v = vals[j]
159+ if (v != null && isFinite(v) && v > 0) {
160+ return (
161+ <text x={f.sx(out.ns[j]) + 8} y={f.sy(v) + 4} fill={color} fontSize={11} fontWeight={600}>
162+ {text}
163+ </text>
164+ )
165+ }
166+ }
167+ return null
168+ }
169+
170+ return (
171+ <>
172+ {stale && <div className="stale-note">Showing the previous study — the settings have changed since.</div>}
173+ <Legend items={legend} />
174+ <Plot
175+ height={360}
176+ xDomain={xd}
177+ yDomain={yd}
178+ xLog
179+ yLog
180+ xLabel="n"
181+ yLabel="max |r - f|"
182+ margin={{ l: 62, r: 58 }}
183+ xTicks={out.ns}
184+ formatX={(v) => String(Math.round(v))}
185+ >
186+ {(f) => (
187+ <>
188+ {showPoly && out.polyErr && (
189+ <>
190+ <path d={linePath(out.ns, out.polyErr, f)} fill="none" stroke={series.poly} strokeWidth={2} />
191+ {dots(f, out.polyErr, series.poly)}
192+ {endLabel(f, out.polyErr, series.poly, 'poly')}
193+ </>
194+ )}
195+ {showSpline && out.splineErr && (
196+ <>
197+ <path d={linePath(out.ns, out.splineErr, f)} fill="none" stroke={series.spline} strokeWidth={2} />
198+ {dots(f, out.splineErr, series.spline)}
199+ {endLabel(f, out.splineErr, series.spline, 'spline')}
200+ </>
201+ )}
202+ {out.E.map((row, i) => {
203+ const c = dColor(i, out.ds.length)
204+ return (
205+ <g key={i}>
206+ <path d={linePath(out.ns, row, f)} fill="none" stroke={c} strokeWidth={2.5} />
207+ {dots(f, row, c)}
208+ {endLabel(f, row, c, `d = ${out.ds[i]}`)}
209+ </g>
210+ )
211+ })}
212+ </>
213+ )}
214+ </Plot>
215+
216+ <h4 className="sub">The same numbers</h4>
217+ <div className="table-wrap">
218+ <table className="conv-table">
219+ <thead>
220+ <tr>
221+ <th>n</th>
222+ {out.ds.map((d) => (
223+ <th key={d} colSpan={2}>
224+ d = {d}
225+ </th>
226+ ))}
227+ {showSpline && out.splineErr && <th colSpan={2}>cubic spline</th>}
228+ {showPoly && out.polyErr && <th>polynomial</th>}
229+ </tr>
230+ <tr className="sub-head">
231+ <th />
232+ {out.ds.map((d) => [
233+ <th key={`e${d}`}>error</th>,
234+ <th key={`o${d}`}>order</th>,
235+ ])}
236+ {showSpline && out.splineErr && [<th key="se">error</th>, <th key="so">order</th>]}
237+ {showPoly && out.polyErr && <th>error</th>}
238+ </tr>
239+ </thead>
240+ <tbody>
241+ {out.ns.map((n, j) => (
242+ <tr key={n}>
243+ <td className="n-cell">{n}</td>
244+ {out.ds.map((d, i) => [
245+ <td key={`e${d}`}>{fmt(out.E[i][j])}</td>,
246+ <td key={`o${d}`} className="ord">
247+ {fmtOrd(out.orders[i][j])}
248+ </td>,
249+ ])}
250+ {showSpline && out.splineErr && [
251+ <td key="se">{fmt(out.splineErr[j])}</td>,
252+ <td key="so" className="ord">
253+ {fmtOrd(out.splineOrders?.[j] ?? null)}
254+ </td>,
255+ ]}
256+ {showPoly && out.polyErr && <td>{fmt(out.polyErr[j])}</td>}
257+ </tr>
258+ ))}
259+ </tbody>
260+ </table>
261+ </div>
262+ <p className="panel-note" style={{ color: ink.muted }}>
263+ The order column is log(e<sub>prev</sub> / e) / log(n / n<sub>prev</sub>), so d + 1 is what Theorem 2
264+ predicts for d ≥ 1. Where a row of errors stops falling, it has reached the point at which the
265+ weights themselves, which grow like h<sup>−d</sup>, cost more accuracy than the higher order
266+ buys.
267+ </p>
268+ </>
269+ )
270+}
src/panels/InterpolantPanel.tsxadded+185−0View file
@@ -0,0 +1,185 @@
1+import { useState } from 'react'
2+import Plot from '../plot/Plot.tsx'
3+import Legend, { type LegendItem } from '../plot/Legend.tsx'
4+import { ink, series } from '../plot/palette.ts'
5+import { extent, linePath, nearestIndex, padDomain, type Frame } from '../plot/scales.ts'
6+import type { ExploreOut, Num } from '../engine/types.ts'
7+
8+interface Props {
9+ out: ExploreOut
10+ showPoly: boolean
11+ showSpline: boolean
12+ onToggle: (which: 'poly' | 'spline', on: boolean) => void
13+}
14+
15+function fmtErr(v: Num): string {
16+ if (v == null || !isFinite(v)) return 'off scale'
17+ if (v === 0) return '0'
18+ return v.toExponential(1)
19+}
20+
21+/** Largest |value| of a series, ignoring the parts that ran off to infinity. */
22+function maxAbs(v: readonly Num[] | undefined, ref: readonly number[]): Num {
23+ if (!v) return null
24+ let m = 0
25+ let any = false
26+ for (let i = 0; i < v.length; i++) {
27+ const d = v[i]
28+ if (d == null || !isFinite(d)) continue
29+ any = true
30+ m = Math.max(m, Math.abs(d - ref[i]))
31+ }
32+ return any ? m : null
33+}
34+
35+export default function InterpolantPanel({ out, showPoly, showSpline, onToggle }: Props) {
36+ const [hoverX, setHoverX] = useState<number | null>(null)
37+
38+ // The polynomial interpolant at equispaced nodes is the whole point of the
39+ // paper's opening paragraph, and at n = 40 it is off scale by a factor of
40+ // 10^8. Scale the axis to f and to r, and let the rest leave the frame.
41+ const yd = padDomain(extent(out.ft, out.r, showSpline ? out.rspline : undefined), 0.1)
42+ const errOf = (v: readonly Num[] | undefined): Num[] | undefined =>
43+ v ? v.map((d, i) => (d == null ? null : d - out.ft[i])) : undefined
44+ const errPoly = showPoly ? errOf(out.rpoly) : undefined
45+ const errSpline = showSpline ? errOf(out.rspline) : undefined
46+ const eAll = extent(out.err, errPoly, errSpline)
47+ const eMax = Math.max(Math.abs(eAll[0]), Math.abs(eAll[1]), 1e-16)
48+ const ed: [number, number] = [-eMax * 1.1, eMax * 1.1]
49+
50+ const hi = hoverX == null ? -1 : nearestIndex(out.t, hoverX)
51+ const at = (v: readonly Num[] | undefined) => (v && hi >= 0 ? v[hi] : null)
52+
53+ const legend: LegendItem[] = [
54+ { label: 'f(x)', color: ink.reference, dash: '5 4', value: undefined },
55+ { label: 'rational r(x)', color: series.r, value: fmtErr(out.maxerr) },
56+ ]
57+ if (showPoly) {
58+ legend.push({
59+ label: `polynomial (degree ${out.n})`,
60+ color: series.poly,
61+ value: fmtErr(maxAbs(out.rpoly, out.ft)),
62+ })
63+ }
64+ if (showSpline) {
65+ legend.push({ label: 'cubic spline', color: series.spline, value: fmtErr(maxAbs(out.rspline, out.ft)) })
66+ }
67+ legend.push({ label: `${out.n + 1} nodes`, color: ink.node, shape: 'dot' })
68+
69+ const nodes = (f: Frame) =>
70+ out.x.map((xv, i) => (
71+ <circle
72+ key={i}
73+ cx={f.sx(xv)}
74+ cy={f.sy(out.y[i])}
75+ r={3.2}
76+ fill={ink.node}
77+ stroke="#151a21"
78+ strokeWidth={1.5}
79+ />
80+ ))
81+
82+ return (
83+ <div className="panel">
84+ <p className="panel-lede">
85+ The rational interpolant r of equation (1) through {out.n + 1} nodes, with blend degree d = {out.d}.
86+ Turn on the degree-{out.n} polynomial to see what the paper's first page is about, and the clamped
87+ C<sup>2</sup> cubic spline for the comparison of Tables 3 and 4.
88+ </p>
89+
90+ <div className="row-controls">
91+ <label className="check">
92+ <input type="checkbox" checked={showPoly} onChange={(e) => onToggle('poly', e.target.checked)} />
93+ polynomial interpolant
94+ </label>
95+ <label className="check">
96+ <input type="checkbox" checked={showSpline} onChange={(e) => onToggle('spline', e.target.checked)} />
97+ cubic spline
98+ </label>
99+ </div>
100+
101+ <Legend items={legend} />
102+
103+ <Plot
104+ height={330}
105+ xDomain={[out.t[0], out.t[out.t.length - 1]]}
106+ yDomain={yd}
107+ xLabel="x"
108+ yLabel="f, r"
109+ onHoverX={setHoverX}
110+ hoverX={hoverX}
111+ >
112+ {(f) => (
113+ <>
114+ <path d={linePath(out.t, out.ft, f)} fill="none" stroke={ink.reference} strokeWidth={2} strokeDasharray="5 4" />
115+ {showSpline && out.rspline && (
116+ <path d={linePath(out.t, out.rspline, f)} fill="none" stroke={series.spline} strokeWidth={2} />
117+ )}
118+ {showPoly && out.rpoly && (
119+ <path d={linePath(out.t, out.rpoly, f)} fill="none" stroke={series.poly} strokeWidth={2} />
120+ )}
121+ <path d={linePath(out.t, out.r, f)} fill="none" stroke={series.r} strokeWidth={2.5} />
122+ {nodes(f)}
123+ </>
124+ )}
125+ </Plot>
126+
127+ {hi >= 0 && (
128+ <div className="readout">
129+ <span>
130+ x = <b>{out.t[hi].toFixed(3)}</b>
131+ </span>
132+ <span style={{ color: ink.reference }}>
133+ f = <b>{out.ft[hi].toFixed(6)}</b>
134+ </span>
135+ <span style={{ color: series.r }}>
136+ r = <b>{at(out.r)?.toFixed(6) ?? '-'}</b>
137+ </span>
138+ {showPoly && (
139+ <span style={{ color: series.poly }}>
140+ poly = <b>{fmtSigned(at(out.rpoly))}</b>
141+ </span>
142+ )}
143+ {showSpline && (
144+ <span style={{ color: series.spline }}>
145+ spline = <b>{at(out.rspline)?.toFixed(6) ?? '-'}</b>
146+ </span>
147+ )}
148+ </div>
149+ )}
150+
151+ <h4 className="sub">Error</h4>
152+ <p className="panel-note">
153+ r(x) − f(x) on the same grid. It vanishes at every node, by construction, and the largest of
154+ the bumps between them is the number Tables 1 to 4 tabulate.
155+ </p>
156+ <Plot
157+ height={190}
158+ xDomain={[out.t[0], out.t[out.t.length - 1]]}
159+ yDomain={ed}
160+ xLabel="x"
161+ yLabel="r - f"
162+ onHoverX={setHoverX}
163+ hoverX={hoverX}
164+ under={(f) => <line x1={0} x2={f.iw} y1={f.sy(0)} y2={f.sy(0)} stroke={ink.axis} strokeWidth={1} />}
165+ >
166+ {(f) => (
167+ <>
168+ {errSpline && <path d={linePath(out.t, errSpline, f)} fill="none" stroke={series.spline} strokeWidth={1.5} />}
169+ {errPoly && <path d={linePath(out.t, errPoly, f)} fill="none" stroke={series.poly} strokeWidth={1.5} />}
170+ <path d={linePath(out.t, out.err, f)} fill="none" stroke={series.r} strokeWidth={2} />
171+ {out.x.map((xv, i) => (
172+ <circle key={i} cx={f.sx(xv)} cy={f.sy(0)} r={2.2} fill={ink.node} stroke="#151a21" strokeWidth={1} />
173+ ))}
174+ </>
175+ )}
176+ </Plot>
177+ </div>
178+ )
179+}
180+
181+function fmtSigned(v: Num): string {
182+ if (v == null || !isFinite(v)) return 'off scale'
183+ if (Math.abs(v) >= 1e5) return v.toExponential(2)
184+ return v.toFixed(6)
185+}
src/panels/PolesPanel.tsxadded+331−0View file
@@ -0,0 +1,331 @@
1+import { useState } from 'react'
2+import Plot from '../plot/Plot.tsx'
3+import Legend from '../plot/Legend.tsx'
4+import { ink, series, status } from '../plot/palette.ts'
5+import { extent, linePath, padDomain, type Frame } from '../plot/scales.ts'
6+import type { ExploreOut } from '../engine/types.ts'
7+
8+interface Props {
9+ out: ExploreOut
10+ showClassical: boolean
11+ onToggleClassical: (on: boolean) => void
12+}
13+
14+export default function PolesPanel({ out, showClassical, onToggleClassical }: Props) {
15+ const [hoverRoot, setHoverRoot] = useState<number | null>(null)
16+ const p = out.poles
17+ if (!p) return <div className="panel">Loading.</div>
18+
19+ const xd: [number, number] = [p.t[0], p.t[p.t.length - 1]]
20+ const clean = p.realPoles.length === 0
21+ // a root counts as sitting on the real axis if its imaginary part is a
22+ // rounding error next to the interval it lives on
23+ const axisTol = (out.x[out.x.length - 1] - out.x[0]) * 1e-7
24+ const nOnAxis = p.rootsIm.filter((v) => Math.abs(v) <= axisTol).length
25+
26+ const rootRange = Math.max(
27+ ...p.rootsRe.map((v) => Math.abs(v - (out.x[0] + out.x[out.x.length - 1]) / 2)),
28+ ...p.rootsIm.map(Math.abs),
29+ (out.x[out.x.length - 1] - out.x[0]) / 2,
30+ )
31+ const cx0 = (out.x[0] + out.x[out.x.length - 1]) / 2
32+ const cd: [number, number] = [cx0 - rootRange * 1.12, cx0 + rootRange * 1.12]
33+ const cyd: [number, number] = [-rootRange * 1.12, rootRange * 1.12]
34+
35+ return (
36+ <div className="panel">
37+ <div className={`verdict ${clean ? 'good' : 'bad'}`}>
38+ <span className="verdict-mark" aria-hidden="true">
39+ {clean ? '✓' : '✕'}
40+ </span>
41+ <span>
42+ {clean ? (
43+ <>
44+ <b>No real poles.</b> The denominator keeps one sign across the whole real line, which is
45+ Theorem 1, and every one of its {p.rootsShown ? p.rootsRe.length : 'roots'} roots sits off the
46+ real axis.
47+ </>
48+ ) : (
49+ <>
50+ <b>
51+ {p.realPoles.length} real pole{p.realPoles.length === 1 ? '' : 's'}.
52+ </b>{' '}
53+ The denominator changes sign, so r blows up inside the interval. These weights are not the ones
54+ of equation (18).
55+ </>
56+ )}
57+ </span>
58+ </div>
59+
60+ <p className="panel-lede">
61+ Writing r as a quotient of polynomials (equation 7) puts everything on the denominator s of equation
62+ (10). Its zeros are exactly the poles of r, so Theorem 1 amounts to the claim that s never crosses
63+ zero.
64+ </p>
65+
66+ <h4 className="sub">The denominator on the real line</h4>
67+ <p className="panel-note">
68+ s(x) spans many orders of magnitude, so what is drawn is the signed n-th root of it. That leaves every
69+ sign and every zero exactly where it was, and brings the rest into a range that fits on a page. The
70+ range shown runs well past both ends of the interpolation interval, since Theorem 1 is a statement
71+ about all of <b>R</b>, not just [a, b].
72+ </p>
73+ <Plot
74+ height={200}
75+ xDomain={xd}
76+ yDomain={[-1.15, 1.15]}
77+ xLabel="x"
78+ yLabel="signed n-th root of s"
79+ margin={{ l: 74 }}
80+ under={(f) => (
81+ <>
82+ <rect
83+ x={f.sx(out.x[0])}
84+ width={f.sx(out.x[out.x.length - 1]) - f.sx(out.x[0])}
85+ y={0}
86+ height={f.ih}
87+ fill={ink.node}
88+ opacity={0.05}
89+ />
90+ <line x1={0} x2={f.iw} y1={f.sy(0)} y2={f.sy(0)} stroke={ink.axis} strokeWidth={1.5} />
91+ </>
92+ )}
93+ >
94+ {(f) => (
95+ <>
96+ <path d={linePath(p.t, p.u, f)} fill="none" stroke={series.r} strokeWidth={2.5} />
97+ {p.realPoles.map((xv, i) => (
98+ <g key={i}>
99+ <line
100+ x1={f.sx(xv)}
101+ x2={f.sx(xv)}
102+ y1={0}
103+ y2={f.ih}
104+ stroke={status.critical}
105+ strokeWidth={1.5}
106+ strokeDasharray="4 3"
107+ />
108+ <path
109+ d={`M${f.sx(xv) - 4} ${f.sy(0) - 4} L${f.sx(xv) + 4} ${f.sy(0) + 4} M${f.sx(xv) + 4} ${
110+ f.sy(0) - 4
111+ } L${f.sx(xv) - 4} ${f.sy(0) + 4}`}
112+ stroke={status.critical}
113+ strokeWidth={2}
114+ />
115+ </g>
116+ ))}
117+ {out.x.map((xv, i) => (
118+ <line
119+ key={i}
120+ x1={f.sx(xv)}
121+ x2={f.sx(xv)}
122+ y1={f.ih}
123+ y2={f.ih - 5}
124+ stroke={ink.node}
125+ strokeWidth={1.5}
126+ />
127+ ))}
128+ </>
129+ )}
130+ </Plot>
131+
132+ <h4 className="sub">Where the roots actually are</h4>
133+ {p.rootsShown ? (
134+ <>
135+ <p className="panel-note">
136+ The {p.rootsRe.length} roots of s in the complex plane. Theorem 1 says none of them are real, so
137+ none of them touch the horizontal axis. The grey band is the interpolation interval; ticks on the
138+ axis are the nodes.
139+ {nOnAxis > 0 && ' Roots that have landed on the axis are marked with a cross.'}
140+ </p>
141+ <Legend
142+ items={[
143+ { label: 'root of s', color: series.r, shape: 'dot' },
144+ ...(nOnAxis > 0 ? [{ label: 'root on the real axis: a pole', color: status.critical, shape: 'cross' as const }] : []),
145+ ]}
146+ />
147+ <Plot
148+ height={300}
149+ xDomain={cd}
150+ yDomain={cyd}
151+ xLabel="Re"
152+ yLabel="Im"
153+ under={(f: Frame) => (
154+ <>
155+ <rect
156+ x={f.sx(out.x[0])}
157+ width={f.sx(out.x[out.x.length - 1]) - f.sx(out.x[0])}
158+ y={0}
159+ height={f.ih}
160+ fill={ink.node}
161+ opacity={0.05}
162+ />
163+ <line x1={0} x2={f.iw} y1={f.sy(0)} y2={f.sy(0)} stroke={ink.axis} strokeWidth={1.5} />
164+ <line x1={f.sx(0)} x2={f.sx(0)} y1={0} y2={f.ih} stroke={ink.grid} strokeWidth={1} />
165+ </>
166+ )}
167+ overlay={
168+ hoverRoot != null && (
169+ <div className="tooltip">
170+ {p.rootsRe[hoverRoot].toPrecision(4)}
171+ {p.rootsIm[hoverRoot] >= 0 ? ' + ' : ' − '}
172+ {Math.abs(p.rootsIm[hoverRoot]).toPrecision(4)} i
173+ </div>
174+ )
175+ }
176+ >
177+ {(f) => (
178+ <>
179+ {out.x.map((xv, i) => (
180+ <line
181+ key={i}
182+ x1={f.sx(xv)}
183+ x2={f.sx(xv)}
184+ y1={f.sy(0) - 4}
185+ y2={f.sy(0) + 4}
186+ stroke={ink.node}
187+ strokeWidth={1.5}
188+ />
189+ ))}
190+ {p.rootsRe.map((re, i) => {
191+ const onAxis = Math.abs(p.rootsIm[i]) <= axisTol
192+ return onAxis ? (
193+ <path
194+ key={i}
195+ d={`M${f.sx(re) - 5} ${f.sy(p.rootsIm[i]) - 5} L${f.sx(re) + 5} ${
196+ f.sy(p.rootsIm[i]) + 5
197+ } M${f.sx(re) + 5} ${f.sy(p.rootsIm[i]) - 5} L${f.sx(re) - 5} ${f.sy(p.rootsIm[i]) + 5}`}
198+ stroke={status.critical}
199+ strokeWidth={2.5}
200+ strokeLinecap="round"
201+ onPointerEnter={() => setHoverRoot(i)}
202+ onPointerLeave={() => setHoverRoot(null)}
203+ />
204+ ) : (
205+ <circle
206+ key={i}
207+ cx={f.sx(re)}
208+ cy={f.sy(p.rootsIm[i])}
209+ r={hoverRoot === i ? 7 : 5}
210+ fill={series.r}
211+ stroke="#151a21"
212+ strokeWidth={2}
213+ onPointerEnter={() => setHoverRoot(i)}
214+ onPointerLeave={() => setHoverRoot(null)}
215+ />
216+ )
217+ })}
218+ </>
219+ )}
220+ </Plot>
221+ </>
222+ ) : (
223+ <p className="panel-note">
224+ {out.n > 40
225+ ? 'Not drawn above n = 40: recovering roots from the coefficients of a polynomial of that degree is not reliable enough to be worth showing. The sign test above still holds at any n.'
226+ : 'The denominator reduces to a constant, so it has no roots at all — which is one way for an interpolant to have no poles. That is what the Lagrange weights of equation (2) do.'}
227+ </p>
228+ )}
229+
230+ <h4 className="sub">The classical alternative</h4>
231+ <div className="row-controls">
232+ <label className="check">
233+ <input
234+ type="checkbox"
235+ checked={showClassical}
236+ onChange={(e) => onToggleClassical(e.target.checked)}
237+ />
238+ fit p<sub>M</sub> / q<sub>N</sub> with M + N = n
239+ </label>
240+ </div>
241+ <p className="panel-note">
242+ The construction the paper's introduction rejects: fit the same data with a quotient of polynomials of
243+ degrees M and N summing to n. It is a good approximation where it is finite, and there is no way to
244+ stop it putting poles wherever it likes.
245+ </p>
246+ {showClassical && p.classical ? (
247+ <>
248+ <Legend
249+ items={[
250+ { label: 'f(x)', color: ink.reference, dash: '5 4' },
251+ { label: 'rational r(x), this page', color: series.r },
252+ { label: 'classical p_M / q_N', color: series.classical },
253+ ...(p.classicalPoles && p.classicalPoles.length
254+ ? [{ label: 'its poles', color: status.critical, shape: 'cross' as const }]
255+ : []),
256+ ]}
257+ />
258+ <ClassicalPlot out={out} />
259+ </>
260+ ) : (
261+ <p className="panel-note muted">Turn it on to draw it.</p>
262+ )}
263+ </div>
264+ )
265+}
266+
267+function ClassicalPlot({ out }: { out: ExploreOut }) {
268+ const p = out.poles!
269+ const xd: [number, number] = [p.t[0], p.t[p.t.length - 1]]
270+ const yd = padDomain(extent(out.ft, out.r), 0.45)
271+ const inView = (p.classicalPoles ?? []).filter((v) => v >= xd[0] && v <= xd[1])
272+ return (
273+ <>
274+ <Plot
275+ height={260}
276+ xDomain={xd}
277+ yDomain={yd}
278+ xLabel="x"
279+ yLabel="f, r, classical"
280+ under={(f) => (
281+ <rect
282+ x={f.sx(out.x[0])}
283+ width={f.sx(out.x[out.x.length - 1]) - f.sx(out.x[0])}
284+ y={0}
285+ height={f.ih}
286+ fill={ink.node}
287+ opacity={0.05}
288+ />
289+ )}
290+ >
291+ {(f) => (
292+ <>
293+ {inView.map((xv, i) => (
294+ <line
295+ key={i}
296+ x1={f.sx(xv)}
297+ x2={f.sx(xv)}
298+ y1={0}
299+ y2={f.ih}
300+ stroke={status.critical}
301+ strokeWidth={1.5}
302+ strokeDasharray="4 3"
303+ />
304+ ))}
305+ <path d={linePath(out.t, out.ft, f)} fill="none" stroke={ink.reference} strokeWidth={2} strokeDasharray="5 4" />
306+ <path d={linePath(p.t, p.classical!, f)} fill="none" stroke={series.classical} strokeWidth={2} />
307+ <path d={linePath(out.t, out.r, f)} fill="none" stroke={series.r} strokeWidth={2.5} />
308+ {out.x.map((xv, i) => (
309+ <circle
310+ key={i}
311+ cx={f.sx(xv)}
312+ cy={f.sy(out.y[i])}
313+ r={3}
314+ fill={ink.node}
315+ stroke="#151a21"
316+ strokeWidth={1.5}
317+ />
318+ ))}
319+ </>
320+ )}
321+ </Plot>
322+ <p className="panel-note">
323+ {inView.length === 0
324+ ? 'On this data it happens to have no poles in the window shown. Change n, or the function, or the nodes, and that is not something you can rely on.'
325+ : `${inView.length} real pole${inView.length === 1 ? '' : 's'} in the window, at ${inView
326+ .map((v) => v.toFixed(3))
327+ .join(', ')}.`}
328+ </p>
329+ </>
330+ )
331+}
src/plot/Legend.tsxadded+47−0View file
@@ -0,0 +1,47 @@
1+import { ink } from './palette.ts'
2+
3+export interface LegendItem {
4+ label: string
5+ color: string
6+ /** stroke-dasharray for a line swatch */
7+ dash?: string
8+ shape?: 'line' | 'dot' | 'cross'
9+ /** e.g. the max error for this series, shown after the label */
10+ value?: string
11+}
12+
13+/**
14+ * Identity is never carried by colour alone: every series that appears in a
15+ * plot appears here too, and the ones the reader most needs are also labelled
16+ * directly on the plot.
17+ */
18+export default function Legend({ items }: { items: LegendItem[] }) {
19+ return (
20+ <div className="legend">
21+ {items.map((it) => (
22+ <span className="legend-item" key={it.label}>
23+ <svg width={18} height={10} aria-hidden="true">
24+ {it.shape === 'dot' ? (
25+ <circle cx={9} cy={5} r={3.5} fill={it.color} stroke={ink.grid} strokeWidth={1} />
26+ ) : it.shape === 'cross' ? (
27+ <path d="M5 1 L13 9 M13 1 L5 9" stroke={it.color} strokeWidth={2} strokeLinecap="round" />
28+ ) : (
29+ <line
30+ x1={1}
31+ x2={17}
32+ y1={5}
33+ y2={5}
34+ stroke={it.color}
35+ strokeWidth={2}
36+ strokeDasharray={it.dash}
37+ strokeLinecap="round"
38+ />
39+ )}
40+ </svg>
41+ <span>{it.label}</span>
42+ {it.value && <span className="legend-value">{it.value}</span>}
43+ </span>
44+ ))}
45+ </div>
46+ )
47+}
src/plot/Plot.tsxadded+194−0View file
@@ -0,0 +1,194 @@
1+import { useCallback, useEffect, useId, useRef, useState, type ReactNode } from 'react'
2+import { ink } from './palette.ts'
3+import { decadeTicks, formatTick, niceTicks, type Frame } from './scales.ts'
4+
5+export interface PlotProps {
6+ height: number
7+ xDomain: [number, number]
8+ yDomain: [number, number]
9+ xLog?: boolean
10+ yLog?: boolean
11+ xLabel?: string
12+ yLabel?: string
13+ xTicks?: number[]
14+ yTicks?: number[]
15+ formatX?: (v: number) => string
16+ formatY?: (v: number) => string
17+ margin?: Partial<{ l: number; r: number; t: number; b: number }>
18+ /** drawn under the data, inside the frame but outside the clip */
19+ under?: (f: Frame) => ReactNode
20+ children: (f: Frame) => ReactNode
21+ /** enables the crosshair; called with the data x under the pointer, or null */
22+ onHoverX?: (x: number | null) => void
23+ /** data x at which to draw the crosshair (usually what onHoverX last gave) */
24+ hoverX?: number | null
25+ /** rendered as an absolutely positioned box over the plot */
26+ overlay?: ReactNode
27+}
28+
29+const DEFAULT_MARGIN = { l: 54, r: 14, t: 10, b: 30 }
30+
31+/** Measures its own width so plots reflow with the panel. */
32+function useWidth(): [(el: HTMLDivElement | null) => void, number] {
33+ const [width, setWidth] = useState(640)
34+ const obs = useRef<ResizeObserver | null>(null)
35+ const ref = useCallback((el: HTMLDivElement | null) => {
36+ obs.current?.disconnect()
37+ if (!el) return
38+ const ro = new ResizeObserver((entries) => {
39+ const w = entries[0]?.contentRect.width
40+ if (w && w > 0) setWidth(w)
41+ })
42+ ro.observe(el)
43+ obs.current = ro
44+ setWidth(el.clientWidth || 640)
45+ }, [])
46+ useEffect(() => () => obs.current?.disconnect(), [])
47+ return [ref, width]
48+}
49+
50+export default function Plot(props: PlotProps) {
51+ const {
52+ height,
53+ xDomain,
54+ yDomain,
55+ xLog = false,
56+ yLog = false,
57+ xLabel,
58+ yLabel,
59+ formatX = formatTick,
60+ formatY = formatTick,
61+ under,
62+ children,
63+ onHoverX,
64+ hoverX,
65+ overlay,
66+ } = props
67+ const m = { ...DEFAULT_MARGIN, ...props.margin }
68+ const [hostRef, width] = useWidth()
69+ const clipId = useId().replace(/:/g, '')
70+ const svgRef = useRef<SVGSVGElement>(null)
71+
72+ const iw = Math.max(10, width - m.l - m.r)
73+ const ih = Math.max(10, height - m.t - m.b)
74+
75+ const fwd = (v: number, [lo, hi]: [number, number], log: boolean, span: number, flip: boolean) => {
76+ const t = log
77+ ? (Math.log10(Math.max(v, Number.MIN_VALUE)) - Math.log10(lo)) / (Math.log10(hi) - Math.log10(lo))
78+ : (v - lo) / (hi - lo)
79+ return flip ? span * (1 - t) : span * t
80+ }
81+
82+ const f: Frame = {
83+ sx: (v) => fwd(v, xDomain, xLog, iw, false),
84+ sy: (v) => fwd(v, yDomain, yLog, ih, true),
85+ ix: (px) => {
86+ const t = px / iw
87+ return xLog
88+ ? Math.pow(10, Math.log10(xDomain[0]) + t * (Math.log10(xDomain[1]) - Math.log10(xDomain[0])))
89+ : xDomain[0] + t * (xDomain[1] - xDomain[0])
90+ },
91+ iw,
92+ ih,
93+ xDomain,
94+ yDomain,
95+ }
96+
97+ const xt = props.xTicks ?? (xLog ? decadeTicks(xDomain[0], xDomain[1]) : niceTicks(xDomain[0], xDomain[1], 7))
98+ const yt = props.yTicks ?? (yLog ? decadeTicks(yDomain[0], yDomain[1]) : niceTicks(yDomain[0], yDomain[1], 5))
99+
100+ const pointer = (e: React.PointerEvent) => {
101+ if (!onHoverX) return
102+ const rect = svgRef.current?.getBoundingClientRect()
103+ if (!rect) return
104+ const px = e.clientX - rect.left - m.l
105+ onHoverX(px < -4 || px > iw + 4 ? null : f.ix(Math.min(iw, Math.max(0, px))))
106+ }
107+
108+ return (
109+ <div className="plot-host" ref={hostRef}>
110+ <svg
111+ ref={svgRef}
112+ width={width}
113+ height={height}
114+ role="img"
115+ onPointerMove={pointer}
116+ onPointerLeave={() => onHoverX?.(null)}
117+ >
118+ <defs>
119+ <clipPath id={clipId}>
120+ <rect x={0} y={0} width={iw} height={ih} />
121+ </clipPath>
122+ </defs>
123+ <g transform={`translate(${m.l},${m.t})`}>
124+ {/* gridlines, recessive */}
125+ {xt.map((v) => (
126+ <line key={`gx${v}`} x1={f.sx(v)} x2={f.sx(v)} y1={0} y2={ih} stroke={ink.grid} strokeWidth={1} />
127+ ))}
128+ {yt.map((v) => (
129+ <line key={`gy${v}`} x1={0} x2={iw} y1={f.sy(v)} y2={f.sy(v)} stroke={ink.grid} strokeWidth={1} />
130+ ))}
131+
132+ {under?.(f)}
133+
134+ <g clipPath={`url(#${clipId})`}>{children(f)}</g>
135+
136+ {hoverX != null && hoverX >= Math.min(...xDomain) && hoverX <= Math.max(...xDomain) && (
137+ <line
138+ x1={f.sx(hoverX)}
139+ x2={f.sx(hoverX)}
140+ y1={0}
141+ y2={ih}
142+ stroke={ink.secondary}
143+ strokeWidth={1}
144+ strokeDasharray="3 3"
145+ opacity={0.7}
146+ pointerEvents="none"
147+ />
148+ )}
149+
150+ {/* axes */}
151+ <line x1={0} x2={iw} y1={ih} y2={ih} stroke={ink.axis} strokeWidth={1} />
152+ <line x1={0} x2={0} y1={0} y2={ih} stroke={ink.axis} strokeWidth={1} />
153+ {xt.map((v) => (
154+ <text
155+ key={`tx${v}`}
156+ x={f.sx(v)}
157+ y={ih + 15}
158+ fill={ink.muted}
159+ fontSize={11}
160+ textAnchor="middle"
161+ style={{ fontVariantNumeric: 'tabular-nums' }}
162+ >
163+ {formatX(v)}
164+ </text>
165+ ))}
166+ {yt.map((v) => (
167+ <text
168+ key={`ty${v}`}
169+ x={-7}
170+ y={f.sy(v) + 4}
171+ fill={ink.muted}
172+ fontSize={11}
173+ textAnchor="end"
174+ style={{ fontVariantNumeric: 'tabular-nums' }}
175+ >
176+ {formatY(v)}
177+ </text>
178+ ))}
179+ {xLabel && (
180+ <text x={iw} y={ih + 27} fill={ink.muted} fontSize={11} textAnchor="end">
181+ {xLabel}
182+ </text>
183+ )}
184+ {yLabel && (
185+ <text x={-m.l + 4} y={-1} fill={ink.muted} fontSize={11} textAnchor="start">
186+ {yLabel}
187+ </text>
188+ )}
189+ </g>
190+ </svg>
191+ {overlay}
192+ </div>
193+ )
194+}
src/plot/palette.tsadded+70−0View file
@@ -0,0 +1,70 @@
1+// Colours, assigned by the job each one does. Validated against the dark chart
2+// surface #151a21 with the data-viz validator:
3+//
4+// categorical, all pairs r / polynomial / spline worst CVD dE 9.4
5+// categorical, all pairs r / classical rational worst CVD dE 27.4
6+// ordinal ramp the six steps of dRamp monotone L, gaps ok
7+// diverging pair positive / negative weights worst CVD dE 19.2
8+//
9+// Every colour is a documented step from the reference palette; none are
10+// eyeballed.
11+
12+export const surface = '#151a21'
13+
14+/** Identity: which interpolant a curve is. Fixed slots, never reassigned. */
15+export const series = {
16+ /** the Floater-Hormann rational interpolant */
17+ r: '#3987e5',
18+ /** the degree-n polynomial interpolant */
19+ poly: '#d95926',
20+ /** the clamped C^2 cubic spline */
21+ spline: '#199e70',
22+ /** the classical rational interpolant p_M / q_N */
23+ classical: '#c98500',
24+} as const
25+
26+/** Polarity: the sign of a barycentric weight. */
27+export const diverging = {
28+ pos: '#3987e5',
29+ neg: '#e66767',
30+ mid: '#383835',
31+} as const
32+
33+/** State. A real pole is a failure, and always ships with a label and a
34+ * different marker shape, never colour alone. */
35+export const status = {
36+ critical: '#d03b3b',
37+ good: '#0ca30c',
38+ warning: '#fab219',
39+} as const
40+
41+/**
42+ * Ordinal: position in a sequence. Used for the blend degree d, where the
43+ * order is the meaning, so the reader should see it in the colour. Six steps
44+ * of the blue ramp, light end kept clear of the surface.
45+ */
46+export const dRamp = ['#184f95', '#256abf', '#3987e5', '#6da7ec', '#9ec5f4', '#cde2fb'] as const
47+
48+export function dColor(index: number, count: number): string {
49+ if (count <= 1) return dRamp[2]
50+ const k = Math.round((index / (count - 1)) * (dRamp.length - 1))
51+ return dRamp[Math.min(dRamp.length - 1, Math.max(0, k))]
52+}
53+
54+/** Chart chrome and ink. */
55+export const ink = {
56+ primary: '#e6e9ef',
57+ secondary: '#c3c2b7',
58+ muted: '#898781',
59+ grid: '#232833',
60+ axis: '#38414f',
61+ /** the exact function f: the reference the interpolants are measured against,
62+ * deliberately not given a series slot */
63+ reference: '#8b95a5',
64+ /** the interpolation nodes: they belong to the data, not to any one method */
65+ node: '#e8e6df',
66+ /** the family of local polynomials and blending functions, drawn as a mass */
67+ family: '#46536a',
68+ /** the one member of that family the reader is following */
69+ familyHi: '#d95926',
70+} as const
src/plot/scales.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+import type { Num } from '../engine/types.ts'
2+
3+export interface Frame {
4+ /** data x to pixel x, within the plotting area */
5+ sx: (v: number) => number
6+ /** data y to pixel y */
7+ sy: (v: number) => number
8+ /** pixel x back to data x, for hover */
9+ ix: (px: number) => number
10+ iw: number
11+ ih: number
12+ xDomain: [number, number]
13+ yDomain: [number, number]
14+}
15+
16+export function niceTicks(lo: number, hi: number, target = 6): number[] {
17+ if (!isFinite(lo) || !isFinite(hi) || hi <= lo) return []
18+ const raw = (hi - lo) / target
19+ const mag = Math.pow(10, Math.floor(Math.log10(raw)))
20+ const norm = raw / mag
21+ const step = (norm >= 5 ? 10 : norm >= 2 ? 5 : norm >= 1 ? 2 : 1) * mag
22+ const out: number[] = []
23+ for (let v = Math.ceil(lo / step) * step; v <= hi + step * 1e-9; v += step) {
24+ out.push(Math.abs(v) < step * 1e-9 ? 0 : v)
25+ }
26+ return out
27+}
28+
29+/** Decade ticks, thinned so the labels do not collide over a wide range. */
30+export function decadeTicks(lo: number, hi: number, maxCount = 9): number[] {
31+ if (!(lo > 0) || !(hi > 0)) return []
32+ const a = Math.floor(Math.log10(lo))
33+ const b = Math.ceil(Math.log10(hi))
34+ const every = Math.max(1, Math.ceil((b - a + 1) / maxCount))
35+ const out: number[] = []
36+ for (let e = a; e <= b; e += every) {
37+ const v = Math.pow(10, e)
38+ if (v >= lo * 0.999 && v <= hi * 1.001) out.push(v)
39+ }
40+ return out
41+}
42+
43+export function formatTick(v: number): string {
44+ if (v === 0) return '0'
45+ const a = Math.abs(v)
46+ if (a >= 1e5 || a < 1e-3) {
47+ const e = Math.round(Math.log10(a))
48+ if (Math.abs(a - Math.pow(10, e)) < Math.pow(10, e) * 1e-6) return `1e${e}`
49+ return v.toExponential(0)
50+ }
51+ if (Number.isInteger(v)) return String(v)
52+ return String(Number(v.toPrecision(3)))
53+}
54+
55+export function padDomain([lo, hi]: [number, number], frac = 0.06): [number, number] {
56+ if (!isFinite(lo) || !isFinite(hi)) return [0, 1]
57+ if (hi === lo) return [lo - 0.5, hi + 0.5]
58+ const p = (hi - lo) * frac
59+ return [lo - p, hi + p]
60+}
61+
62+/** Range of the finite values in one or more series. */
63+export function extent(...series: (readonly Num[] | undefined)[]): [number, number] {
64+ let lo = Infinity
65+ let hi = -Infinity
66+ for (const s of series) {
67+ if (!s) continue
68+ for (const v of s) {
69+ if (v == null || !isFinite(v)) continue
70+ if (v < lo) lo = v
71+ if (v > hi) hi = v
72+ }
73+ }
74+ return isFinite(lo) ? [lo, hi] : [0, 1]
75+}
76+
77+/**
78+ * A polyline through (xs, ys). Runs of missing values break the path rather
79+ * than being bridged, and values far outside the frame are clamped to a band
80+ * just off-screen so that a curve running off to infinity still leaves the
81+ * frame in the right direction instead of producing unusable path data.
82+ */
83+export function linePath(xs: readonly number[], ys: readonly Num[], f: Frame): string {
84+ const bound = f.ih * 12
85+ let out = ''
86+ let pen = false
87+ for (let i = 0; i < xs.length && i < ys.length; i++) {
88+ const y = ys[i]
89+ if (y == null || !isFinite(y)) {
90+ pen = false
91+ continue
92+ }
93+ const px = f.sx(xs[i])
94+ const py = Math.min(bound, Math.max(-bound, f.sy(y)))
95+ out += `${pen ? 'L' : 'M'}${px.toFixed(2)} ${py.toFixed(2)}`
96+ pen = true
97+ }
98+ return out
99+}
100+
101+/** Index of the sample nearest a data-space x. */
102+export function nearestIndex(xs: readonly number[], v: number): number {
103+ if (xs.length === 0) return -1
104+ let lo = 0
105+ let hi = xs.length - 1
106+ if (v <= xs[lo]) return lo
107+ if (v >= xs[hi]) return hi
108+ while (hi - lo > 1) {
109+ const mid = (lo + hi) >> 1
110+ if (xs[mid] <= v) lo = mid
111+ else hi = mid
112+ }
113+ return v - xs[lo] <= xs[hi] - v ? lo : hi
114+}
tsconfig.app.jsonadded+25−0View file
@@ -0,0 +1,25 @@
1+{
2+ "compilerOptions": {
3+ "tsBuildInfoFile": "./node_modules/.tmp/tsconfig.app.tsbuildinfo",
4+ "target": "es2023",
5+ "lib": ["ES2023", "DOM", "WebWorker"],
6+ "module": "esnext",
7+ "types": ["vite/client"],
8+ "allowArbitraryExtensions": true,
9+ "skipLibCheck": true,
10+
11+ "moduleResolution": "bundler",
12+ "allowImportingTsExtensions": true,
13+ "verbatimModuleSyntax": true,
14+ "moduleDetection": "force",
15+ "noEmit": true,
16+ "jsx": "react-jsx",
17+
18+ "strict": true,
19+ "noUnusedLocals": true,
20+ "noUnusedParameters": true,
21+ "erasableSyntaxOnly": true,
22+ "noFallthroughCasesInSwitch": true
23+ },
24+ "include": ["src"]
25+}
tsconfig.jsonadded+4−0View file
@@ -0,0 +1,4 @@
1+{
2+ "files": [],
3+ "references": [{ "path": "./tsconfig.app.json" }, { "path": "./tsconfig.node.json" }]
4+}
tsconfig.node.jsonadded+21−0View file
@@ -0,0 +1,21 @@
1+{
2+ "compilerOptions": {
3+ "tsBuildInfoFile": "./node_modules/.tmp/tsconfig.node.tsbuildinfo",
4+ "target": "es2023",
5+ "lib": ["ES2023"],
6+ "module": "esnext",
7+ "types": ["node"],
8+ "skipLibCheck": true,
9+
10+ "moduleResolution": "bundler",
11+ "allowImportingTsExtensions": true,
12+ "verbatimModuleSyntax": true,
13+ "moduleDetection": "force",
14+ "noEmit": true,
15+
16+ "strict": true,
17+ "noUnusedLocals": true,
18+ "noUnusedParameters": true
19+ },
20+ "include": ["vite.config.ts"]
21+}
vite.config.tsadded+11−0View file
@@ -0,0 +1,11 @@
1+import { defineConfig } from 'vite'
2+import react from '@vitejs/plugin-react'
3+
4+export default defineConfig({
5+ plugins: [react()],
6+ // relative base so the same build works at the repository root and under
7+ // /barycentric-rational/ on GitHub Pages
8+ base: './',
9+ // .m files are pulled in with ?raw; tell vite they are assets, not modules
10+ assetsInclude: ['**/*.m'],
11+})