Unify hankel transforms by bclevine · Pull Request #40 · LSSTDESC/dsf

bclevine · 2026-06-12T13:47:18Z

Description

This PR addresses #33, combining the Hankel functionality into a single public module. Changes include the following:

Created a new HankelTransformBase parent class for all the Hankel transform methods.
Moved the matrix method currently implemented in the covariance module to HankelTransformMatrixZeros, and the FFTLog method currently implemented in the utils module to HankelTransformFFTLog.
Added a new matrix Hankel method, HankelTransformMatrixDirect, which in some (but not all!) instances can be an improvement over HankelTransformMatrixZeros.
Added a new public-facing HankelTransform class which uses one of the three aforementioned methods as a backend and provides consistent, public access to the major Hankel transform operations.
Added associated tests and benchmarks.

…relation from matrix methods

…at don't cover the full matrix k-grid

nikosarcevic

first of all great work, ben!

i am goign to split this in several chunks as i am afrid ill loose my comments

nikosarcevic · 2026-06-12T15:24:08Z

+            order: Bessel order to use.
+
+        Returns:
+            Radial grid and projected radial statistic.


add Raises section after Returns.

nikosarcevic · 2026-06-12T15:25:05Z

+from dsf.utils.validators import validate_positive_scalar
+
+
+class HankelTransformMatrixDirect(HankelTransformMatrixZeros):


do we need to call it MatrixZeros? WDYT? Perhpas just calling it HankelTransformMatrix woudl be ok?

nikosarcevic · 2026-06-12T15:26:10Z

+        k_max: Maximum wavenumber to cover.
+        n_r: Number of radial grid points to use for the Hankel transform.
+        n_k: Number of wavenumber grid points to use for Hankel transform.
+        orders: Bessel orders to precompute.


sohuld this arg be order instead of orders?

nikosarcevic · 2026-06-12T15:27:41Z

-This module provides the ``HankelTransform`` class, which converts
-Fourier-space spectra into projected radial-space quantities. It is useful for
-computing projected correlation functions, covariance matrices, and higher-order
+This module provides the ``HankelTransformMatrixZeros`` class, which converts


if we change to HankelTransformMatrix, do not forget to change it here

nikosarcevic · 2026-06-12T15:45:35Z

-        )
+        return self._project_spectra_to_radial([c_ell_eval], order)
+
+    # def spherical_correlation(


if obsolete now, delete everything. if we ever need it, we can go back in git history.

nikosarcevic · 2026-07-07T02:00:19Z

+    dln_ell = float(np.log(ell_arr[1] / ell_arr[0]))
+    offset = fhtoffset(dln=dln_ell, mu=order) if use_offset else 0.0
+
+    transformed_power = ifht(c_ell_arr * ell_arr, dln=dln_ell, mu=order, offset=offset)


wrap it if you can multiline

secondlycan we add or point to an analytic/numerical regression test for this normalization? The formula is easy to get wrong by a factor of ell, theta, or 2*pi, so it wold be good to compare against direct quadrature for a smooth test spectrum

nikosarcevic · 2026-07-07T02:01:23Z

+    """Convert power spectrum to 3D correlation function using FFTLog:
+
+    :math:`\\xi(r) = \\int \\frac{k^2 dk}{2\\pi^2} P(k) j_\\mu(kr)`.


snice only order=0 is implemented below, can we make this docstring slightly more specific?

the formula currently writes j_mu, but for spherical bessel transforms the order is usually written as j_l, with the fftlog cylindrical bessel order being mu = l + 1/2.

nikosarcevic · 2026-07-07T02:02:06Z

+    Returns:
+        Radial grid and 3D correlation function.
+    """
+    if order != 0:


can this check happen before the docstring claims a general j_mu spherical transform?

or perhpas even simpler--- can the docstring say this helper currently implements only the l=0 spherical correlation?

nikosarcevic · 2026-07-07T02:02:37Z

