CLAUDE.md#

Notes for anyone, human or agent, working in this repository. Read README.md first for what the tutorial contains.

The shape of the thing#

A Vite + React app that embeds numbl as a library (numbl/browser, createNumblSession). This is the seqlab / barycentric-rational pattern, not the numbl-project.json + site-viewer pattern: there is no project file and the deploy is a plain vite build.

src/
  App.tsx              sidebar, hash routing, chapter switching
  ui/                  Lab (one runnable section), ChapterView, Controls, Markdown
  engine/              the numbl session wrapper
  plot/                a small hand-rolled SVG layer, plus a canvas field plot
  matlab/              the shared .m library, shown to the reader on its own page
  chapters/
    onedim/driver.m    the driver for chapters 1-4, plus their shared panels
    twodim/driver.m    the driver for chapters 5-7, plus their shared panels
    chNN/              per-chapter text.md, editable .m script, chapter.ts

A chapter is prose in one text.md with <!--lab:id--> markers where its labs go. A lab is a driver (fixed, hidden), an editable script (the method), a list of declarative controls, and a render function. Each lab owns one numbl session for as long as it is on screen: moving a control rewrites params.json and re-runs, while editing the script reboots the session, because the script is compiled into main.m at boot.

main.m is driver.m + "\n" + the editable script, so the script's functions become local functions of the driver and the driver can call them. That is the whole of the contract; each driver's header comment states which functions its script must define.

Gotchas in numbl 0.4.18#

These were all found the hard way, by measurement. Most of them cost a factor of twenty in speed rather than producing a wrong answer, which makes them easy to miss in a project whose subject is speed.

Invoke the script with run('main.m'), not main;#

numbl mis-binds the arguments of a script's local functions when the script is invoked by name from the REPL, which is what session.execute gives us: the callee sees its parameters as undefined. engine/runner.ts goes through run() for this reason.

The compiler's fast path, and the eight ways to fall off it#

numbl compiles MATLAB to JavaScript and infers types as it goes. A loop whose body it cannot type is run through a generic path that is 20 to 70 times slower. The following all cause that, and none of them are errors:

1. Cell or struct indexing inside a loop. C{k} or p.tol in a loop body deoptimises the whole loop and every function it calls. Hoist the read above the loop.
2. A multiple-output call inside a loop. [a, b] = f(...). Return one matrix and unpack it afterwards.
3. A slice assignment inside a loop. A(k, :) = row. Do the copy inside a called function, or with a scalar loop.
4. Bare tic / toc inside a loop. Use handles: t = tic; ...; toc(t).
5. Growing an array inside a loop, that is, v(k) = ... for a v that was never preallocated.
6. A loop variable with more than one kind of assignment. This is the subtle one. Y = V; for ... Y = f(Y); end is slow, and Y = f(V); for ... Y = f(Y); end is fast, because in the first case the type of Y at the top of the loop is the merge of "a parameter" and "whatever f returns". It applies to dead branches too: an if arm that assigns Y = reshape(...) slows the loop down even when it never runs. Give each loop variable one source, or use separate names.
7. Passing a loop-written array into a function call, even after the loop. for ... res(k) = ...; end; H = pack(res, ...) is slow, while H = [res; ...] is fast. This is why the iterative chapter hands the page its squared residuals and lets JavaScript take the square roots.
8. A call inside a loop, to a function that itself loops and calls further functions. The chain "loop → solver → kernel with loops" is compiled generically even when every individual piece is fine. Recursion in place of the outer loop is compiled normally, which is what drv_res, drv_methods and their relatives in the drivers are for. This one was worth a factor of 30 on the L-shaped conjugate gradient solve.

The rule of thumb that follows: hot loops should touch only preallocated numeric arrays and scalars, nothing they write should later be handed to a function, and any loop that drives a timed solve should be a recursion. When in doubt, measure: time the same call inside and outside a loop and compare.

