[lachlangrose]

Draft, stacked on #372 (and so #365). Its diff against master includes those commits until they merge; the change here is the last commit.

New examples (two batches) planned from the interpolation golden tests (packages/loop_interpolation/tests/golden). They use the same idea: synthetic geological models whose scalar field is known exactly, sampled for data, so every result is measured against the truth. They are self-contained (no shared helper, so the downloadable .py files work on their own) and use the interpolators directly. New gallery section examples/5_interpolation/, listed in examples/README.rst.

Examples

• 5a Choosing an interpolator: FDI, PLI and P2 on a gaussian dome. No interpolator wins on every measure: P2 has the lowest interior error and gradient-angle error, FDI the lowest crest error, P2 passes through the data while FDI has the largest misfit, FDI is fastest.

• 5b Regularisation: sweeps the regularisation weight on exact and noisy data (five noise draws averaged). With exact data the error rises with the weight; with noisy data a very small weight fits the noise (error 0.35 against 0.09), a very large one flattens the dome, and the best is near the default of 1.

• 5c Fold constraints with sparse data: two contacts and one orientation on a cylindrical fold. The three plain interpolators give straight tilted horizons about 45 degrees out; FDFold, DFI and P2Fold bring the direction to within about 5 degrees. Only P2Fold is also close in the scalar field; FDFold and DFI recover the shape but exaggerate the amplitude.

• 5d Refolded folds: a dome-and-basin interference pattern. With the axis rotation in the fold event the fold interpolators recover it (median direction error 2-10 degrees); with a constant axis the error is 18-19 degrees and the horizon collapses to straight stripes. The scalar-field error with a constant axis (0.065) is as small as with the rotating axis (0.062), so one error number hides the failure; the example says so. It also reports that FDFold ends 40 degrees from one of its orientation data, without explaining why.

• 5e Kinds of orientation data (InterpolatorBuilder): with six contacts the error falls from contacts only (0.172) through one tangent (0.158), two tangents (0.127) or a gradient direction (0.129), to a full normal (0.110).

• 5f The native solver: about 55 000 unknowns. The native field agrees with the Python one to 9e-9 and is faster (2.2 s against 3.9 s); Python-side peak memory is about 1 MB against 102 MB when the solver is chosen before setup (28 MB when chosen after); warm_start cuts a re-solve from 5062 to 3354 iterations. It also shows the trap: the native solver stops at its default maxiter=5000 without converging, logs a warning, leaves up_to_date False and returns a field that is 0.9 off, even though the residual in the warning (3e-8) looks close to the tolerance (1e-8). Raising maxiter fixes it. The user guide gets a short "checking convergence" section.

How I checked

• Each example was run before its text was written, and the conclusions were rewritten from the printed tables and the rendered figures. Several first drafts were wrong and were corrected this way: I had assumed FDI smooths the crest (it has the lowest crest error), that a plane would show differences (every interpolator reproduces it exactly, so the examples use curved models), and that the weakest weight stops the solver converging (it does not with the default solve).
• Results are stable run to run (fold example run twice; the fold interpolators shuffle rows with a shared generator).
• Found while writing the first batch (5a): the 3D supports take nsteps as a cell count, so my first grid was 17 cells (18 nodes) and extended past the unit cube. The examples pass nsteps_cells explicitly.
• Gallery build: succeeds with the same 5 warnings as docs: rewrite the visualisation guide, fix the API reference, add a docs check #372, none from the new section. The six examples add about 55 s (23 s the regularisation sweep, 21 s the native solver). ruff clean.

CI / requirements

The native example imports the compiled loop_interpolation_core. It has wheels on PyPI for Linux, macOS and Windows, so docs/requirements.txt, docs-check.yml and documentation.yml now install it. I could not run GitHub Actions to confirm the workflows.

Deliberately not included (findings, not examples)

• Saving and loading (recipes): (pinned by xfail tests in the PR below) to_recipe_json / from_recipe_json do not store the build parameters. A feature entry has only name, kind, faults, regions, so a model built with nelements=2000, interpolatortype="PLI" or regularisation=10 rebuilds with the defaults and differs (max difference 27, 12 and 62 on a field spanning about 560). With the defaults the round trip is exact (difference 0.0). An example would mislead until this is fixed or documented.
• Inequality constraints: with the default cg solver they are dropped with no warning, and with admm they are enforced only for some argument choices (default arguments leave the bounds violated; the arguments used in the existing builder test enforce them). Cause not diagnosed, so no example. (An earlier version of this note said they never changed the solution with admm; that was wrong, and is corrected in the xfail-test PR.)
• Intrusions: INTRUSIONS.md lists known regressions and unimplemented parts (intrusion_steps broken against the current stratigraphic column, geometric_scaling_parameters raising, no gradients). I would not put it in the gallery as it stands.

🤖 Generated with Claude Code