+    r = 1.0 / k_arr[::-1]
+    dln_k = float(np.log(k_arr[1] / k_arr[0]))
+    offset = fhtoffset(dln=dln_k, mu=0.5) if use_offset else 0.0


same offset-grid question as above: if use_offset=True, should the returned r grid include the ftlog offset factor rather than always using 1.0 / k_arr[::-1]?

nikosarcevic · 2026-07-07T02:03:18Z

+    transformed_power = ifht(
+        k_arr**1.5 * pk_arr,
+        dln=dln_k,
+        mu=0.5,
+        offset=offset,
+    )
+
+    prefactor = 1.0 / (2.0 * np.pi * r) ** 1.5
+    xi = np.asarray(prefactor * transformed_power, dtype=float)


can we maybe add a direct-quadrature regression test for the spherical normalization too?
mapping from the spherical bessel transform to scipy's cylindrical fftlog convention is subtle, so a test would give us confidence in the k_arr**1.5 input factor and (2*pi*r)**(-1.5) prefactor

nikosarcevic · 2026-07-07T02:06:38Z

@@ -0,0 +1,148 @@
+"""Hankel transforms for projected radial statistics.


we have "projected radial statistics", which is fine, but maybe worth saying this is the direct-matrix implementation specifically.

so perhaps we can say

"""Direct matrix Hankel transforms for projected radial statistics."""

nikosarcevic · 2026-07-07T02:06:58Z

@@ -0,0 +1,148 @@
+"""Hankel transforms for projected radial statistics.
+
+This module provides the ``HankelTransformMatrix`` class, which converts


The docstring says this module provides HankelTransformMatrix, but the class below is HankelTransformMatrixDirect. Could we update the class name in the docstring?

nikosarcevic · 2026-07-07T02:07:31Z

+Hankel operator grids along a predefined radial grid. It is useful for computing
+projected correlation functions, covariance matrices, and higher-order
+radial tensors that appear in weak-lensing and Delta Sigma calculations.


can we clarify the covnention here: this direct matrix backend appears to compute projected transforms of the form int k dk / (2 pi) P(k) J_n(k r).

since this is now a standalone Hankel module, it would help to state the normalization and reciprocal units explicitly

nikosarcevic · 2026-07-07T02:09:25Z

+The class owns the public validation layer and delegates low-level radial and
+Bessel operations to ``hankel_utils``.


this file does not directly delegate bessel operations to hankel_utils; it builds the bessel matrix directly with scipy.special.jv. can you udpate this sentence or just remove it?

nikosarcevic · 2026-07-07T02:10:10Z

+import numpy as np
+from scipy.special import jv
+
+from dsf.hankel.hankel_transform_matrix_zeros import HankelTransformMatrixZeros


this class inherits from HankelTransformMatrixZeros, which is useful for sharing the public matrix API, but it is a little surprising because the numerical grids are not zero-crossing grids.

can we add a short comment/docstring note that this subclass intentionally reuses the matirx interface/binning helpers from the zeros backend while overriding the grid construction?

nikosarcevic · 2026-07-07T02:15:41Z

+            order: Bessel order used by the projected statistic.
+        """
+
+        k = np.geomspace(self.k_min, self.k_max, self.n_k)


np.gradien(np.log(k)) gives endpoint weights that are not exactly trapezoidal integration weights

is this intentional? for direct quadrature, should these weights instead use explicit trapezoidal/log-space integration weights to better control endpoint behavior?

nikosarcevic · 2026-07-07T02:16:21Z

+        dlnk = np.gradient(np.log(k))
+        weight = k**2 * dlnk / (2.0 * np.pi)
+
+        r = np.geomspace(self.r_min, self.r_max, self.n_r)


can u confirm that all allowed order values are non-negative? jv can evaluate more general orders, but the physical bessel orders used here should probably be validated epxlicitly

nikosarcevic · 2026-07-07T02:16:56Z

+        """
+        self._check_order(order)
+
+        product = np.ones_like(self.k[order])
+        for spectrum in spectra:


