remove_trend_linear: replace SVD lstsq with an equivalent separable plane fit (faster; fixes non-square grids)

[brmather]

Summary

CurieGrid.remove_trend_linear fits the best linear trend with a full
np.linalg.lstsq (SVD) over an (N**2, 3) design matrix. For a regular grid
this can be replaced by a separable closed-form plane fit that is numerically
identical to the SVD result (to ~1e-12 on square windows), roughly an order of
magnitude faster, and which additionally fixes a latent bug on non-square grids.

Current implementation (pycurious/grid.py)

def remove_trend_linear(self, data):
    nr, nc = data.shape
    yq, xq = np.mgrid[0:nc, 0:nr]
    A = np.c_[xq.ravel(), yq.ravel(), np.ones(xq.size)]
    c, resid, rank, sigma = np.linalg.lstsq(A, data.ravel(), rcond=None)
    return data - np.dot(A, c).reshape(data.shape)

Problem 1 — performance

remove_trend_linear is the usual process_subgrid callback, so it runs once
per window whenever a spectrum is computed. The SVD over the (N**2, 3) matrix
dominates: on a 2001×2001 window it is ~350 ms, roughly half the cost of the
entire window_spectrum call. In a workload that fits many windows (a global
mesh across several window sizes) this becomes a large fraction of total runtime.

Problem 2 — correctness on non-square grids

np.mgrid[0:nc, 0:nr] produces index arrays of shape (nc, nr), but data has
shape (nr, nc). When nr != nc the raveled indices no longer line up with
data.ravel(), so the fitted plane is effectively transposed and the trend is
not removed. On a 40×25 pure plane:

base() residual on a pure plane = 3.361e+00     # expected ~0

It is correct when nr == nc (the two ravel orders coincide), which is why it
works fine for the usual square windows and has gone unnoticed.

Proposed implementation

Over a regular grid the centred row and column indices are mutually orthogonal
and both orthogonal to the constant, so the normal equations decouple: the three
plane coefficients reduce to one mean and two 1-D inner products — no design
matrix, no SVD.

def remove_trend_linear(self, data):
    nr, nc = data.shape
    i = np.arange(nr) - (nr - 1) / 2.0      # centred row index
    j = np.arange(nc) - (nc - 1) / 2.0      # centred column index
    mean = data.mean()
    ci = (i * (data.mean(axis=1) - mean)).sum() / (i * i).sum()
    cj = (j * (data.mean(axis=0) - mean)).sum() / (j * j).sum()
    return data - (mean + ci * i[:, None] + cj * j[None, :])

(Undefined only for a 1×N or N×1 grid, where a trend along the singleton axis
is not identifiable anyway — trivial to guard if wanted.)

Verification

• Square grids (200², 411², 2001²): matches the current SVD fit to ~2e-12
  relative; removes a pure plane exactly.
• Non-square grids: removes a pure plane exactly (0.0) where the current
  version leaves a finite residual (see above).
• Timing (2001², single-threaded): ~350 ms → ~49 ms; the gap widens further
  with multithreaded BLAS.

Happy to open a PR with this change plus a small regression test: pure-plane
removal on square and non-square inputs, and equivalence to the SVD fit on
random square data.