[roytsmart]

method="multilinear" used to work only along one axis. It now interpolates on logically-rectangular grids of any number of dimensions, rectilinear or curvilinear, at output points arranged along any number of axes. This is the first step toward evaluating a sampled na.FunctionArray (such as a map of a mirror's surface error) at scattered points.

What changes

• Any number of dimensions. Each output point gets the $2^N$ corners of the input cell which contains it, weighted by its local coordinates in that cell.
• Output points on any axes. For this method, its outputs are points rather than cells, so an explicit axis_output may have any number of axes. A 2D grid can be interpolated at one axis of scattered points (axis_output=0), or at a single point (axis_output=()). Without axis_output, the points must still have as many axes as the grid, so a forgotten axis_output raises instead of turning the orthogonal axes into points. The conservative method still requires the two to match.
• Rectilinear grids are searched along each axis by binary search, which is exact and needs no setup. A grid counts as rectilinear if each coordinate varies along one axis only (a different axis for each) and is strictly monotonic along it. The check uses strides first, so a grid broadcast from 1D coordinates, as named-arrays passes them, is recognized without reading its values. Axes which run backwards work, including unsigned integer ones. A 1D grid which ran backwards used to give wrong answers silently.
• Curvilinear grids are searched through a hierarchy of bounding boxes. Each level groups the level below in blocks of two along every axis. A point descends through the boxes that contain it, and in each candidate cell Newton's method finds its local coordinates. The point is in the first cell where they fall inside the unit cube. The search takes about log(N) steps and doesn't depend on the grid being convex. Newton's method has converged once a step is short, or once the distance to the point is down to the rounding error of the position. The second matters in thin cells at an angle to the axes, where rounding along a cell's length moves its local coordinate across its width.
• Many elements at once. Every element along the orthogonal axes is searched in one compiled call. That includes rectilinear grids which differ from element to element, such as a spectrum for each pixel, and the points of every element on a curvilinear grid they share.
• bounds. On a rectilinear grid, "extrapolate" behaves as before, along each axis. On a curvilinear grid, the same Newton routine runs across the whole grid in index space, with the edge cells continued beyond it, starting from the cell whose box is nearest. On a grid of parallelograms that is the same as clamping each index. A point the continuation can't reach (beyond a strongly curved edge) gets NaN, as does a point with a non-finite coordinate.
• regrid_from_weights flattens the resampled axes of each element, giving the number of elements explicitly. That lets through outputs with no resampled axes (a single point), no points, or no elements. The CUDA apply path is untouched; multilinear weights are still built on the host only.

Why not index_of_point_secant

The conservative sweep already has a 2D Newton locator, so I tested it against index_of_point_brute first (300 random interior points each):

Grid	Wrong	Slow
uniform 101×101, spacing 1e-12	300 of 300 reported outside
annulus, 0.9 of a turn	109 of 300 reported outside
disk-like annulus, r from 0.05 to 1	117 of 300 reported outside
cell widths alternating 0.4 / 1	0	20 fall back to brute force, about 100× slower

The causes:

• its convergence test is an absolute 1e-10 in coordinate units;
• on an annulus it follows the continuation of a cell across the hole to a root outside the grid;
• on unequal neighboring cells it cycles until it gives up.

The first two faults also reached weights(..., method="conservative") through the sweep's row seeds; #74 fixes that separately. Here the boxes do the global search and Newton's method only runs inside one cell, where it is well behaved. Its tests are in local coordinates, and positions are measured from the cell's first corner. That means a grid's scale, and its distance from the origin, don't matter.

Tests

• Rectilinear grids in 1D, 2D and 3D, one axis reversed, against scipy.interpolate.RegularGridInterpolator, with and without extrapolation. They agree to 1e-14, and the NaNs fall in the same places. The same for 50 elements on grids of their own, some reversed.
• Curvilinear grids (two annuli, a distorted grid at scales of 1e-12 and 1e9, cells of alternating width, thin rotated cells of aspect ratio 1e4 and 1e6 in 2D and 1e4 in 3D, a distorted 3D grid, a 1D grid that folds back, a single cell): which points are inside matches brute force over every cell. The brute force shares no code with this change. In 2D it uses point_is_inside_polygon. In other dimensions it uses scipy.optimize.root on each cell, judged by its residual.
• Any affine function of the coordinates is reproduced exactly on any grid, inside and extrapolated, since the corner weights sum to one and average the corners to the point.
• Output points along one axis, along two axes from a 3D grid, a single point, no points, and no elements. A forgotten axis_output raises.
• Every vertex of a grid, boundary included, under bounds="raise"; grids that vary along an orthogonal axis, and grids broadcast across one; weights_input, shared or per element; units; non-finite points; cells of zero size; Newton's halved steps across cells of very different widths.
• Putting back each of the four bugs the review found (below) makes these tests fail.
• The existing 1D tests pass unchanged. The new package has full coverage with NUMBA_DISABLE_JIT=1.

Speed

On a 1001×1000 grid, building the weights for scattered points with coalesce=False (the new asv benchmark time_weights_multilinear_2d, best of three):

Grid	1 point	10⁴ points	10⁶ points
rectilinear	3 ms	3 ms	52 ms
curvilinear (rotated and distorted)	23 ms	22 ms	213 ms

The single-point cost is the setup: comparing a materialized meshgrid's values, or building the boxes. For a rectilinear grid broadcast from 1D coordinates, the setup takes 0.4 ms.

With many elements, a grid of 50 points for each pixel and 40 points to interpolate at, weights() takes about as long as on main, where most of the time goes to merging repeated pairs:

Pixels	main	this branch
10⁴, each its own grid	0.39 s	0.35 s
10⁵, each its own grid	3.6 s	3.3 s
10⁴, one grid tiled	0.39 s	0.36 s

Review

An xhigh review of the first version found ten issues, all fixed here:

• points in thin rotated cells reported outside;
• a forgotten axis_output returning a wrongly shaped result instead of raising;
• descending unsigned grids taken as ascending;
• a slowdown of about 5× for many orthogonal elements;
• a crash with no output points;
• a 3D test oracle which reused the code under test, and no thin cells or boundary points in the tests;
• a stale regrid_from_weights docstring;
• a Newton loop copied between the in-cell search and the extrapolation walk;
• axis bookkeeping duplicated from the conservative path;
• double copies of the coordinates.

find_indices is unchanged. Its brute and searchsorted methods could use these locators to work in more than 1D in a later PR.

🤖 Generated with Claude Code

https://claude.ai/code/session_01BYjDL98znSud1yFh9chnkP

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.82%. Comparing base (e3acd29) to head (74934fe).

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #73      +/-   ##
==========================================
+ Coverage   99.72%   99.82%   +0.09%     
==========================================
  Files          65       72       +7     
  Lines        4346     5612    +1266     
==========================================
+ Hits         4334     5602    +1268     
+ Misses         12       10       -2     

Flag	Coverage Δ
unittests	99.82% <100.00%> (+0.09%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.