can we validate that each spectrum has the same shape as self.k[order] before multiplying?

this may already happen upstream, but checking here would give a clearer error if a backend method passes a weird array

so smth like

for spectrum in spectra: if spectrum.shape != self.k[order].shape: raise ValueError( "Each spectrum must be evaluated on the internal k grid." ) product *= spectrum

nikosarcevic · 2026-07-07T02:17:49Z

+        j = self.j[order]
+        j_t = self.j_t[order]
+
+        ndim = len(spectra)


the one- and two-spectrum branches look consistent with the projected convention: j @ (weight * product) for a radial vector and (j * weight * product) @ j_t for a radial covariance matrix

would be great to have a direct-quadrature regression test covering both of these normalizations

nikosarcevic · 2026-07-07T02:19:25Z

+            transformed = j @ (weight * product)
+        elif ndim == 2:
+            transformed = (j * weight * product) @ j_t
+        elif ndim == 3:


so hmm i think this 3-spectrum branch is not computing a third-order radial tensor correctly

shape-wise, j @ (weight * product) gives a 1D radial vector, and then multiplying by j_t @ j_t only works accidentally when n_r == n_k; it does not produce an (n_r, n_r, n_r) tensor

If the intended statistic is

T(r1, r2, r3) = int k dk / (2 pi) product(k) J(r1 k) J(r2 k) J(r3 k),

then this should probably be an einsum

transformed = np.einsum( "ak,bk,ck,k->abc", j, j, j, weight * product, )

can u fix this and add a test with n_r != n_k so the shape bugzzzz cannot hide?

nikosarcevic · 2026-07-07T02:36:49Z

+Fourier-space spectra into projected radial-space quantities using precomputed
+Hankel operator grids along the Bessel zeros. It is useful for computing
+projected correlation functions, covariance matrices, and higher-order
 radial tensors that appear in weak-lensing and Delta Sigma calculations.


Tiny terminology nit: elsewhere we usually write DeltaSigma I think. Could we make this consistent with the rest of DSF docs?

nikosarcevic · 2026-07-07T02:37:48Z

 radial tensors that appear in weak-lensing and Delta Sigma calculations.

 The class owns the public validation layer and delegates low-level radial and
 Bessel operations to ``hankel_utils``.


this says low-level radial and bessel operations are delegated to hankel_utils, but this class directly builds the bessel matrix with scipy.special.jv.

so maybe rephrase to say that shared helper operations, such as bessel zeros, binning, and covariance utilities, are delegated to hankel_utils.

nikosarcevic · 2026-07-07T02:38:50Z

        r_min: Minimum radial scale to cover.
        r_max: Maximum radial scale to cover.
        k_min: Minimum wavenumber to cover.
        k_max: Maximum wavenumber to cover.


can u state the reciprocal unit reqs here? r_min/r_max and k_min/k_max need to use reciprocal units so that k * r is dimensionless. this is mui importante in DSF because user-facing radial quantities may be Mpc/h while some internal calculations may use Mpc.

nikosarcevic · 2026-07-07T02:39:58Z

        self.n_zeros = int(n_zeros)
        self.n_zeros_step = int(n_zeros_step)
        self.prune_r = prune_r
        self.prune_log_space = bool(prune_log_space)
        self.verbose = bool(verbose)
        self.max_iterations = int(max_iterations)


can we avoid silently truncating integer-like inputs? int(n_zeros), int(n_zeros_step), and int(max_iterations) will silently turn e.g. 1000.9 into 1000.

it would be safer to validate that these are integer-like before casting

nikosarcevic · 2026-07-07T02:40:52Z

        if self.n_zeros <= 0:
            raise ValueError("n_zeros must be positive.")
        if self.n_zeros_step <= 0:
            raise ValueError("n_zeros_step must be positive.")
        if self.max_iterations <= 0:
            raise ValueError("max_iterations must be positive.")


can we vaildate that n_zeros, n_zeros_step, and max_iterations are at least 1 as int-like values before casting?

right now non-integer floats can be silently truncated in __init__.

nikosarcevic · 2026-07-07T02:49:47Z

        product = np.ones_like(self.k[order])
        for spectrum in spectra:
            product *= spectrum


can we validate each spectrum shape before multiplying? _evaluate_spectrum() should already enforce this, but this private helper accepts list[FloatArray] directly, so a local check would give a clearer error if it is reused internally

so smth like

for spectrum in spectra: if spectrum.shape != self.k[order].shape: raise ValueError( "Each spectrum must match the internal k-grid shape." ) product *= spectrum

nikosarcevic · 2026-07-07T02:51:18Z

        weighted = product / self.j_next_at_zeros[order] ** 2
        j_matrix = self.j[order]
        norm = self.normalization[order]
        ndim = len(spectra)

        if ndim == 1:
            transformed = np.dot(j_matrix, weighted) * norm
        elif ndim == 2:
            transformed = np.dot(j_matrix, (j_matrix * weighted).T) * norm
        elif ndim == 3:
            transformed = (
                np.einsum(
                    "az,bz,cz,z->abc",
                    j_matrix,
                    j_matrix,
                    j_matrix,
                    weighted,
                )
                * norm
            )


projection code is nicely centralized :)

can we document the output shapes here? For one, two, and three spectra I think these are (n_r,), (n_r, n_r), and (n_r, n_r, n_r).

the 3-spectr branch looks much better here than in the direct backend because it uses an explicit einsum and should produce an (n_r, n_r, n_r) tensor

can we add a test that confirms the shape and symmetry under swapping the three input spectra?

nikosarcevic · 2026-07-07T02:53:23Z

        k_pk: ArrayLike | None = None,
        pk1: SpectrumInput | None = None,
        pk2: SpectrumInput | None = None,
        order: float | int = 0,
        taper: bool = False,
        taper_kwargs: dict | None = None,
        **kwargs,
    ) -> tuple[FloatArray, FloatArray]:
        """Compute a projected covariance matrix from two spectra.

        Args:
            k_pk: Wavenumber grid for tabulated spectra.


clarify whether k_pk is actually a 3D wavenumber grid here or an angular multipole grid?
In projected_correlation() the analogous input is ell, but covariance uses k_pk

If both are using the same projected hankel convention, consistent naming would reduce confusion

nikosarcevic · 2026-07-07T02:53:50Z

        expected_shape = tuple([r_arr.size] * matrix_arr.ndim)
        if matrix_arr.shape != expected_shape:
            raise ValueError(f"matrix shape must be {expected_shape}. Got {matrix_arr.shape}.")


ruff

raise ValueError( f"matrix shape must be {expected_shape}. Got {matrix_arr.shape}." )

nikosarcevic · 2026-07-07T02:54:28Z

    def correlation_matrix(self, covariance: FloatArray) -> FloatArray:
        """Return the correlation matrix associated with a covariance matrix.

        Args:
            covariance: Covariance matrix.

        Returns:
            Dimensionless correlation matrix.
        """
        covariance_arr = np.asarray(covariance, dtype=float)

        if covariance_arr.ndim != 2:
            raise ValueError("covariance must be two-dimensional.")
        if covariance_arr.shape[0] != covariance_arr.shape[1]:
            raise ValueError("covariance must be square.")
        if np.any(~np.isfinite(covariance_arr)):
            raise ValueError("covariance must contain only finite values.")

        return compute_correlation_matrix(covariance_arr)

    def diagonal_error(self, covariance: FloatArray) -> FloatArray:
        """Return one-sigma errors from a covariance matrix.

        Args:
            covariance: Covariance matrix.

        Returns:
            Square root of the covariance diagonal.
        """
        covariance_arr = np.asarray(covariance, dtype=float)

        if covariance_arr.ndim != 2:
            raise ValueError("covariance must be two-dimensional.")
        if covariance_arr.shape[0] != covariance_arr.shape[1]:
            raise ValueError("covariance must be square.")
        if np.any(~np.isfinite(covariance_arr)):
            raise ValueError("covariance must contain only finite values.")

        return compute_diagonal_error(covariance_arr)


these covariance utility wrappers are clean

can we add or confirm tests for zero or negative diagonal entries? correlation_matrix and diagonal_error should fail clearly if the covariance diagonal is not physically valid, unless hankel_utils already handles that

nikosarcevic · 2026-07-07T03:04:25Z

@@ -1 +1 @@
 """Numerical scripts for projected radial statistics.


