Skip to content

feat(loop-interpolation-core): add a low-memory native solver for Python and the browser - #354

Merged
lachlangrose merged 54 commits into
Loop3D:masterfrom
lachlangrose:feat/001-low-memory-interpolator
Oct 28, 2026
Merged

lachlangrose merged 54 commits into
Loop3D:masterfrom
lachlangrose:feat/001-low-memory-interpolator

Conversation

@lachlangrose

Copy link
Copy Markdown
Member

Summary

Adds solver="native" for the finite-difference interpolator. A new Rust package, packages/loop_interpolation_core, solves the system without building a sparse matrix. It runs as a Python wheel (PyO3), a Rust library, and a browser module (WebAssembly in a Worker). Spec: specs/001-low-memory-interpolator/ (spec.md, plan.md, tasks.md, measurements.md).

  • Memory: about 63 bytes per node. The 100 x 100 x 100 reference problem (R1) peaks at 63 MB. The Python matrix solver needs about 2.4 GB for the same problem (29 times more).
  • Same answer as the matrix solver: the native field equals the pinned FD golden fields to about 1e-9, and the browser result equals the Rust result bit for bit.
  • Python is unchanged unless you select it. solver="native" needs pip install loop-interpolation[native] (the core package). Other solvers behave as before. The core is imported only inside the adapter.
  • Inequality constraints: ADMM only, with a feasibility pre-solve and a stop at violation 1e-4 of the field range (SC-004). No active-set or hybrid method.
  • Single precision: 0.56 times the memory of double precision. The final residual is recomputed in double precision.
  • Browser: the module is 92.5 KB after gzip (limit 2 MB). The Worker reports progress, can cancel, refuses a second solve, and refuses a solve over a memory budget.

What changes in existing code

  • packages/loop_interpolation: new _native.py (adapter), a native branch in solve_system, a branch in setup_interpolator that skips the sparse rows when the solver is native, and "native" accepted by resolve_solver_choice. New optional extra native.
  • Release and workspace files: release-please entries for the new package (and exclude-paths for the root package), the uv workspace excludes the Rust package, .gitignore entries.
  • .github/workflows/native.yml: new CI with path filters (Rust, wheel on three systems, WebAssembly, and a manual reference job).
  • docs/source/user_guide/native_solver.rst: user guide page.

Decisions to review

  • The Rust numerics follow the Python matrix solver where loop-rs differed: curvature rows only at nodes that are interior on all three axes, centred second-derivative border rows, and the horizontal strike/dip basis for gradient data. Reason: the native field must equal the matrix solver (FR-025, SC-003), and existing Python results must not change (FR-019). See "Answers" in plan.md.
  • The native adapter ignores column scaling (apply_scaling_matrix), because it is on by default and native cannot honour it. It raises NotImplementedError for the other unsupported settings.
  • The Python wheel is built with the parallel feature (about 3 times faster than the matrix solver at 100^3). The browser build is single thread.

Tests

  • Rust: cargo test --release (121 passed, and 122 with --features parallel), cargo clippy --all-targets -- -D warnings clean for no features, parallel, python and wasm.
  • Python: 607 passed with the core installed, 588 passed without it. The existing golden suite is unchanged (143 passed). Native against the FD golden cases: all 7 pass.
  • Browser (Chromium): 8 passed, plus the R1 reference run (RUN_REFERENCE=1).
  • Memory checks: SC-001, SC-002, SC-008, SC-009 (10^3 to 200^3), SC-011 and SC-012 pass. Numbers are in specs/001-low-memory-interpolator/measurements.md.

Not done / known issues

  • SC-007 (R1 in 60 s on the reference laptop) is not met on the development machine: R1 takes 93 to 107 s in Chromium and 53 to 78 s native on one thread. It must be checked on the reference machine. The Firefox and Safari smoke test is manual.
  • The Python ADMM solver does not scale its bound matrix when column scaling is on (the default), so bounds apply to scaled values. This is an existing bug. The native tests set apply_scaling_matrix = False for the Python reference. It is not fixed here.
  • ruff reports 2 old findings in _interpolator_builder.py. This PR does not change that file.
  • The version of the core is 0.1.0 in pyproject.toml and 0.0.0 in Cargo.toml and the release-please manifest. This needs a decision.
  • The CI workflow was checked as YAML only. The runner steps (macOS, Windows, WebAssembly) have not run yet.

🤖 Generated with Claude Code

…tive extra

