Performance

The design rules, and the measurements behind them.

Rules

1. No superlinear algorithm where a linear one exists.
2. Intrinsic data is O(∑ Nᵈ); anything O(∏ Nᵈ) is a materialization — make it opt-in.
3. Compute nothing twice, across calls as well as within one.
4. One algorithm, one implementation.
5. No transcendental in an inner loop that a recurrence, sincos, or a hoist can remove.
6. Never compute at runtime what is a mathematical constant.
7. An API must not be able to request a terabyte by accident.

Scaling

Every public entry point is linear or better in its problem size. Measured as cost ~ nᵖ across a size step:

Allocation counts are flat in n everywhere. build_connectivity's 91 MiB is the CSR output itself for ~2M nodes, not overhead.

Ball queries: nothing per query that belongs to the grid

A ball query has costs that are properties of the grid, not of the query: the smallest step per direction, which bounds the candidate window, the span of each axis, and — where there are no separable axes to bound with — a spatial index. The first two are reduced once when the grid is built and stored as Grids.AxisStats; Connectivity.MetricTopology then reads them in O(1).

The effect is that the per-query default costs nothing, so there is no hoisting to get right. On a stretched axis, with the radius and the window's cell count held fixed:

Both columns allocate nothing. Before the reductions were cached the default column was a per-query O(N_d) rescan — 15.2 µs at 16 384, growing without bound in the axis length — and every grid with a non-uniform axis paid it, which includes every Gaussian-latitude grid.

A heterogeneous coordinate tuple

Each direction keeps whatever AbstractVector type it was built from, so a grid with a range for longitude and a vector for latitude has a heterogeneous coordinate tuple. Indexing that tuple with a loop variable is a dynamic lookup, and one such loop — a bounds check in _raw_coords — used to defeat inference for every coordinate read:

Every per-cell entry point is now allocation-free on every grid shape, and the suite checks that against a matrix of shapes — including mixed-axis ones, which is the case it did not previously have.

Which index

Two are available. Grids.cell_list needs no package, and on a 65 536-cell curvilinear grid it builds faster, queries faster and allocates less than the k-d tree, so it is what the sweeps build by default:

The cell list also enumerates through a fold rather than a returned list, so a query holds no buffer and runs inside a kernel; a tree cannot, because it searches replicated points and has to deduplicate.

Its build is O(n), but only once the dimension reaches the type: taking it from size(pts, 1) — a runtime value — left the construction loop dynamically dispatched at 196 ms and 34 MiB for those 65 536 points, against 2.75 ms and 3 MiB behind a function barrier.

Indexing the architectures with no axes

Curvilinear and node grids have no window to bound, so a query without an index tests every cell. With Connectivity.indexed (needs NearestNeighbors), one query on a 2-D curvilinear grid, radius fixed at 2.5 cells:

and the whole-grid build, which is n of those queries:

O(n log n) against O(n²). The index returns a superset and the exact distance gate still decides membership, so both columns give the same graph — the index buys speed and changes nothing else. Rows come out in whatever order enumerated them, as everywhere else here; sort_neighbors! if you need them ordered.

Because the index cannot be built per query, Connectivity.foreach_within and Connectivity.mapreduce_within exist to build it once for a whole sweep: 800.7 ms → 88.8 ms (9.0×) against the same sweep written as a hand loop over the per-cell entry points, on 9 216 cells.

Connectivity.ball_scratch supplies the candidate buffer, which takes an indexed query to zero allocation whatever the grid size, against 480 bytes on a small ball and 6.1 KB on a 310-candidate one without it.

Quadrature

Gauss–Legendre holds machine precision past degree 2N−1; Driscoll–Healy and Clenshaw–Curtis lose exactness just after N−1, which is the distinction admits_exact_bandlimited_quadrature encodes. The fitted exponents on the right are measured per run, not asserted.

Gauss–Legendre uses asymptotic expansions (Bogaert 2014; Hale & Townsend 2013) above n = 60: each node sits near j_k/(n+½) for j_k a zero of J₀, with corrections in powers of (n+½)⁻². Nothing iterates and no root depends on its neighbours, so it is O(1) per node and allocation-free.

Accuracy against a 256-bit reference is ~6e-16 relative weight error at every n from 64 to 4096.

Below n = 60, and for element types wider than Float64, the solve falls back to Newton on the Bonnet recurrence: the expansion's coefficients are a fixed Float64 set and cannot exceed that precision, while Newton converges to eps(T) — BigFloat at 256 bits gives |Σw − 2| = 1.7e-77.

Equiangular (DH/CC) weights are O(n log n) with an FFT loaded and O(n²) without, with the two agreeing to 1.4e-14.

Memory

The two big materializations are gone:

• Cell measure — stored as per-axis factors. At 2000²: 61.0 MiB → 0.046 MiB.
• Curvilinear corner directions — two rows live at a time instead of the whole field. At 1000²: 61.3 MiB → 23.1 MiB, with bit-identical areas.

The factored measure is also faster to read, which was not the expectation. At 2000² the dense array is 61 MiB and DRAM-bound while the factors stay in cache:

Threading

Opt-in through ComputationalBackends tags, serial by default, bit-identical results:

using ComputationalBackends: ThreadedBackend
FG.Connectivity.build_connectivity(grid; backend = ThreadedBackend())

Chunks are contiguous, so each thread touches one span of every array rather than striding across all of them. Passing nothing (the default) hands the loop body the whole range in one call, so the serial path adds no partitioning at all.

Two passes are deliberately not threaded, and the reasons are structural rather than pending work. The CSR compacting move can have row j's destination fall inside row i's source block for i < j, so concurrent rows would overwrite unread candidates. And the prefix scan between the count and fill passes is inherently sequential — it is O(n) against the O(n·stencil) passes it separates.

Things measured and found not to be problems

Recorded so they are not re-litigated:

• getproperty coordinate-name lookup is free — 0.065 vs 0.077 ms per 10⁶ reads; it const-folds.
• The dense-mask branch is free in a predictable loop.
• Cross-point SIMD is unavailable — @simd over 10⁶ haversines is 1.01×; LLVM cannot vectorize sin/cos/atan. Eliminating trig is the only lever, so the geometry kernels use sincos and hoist atan out of loops.
• Vector-axis construction is not superlinear — the gap against a uniform axis is a constant factor, not a growing one.
• The ball-query candidate buffer is mostly an allocation win — 1.43× in time on a 20-candidate ball and 1.06× on a 310-candidate one, against zero bytes per query rather than 480 B and 6.1 KB. Worth passing in a sweep; the time it saves is only visible where the ball is small enough that the allocation is most of the query.