"scripts" sounds like executable scripts, but this is a utility module. use "Numerical utilities for projected radial statistics" instead?

nikosarcevic · 2026-07-07T03:04:44Z

 This module provides small utilities used when converting Fourier-space
 power spectra into radial-space projected statistics. The scripts cover
 Bessel-function roots, smooth spectrum tapering, radial bin centers,
 radial integration weights, covariance-to-correlation conversion, and
 bin-averaging of radial matrices or tensors.


mention that these helpers are intentionally low-level and assume validated inputs? That is stated below, but I would avoid calling them "scripts" here too. Maybe "utilities" throughout.

nikosarcevic · 2026-07-07T03:05:06Z

 def _is_non_negative_integer(value: float | int) -> bool:
    """Return whether a value is a valid integer Bessel order."""
    value_float = float(value)
    return value_float >= 0.0 and value_float.is_integer()


Could this helper reject non-finite values explicitly? float(np.nan) works and then returns false, but for a validation-style helper it may be clearer to say finite non-negative integer orders only

nikosarcevic · 2026-07-07T03:05:59Z

    if n_zeros <= 0:
        raise ValueError("n_zeros must be positive.")

    if _is_non_negative_integer(order):
        return np.asarray(jn_zeros(int(order), n_zeros), dtype=float)

    order_float = float(order)
    roots = np.empty(n_zeros, dtype=float)


can we validate that order is finite and non-negative before computing zeros? the int path handles negative integers by falling through, but then non-integer negative orders can reach the root bracketing logic. Since the backends only allow non-negative Bessel orders, it would be safer for this utility to enforce the same contract

smth like

order_float = float(order) if not np.isfinite(order_float): raise ValueError("order must be finite.") if order_float < 0: raise ValueError("order must be non-negative.")

then use oder_float below!!

nikosarcevic · 2026-07-07T03:06:39Z

    if n_zeros <= 0:
        raise ValueError("n_zeros must be positive.")

    if _is_non_negative_integer(order):
        return np.asarray(jn_zeros(int(order), n_zeros), dtype=float)

    order_float = float(order)
    roots = np.empty(n_zeros, dtype=float)


validate that n_zeros is integer-like. The type hint says int, but a float can still be passed at runtime and np.empty(n_zeros, ...) will fail with a less clear error.

nikosarcevic · 2026-07-07T03:15:16Z

    bin_grids = np.meshgrid(*([bin_index] * ndim), indexing="ij")
    valid_grids = np.meshgrid(*([valid] * ndim), indexing="ij")
    valid_matrix = np.logical_and.reduce(valid_grids)

    output_indices = tuple(grid[valid_matrix] for grid in bin_grids)
    np.add.at(binned_sum, output_indices, weighted_matrix[valid_matrix])


for high-rank tensors, this builds full meshgrid arrays with shape (n_r,) * ndim.
That is fine for 2D/3D with moderate grids, but it can get very memory-heavy

we can document that this helper is intended for small binned tensors, or use an einsum/loop approach if large tensors are expected?

