Source code for wale.CommonUtils

import numpy as np
from scipy.integrate import simpson
from astropy import units as u

from .FilterFunctions import *


[docs] def apply_pixel_window(ells, theta_deg=10.0, npix=1200): """ Compute the pixel window function for a square map and apply it to multipoles. Parameters ---------- ells : array_like Multipole values (ℓ) at which the window function is evaluated. theta_deg : float, optional Total angular size of the map in degrees (default is 10.0). npix : int, optional Number of pixels per side of the square map (default is 1200). Returns ------- W_ell : ndarray The pixel window function evaluated at each ℓ. """ theta_pix_rad = np.deg2rad(theta_deg / npix) arg = ells * theta_pix_rad / 2 W_ell = np.sinc(arg / np.pi) ** 2 return W_ell
[docs] def fourier_coordinate(x, y, map_size): """ Return the 1D Fourier coordinate index corresponding to 2D (x, y) on a square map. Parameters ---------- x : int X-coordinate (horizontal index). y : int Y-coordinate (vertical index). map_size : int Size of one side of the square map. Returns ------- idx : int Flattened Fourier-space index. """ return (((map_size // 2) + 1) * x) + y
[docs] def get_moments(kappa_values, pdf_values): """ Compute the first four moments of a given 1D probability distribution. Parameters ---------- kappa_values : ndarray Bin centers or sample points along the kappa (x-axis). pdf_values : ndarray Corresponding PDF values at each kappa. Returns ------- mean_kappa : float Mean of the distribution. variance : float Variance of the distribution. S_3 : float Skewness (third standardized moment). K : float Kurtosis minus 3 (excess kurtosis). norm : float Normalization constant of the input PDF. """ norm = np.trapz(pdf_values, kappa_values) normalized_pdf_values = pdf_values / norm mean_kappa = np.trapz(kappa_values * normalized_pdf_values, kappa_values) variance = np.trapz( (kappa_values - mean_kappa) ** 2 * normalized_pdf_values, kappa_values ) third_moment = np.trapz( (kappa_values - mean_kappa) ** 3 * normalized_pdf_values, kappa_values ) fourth_moment = np.trapz( (kappa_values - mean_kappa) ** 4 * normalized_pdf_values, kappa_values ) S_3 = third_moment / (variance**2.0) K = fourth_moment / variance**2 - 3 return mean_kappa, variance, S_3, K, norm
[docs] def get_l1_from_pdf(counts, bins): """ Compute the L1 norm (∫|x|P(x)dx) from a histogram representation of a PDF. Parameters ---------- counts : ndarray Histogram bin counts or PDF values (P(x)). bins : ndarray Bin centers or values corresponding to the counts. Returns ------- l1_norm : ndarray L1 norm approximation (P(x) * |x| per bin). """ return counts * np.abs(bins)
[docs] def compute_sigma_kappa_squared( theta_arcmin, chis, lensingweights, redshifts, k, pnl, filter_type, h ): """ Compute the smoothed convergence variance σ²_κ(θ) at a given angular scale using a filter. This function computes the convergence power spectrum P_κ(ℓ) from a 3D P(k, z) and integrates over ℓ using a top-hat or starlet filter. Parameters ---------- theta_arcmin : float Angular smoothing scale θ in arcminutes. chis : ndarray Comoving distances χ (in Mpc) corresponding to redshifts. lensingweights : ndarray Lensing kernel W(χ) evaluated at each χ. redshifts : ndarray Redshifts corresponding to chis. k : ndarray Wavenumber grid (in h/Mpc). pnl : 2D ndarray Nonlinear power spectrum P(k, z), shape (n_z, len(k)). filter_type : str Type of filter to apply ("tophat" or "starlet"). h : float Reduced Hubble constant (H0 / 100). Returns ------- sigma2 : float Smoothed convergence variance σ²_κ(θ). """ theta_rad = (theta_arcmin * u.arcmin).to(u.rad).value ell = np.logspace(1, 5, 200) P_kappa = np.zeros_like(ell) pnl_array = np.array([pnl[z_] for z_ in redshifts]) # shape: (n_chi, nk) # plt.loglog(k, pnl_array.T) # plt.show() for i, l in enumerate(ell): chis_h_inv = chis * h # now in h⁻¹ Mpc k_l = l / chis_h_inv # Interpolate P(k) at each redshift slice pk_vals = np.array( [ np.interp(k_l[j], k, pnl_array[j], left=0, right=0) for j in range(len(chis)) ] ) integrand = (lensingweights / chis) ** 2 * pk_vals P_kappa[i] = simpson(integrand, chis) # Apply top-hat filter window in Fourier space if filter_type == "tophat": W = top_hat_filter(ell, 2.0 * theta_rad) - top_hat_filter(ell, theta_rad) elif filter_type == "starlet": W = starlet_filter(ell, 2.0 * theta_rad) - starlet_filter(ell, theta_rad) pixel_window = apply_pixel_window(ell, theta_deg=theta_rad * u.rad.to(u.deg)) integrand = ell * P_kappa * (W**2) sigma2 = simpson(integrand, ell) / (2.0 * np.pi) return sigma2