quadsv.kernels.nufft ==================== .. py:module:: quadsv.kernels.nufft .. autoapi-nested-parse:: Non-uniform FFT (NUFFT) spectra, kernel and spatial tests for irregular data When data sit on a regular grid (e.g., a rasterized Visium slide), :func:`quadsv.power_spectrum_2d` computes :math:`|\hat{x}(k)|^2` with a plain 2D FFT. For data whose spatial coordinates are **irregular** — e.g., imaging-based in situ platforms, Slide-seq, or a Visium slide read straight from ``adata.obsm['spatial']`` without rasterization — :func:`power_spectrum_2d_nufft` evaluates the type-1 NUFFT .. math:: \hat c(k_y, k_x) = \sum_{j=1}^{n} c_j \, \exp\!\bigl[-i(k_y\,y_j + k_x\,x_j)\bigr] on the same uniform ``(ny, nx)`` k-space grid that :func:`power_spectrum_2d` would produce for a rasterized input of the same physical extent, and returns :math:`|\hat c|^2` in the scipy FFT layout (DC at ``[0, 0]``). Anything downstream — :func:`quadsv.comparators.multisample.radial_bin_spectrum`, :class:`quadsv.ComparatorIrregular` — works identically. Notation (shared across this module) ------------------------------------ Dimensions: - ``n``: number of spots (on the irregular grid). - ``(ny, nx)``: internal uniform k-grid dimensions; ``n' = ny · nx``. - ``(dy, dx)``: **physical** spacing per k-grid cell, same unit as the spatial coordinates after multiplying ``unit_scale``. - ``unit_scale``: multiplier that converts the input coordinates ``S`` to the same unit as ``(dy, dx)`` (e.g., 0.35 if ``S`` are in pixels at 0.35 μm/pixel). Samples from different slides and platforms may ship coordinates in different units; this parameter harmonizes them onto the same **physical** unit for the internal k-grid and all downstream spectra and tests. Vectors and matrices: - ``S``: the ``n × 2`` spatial coordinate matrix of the irregular points, ordered as ``(y, x)``. - ``K``: the ``n × n`` translation-invariant kernel at the irregular points. - ``K'``: the ``n' × n'`` grid kernel with FFT eigenvalues ``λ(k) = F(K')(k)``. - ``U``: the ``n × n'`` type-2 NUFFT evaluation matrix; the band-limited approximation is ``K ≈ (1/n') · U · diag(λ) · Uᴴ``. - ``x̂ = Uᴴ x``: type-1 NUFFT of a length-``N`` signal onto the k-grid ``(ny, nx)`` (vectorized). Classes ------- .. autoapisummary:: quadsv.kernels.nufft.NUFFTKernel Functions --------- .. autoapisummary:: quadsv.kernels.nufft.power_spectrum_2d_nufft Module Contents --------------- .. py:class:: NUFFTKernel(coords, grid_shape = None, spacing = None, method = 'matern', unit_scale = 1.0, oversample = 2.0, eps = 1e-06, workers = None, *, centering = True, **kwargs) Bases: :py:obj:`quadsv.kernels.base.Kernel` Spatial kernel over **irregular** 2D coordinates evaluated via NUFFTs. Parallels :class:`quadsv.fft.FFTKernel` (which requires a regular grid) and implements the :class:`~quadsv.kernels.Kernel` interface so it plugs into :func:`quadsv.statistics.spatial_q_test` / :func:`quadsv.statistics.spatial_r_test` the same way. The band-limited approximation of the ``n × n`` irregular-point operator is ``K ≈ (1/n') · U · diag(λ) · Uᴴ``, where ``U`` is the ``n × n'`` type-2 NUFFT matrix and ``λ = F(K')`` is the grid kernel's spectrum. Under this approximation, Parseval's identity gives the fast quadratic form ``xᵀ K x = (1/n') Σ_k λ(k) |x̂(k)|²`` with ``x̂ = Uᴴ x`` (a single type-1 NUFFT). The matrix-vector primitive :meth:`Kx` uses the companion two-shot NUFFT ``K z = (1/n') · U · (λ ⊙ Uᴴ z)`` and backs the Hutchinson-based cumulant estimator (:func:`quadsv.statistics._hutchinson_cumulants`) and the bipartite R-test cross matrix. :func:`quadsv.spatial_q_test` always uses the k-space Parseval path (:meth:`xtKx`); :func:`quadsv.spatial_r_test` dispatches on shape — paired diagonal for ``M_x == M_y`` via :meth:`xtKy`, full ``(M_x, M_y)`` cross matrix via :meth:`Kx` otherwise. The matmul counterparts (:meth:`xtKx_matmul` / :meth:`xtKy_matmul`) are exposed for callers that want to compute the same form directly at the ``n`` irregular points — they agree with the spectral path to NUFFT precision (``eps``) on both regular and irregular coordinates. :param coords: Spot coordinates of shape ``(n, 2)`` in order ``(y, x)``. :type coords: np.ndarray :param grid_shape: ``(ny, nx)`` of the internal uniform k-grid. If ``None`` (default), auto-inferred from ``coords``: the grid is sized to cover the bounding box and to resolve the sampling Nyquist set by the median nearest-neighbor distance (fully coordinate-driven, kernel-agnostic). Override only when you know you need a finer or coarser grid. :type grid_shape: tuple[int, int], optional :param spacing: ``(dy, dx)`` physical spacing per k-grid cell (same unit as ``coords`` after ``unit_scale``). If ``None`` (default), auto-inferred alongside ``grid_shape``. When both are supplied, users are responsible for ensuring ``ny · dy``, ``nx · dx`` covers the coordinate extent. :type spacing: tuple[float, float], optional :param method: Kernel method forwarded to :class:`FFTKernel`. One of ``'gaussian'``, ``'matern'``, ``'moran'``, ``'graph_laplacian'``, ``'car'``. :type method: str, default ``'matern'`` :param unit_scale: Multiplier applied to ``coords`` so they share the same unit as ``spacing`` (e.g. ``0.35`` if coords are in pixels at 0.35 μm/pixel). :type unit_scale: float, default 1.0 :param oversample: Auto-grid oversampling factor above the sampling Nyquist. Used only when ``grid_shape`` / ``spacing`` are auto-derived. Larger values give a finer k-grid (more accurate, slower); 2.0 is safe for all tested kernels. :type oversample: float, default 2.0 :param eps: NUFFT tolerance forwarded to finufft. :type eps: float, default 1e-6 :param workers: Forwarded to :mod:`scipy.fft` (used by :meth:`Kx_grid`) and reserved for future finufft parallelism. ``None`` uses the SciPy default. :type workers: int, optional :param \*\*kwargs: Method-specific kernel hyperparameters forwarded to the internal :class:`FFTKernel`. ``bandwidth`` / ``nu`` for ``gaussian`` / ``matern``; ``rho`` for ``car``; plus a *coord-aware* ``k_neighbors`` for the graph methods (``moran`` / ``graph_laplacian`` / ``car``). Note the ``k_neighbors``-to-``neighbor_degree`` mapping differs from that of :class:`FFTKernel` due to an oversampled internal grid. ``neighbor_degree`` is chosen so the internal-grid-ring cutoff matches the median k-th nearest-neighbor distance among the coords — matching :class:`~quadsv.kernels.MatrixKernel`'s mutual-k-NN graph semantic up to the band-limit of the internal grid. See :func:`_resolve_k_neighbors_on_coords` for the mapping detail. Pass ``neighbor_degree`` directly to bypass this and use the raw grid-ring semantic as in :class:`FFTKernel`. :ivar coords: Original ``(n, 2)`` coordinates. :vartype coords: np.ndarray :ivar n: Number of spots ``n``. :vartype n: int :ivar grid_shape: Internal k-grid shape ``(ny, nx)``. :vartype grid_shape: tuple[int, int] :ivar spacing: Physical spacing per k-grid cell ``(dy, dx)``. :vartype spacing: tuple[float, float] :ivar method: Kernel method name. :vartype method: str :ivar params: Resolved kernel hyperparameters (snapshot of the internal FFT kernel). :vartype params: dict :ivar workers: scipy.fft worker count used by :meth:`Kx_grid`. :vartype workers: int or None :ivar stores_precision: Always ``False`` — NUFFTKernel never holds an ``n × n`` matrix. :vartype stores_precision: bool .. py:method:: Kx(z) Matrix–vector product ``K z`` at the ``n`` irregular coordinates. Implements the band-limited apply .. math:: K z \;\approx\; \tfrac{1}{n'} \, U \bigl(\lambda \odot U^{\mathsf H} z\bigr), evaluated as type-1 NUFFT → elementwise multiply by ``λ(k) / n'`` → type-2 NUFFT. Output length ``n``, same shape as ``z``. Base primitive for :meth:`xtKx_matmul`, :meth:`xtKy_matmul`, the Hutchinson cumulant estimator used by Liu's null approximation, and the bipartite R-test in :class:`DetectorIrregular`. :param z: ``(n,)`` or ``(n, M)``. :type z: np.ndarray :returns: Same shape as ``z``. :rtype: np.ndarray .. py:method:: Kx_grid(x) Grid-domain companion of :meth:`Kx` — ``(ny, nx)`` spatial output. Whereas :meth:`Kx` returns the length-``n`` apply at the irregular coordinates, :meth:`Kx_grid` returns the apply evaluated on the internal uniform grid. Pipeline: type-1 NUFFT → undo the coordinate-centering phase (needed here because we keep complex coefficients; the square-magnitude and adjoint-round-trip paths of :meth:`xtKx` / :meth:`Kx` absorb it automatically) → multiply by ``λ(k)`` → ``ifftshift`` → ``ifft2`` → real. :param x: ``(n,)`` or ``(n, M)``. :type x: np.ndarray :returns: Real ``(ny, nx)`` or ``(ny, nx, M)`` in the scipy FFT layout (DC at ``[0, 0]``). :rtype: np.ndarray .. py:method:: eigenvalues(k = None, return_full_layout = False) Eigenvalues of the ``n × n`` irregular-point operator. Returns the spectrum of the realized operator ``K_n`` (raw) or ``H K_n H`` (centered) at the irregular coordinates — **not** the internal FFT-grid spectrum. Purely matrix-free: no dense ``n × n`` construction. Two paths: - **``k`` given — Lanczos top-k.** Wrap :meth:`Kx` as a :class:`~scipy.sparse.linalg.LinearOperator` and call :func:`~scipy.sparse.linalg.eigsh(which='LM')` for the top ``k`` eigenvalues by magnitude. Cost: ``O(k × nufft)``. Not cached. - **``k=None`` — Toeplitz-M (analytic).** See :meth:`_eigvals_toeplitz_M`. Applies when ``Λ`` is PSD and the support size ``r = #{k : |λ(k)| > 10⁻⁵ · max|λ|}`` is below ``_TOEPLITZ_R_THRESHOLD`` (default 2000). Exact to NUFFT ``eps``; cost dominated by ``O(r³)`` eigvalsh on a dense ``r × r`` reduced matrix. Cached per ``centering`` mode. Full-spectrum requests that fall outside Toeplitz-M's reach (indefinite ``Λ`` like Moran, or broad-support kernels like CAR at strong coupling) raise ``NotImplementedError`` rather than falling back to an approximate density reconstruction. For Liu's method, use :func:`compute_null_params(..., method='liu', liu_n_probes=...)` to get cumulant-based Liu directly from :math:`2m` matvecs. :param k: Top-``k`` eigenvalues by magnitude (Lanczos). ``None`` returns the full spectrum via Toeplitz-M. :type k: int, optional :param return_full_layout: Kept for API compatibility with :meth:`FFTKernel.eigenvalues`; ignored. :type return_full_layout: bool, default False :returns: Eigenvalues in descending order. Length ``k`` when ``k`` is given, ``n`` for the full-spectrum path. :rtype: np.ndarray .. py:method:: square_trace() ``trace(K²)`` (raw) or ``trace((HKH)²)`` (centered). Closed-form ``(n²/n'²) · λᵀ Ψ λ`` with Toeplitz ``Ψ_{k,k'} = |φ(k'-k)|²``, evaluated as a *linear* (non-circular) 2D convolution of ``|φ|²`` with ``λ`` via a doubled-grid FFT in ``O(n' log n')``. ``φ(j) = (1/n) Σ_i exp(-ij·y_i)`` is evaluated on a ``(2·ny, 2·nx)`` mode grid by a separate type-1 NUFFT of the ones vector (see :meth:`_coord_power_spectrum_doubled`), and ``λ`` is zero-padded to the same doubled layout so the FFT convolution does not wrap values of ``|φ|²`` across the ``n'``-period — a silent bias of up to ~45% on broad-spectrum kernels (CAR, graph_laplacian) on the typical oversampled NUFFT grid. On a regular grid where coords coincide with the k-grid, ``|φ|² = δ`` and the formula collapses to ``(n/n')² · Σ_k λ(k)²``. Adjusts by ``-2·s₂/n + s₁²/n²`` when ``centering=True``. Observed band-limit residuals vs. explicit ``Kx(I)`` truth: ≲ ``1e-7`` on Gaussian / Matern, ``~0.1 %`` on CAR, ``~1 %`` on graph_laplacian, and ``~0.05–1.2 %`` on Moran (indefinite ``Λ``) — accurate across regular, irregular, and clustered coord layouts. .. py:method:: trace() ``trace(K)`` (raw) or ``trace(HKH)`` (centered). Closed-form ``(n/n') · Σ_k λ(k)`` — exact because the diagonal of ``G = UᴴU`` is ``n`` regardless of coord arrangement, so ``trace(K_n) = (n/n') · trace(K_grid)`` is independent of the coord layout. Adjusts by ``-s₁/n`` when ``centering=True`` (``s₁ = 𝟏ᵀ K 𝟏`` via a single ``K·𝟏`` apply in :meth:`Kernel._ones_stats`). .. py:method:: xtKx(x) Quadratic form ``xᵀ K x`` via **k-space Parseval**. Implements the default path: .. math:: x^T K x \;=\; \frac{1}{n'} \sum_k \lambda(k) \, |\hat x(k)|^{2} using one type-1 NUFFT of ``x`` and an elementwise Parseval sum. Only the real power spectrum ``|x̂|²`` of shape ``(ny, nx)`` is materialized — no ``ifft2``, no spatial-grid copy of ``x``. :param x: ``(n,)`` for one feature or ``(n, M)`` for ``M`` features. :type x: np.ndarray :returns: Scalar if ``x`` is 1-D, shape ``(M,)`` otherwise. :rtype: float or np.ndarray .. seealso:: :obj:`xtKx_matmul` Compute ``xᵀ · Kx`` via the length-``n`` matrix product. .. py:method:: xtKx_matmul(x) Quadratic form ``xᵀ K x`` via **direct matmul**. Computes ``Q_B = xᵀ · self.Kx(x)`` end-to-end at the ``n`` irregular points. Sparse-aware on ``x`` (``x.multiply(Kx).sum``). ~2× the NUFFT work of :meth:`xtKx` per feature; agrees with it to NUFFT precision on regular grids and to the torus-BC band (~1–2 %) on irregular ones. :param x: ``(n,)`` or ``(n, M)``. :type x: np.ndarray or scipy.sparse matrix :returns: Scalar for 1-D input; ``(M,)`` for batched. :rtype: float or np.ndarray .. py:method:: xtKy(x, y) Bilinear form ``xᵀ K y`` via **cross Parseval**. Implements the default path: .. math:: x^T K y \;=\; \frac{1}{n'} \sum_k \lambda(k) \, \overline{\hat x(k)}\, \hat y(k). Paired same-``M`` convention — returns the diagonal of ``Xᵀ K Y`` (shape ``(M,)``) for batched inputs, scalar for 1-D inputs. For the bipartite ``(M_x, M_y)`` cross matrix use :meth:`xtKy_matmul` (or build it explicitly via ``X.T @ self.Kx(Y)``). :param x: ``(n,)`` or ``(n, M)`` — must share shape. :type x: np.ndarray :param y: ``(n,)`` or ``(n, M)`` — must share shape. :type y: np.ndarray :returns: Scalar for 1-D inputs; ``(M,)`` for batched. :rtype: float or np.ndarray .. py:method:: xtKy_matmul(x, y) Bilinear form ``xᵀ K y`` via **direct matmul**. Returns the paired ``(M,)`` diagonal of ``Xᵀ K Y`` (sparse-aware on ``x``). For the full ``(M_x, M_y)`` bipartite cross matrix build it explicitly as ``X.T @ self.Kx(Y)`` — that's what :class:`DetectorIrregular` does for ``compute_rstat`` when the two feature blocks have different widths. :param x: ``(n,)`` or ``(n, M)``. :type x: np.ndarray or scipy.sparse matrix :param y: ``(n,)`` or ``(n, M)``. :type y: np.ndarray or scipy.sparse matrix :returns: Scalar for 1-D inputs; ``(M,)`` for batched. :rtype: float or np.ndarray .. py:attribute:: coords :type: numpy.ndarray .. py:attribute:: grid_shape :type: tuple[int, int] .. py:attribute:: method :type: str :value: 'matern' .. py:attribute:: n :type: int .. py:attribute:: params :type: dict .. py:attribute:: spacing :type: tuple[float, float] .. py:attribute:: stores_precision :type: bool :value: False .. py:attribute:: workers :type: int | None :value: None .. py:function:: power_spectrum_2d_nufft(coords, values, grid_shape, spacing, unit_scale = 1.0, eps = 1e-06, center_coords = True) Compute the 2D power spectrum via type-1 NUFFT This function computes the power spectrum :math:`P(k) = |\hat{c}(k)|^2` of one or more non-uniform spatial signals via a type-1 NUFFT. The output has the same ``(ny, nx)`` layout as :func:`quadsv.kernels.fft.power_spectrum_2d` with ``fft_solver='fft2'``: DC at ``[0, 0]``, Nyquist at ``[ny/2, nx/2]`` (when dimensions are even). :param coords: Non-uniform spatial coordinates, shape ``(n, 2)`` in the order ``(y, x)``. Values outside the physical domain implied by ``grid_shape`` and ``spacing`` are folded into ``[-π, π)`` by finufft. :type coords: np.ndarray :param values: Signal strengths at each coordinate. Shape ``(n,)`` for a single feature, or ``(n, M)`` for ``M`` stacked features (e.g., genes) on the same coordinates. Real-valued; promoted to complex internally. :type values: np.ndarray :param grid_shape: ``(ny, nx)`` of the target uniform k-space grid. Match whatever grid you use for rasterized samples so the two paths produce comparable spectra. :type grid_shape: tuple[int, int] :param spacing: ``(dy, dx)`` physical spacing per cell of the target grid, in the same unit as ``unit_scale * coords``. Together with ``grid_shape`` this defines the physical domain extent ``(ny · dy, nx · dx)``. :type spacing: tuple[float, float] :param unit_scale: Multiplier applied to ``coords`` before scaling into ``[-π, π)``. Use this to convert per-sample coordinate units into the common unit of ``spacing`` (e.g., 0.35 if ``coords`` are in Visium full-res pixels at 0.35 μm/pixel and ``spacing`` is in μm). :type unit_scale: float, default 1.0 :param eps: NUFFT tolerance forwarded to finufft. :type eps: float, default 1e-6 :param center_coords: If True, subtract the mean of ``coords`` before scaling — avoids wrapping artefacts when coordinates are stored with an arbitrary origin offset (e.g., Visium pixel coordinates start at a few thousand). Power spectra are translation-invariant so recentering does not change the result. :type center_coords: bool, default True :returns: Power spectrum. Shape ``(ny, nx)`` for 1D ``values`` or ``(ny, nx, M)`` for 2D ``values``, with DC at index ``[0, 0]``. :rtype: np.ndarray :raises ImportError: If :mod:`finufft` is not installed. :raises ValueError: If input shapes are inconsistent. .. rubric:: Examples >>> import numpy as np >>> coords = np.random.default_rng(0).uniform(0, 100, size=(500, 2)) >>> vals = np.random.default_rng(1).standard_normal(500) >>> P = power_spectrum_2d_nufft(coords, vals, grid_shape=(32, 32), spacing=(4.0, 4.0)) >>> P.shape (32, 32)