Wave 1 of spec 001: release-please entries, uv workspace exclude, Cargo and
maturin scaffolds, 20x20x20 reference fixture, native extra, and a failing
test for the native solver choice.
…native solver tests

Wave 2 of spec 001: copy cg, grid, stencil, interpolator and edge-case tests
from loop-rs, accept the native solver name, and add Python tests that skip
until the core is installed.
…y and D8 stubs

Baseline for the 100^3 reference problem R1 is 84.6 bytes per dof.
…nput and memory tracker tests

After table cuts, R1 uses 79.1 bytes per dof (baseline 84.6).
T013 also removes the Clone derive from SolveSession (approved by the user)
because Problem now holds the interpolator.
…lazy grid tables

Constraint record sizes follow the real 144-byte stored row. The scaling test
now keeps the row density constant (approved by the user). The table size test
builds the tables before it measures them.
R1 uses 63.3 bytes per dof through the session (target 75). Also adds the
session_r1 measurement test, makes every session_cg test take the allocator
lock, and implements Field::values_f64 for double precision.
… tests and the fixture parity test

The parity test fails with a difference of 24 percent of the range. T026
stays open until the mapping error is found.
…ure, border and gradient basis

Curvature rows are interior on all three axes, border rows are centred
second derivatives on layers 1 and n-2, and the gradient strike and dip use the
horizontal strike of the Python solver. The fixture parity test now passes.
Four ported unit tests that depended on the old one-sided border are updated.
…y pre-solve

ADMM stops only when the largest violation is at most 1e-4 of the field range
(SC-004). A pre-solve without bounds ends the solve when no bound is active.
Infeasible problems return the 20 worst rows after a stall test.
R3 uses 1.10 times the peak of R1 (limit 1.2). ADMM rebuilds the data
right-hand side each outer iteration and the constraint rows are reserved
before they are added. Adds the inequality scale tests.
…neric over Scalar

The native solver ignores column scaling (plan answer), so the column
scaling row is removed from the unsupported settings test.
The native tests use the same tolerance and iteration limit as the reference
solver. The inequality test uses a well-posed problem and turns column
scaling off for the Python ADMM reference, which does not scale the bounds.
…ive solver guide

Single precision uses 0.56 times the peak of double precision on the 50^3
test. The final residual is recomputed in double precision.
…e memory estimate

The estimate follows the measured peak within 12 percent from 10^3 to 200^3
(SC-009). R1 uses 63.2 MB (SC-001) and 250^3 in single precision uses 456 MB
(SC-008). The scaling test range is 5.5 to 8.5 because the shell scratch grows
with the surface (approved by the user).
…el feature

The wheel is about 3 times faster than the matrix solver on a 100^3 grid
(single thread: about 10 percent slower). The browser build stays single thread.
…ild script

The module is 92.5 KB after gzip (limit 2 MB). The fixture test and the
responsiveness tests pass in Chromium.
…b README

TypeScript is a dev dependency of the browser test package, so the type check
command of the task works.
… the solver comparison

The browser tests pass in Chromium.
The browser result equals the Rust result bit for bit. Memory 65.8 MB.
Set the version to 0.1.0 in the release manifest, Cargo.toml and
Cargo.lock. Add Cargo.toml to the extra files of the package, so that
release-please bumps it with pyproject.toml.

Review finding 9 (R010).
…ssion

When finalize returns an error, store it in the session. Later calls to
step, best_so_far and finish return the stored error. Before, a second
step could report a zero-length field as converged.

The f32 buffer is not reserved before the f64 residual check. That
would raise the f32 peak by 4 bytes per dof and break the 60 % peak
limit. The stored error is the guard for that failure point.

Review finding 7 (R006).
…rted input

- Reject interpolator classes that the core cannot represent. Only
  FiniteDifferenceInterpolator, or a subclass that sets
  supports_native = True in its own body, is accepted (finding 1).
- Reject sparse rows that setup did not make. Setup records the names of
  its rows; any other name, and add_inequality_feature, raise
  NotImplementedError (finding 1).
- A native setup sets _native_setup. A later solve with a matrix solver
  raises ValueError before the field is changed (finding 2).
- Clear solver_history and reset the timing on the native path
  (finding 10).
- Keep the old field when the native solve does not converge. The
  partial field is in latest_native_report["field"] (finding 8).
- Cast the field to float64 and check its length (finding 6).

Review findings 1, 2, 6, 8, 10 (R001-R005).
…ative adapter

- Drop data rows like the matrix path: a non-finite position, value or
  vector, a weight of 0 or NaN, and a gradient or tangent of length 0.
  A negative weight is sent as its absolute value, because least
  squares squares the weight (finding 3).