sqrt, log and abs widen a scalar's type to complex#

v = v + r*r in a loop followed by v = sqrt(v) after it makes the loop 20 times slower, because sqrt in numbl can return a complex number and the inferred type of v is widened accordingly. Accumulate into one variable and take the root into another, or better, return the sum of squares and take the root in the caller. Every res_sumsq in this repository exists for that reason.

sqrt(NaN) is 0 + NaN i, and log(NaN) is complex#

MATLAB returns NaN. numbl returns a complex value whose real part is zero, so real(sqrt(NaN)) is 0, not NaN. Square-rooting an array with a NaN tail silently turns the tail into zeros, which in an error history reads as the iteration having converged perfectly. The iterative chapter therefore never takes a square root in MATLAB at all, and lp_orders guards its ratio before taking a logarithm.

jsondecode of a JSON array of arrays gives a cell of column vectors#

MATLAB gives a matrix when the rows are the same length. The L-shaped driver unpacks p.nsm by hand.

jsonencode refuses complex values, and flattens single-row matrices#

jsonencode([1 2 3]) is [1,2,3], not [[1,2,3]], so anything matrix-shaped goes through lp_rows, which converts it to a cell array of rows first. NaN and Inf cross as null, which is why the TypeScript side uses Num = number | null everywhere and every plot helper skips nulls.

Backslash densifies sparse matrices#

A \ b converts a sparse A to dense before calling LAPACK, so there is no sparse direct solver here. That is why Chapter 5 stops at $N = 48$ and why Chapters 6 to 10 are about methods that do not need one. The text says so where it matters.

Scalar loops are often faster than vectorised slices#

A hand-written double loop over a stencil beats the equivalent A(2:end-1,...) expression by about a factor of three, because the vectorised form allocates several temporaries per sweep. The kernels in ch08/iterative.m are written as loops for that reason, which also happens to make them readable.

The browser is slower than the CLI#

The CLI uses a native LAPACK addon; the browser falls back to a TypeScript implementation, and it announces which one on console.error. Dense solves are noticeably slower in the page than in npm run test:matlab. Nothing in the tutorial depends on the absolute numbers, but do not quote a CLI timing as if it were what a reader will see.

Measuring time honestly#

Timings come from tic/toc inside the MATLAB layer, around the solve only, so they exclude session boot, JSON encoding and rendering. Each method gets a warm-up call before its first timed solve, because the first call also pays for compilation, and any solve under 5 ms is repeated a few times with the minimum kept. Below about 0.1 ms the measurement is dominated by fixed overhead and the work-precision curves flatten out on the left; that is the clock, not the method.

Testing#

npm run test:matlab     # ~2 min: the .m layer through the numbl CLI
npm run build && npm run test:browser [--shots]

Every assertion in test:matlab is a number stated in the prose: an observed order, an error at a given resolution, a ratio between two methods. If one fails, either the code or the text is wrong, and the text is not allowed to win by default.

test:browser opens every chapter, presses run on the labs that wait to be asked, and checks that nothing throws and that the figures contain data. It cannot judge whether a picture is any good; --shots writes full-page screenshots into tmp/ so that a human, or an agent with vision, can look.

Plotting#

src/plot/ is a small SVG layer with no chart library. palette.ts records which colour does which job and the validator run that justifies the set; if you add a series, re-run the validator rather than picking a hex that looks nice. Fields are drawn into a canvas at one pixel per grid cell and scaled, so a 512-by-512 grid costs no more to draw than a 16-by-16 one, and points outside the domain are passed as null and left as background.

Things deliberately not done#

• A sparse direct solver, for the reason above.
• Graded meshes in Chapter 9. The chapter names local refinement and singular enrichment as the standard remedies and then demonstrates a third, because the third is the one that is short enough to read.
• Close-evaluation quadrature in Chapter 10. The failure is shown and quantified; the fix is a chapter of its own.
• Neumann and mixed boundary conditions, anywhere.