nikosarcevic · 2026-07-07T03:15:40Z

    norm = _outer_product(bin_weight_sums, ndim)
    binned = np.zeros_like(binned_sum)

    nonzero = norm != 0.0
    binned[nonzero] = binned_sum[nonzero] / norm[nonzero]


same empty-bin concern: leaving binned values at zero when norm == 0 can look like a real measurement. I would prefer raising if any requested bin has zero radial weight

nikosarcevic · 2026-07-07T03:16:11Z

+def apply_taper_spectrum(
+    k: FloatArray,
+    pk: FloatArray,
+    large_k_lower: float = 10.0,
+    large_k_upper: float = 100.0,
+    low_k_lower: float = 0.0,
+    low_k_upper: float = 1.0e-5,
+) -> FloatArray:


can we validate the taper bounds? we need large_k_upper > large_k_lower, low_k_upper > low_k_lower, and probably low_k_upper <= large_k_lower. otherwise this can divide by zero or create overlapping low/high tapers

nikosarcevic · 2026-07-07T03:16:48Z

+    high = k > large_k_lower
+    pk_out[high] *= np.cos(
+        (k[high] - large_k_lower) / (large_k_upper - large_k_lower) * np.pi / 2.0
+    )
+    pk_out[k > large_k_upper] = 0.0
+
+    low = k < low_k_upper
+    pk_out[low] *= np.cos((k[low] - low_k_upper) / (low_k_upper - low_k_lower) * np.pi / 2.0)


the taper formula itself looks good, but can we add tests checking continuity at low_k_upper, large_k_lower, and zeroing beyond low_k_lower/large_k_upper?

nikosarcevic · 2026-07-07T03:18:03Z

+    pk_out = np.copy(pk)
+
+    high = k > large_k_lower
+    pk_out[high] *= np.cos(


ruff

phase = ( (k[high] - large_k_lower) / (large_k_upper - large_k_lower) * np.pi / 2.0 ) pk_out[high] *= np.cos(phase)

check all lines L261, low k taper etc

bclevine added 11 commits June 12, 2026 09:30

Merged hankel backends into a single Hankel module

0a15426

Integrated new hankel framework into the mag bias and covariance code

7ed8cf6

updated tests for hankel module

00bc47a

deleted empty folder

d5f4f41

Added parent handling of unimplemented methods; removed spherical cor…

7382619

…relation from matrix methods

Fixed public api for correlation_matrix and diagonal_error

065b6bd

Added NotImplemented fallbacks for correlation_matrix and diagonal_error

682b876

Updated docs

0a8ac95

Updated docs formatting

45c545f

Added benchmarking tests for hankel module

ca1e99f

Fixed linting mistakes

7e5374b

bclevine added enhancement New feature or request modelling labels Jun 12, 2026

bclevine added 2 commits June 12, 2026 14:49

Changed behavior of _evaluate_tabulated_spectrum to reject spectra th…

14bc563

…at don't cover the full matrix k-grid

Updated docs with better info on hankel module

d6e1a24

nikosarcevic reviewed Jul 7, 2026

View reviewed changes

		from dsf.utils.validators import validate_positive_scalar


		class HankelTransformMatrixDirect(HankelTransformMatrixZeros):

		"""Convert power spectrum to 3D correlation function using FFTLog:

		:math:`\\xi(r) = \\int \\frac{k^2 dk}{2\\pi^2} P(k) j_\\mu(kr)`.

		@@ -0,0 +1,148 @@
		"""Hankel transforms for projected radial statistics.

		@@ -0,0 +1,148 @@
		"""Hankel transforms for projected radial statistics.

		This module provides the ``HankelTransformMatrix`` class, which converts

		The class owns the public validation layer and delegates low-level radial and
		Bessel operations to ``hankel_utils``.

		@@ -1 +1 @@
		"""Numerical scripts for projected radial statistics.

Uh oh!

Conversation

bclevine commented Jun 12, 2026

Description

Uh oh!

nikosarcevic left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!