- Interface family: drop NaN ids and non-finite positions, then groups
  with fewer than two points. A negative ipw is used as abs(ipw)
  (finding 4).
- Inequality pairs: the new setup keyword inequality_pair_max_rows
  (default 100 000 for each pair) caps the rows. Above the cap, the
  adapter samples the pairs at random with a fixed seed
  (inequality_pair_seed, default 0) and logs a warning. Rows are built
  in chunks of 100 000 with index arithmetic, never the full cross
  product. Rows with a non-finite position are dropped (finding 5).

Review findings 3, 4, 5 (R007-R009).
…o solve_system

solve_system reads ridge_factor, add_ridge_regularisation, warm_start and
x0 from solver_kwargs. cg, lsmr and admm accept a start value. The start
value is read before _pre_solve clears the previous field. With column
scaling the start value is divided by the scaling diagonal.
InterpolatorBuilder.solve has ridge_factor and warm_start parameters.

The native path still receives the raw kwargs, so check_supported
rejects these keys until the core supports them.
…oblem

Problem.set_ridge_factor sets the ridge term of the normal operator. 0
removes it. Problem.set_initial_guess sets the start value of the solve.
CG starts with r = rhs - N x0, and the ADMM pre-solve starts from x0.
The residual stays relative to |rhs|. The memory estimate adds 8 bytes
for each dof when a start value is set.

The Python and wasm bindings expose both methods.
…tive solver

solve_system gives x0 (from x0 or warm_start) and the effective ridge
factor to _native.solve, which sends them to the core. The native
solver accepts ridge_factor, add_ridge_regularisation, warm_start and
x0. The user guide lists the new options.
Use as_chunks for fixed-size chunks (chunks_exact_to_as_chunks), and put
getter_with_clone only on the String and Vec fields of StepResult
(clone_on_copy). The JavaScript API does not change.
The core is not a workspace member, so the extra loop-interpolation[native]
had no source and uv could not resolve the workspace. A path source lets
uv read its static metadata. The core is built only when the extra is used.
…ields

macOS arm64 differs by about 1e-6 from the golden field because of the
summation order. The solver tolerance does not resolve 1e-6.
A fork pull request gets a read-only GITHUB_TOKEN, so the push of the
lint/style-fixes branch failed with 403. Push and master runs are not
changed.
…ance

The parallel sums of the core do not add in a fixed order, so two runs
can differ in the last bits (SC-003). On macOS the difference was 1.4e-10
and the bit-equal check failed. The sampled rows are still compared bit
for bit.
…s stencils

The row of the normal operator at a node depends only on the class of
each coordinate (first, second, third, interior, and the mirrored three),
so 343 kernels cover the curvature and border families. Each kernel is
built from the family stencils and weights. apply_normal_operator_into
does one sweep: edge nodes with their class kernel, interior runs as
contiguous multiply-add passes with the 25-tap kernel, then the data
rows. This replaces the fused-interior, shell-correction and border
passes and their scratch buffers. Grids with an axis of fewer than 7
nodes, and directional regularisation, keep the old path. The unfused
path stays the reference in the tests.

One apply on 100^3: 46.4 ms -> 7.7 ms on one thread, 40.8 ms -> 8.1 ms
with rayon. The peak is 40 bytes per dof (was about 53), and the memory
estimate has a class-stencil model (ratio 1.09 to 1.14 of the peak).
The f32 residual check needs no rhs vector on this path.

Rename the start value API to the spec names: Problem::set_start,
CgState::restart_from, Python set_start, wasm setStart. Add tests: the
class stencil equals the reference to 2e-15 on 9 grid shapes, a start
at the solution needs 0 iterations, and nested 10^3 -> 20^3 equals a
cold 20^3 solve within SC-003.
…MM solves

ADMM stops at its tolerance, so a different summation order can change
the field by about that tolerance. On macOS two runs differed by 2.5e-8,
above the 1e-9 limit of the previous commit. SC-003 allows 1e-6 of the
field range.
core-wheels.yml builds abi3 wheels for Linux x86_64 and aarch64, macOS
arm64 and x86_64 and Windows x64, and an sdist, with maturin-action. Each
wheel that the runner can run is installed and solves a small problem.
It runs on pull requests that change the crate, and pypi.yml calls it.

pypi.yml uploads the core with PyPI trusted publishing. No other package
waits for it. The core pyproject gets authors, licence, classifiers and
project URLs.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant