Refactor DeltaSigma into a standalone class

benchmarks/covariance/benchmark_custom_power_spectra.py

@@ -129,11 +129,12 @@ def compare_covariances(reference, equivalent, modified):
 
         np.testing.assert_allclose(custom, baseline, rtol=1e-6, atol=0.0)
 
+        equivalent_difference = np.linalg.norm(custom - baseline) / np.linalg.norm(baseline)
         difference = np.linalg.norm(changed - baseline)
         relative_change = difference / np.linalg.norm(baseline)
 
         print(f"\n{block}")
-        print("  Equivalent custom spectra: PASS")
+        print(f"  Equivalent custom spectra: PASS ({equivalent_difference:.3e})")
         print(f"  Relative Frobenius difference: {relative_change:.3e}")
 
         if block != "cov_gm_gg":

benchmarks/data_vector/benchmark_delta_sigma.py

@@ -0,0 +1,162 @@
+"""Benchmark standalone Delta Sigma against the legacy CCL projection."""
+
+from __future__ import annotations
+
+import time
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pyccl as ccl
+
+from dsf.data_vector.ccl_integration import DeltaSigma
+
+
+__all__ = ["main"]
+
+OUTPUT_DIR = Path(__file__).resolve().parent / "figures" / "delta_sigma"
+
+SCALE_FACTORS = (0.5, 0.7, 1.0)
+RADII = np.geomspace(0.6, 110.0, 50)
+
+RTOL = 1e-3
+
+
+class LegacyHaloProfile(ccl.halos.profiles.profile_base.HaloProfile):
+    """Reproduce the original CCL-based Delta Sigma projection."""
+
+    def __init__(self, *, pk2d: ccl.Pk2D) -> None:
+        """Initialize the legacy profile with its original FFTLog settings."""
+        self.pk2d = pk2d
+
+        self.fourier_analytic = True
+        self._fourier = self._fourier_analytic
+        self.projected_analytic = False
+        self.cumul2d_analytic = False
+
+        super().__init__(
+            mass_def="200m",
+            concentration=ccl.halos.concentration.ConcentrationConstant(
+                5.0, mass_def="200m"
+            ),
+        )
+
+        self.update_precision_fftlog(
+            padding_lo_fftlog=1e-6,
+            padding_hi_fftlog=1e6,
+            n_per_decade=5000,
+        )
+
+    def _fourier_analytic(self, cosmo, k, M, a):
+        """Evaluate the Fourier-space profile using the legacy CCL interface."""
+        if np.size(M) != 1:
+            raise ValueError("M must be a single value.")
+
+        k_arr = np.atleast_1d(np.asarray(k, dtype=float))
+        values = np.asarray(self.pk2d(k_arr, a), dtype=float)
+
+        return values.reshape((1, k_arr.size))
+
+
+def build_cosmology():
+    """Construct the cosmology used in the covariance benchmarks."""
+    h = 0.69
+    omega_b = 0.022 / h**2
+
+    return ccl.Cosmology(
+        Omega_c=0.292 - omega_b,
+        Omega_b=omega_b,
+        h=h,
+        A_s=2.115e-9,
+        n_s=0.965,
+        transfer_function="boltzmann_camb",
+        matter_power_spectrum="halofit",
+    )
+
+
+def legacy_delta_sigma(profile, cosmo, r, a):
+    """Evaluate Delta Sigma using the original CCL implementation."""
+    sigma = profile.projected(cosmo, r, 1.0, a)
+    mean_sigma = profile.cumul2d(cosmo, r, 1.0, a)
+
+    return np.asarray(mean_sigma - sigma, dtype=float) / 1e12
+
+
+def benchmark_projection(legacy, standalone, cosmo, r, a):
+    """Compare numerical accuracy and execution time at one scale factor."""
+    start = time.perf_counter()
+    reference = legacy_delta_sigma(legacy, cosmo, r, a)
+    legacy_time = time.perf_counter() - start
+
+    start = time.perf_counter()
+    result = standalone.delta_sigma(cosmo=cosmo, r=r, a=a) / 1e12
+    standalone_time = time.perf_counter() - start
+
+    np.testing.assert_allclose(result, reference, rtol=RTOL, atol=0.0)
+
+    fractional_difference = (result - reference) / reference
+
+    print(f"\na = {a:.2f}")
+    print(f"  Legacy CCL:       {legacy_time:.3f} s")
+    print(f"  Standalone FFTLog: {standalone_time:.3f} s")
+    print(
+        "  Maximum fractional difference: "
+        f"{np.max(np.abs(fractional_difference)):.3e}"
+    )
+
+    return reference, result, fractional_difference
+
+
+def plot_comparison(results):
+    """Plot Delta Sigma predictions and their fractional differences."""
+    fig, axes = plt.subplots(2, 1, figsize=(7, 8), sharex=True)
+
+    for a, (reference, result, fractional_difference) in results.items():
+        axes[0].plot(RADII, reference, label=f"Legacy, a={a:.2f}")
+        axes[0].plot(RADII, result, linestyle="--", label=f"Standalone, a={a:.2f}")
+
+        axes[1].plot(RADII, fractional_difference, label=f"a={a:.2f}")
+
+    axes[0].set_xscale("log")
+    axes[0].set_yscale("symlog")
+    axes[0].set_ylabel(r"$\Delta\Sigma\ [M_\odot/\mathrm{pc}^2]$")
+    axes[0].legend(fontsize=8)
+
+    axes[1].axhline(0.0, color="gray", linestyle="--")
+    axes[1].axhline(RTOL, color="gray", linestyle=":")
+    axes[1].axhline(-RTOL, color="gray", linestyle=":")
+    axes[1].set_xscale("log")
+    axes[1].set_xlabel("Comoving projected radius [Mpc]")
+    axes[1].set_ylabel("Fractional difference")
+    axes[1].legend()
+
+    fig.tight_layout()
+    fig.savefig(OUTPUT_DIR / "delta_sigma_comparison.png", dpi=200)
+    plt.close(fig)
+
+
+def main():
+    """Run the legacy-versus-standalone Delta Sigma benchmark."""
+    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
+
+    cosmo = build_cosmology()
+    pk2d = cosmo.get_nonlin_power()
+
+    legacy = LegacyHaloProfile(pk2d=pk2d)
+    standalone = DeltaSigma(pk2d=pk2d)
+
+    results = {}
+
+    for a in SCALE_FACTORS:
+        results[a] = benchmark_projection(
+            legacy, standalone, cosmo, RADII, a
+        )
+
+    plot_comparison(results)
+
+    print(f"\nFigures saved to: {OUTPUT_DIR}")
+    print("All Delta Sigma benchmark checks passed.")
+
+
+if __name__ == "__main__":
+    main()

docs/api/dsf.hankel.hankel_transform.rst

@@ -0,0 +1,7 @@
+dsf.hankel.hankel\_transform module
+===================================
+
+.. automodule:: dsf.hankel.hankel_transform
+   :members:
+   :show-inheritance:
+   :undoc-members:

docs/api/dsf.hankel.rst

@@ -7,7 +7,7 @@ Submodules
 .. toctree::
    :maxdepth: 2
 
-   dsf.hankel.hankel
+   dsf.hankel.hankel_transform
    dsf.hankel.hankel_transform_base
    dsf.hankel.hankel_transform_fftlog
    dsf.hankel.hankel_transform_matrix_direct

src/dsf/covariance/cov_builder.py

@@ -40,7 +40,7 @@
     lens_averaged_power_spectrum,
 )
 from dsf.covariance.ingredients.sigma_crit import effective_squared_sigma_crit
-from dsf.hankel.hankel import HankelTransform
+from dsf.hankel.hankel_transform import HankelTransform
 from dsf.utils.converters import (
     resolve_h,
     resolve_omega_m,

src/dsf/covariance/ingredients/cov_blocks.py

@@ -22,7 +22,7 @@
 
 import numpy as np
 
-from dsf.hankel.hankel import HankelTransform
+from dsf.hankel.hankel_transform import HankelTransform
 from dsf.utils.types import ArrayLike
 from dsf.utils.validators import (
     as_1d_float_array,

src/dsf/data_vector/ccl_integration.py

@@ -1,85 +1,100 @@
-"""CCL integration scripts for generic Delta Sigma profiles."""
+r"""Excess surface density projection from CCL power spectra.
+
+This module provides the standalone ``DeltaSigma`` class for calculating
+excess surface density profiles from CCL ``Pk2D`` objects.
+
+The calculation uses a second-order projected Hankel transform,
+
+.. math::
+
+    \Delta\Sigma(R, a) =
+    \int_0^\infty \frac{k\,\mathrm{d}k}{2\pi}
+    \tilde{\rho}(k, a) J_2(kR),
+
+where the input ``Pk2D`` supplies the Fourier-space density profile.
+The transform is evaluated using DSF's FFTLog implementation.
+"""
 
 from __future__ import annotations
 
 import numpy as np
 import pyccl as ccl
 from numpy.typing import NDArray
 
+from dsf.hankel.hankel_transform import HankelTransform
+from dsf.utils.validators import validate_positive_1d_array, validate_scale_factor
+
+__all__ = ["DeltaSigma"]
+
 
-class HaloProfileGeneric(ccl.halos.profiles.profile_base.HaloProfile):
-    """CCL halo-profile wrapper around a generic Pk2D.
+class DeltaSigma:
+    r"""Calculate excess surface density from a CCL ``Pk2D`` object.
 
-    This class uses CCL's HaloProfile API only as a projection backend.
-    The mass definition, concentration, and input mass are dummy values
-    required by the API and are not physically used by this wrapper.
+    The class evaluates the Fourier-space profile on a logarithmic wavenumber
+    grid and applies a second-order Hankel transform to obtain
+    :math:`\Delta\Sigma(R, a)`.
+
+    The input ``Pk2D`` must include the required density normalization.
+    Projected radii are comoving and expressed in Mpc.
     """
 
     def __init__(
         self,
         *,
         pk2d: ccl.Pk2D,
-        padding_hi_fftlog: float = 1.0e3,
-        padding_lo_fftlog: float = 1.0e-2,
-        n_per_decade: int = 5000,
+        k_min: float = 1.0e-6,
+        k_max: float = 1.0e6,
+        n_k: int = 5000,
     ) -> None:
-        """Initialize the generic CCL profile.
+        """Initialize the excess surface density projection.
 
         Args:
-            pk2d: CCL ``Pk2D`` object used as the Fourier-space profile.
-            padding_hi_fftlog: High-k FFTLog padding used by CCL.
-            padding_lo_fftlog: Low-k FFTLog padding used by CCL.
-            n_per_decade: FFTLog sampling density.
+            pk2d: CCL ``Pk2D`` object defining the Fourier-space density profile.
+            k_min: Minimum wavenumber in 1/Mpc.
+            k_max: Maximum wavenumber in 1/Mpc.
+            n_k: Number of logarithmically spaced wavenumbers.
         """
-        self.pk2d = pk2d
+        if not np.isfinite(k_min) or not np.isfinite(k_max):
+            raise ValueError("k_min and k_max must be finite.")
 
-        self.fourier_analytic = True
-        self._fourier = self._fourier_analytic
-        self.projected_analytic = False
-        self.cumul2d_analytic = False
-
-        super().__init__(
-            mass_def="200m",
-            concentration=ccl.halos.concentration.ConcentrationConstant(
-                5.0,
-                mass_def="200m",
-            ),
-        )
+        if not 0.0 < k_min < k_max:
+            raise ValueError("Require 0 < k_min < k_max.")
 
-        self.update_precision_fftlog(
-            padding_hi_fftlog=padding_hi_fftlog,
-            padding_lo_fftlog=padding_lo_fftlog,
-            n_per_decade=n_per_decade,
-        )
+        if n_k < 2:
+            raise ValueError("n_k must be at least 2.")
 
-    def _fourier_analytic(
+        self.pk2d = pk2d
+        self.k = np.geomspace(k_min, k_max, n_k)
+        self.hankel = HankelTransform(backend="fftlog")
+
+    def delta_sigma(
         self,
         cosmo: ccl.Cosmology,
-        k: float | NDArray[np.float64],
-        M: float | NDArray[np.float64],
+        r: NDArray[np.float64],
         a: float,
     ) -> NDArray[np.float64]:
-        """Evaluate the stored ``Pk2D`` object as a Fourier-space profile.
+        r"""Calculate :math:`\Delta\Sigma(R, a)` at projected radii.
 
         Args:
-            cosmo: CCL cosmology object. Included for compatibility with the
-                CCL ``HaloProfile`` API.
-            k: Wavenumber in 1/Mpc.
-            M: Halo mass in solar masses. This wrapper only supports a single
-                dummy mass value.
-            a: Scale factor.
+            cosmo: CCL cosmology object.
+            r: Projected comoving radii in Mpc.
+            a: Scale factor at which the profile is evaluated.
 
         Returns:
-            Fourier-space profile values with shape ``(1, n_k)``.
+            Excess surface density in the projected-density units implied
+            by the input ``Pk2D``.
         """
-        _ = cosmo  # Required by the CCL HaloProfile API, unused here.
+        r_arr = validate_positive_1d_array(r, "r")
+        validate_scale_factor(a)
 
-        if np.size(M) != 1:
-            raise ValueError(
-                "M must be a single value for HaloProfileGeneric."
-            )
+        pk = np.asarray(self.pk2d(self.k, a), dtype=float)
 
-        k_arr = np.atleast_1d(np.asarray(k, dtype=float))
-        values = np.asarray(self.pk2d(k_arr, a), dtype=float)
+        _, delta_sigma = self.hankel.projected_correlation_interpolated(
+            theta=r_arr,
+            radial_input=self.k,
+            spectrum=pk,
+            order=2,
+            grid_spacing="linear",
+        )
 
-        return values.reshape((1, k_arr.size))
+        return np.asarray(delta_sigma, dtype=float)