Codebase audit: fix wrong-result, correctness, cleanup and performance findings

[pastewka]

This PR works through a full audit of the codebase, tier by tier: high-severity wrong-result bugs, medium-severity correctness and metadata issues, low-severity cleanups, and the performance findings. All fixes include regression tests. What was deliberately left unchanged, and why, is recorded in the sections below.

Fixes #282 — the +q/-q folding in the power spectrum halved the q = 0 entry, which receives no mirrored contribution. fold_fft_half now halves only [1:n//2], i.e. exactly the C_all[1:] /= 2 the issue asks for. The Nyquist entry is structurally excluded, since the folded output stops at n//2.

1. High severity — silent data corruption / wrong results

Geometry and derivatives

• TransposedUniformTopography now swaps pixel_size, nb_subdomain_grid_pts and subdomain_locations
• DownsampledUniformTopography reports correct pixel_size and area_per_pt
• DecoratedUniformTopography.pixel_size derives from the decorator's own geometry instead of the parent's
• Derivative normalization infers the operator direction from the stencil's nonzero coefficients, so anisotropic pixels scale the right axis

File I/O readers

• Reshape order fixed for non-square scans in the DI, EZD, MI and PS readers
• MI reader converts lateral sizes from meters to the requested unit
• BCR reader interprets voidpixels as a count and masks against format markers
• Channel selection fixed in the NMM and OIR readers (late-binding closures, index tracking)
• EZD magic check fixed (bytes vs. str, inverted logic)
• ZAG container reader re-opens the ZIP per read

C++ extensions

• _TriangleMoment and _LineScanMoment rewritten in exact polynomial form (complete homogeneous symmetric polynomials), removing sign errors, missing prefactors and division-by-zero crashes
• Bicubic.__call__ enforces C-contiguity and respects strides; scalar derivative returns match the requested order

Statistics and analysis

• Periodic 2D bearing-area bounds use explicit axes in np.roll
• to_uniform() interpolates at the correct x-positions for scans whose origin is not zero
• _slope_detrend_coeffs uses a length-weighted slope and the matching offset
• Hann window evaluation in the nonuniform power spectrum is origin-invariant
• Scale-dependent curvature stencil handles log-spaced and reliability-trimmed ACF grids
• fourier_synthesis zeroes the DC mode instead of setting it to C(q=1)
• scale_dependent_statistical_property distance calculation fixed; the container variant no longer shadows its outer loop variable and accepts a scalar distances
• make_sphere(..., periodic=True) passes the periodic flag through
• Scan-line alignment direction semantics corrected (transpose for 'x', not 'y')
• Nonuniform.moment() handles equal adjacent heights analytically
• interp_derxy was assigned interp_derxx in the bicubic interpolator
• Fourier interpolation only splits the Nyquist entry along directions that are actually enlarged
• Linear interpolation uses floor rather than truncation, so negative positions on periodic data interpolate instead of extrapolating

2. Medium severity

• Undefined (masked) data is normalized by the number of defined points throughout, instead of by the full grid size. This changed the reference values of several IO tests whose fixtures contain undefined data — the old values were dominated by a mean-normalization artifact, and each updated value carries a comment saying so.
• Analysis-function registration is copy-on-write per class with MRO-merged lookup, so registering on a subclass no longer mutates the base class's table.
• Correctness under non-default options: reliability cutoffs, resampling methods and window arguments.
• Metadata and unit handling: unit-less channels no longer fail info validation.

3. Low severity

Reader edge cases (AL3D row stride, WSXM and DZI flat-scan guards, ZON strides and a non-float32 rejection, Text channel bounds, JPK IndexError, FRT explicit byte order), numerical conditioning (the nonuniform polyfit normal equations are solved in scan-centered coordinates, so detrending no longer depends on where the scan sits on the x-axis), a corrected Nyquist frequency for anisotropic pixels, stacklevel on deprecation warnings, and a batch of docstring corrections (two-sided integrals, polyfit2 coefficient order, moment prefactors, derivative trimming semantics).

Two items were left deliberately unchanged: the NMS height-scale divisor (2**16 - 2 vs 2**16 - 1), because no reference implementation or fixture can decide it, and the latent non-square axis pairing in DATX/LEXT, which cannot be proven without a non-square fixture.

4. Architecture

test/test_invariants.py adds property-based invariant tests: transpose/translate/scale invariance of scalar parameters with anisotropic pixels, decorator geometry consistency (pixel_size == physical_sizes / nb_grid_pts for every decorator), bearing-area CDF properties, and origin invariance of nonuniform detrending.

5. Performance, under the library's memory constraints

The performance findings were re-assessed against the memory boundary conditions of the library: a container streams its topographies from disk and must never hold more than one of them in memory, and readers are long-lived metadata handles that must not pin file data. Several of the apparent performance penalties exist precisely because of those constraints, so these changes either reduce memory use or leave it unchanged.

Fixed, and memory use goes down:

• Bicubic spline coefficients are computed per cell with a single-cell cache instead of being tabulated for every pixel, which cost 128 bytes/pixel before the first evaluation. A 16384² map now constructs and evaluates in 2.1 GB peak RSS; the eager table alone would have needed 34.4 GB.
• The IBW, MI and MNT readers no longer retain wave data, the whole file, or all decompressed blocks past __init__.
• The flood-fill stack is allocated once per call instead of 16 MB per patch.
• Gaussian process regression reuses one Cholesky factorization and computes only the variance diagonal instead of the full predictive covariance.

Fixed, memory-neutral:

• Container PSD integration precomputes each topography's bandwidth interval in a single streaming pass (two scalars each; heights are never cached) and answers counting queries with searchsorted, turning N² file reads into 2N as the pre-existing TODO proposed.
• The MNT zlib scan measures candidate streams through a chunked decompressobj over a memoryview and skips past accepted streams, instead of slicing the remaining buffer and fully decompressing it at every candidate byte.
• The nonuniform height-height autocorrelation moved into a C++ kernel that runs in constant additional memory, rather than being vectorized into an O(N·M) intermediate.
• The rigid-sphere scan is vectorized over blocks sized so the temporary stays bounded irrespective of scan length and tip radius.

Deliberately unchanged, because the speedup would cost memory: the correlation function's pixel sweep (its FFT alternative needs several full-map complex arrays, ~1 GB at 4096², versus O(1) today) and the bearing-area full-grid sweep per query height (a sorted cache adds another O(N) array and only pays off when many heights are queried).

Related issues not closed by this PR

• Interpolate bicubic throws segmentation fault on surface topographies of size 32768 x 32768 #318 (bicubic crash on very large maps) is materially improved — the 128 bytes/pixel eager table is gone, so a 16384² map now works in 2.1 GB — but not fixed: n1_ * n2_ is still an int × int product in cpp/bicubic.h, which overflows at and beyond 46341², including the ~65000² case in the issue body. That needs a wider index type throughout and is out of scope here.
• XYZ reader should not parse file twice #339 (XYZ reader parses the file twice) is intentionally not fixed. Caching the parsed columns on the reader would pin the whole point cloud for the reader's lifetime, which is exactly what the lazy-reader rule forbids; only the per-token allocation cost was reduced. Closing it would require deciding whether the Celery-worker argument in the issue outweighs the memory rule.

Testing

Full suite: 2375 passed, 6 skipped locally. flake8 SurfaceTopography test is clean. Every numerical change was checked against the previous implementation directly, not only through the test suite — the new ACF kernel matches the old Python loop to 6e-17, the chunked rigid-sphere scan is bit-identical across four tip radii on uniform and nonuniform scans, GP mean and variance match to 1e-11 and 1e-15, and the MI reader returns identical heights via path, open stream, repeated calls and a stream handed over at a nonzero offset.

The three tests that require downloading from contact.engineering are excluded locally (sandbox has no network access to that host) and are covered by CI.

https://claude.ai/code/session_01TPEpAGDASUXBVpV5BVDnUZ