From 621c521b46c77381cb3ec4c4237d0191d12107de Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Sun, 31 Aug 2025 11:46:21 +0900 Subject: [PATCH 01/10] implement theta bin average --- tjpcov/covariance_builder.py | 21 +- tjpcov/covariance_gaussian_fsky.py | 21 +- tjpcov/legendre.py | 443 +++++++++++++++++++++++++++++ tjpcov/wigner_transform.py | 48 +++- 4 files changed, 504 insertions(+), 29 deletions(-) create mode 100644 tjpcov/legendre.py diff --git a/tjpcov/covariance_builder.py b/tjpcov/covariance_builder.py index 7fbf1ef0..57971d5c 100644 --- a/tjpcov/covariance_builder.py +++ b/tjpcov/covariance_builder.py @@ -6,7 +6,7 @@ import pyccl as ccl import sacc -from .wigner_transform import bin_cov, WignerTransform +from .wigner_transform import WignerTransform from . import tools from .covariance_io import CovarianceIO @@ -1086,13 +1086,14 @@ def get_Wigner_transform(self): :obj:`~tjpcov.wigner_transform.WignerTransform` instance """ if self.WT is None: - # Removing ell <= 1 (following original implementation) - ell = np.arange(2, self.lmax + 1) - theta, _, _ = self.get_binning_info(in_radians=True) + # Removing ell <= 1 is done in legendre.py + ell = np.arange(0, self.lmax + 1) + theta, _, theta_edges = self.get_binning_info(in_radians=True) WT_kwargs = { "ell": ell, "theta": theta, + "theta_edges": theta_edges, "s1_s2": [(2, 2), (2, -2), (0, 2), (2, 0), (0, 0)], } @@ -1127,7 +1128,7 @@ def get_covariance_block( """ # For now we just use the EE block which should be dominant over the # EB, BE and BB pieces - cov = self._get_fourier_block(tracer_comb1, tracer_comb2) + cov, SN = self._get_fourier_block(tracer_comb1, tracer_comb2) WT = self.get_Wigner_transform() @@ -1137,15 +1138,11 @@ def get_covariance_block( s1_s2_1 = s1_s2_1[xi_plus_minus1] if isinstance(s1_s2_2, dict): s1_s2_2 = s1_s2_2[xi_plus_minus2] - # Remove ell <= 1 for WT (following original implementation) - ell = np.arange(2, self.lmax + 1) - cov = cov[2:][:, 2:] + # Removing ell <= 1 is done in legendre.py + ell = np.arange(0, self.lmax + 1) th, cov = WT.projected_covariance( - ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov + ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov, SN=SN ) - if binned: - theta, _, theta_edges = self.get_binning_info(in_radians=False) - thb, cov = bin_cov(r=theta, r_bins=theta_edges, cov=cov) return cov diff --git a/tjpcov/covariance_gaussian_fsky.py b/tjpcov/covariance_gaussian_fsky.py index 6adc2d4b..d0f50342 100644 --- a/tjpcov/covariance_gaussian_fsky.py +++ b/tjpcov/covariance_gaussian_fsky.py @@ -160,16 +160,23 @@ def get_covariance_block( else 0 ) + if for_real: + cov = np.diag( + (cl[13] + SN[13]) * (cl[24] + SN[24]) + + (cl[14] + SN[14]) * (cl[23] + SN[23]) + - np.ones_like(cl[13]) * SN[13] * SN[24] + - np.ones_like(cl[14]) * SN[14] * SN[23] + ) + # If it is to compute the real space covariance, return the + # covariance before binning or normalizing + # pure shot/shape noise term will be computed separately + return cov, SN[13] * SN[24] + SN[14] * SN[23] + cov = np.diag( (cl[13] + SN[13]) * (cl[24] + SN[24]) + (cl[14] + SN[14]) * (cl[23] + SN[23]) ) - if for_real: - # If it is to compute the real space covariance, return the - # covariance before binning or normalizing - return cov - norm = (2 * ell + 1) * np.gradient(ell) * self.fsky cov /= norm @@ -227,9 +234,9 @@ def _get_fourier_block(self, tracer_comb1, tracer_comb2): """ # For now we just use the EE block which should be dominant over the # EB, BE and BB pieces when projecting to real space - cov = self.fourier.get_covariance_block( + cov, SN = self.fourier.get_covariance_block( tracer_comb1, tracer_comb2, for_real=True, lmax=self.lmax ) norm = np.pi * 4 * self.fsky - return cov / norm + return cov / norm, SN / norm diff --git a/tjpcov/legendre.py b/tjpcov/legendre.py new file mode 100644 index 00000000..41c05477 --- /dev/null +++ b/tjpcov/legendre.py @@ -0,0 +1,443 @@ +# copied from +# https://github.com/joezuntz/cosmosis-standard-library/blob/main +# /shear/cl_to_xi_fullsky/legendre.py +from __future__ import print_function +from builtins import range +import numpy as np +from scipy.special import lpn + +PI = np.pi + + +def sin_filter(ell_max, ell_right): + ells = np.arange(ell_max + 1) + y = (ell_max - ells) / (ell_max - ell_right) + return np.where( + ells > ell_right, y - np.sin(2 * PI * y) / 2 / PI, np.ones(len(ells)) + ) + + +def get_N_ell(ell): + """N_ell (eq. 2.16)""" + ell = np.atleast_1d(ell) + N_ell = np.sqrt(2.0 / (ell - 1) / ell / (ell + 1) / (ell + 2)) + N_ell[ell < 2] = 0.0 + N_ell[ell == 2] = np.sqrt(2.0 / (4 * 3 * 2)) + return N_ell + + +def apply_filter(ell_max, high_l_filter, legfacs): + f = sin_filter(ell_max, ell_max * high_l_filter) + return legfacs * np.tile(f, (legfacs.shape[0], 1)) + + +def get_legfactors_00(ells, thetas): + n_ell, n_theta = len(ells), len(thetas) + legfacs = np.zeros((n_theta, n_ell)) + for it, t in enumerate(thetas): + legfacs[it] = (2 * ells + 1) * lpn(ells[-1], np.cos(t))[0] / 4.0 / PI + return legfacs + + +def get_legfactors_02(ells, thetas): + n_ell, n_theta = len(ells), len(thetas) + legfacs = np.zeros((n_theta, n_ell)) + ell_factor = np.zeros(len(ells)) + ell_factor[1:] = (2 * ells[1:] + 1) / 4.0 / PI / ells[1:] / (ells[1:] + 1) + for it, t in enumerate(thetas): + P2l = P2l_rec_norm(ells, np.cos(t)) + legfacs[it] = P2l * ell_factor + return legfacs + + +# def get_legfactors_22(ells, thetas): +# gp, gm = precomp_GpGm(ells, thetas) +# N_ell = get_N_ell(ells) +# prefac = ((2 * ells + 1) / 2.0 / PI) * N_ell * N_ell / 2 +# leg_factors_p = prefac * (gp + gm) +# leg_factors_m = prefac * (gp - gm) +# return (leg_factors_p, leg_factors_m) + + +def P2l_rec(ells, cost): + """Calculate P2l using recurrence relation""" + P22 = 3 * (1 - cost**2) + P23 = 15 * cost * (1 - cost**2) + P2l = np.zeros(len(ells)) + P2l[0] = 0.0 + P2l[1] = 0.0 + P2l[2] = P22 + P2l[3] = P23 + for ell in ells[4:]: + # print ell, P2l[ell-1], P2l[ell-2] + P2l[ell] = ( + (2 * ell - 1) * cost * P2l[ell - 1] - (ell + 2 - 1) * P2l[ell - 2] + ) / (ell - 2) + return P2l + + +def P2l_norm_prefac(ell): + return np.sqrt( + (2.0 * ell + 1) + / ((ell - 1.0) * ell * (ell + 1.0) * (ell + 2.0)) + / 4.0 + / PI + ) + + +def P2l_rec_norm(ells, cost): + """Calculate P2l using recurrence relation for normalised P2l""" + P22 = 3.0 * (1.0 - cost**2) + P23 = 15.0 * cost * (1.0 - cost**2) + P2l = np.zeros(len(ells)) + P2l[0] = 0.0 + P2l[1] = 0.0 + P2l[2] = P22 + P2l[3] = P23 + P2l_norm = np.copy(P2l) + P2l_norm[2] *= P2l_norm_prefac(2) + P2l_norm[3] *= P2l_norm_prefac(3) + for ell in ells[4:]: + # print ell, P2l[ell-1], P2l[ell-2] + a = np.sqrt((4 * ell**2 - 1.0) / (ell**2 - 4)) + b = cost * P2l_norm[ell - 1] + c = ( + np.sqrt(((ell - 1.0) ** 2 - 4) / (4 * (ell - 1.0) ** 2 - 1)) + * P2l_norm[ell - 2] + ) + # print a,b,c + P2l_norm[ell] = a * (b - c) + # print ell, P2l_norm[ell], P2l_norm_prefac(ell) + P2l[ell] = P2l_norm[ell] / P2l_norm_prefac(ell) + return P2l + + +def P2l_rec_binav(ells, cost_min, cost_max): + """Calculate P2l using recurrence relation for normalised P2l""" + P2l_binav = np.zeros(len(ells)) + P2l_binav[0] = 0.0 + P2l_binav[1] = 0.0 + # coefficients that are a function of ell only + ell = ells[2:] + coeff_lm1 = ell + 2.0 / (2.0 * ell + 1.0) + coeff_lp1 = 2.0 / (2.0 * ell + 1.0) + coeff_l = 2.0 - ell + # computation of legendre polynomials + # this computes all polynomials of order 0 to ell_max+1 and for all ell's + lpns_min = lpn(ell[-1] + 1, cost_min)[0][1:] + lpns_max = lpn(ell[-1] + 1, cost_max)[0][1:] + # terms in the numerator of average P2l + term_lm1 = coeff_lm1 * (lpns_max[:-2] - lpns_min[:-2]) + term_lp1 = coeff_lp1 * (lpns_max[2:] - lpns_min[2:]) + term_l = coeff_l * (cost_max * lpns_max[1:-1] - cost_min * lpns_min[1:-1]) + # denominator in average P2l + dcost = cost_max - cost_min + # computation of bin-averaged P2l(ell) + P2l_binav[ell] = (term_lm1 + term_l - term_lp1) / dcost + return P2l_binav + + +def Gp_plus_minus_Gm_binav_dep1(ells, cost_min, cost_max): + """Calculate bin-averaged G_{l,2}^{+/-}""" + Gp_plus_Gm = np.zeros(len(ells)) + Gp_minus_Gm = np.zeros(len(ells)) + + # for ell=0,1 it is 0 + Gp_plus_Gm[0:1] = 0.0 + Gp_minus_Gm[0:1] = 0.0 + + # for the rest of ell's compute equation (5.8) in + # https://arxiv.org/abs/1911.11947 + ell = ells[2:] + # ---coefficients including only P_l + coeff_lm1 = -ell * (ell - 1.0) / 2.0 * (ell + 2.0 / (2.0 * ell + 1)) - ( + ell + 2.0 + ) + coeff_lp1 = ell * (ell - 1.0) / (2.0 * ell + 1.0) + coeff_l = -ell * (ell - 1.0) * (2.0 - ell) / 2.0 + # coeff_l_plus = coeff_l - 2.0 * (ell - 1.0) + # coeff_l_minus = coeff_l + 2.0 * (ell - 1.0) + # ---coefficients including dP_l/dx + coeff_dlm1_plus = -2.0 * (ell + 2.0) + # coeff_dlm1_minus = -coeff_dlm1_plus + coeff_xdl_plus = 2.0 * (ell - 1.0) + # coeff_xdl_minus = -coeff_xdl_plus + coeff_xdlm1 = ell + 2.0 + coeff_dl = 4.0 - ell + + # computation of legendre polynomials + # this computes all polynomials of order 0 to ell_max+1 and for all ell's + lpns_min = lpn(ell[-1] + 1, cost_min)[0][1:] + lpns_max = lpn(ell[-1] + 1, cost_max)[0][1:] + dlpns_min = lpn(ell[-1] + 1, cost_min)[1][1:] + dlpns_max = lpn(ell[-1] + 1, cost_max)[1][1:] + + # denominator in average + dcost = cost_max - cost_min + + # numerator in average + # ---common part in both plus and minus + common_part = coeff_lm1 * (lpns_max[:-2] - lpns_min[:-2]) + common_part += coeff_l * ( + cost_max * lpns_max[1:-1] - cost_min * lpns_min[1:-1] + ) + common_part += coeff_lp1 * (lpns_max[2:] - lpns_min[2:]) + common_part += coeff_dl * (dlpns_max[1:-1] - dlpns_min[1:-1]) + common_part += coeff_xdlm1 * ( + cost_max * dlpns_max[1:-1] - cost_min * dlpns_min[1:-1] + ) + # ---plus + Gp_plus_Gm_extra = coeff_xdl_plus * ( + cost_max * dlpns_max[:-2] - cost_min * dlpns_min[:-2] + ) + Gp_plus_Gm_extra += coeff_dlm1_plus * (dlpns_max[:-2] - dlpns_min[:-2]) + Gp_plus_Gm[2:] = common_part + Gp_plus_Gm_extra + Gp_plus_Gm /= dcost + # ---minus + Gp_minus_Gm_extra = -Gp_plus_Gm_extra + Gp_minus_Gm[2:] = common_part + Gp_minus_Gm_extra + Gp_minus_Gm /= dcost + + return Gp_plus_Gm, Gp_minus_Gm + + +def Pl_rec_binav(ells, cost_min, cost_max): + """Calculate average Pl""" + Pl_binav = np.zeros(len(ells)) + Pl_binav[0] = 1.0 + # coefficients that are a function of ell only + ell = ells[1:] + coeff = 1.0 / (2.0 * ell + 1.0) + # computation of legendre polynomials + # this computes all polynomials of order 0 to ell_max+1 and for all ell's + lpns_min = lpn(ell[-1] + 1, cost_min)[0] + lpns_max = lpn(ell[-1] + 1, cost_max)[0] + # terms in the numerator of average Pl + term_lm1 = lpns_max[:-2] - lpns_min[:-2] + term_lp1 = lpns_max[2:] - lpns_min[2:] + # denominator in average Pl + dcost = cost_max - cost_min + # computation of bin-averaged Pl(ell) + Pl_binav[ell] = coeff * (term_lp1 - term_lm1) / dcost + return Pl_binav + + +def theta_bin_means_to_edges(thetas, binning="log"): + """This function is deprecated now that we are passing in theta_edge + values. Possibly useful for other purposes, so leaving it here.""" + print("Calculating theta bin edges") + print("n_theta_bins=", len(thetas)) + print("thetas = ", thetas) + print("thetas in arcmin = ", thetas / PI * 180 * 60) + # array of theta edges from mean values + tedges = np.zeros(len(thetas) + 1) + for i in range(len(thetas)): + # bin width selection + if binning == "log": + # logarithmic spacing + if i == 0: + dtheta = np.log10(thetas[i + 1]) - np.log10(thetas[i]) + else: + dtheta = np.log10(thetas[i]) - np.log10(thetas[i - 1]) + tedges[i] = 10.0 ** (np.log10(thetas[i]) - dtheta / 2.0) + tedges[i + 1] = 10.0 ** (np.log10(thetas[i]) + dtheta / 2.0) + else: + # linear spacing + if i == 0: + dtheta = thetas[i + 1] - thetas[i] + else: + dtheta = thetas[i] - thetas[i - 1] + tedges[i] = thetas[i] - dtheta / 2.0 + tedges[i + 1] = thetas[i] + dtheta / 2.0 + # if the spacing is large, first value might be negative + if tedges[0] < 0.0: + tedges[0] = 0.0 + # print('theta_edges = ',tedges) + return tedges + + +def get_legfactors_02_binav(ells, theta_edges): + print("getting bin averaged leg factors for 02") + n_ell, n_theta = len(ells), len(theta_edges) - 1 + # theta_edges = theta_bin_means_to_edges(thetas) # this does geometric mean + legfacs = np.zeros((n_theta, n_ell)) + ell_factor = np.zeros(len(ells)) + ell_factor[1:] = (2 * ells[1:] + 1) / 4.0 / PI / ells[1:] / (ells[1:] + 1) + for it, t in enumerate(theta_edges[1:]): + t_min = theta_edges[it] + t_max = t + cost_min = np.cos(t_min) # thetas are already converted to radians + cost_max = np.cos(t_max) + P2l = P2l_rec_binav(ells, cost_min, cost_max) + legfacs[it] = P2l * ell_factor + return legfacs + + +def get_legfactors_00_binav(ells, theta_edges): + print("getting bin averaged leg factors for 00") + n_ell, n_theta = len(ells), len(theta_edges) - 1 + # theta_edges = theta_bin_means_to_edges(thetas) # this does geometric mean + legfacs = np.zeros((n_theta, n_ell)) + ell_factor = np.zeros(len(ells)) + ell_factor[1:] = (2 * ells[1:] + 1) / 4.0 / PI + for it, t in enumerate(theta_edges[1:]): + t_min = theta_edges[it] + t_max = t + cost_min = np.cos(t_min) # thetas are already converted to radians + cost_max = np.cos(t_max) + Pl = Pl_rec_binav(ells, cost_min, cost_max) + legfacs[it] = Pl * ell_factor + return legfacs + + +def Gp_plus_minus_Gm_binav_dep2(ells, cost_min, cost_max): + """Calculate bin-averaged G_{l,2}^{+/-}""" + Gp_plus_Gm = np.zeros(len(ells)) + Gp_minus_Gm = np.zeros(len(ells)) + # for ell=0,1 it is 0 + Gp_plus_Gm[0:1] = 0.0 + Gp_minus_Gm[0:1] = 0.0 + + # for the rest of ell's compute equation (5.8) in + # https://arxiv.org/abs/1911.11947 + ell = ells[2:] + # ---coefficients including only P_l + coeff_lm1 = -ell * (ell - 1.0) / 2.0 * (ell + 2.0 / (2.0 * ell + 1)) - ( + ell + 2.0 + ) + coeff_lp1 = ell * (ell - 1.0) / (2.0 * ell + 1.0) + coeff_l = -ell * (ell - 1.0) * (2.0 - ell) / 2.0 + # coeff_l_plus = coeff_l - 2.0 * (ell - 1.0) + # coeff_l_minus = coeff_l + 2.0 * (ell - 1.0) + # ---coefficients including dP_l/dx + coeff_dlm1_plus = -2.0 * (ell + 2.0) + coeff_xdl_plus = 2.0 * (ell - 1.0) + coeff_xdlm1 = ell + 2.0 + coeff_dl = 4.0 - ell + + # computation of legendre polynomials and derivatives + # this computes all polynomials of order 0 to ell_max+1 and for all ell's + Pl_calc_min = np.asarray(lpn(ell[-1] + 1, cost_min))[:, 1:] + Pl_calc_max = np.asarray(lpn(ell[-1] + 1, cost_max))[:, 1:] + lpns_min = Pl_calc_min[0, :] + dlpns_min = Pl_calc_min[1, :] + lpns_max = Pl_calc_max[0, :] + dlpns_max = Pl_calc_max[1, :] + + # denominator in average + dcost = cost_max - cost_min + + # numerator in average + # ---common part in both plus and minus + common_part = coeff_lm1 * (lpns_max[:-2] - lpns_min[:-2]) + common_part += coeff_l * ( + cost_max * lpns_max[1:-1] - cost_min * lpns_min[1:-1] + ) + common_part += coeff_lp1 * (lpns_max[2:] - lpns_min[2:]) + common_part += coeff_dl * (dlpns_max[1:-1] - dlpns_min[1:-1]) + common_part += coeff_xdlm1 * ( + cost_max * dlpns_max[:-2] - cost_min * dlpns_min[:-2] + ) + # ---plus + Gp_plus_Gm_extra = coeff_xdl_plus * ( + cost_max * dlpns_max[1:-1] - cost_min * dlpns_min[1:-1] + ) + Gp_plus_Gm_extra += -coeff_xdl_plus * (lpns_max[1:-1] - lpns_min[1:-1]) + Gp_plus_Gm_extra += coeff_dlm1_plus * (dlpns_max[:-2] - dlpns_min[:-2]) + Gp_plus_Gm[2:] = common_part + Gp_plus_Gm_extra + Gp_plus_Gm /= dcost + # ---minus + Gp_minus_Gm_extra = -Gp_plus_Gm_extra + Gp_minus_Gm[2:] = common_part + Gp_minus_Gm_extra + Gp_minus_Gm /= dcost + + return Gp_plus_Gm, Gp_minus_Gm + + +def Gp_plus_minus_Gm_binav(ells, cost_min, cost_max): + """Calculate bin-averaged G_{l,2}^{+/-}""" + Gp_plus_Gm = np.zeros(len(ells)) + Gp_minus_Gm = np.zeros(len(ells)) + + # for ell=0,1 it is 0 + Gp_plus_Gm[0:1] = 0.0 + Gp_minus_Gm[0:1] = 0.0 + + # for the rest of ell's compute equation (5.8) in + # https://arxiv.org/abs/1911.11947 + ell = ells[2:] + # ---coefficients including only P_l + coeff_lm1 = -ell * (ell - 1.0) / 2.0 * (ell + 2.0 / (2.0 * ell + 1)) - ( + ell + 2.0 + ) + coeff_lp1 = ell * (ell - 1.0) / (2.0 * ell + 1.0) + coeff_l = -ell * (ell - 1.0) * (2.0 - ell) / 2.0 + coeff_l_plus = -2.0 * (ell - 1.0) + # ---coefficients including dP_l/dx + coeff_dlm1_plus = -2.0 * (ell + 2.0) + coeff_xdl_plus = 2.0 * (ell - 1.0) + coeff_xdlm1 = ell + 2.0 + coeff_dl = 4.0 - ell + + # computation of legendre polynomials + # this computes all polynomials of order 0 to ell_max+1 and for all ell's + lpns_min = lpn(ell[-1] + 1, cost_min)[0][1:] + lpns_max = lpn(ell[-1] + 1, cost_max)[0][1:] + dlpns_min = lpn(ell[-1] + 1, cost_min)[1][1:] + dlpns_max = lpn(ell[-1] + 1, cost_max)[1][1:] + + # denominator in average + dcost = cost_max - cost_min + + # numerator in average + # ---common part in both plus and minus + common_part = coeff_lm1 * (lpns_max[:-2] - lpns_min[:-2]) + common_part += coeff_l * ( + cost_max * lpns_max[1:-1] - cost_min * lpns_min[1:-1] + ) + common_part += coeff_lp1 * (lpns_max[2:] - lpns_min[2:]) + common_part += coeff_dl * (dlpns_max[1:-1] - dlpns_min[1:-1]) + common_part += coeff_xdlm1 * ( + cost_max * dlpns_max[:-2] - cost_min * dlpns_min[:-2] + ) + # ---plus + Gp_plus_Gm_extra = coeff_xdl_plus * ( + cost_max * dlpns_max[1:-1] - cost_min * dlpns_min[1:-1] + ) + Gp_plus_Gm_extra += coeff_l_plus * (lpns_max[1:-1] - lpns_min[1:-1]) + Gp_plus_Gm_extra += coeff_dlm1_plus * (dlpns_max[:-2] - dlpns_min[:-2]) + Gp_plus_Gm[2:] = common_part + Gp_plus_Gm_extra + Gp_plus_Gm /= dcost + # ---minus + Gp_minus_Gm_extra = -Gp_plus_Gm_extra + Gp_minus_Gm[2:] = common_part + Gp_minus_Gm_extra + Gp_minus_Gm /= dcost + + return Gp_plus_Gm, Gp_minus_Gm + + +def get_legfactors_22_binav(ells, theta_edges): + print("getting bin averaged leg factors for 22") + n_ell, n_theta = len(ells), len(theta_edges) - 1 + # theta_edges = theta_bin_means_to_edges(thetas) # this does geometric mean + leg_factors_p = np.zeros((n_theta, n_ell)) + leg_factors_m = np.zeros((n_theta, n_ell)) + ell_factor = np.zeros(len(ells)) + ell_factor[2:] = ( + (2 * ells[2:] + 1) + / 2.0 + / PI + / ells[2:] + / ells[2:] + / (ells[2:] + 1.0) + / (ells[2:] + 1.0) + ) + for it, t in enumerate(theta_edges[1:]): + t_min = theta_edges[it] + t_max = t + cost_min = np.cos(t_min) # thetas are already converted to radians + cost_max = np.cos(t_max) + gp, gm = Gp_plus_minus_Gm_binav(ells, cost_min, cost_max) + leg_factors_p[it] = gp * ell_factor + leg_factors_m[it] = gm * ell_factor + return (leg_factors_p, leg_factors_m) diff --git a/tjpcov/wigner_transform.py b/tjpcov/wigner_transform.py index eb612114..704a2c19 100755 --- a/tjpcov/wigner_transform.py +++ b/tjpcov/wigner_transform.py @@ -10,7 +10,11 @@ from scipy.special import binom from scipy.special import eval_jacobi as jacobi from scipy.special import jn - +from .legendre import ( + get_legfactors_00_binav, + get_legfactors_02_binav, + get_legfactors_22_binav, +) # FIXME: # 1. Do we need to pass logger? @@ -20,7 +24,7 @@ class WignerTransform: """Class to compute curved sky Hankel transforms with wigner-d matrices.""" - def __init__(self, theta, ell, s1_s2, ncpu=None): + def __init__(self, theta, theta_edges, ell, s1_s2, ncpu=None): """Initialize the class for the given angles, scales and spins. Args: @@ -59,12 +63,26 @@ def __init__(self, theta, ell, s1_s2, ncpu=None): self.s1_s2s = s1_s2 self.theta = {} self.theta = theta - - # compute the wigner-d matrices. + self.theta_edges = theta_edges + # compute the bin-averaged legendre polynomials. for s1, s2 in s1_s2: - self.wig_d[(s1, s2)] = wigner_d_parallel( - s1, s2, theta, self.ell, ncpu=ncpu - ) + if s1 == s2 == 0: + self.wig_d[(s1, s2)] = get_legfactors_00_binav( + self.ell, theta_edges + ) + elif s1 == s2 == 2: + self.wig_d[(s1, s2)] = get_legfactors_22_binav( + self.ell, theta_edges + )[0] + elif (np.abs(s1) == np.abs(s2) == 2) and (s1 * s2 < 0): + self.wig_d[(s1, s2)] = get_legfactors_22_binav( + self.ell, theta_edges + )[1] + else: + self.wig_d[(s1, s2)] = get_legfactors_02_binav( + self.ell, theta_edges + ) + self.taper_f = None self.taper_f2 = None @@ -188,7 +206,14 @@ def cl_cov_grid(self, ell_cl, cl_cov, taper=False, **taper_kwargs): # return self.ell, cl def projected_covariance( - self, ell_cl, cl_cov, s1_s2, s1_s2_cross=None, taper=False, **kwargs + self, + ell_cl, + cl_cov, + SN, + s1_s2, + s1_s2_cross=None, + taper=False, + **kwargs, ): """Convert C_ell covariance to correlation function. @@ -240,11 +265,14 @@ def projected_covariance( cov = np.einsum( "rk,kl,sl->rs", - self.wig_d[s1_s2] * np.sqrt(self.norm) * self.grad_ell, + self.wig_d[s1_s2] * self.grad_ell, cl_cov2, - self.wig_d[s1_s2_cross] * np.sqrt(self.norm), + self.wig_d[s1_s2_cross], optimize=True, ) + + # TODO: add SN term + # FIXME: Check normalization return self.theta, cov From c5ebc6fb745cf0b59db6614dea5926bde674d5db Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Fri, 23 Jan 2026 17:00:06 +0900 Subject: [PATCH 02/10] exact N_pair --- tjpcov/covariance_builder.py | 40 +++++++++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/tjpcov/covariance_builder.py b/tjpcov/covariance_builder.py index 57971d5c..fcca678d 100644 --- a/tjpcov/covariance_builder.py +++ b/tjpcov/covariance_builder.py @@ -1017,7 +1017,7 @@ def get_binning_info(self, binning="log", in_radians=True): # TODO: This should be obtained from the sacc file or the input # configuration. Check how it is done in TXPipe: # https://github.com/LSSTDESC/TXPipe/blob/a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/covariance.py#L23 - + theta_eff = self.get_theta_eff() nbpw = theta_eff.size @@ -1138,11 +1138,45 @@ def get_covariance_block( s1_s2_1 = s1_s2_1[xi_plus_minus1] if isinstance(s1_s2_2, dict): s1_s2_2 = s1_s2_2[xi_plus_minus2] + + # Load Npair + # see https://github.com/LSSTDESC/TXPipe/blob/a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/twopoint_plots.py#L215 + if tracer_comb1 == tracer_comb2: + sacc_file = self.io.get_sacc_file() + data_type = self.get_tracer_comb_data_types(tracer_comb1)[0] + D = sacc_file.get_data_points(data_type, (tracer_comb1[0], tracer_comb1[1])) + Npair = np.array([d.get_tag("npair") for d in D]) + + if np.any(Npair == None): + # assuming no survey boundaries. + if np.abs(s1_s2_1[0]) == np.abs(s1_s2_1[1]) == 2: + SN *= 2 + + cov_sn = SN / np.pi / (WT.theta_edges[1:]**2 - WT.theta_edges[:-1]**2) + + else: + # catalog level N_pair from treecorr + T_sn = 1 + if tracer_comb1[0] in self.sigma_e: + T_sn *= self.sigma_e[tracer_comb1[0]] ** 2 + if tracer_comb1[1] in self.sigma_e: + T_sn *= self.sigma_e[tracer_comb1[1]] ** 2 + + if (tracer_comb1[0] in self.sigma_e) and (tracer_comb1[1] in self.sigma_e): + T_sn *= 2 + + # Eq. 64 of https://arxiv.org/abs/2410.06962 + cov_sn = 2 * T_sn / Npair + + # Project sample variance term and mixed term. # Removing ell <= 1 is done in legendre.py ell = np.arange(0, self.lmax + 1) th, cov = WT.projected_covariance( - ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov, SN=SN - ) + ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov) + + # Add pure shot/shape noise contribution. + if tracer_comb1 == tracer_comb2: + cov += np.diag(cov_sn) return cov From 320b30aeace93b60f2dc094f88ffa01d6c2bb15e Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Fri, 23 Jan 2026 17:08:47 +0900 Subject: [PATCH 03/10] lint --- tjpcov/covariance_builder.py | 34 ++++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/tjpcov/covariance_builder.py b/tjpcov/covariance_builder.py index fcca678d..e58a7d2d 100644 --- a/tjpcov/covariance_builder.py +++ b/tjpcov/covariance_builder.py @@ -1017,7 +1017,7 @@ def get_binning_info(self, binning="log", in_radians=True): # TODO: This should be obtained from the sacc file or the input # configuration. Check how it is done in TXPipe: # https://github.com/LSSTDESC/TXPipe/blob/a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/covariance.py#L23 - + theta_eff = self.get_theta_eff() nbpw = theta_eff.size @@ -1140,20 +1140,27 @@ def get_covariance_block( s1_s2_2 = s1_s2_2[xi_plus_minus2] # Load Npair - # see https://github.com/LSSTDESC/TXPipe/blob/a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/twopoint_plots.py#L215 + # see https://github.com/LSSTDESC/TXPipe/blob/ + # a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/twopoint_plots.py#L215 if tracer_comb1 == tracer_comb2: sacc_file = self.io.get_sacc_file() - data_type = self.get_tracer_comb_data_types(tracer_comb1)[0] - D = sacc_file.get_data_points(data_type, (tracer_comb1[0], tracer_comb1[1])) + data_type = self.get_tracer_comb_data_types(tracer_comb1)[0] + D = sacc_file.get_data_points( + data_type, (tracer_comb1[0], tracer_comb1[1]) + ) Npair = np.array([d.get_tag("npair") for d in D]) - if np.any(Npair == None): + if np.any(Npair is None): # assuming no survey boundaries. if np.abs(s1_s2_1[0]) == np.abs(s1_s2_1[1]) == 2: SN *= 2 - - cov_sn = SN / np.pi / (WT.theta_edges[1:]**2 - WT.theta_edges[:-1]**2) - + + cov_sn = ( + SN + / np.pi + / (WT.theta_edges[1:] ** 2 - WT.theta_edges[:-1] ** 2) + ) + else: # catalog level N_pair from treecorr T_sn = 1 @@ -1162,17 +1169,20 @@ def get_covariance_block( if tracer_comb1[1] in self.sigma_e: T_sn *= self.sigma_e[tracer_comb1[1]] ** 2 - if (tracer_comb1[0] in self.sigma_e) and (tracer_comb1[1] in self.sigma_e): + if (tracer_comb1[0] in self.sigma_e) and ( + tracer_comb1[1] in self.sigma_e + ): T_sn *= 2 - + # Eq. 64 of https://arxiv.org/abs/2410.06962 cov_sn = 2 * T_sn / Npair - + # Project sample variance term and mixed term. # Removing ell <= 1 is done in legendre.py ell = np.arange(0, self.lmax + 1) th, cov = WT.projected_covariance( - ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov) + ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov + ) # Add pure shot/shape noise contribution. if tracer_comb1 == tracer_comb2: From 71f8b2f12afbe4b063bd02eb03a7af82d3e41796 Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Fri, 23 Jan 2026 17:32:20 +0900 Subject: [PATCH 04/10] update black --- run_tjpcov.py | 1 - tests/test_covariance_calculator.py | 1 - tests/test_covariance_fourier_gaussian_nmt.py | 1 - tjpcov/covariance_gaussian_fsky.py | 6 ++---- 4 files changed, 2 insertions(+), 7 deletions(-) diff --git a/run_tjpcov.py b/run_tjpcov.py index 0f1ce0bb..322f5b93 100644 --- a/run_tjpcov.py +++ b/run_tjpcov.py @@ -2,7 +2,6 @@ from tjpcov.covariance_calculator import CovarianceCalculator - if __name__ == "__main__": import argparse diff --git a/tests/test_covariance_calculator.py b/tests/test_covariance_calculator.py index 2a30b396..13044aaf 100644 --- a/tests/test_covariance_calculator.py +++ b/tests/test_covariance_calculator.py @@ -8,7 +8,6 @@ import sacc import shutil - INPUT_YML = "./tests/data/conf_covariance_calculator.yml" OUTDIR = "./tests/tmp/" diff --git a/tests/test_covariance_fourier_gaussian_nmt.py b/tests/test_covariance_fourier_gaussian_nmt.py index 8e3016d9..d9e99a33 100644 --- a/tests/test_covariance_fourier_gaussian_nmt.py +++ b/tests/test_covariance_fourier_gaussian_nmt.py @@ -13,7 +13,6 @@ from tjpcov.covariance_fourier_gaussian_nmt import FourierGaussianNmt from tjpcov.covariance_io import CovarianceIO - ROOT = "tests/benchmarks/32_DES_tjpcov_bm/" OUTDIR = "tests/tmp/" INPUT_YML = "tests/data/conf_covariance_gaussian_fourier_nmt.yaml" diff --git a/tjpcov/covariance_gaussian_fsky.py b/tjpcov/covariance_gaussian_fsky.py index d0f50342..766d4449 100644 --- a/tjpcov/covariance_gaussian_fsky.py +++ b/tjpcov/covariance_gaussian_fsky.py @@ -50,10 +50,8 @@ def get_binning_info(self, binning="linear"): ell_edges = out[2] return ell, ell_eff, ell_edges else: - warnings.warn( - "No bandpower windows found, \ - falling back to linear method" - ) + warnings.warn("No bandpower windows found, \ + falling back to linear method") ell_eff = self.get_ell_eff() nbpw = ell_eff.size From 6654942b7ca84a352b14fae591e38573e8a43274 Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Fri, 23 Jan 2026 17:50:07 +0900 Subject: [PATCH 05/10] require python<3.14 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index e5fc2212..17fcf9c0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ authors = [ ] description = "Covariances for LSST DESC" readme = "README.md" -requires-python = ">=3.10" +requires-python = ">=3.10,<3.14" classifiers = [ "Development Status :: 4 - Beta", "Programming Language :: Python :: 3", From a21801eeb41537a7c2e1c1bb11298c336cee1d5d Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Fri, 23 Jan 2026 17:54:14 +0900 Subject: [PATCH 06/10] environment.yml python<3.14 --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 44597528..ee60ea86 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ name: tjpcov channels: - conda-forge dependencies: - - python>=3.10 + - python>=3.10,<3.14 - pip - wheel - black From 28f069e7ca3a76f423d156d3efd622c7a480487b Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Fri, 23 Jan 2026 18:04:34 +0900 Subject: [PATCH 07/10] replace lpn --- tjpcov/legendre.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tjpcov/legendre.py b/tjpcov/legendre.py index 41c05477..6d5a70fb 100644 --- a/tjpcov/legendre.py +++ b/tjpcov/legendre.py @@ -4,7 +4,7 @@ from __future__ import print_function from builtins import range import numpy as np -from scipy.special import lpn +from scipy.special import legendre_p_all as lpn PI = np.pi From 7b300074df7b04792ec8f6ccaed080be541b2b3a Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Sun, 25 Jan 2026 14:20:48 +0900 Subject: [PATCH 08/10] fix tests --- tests/test_covariance_gaussian_fsky.py | 8 ++- tests/test_covariance_real_base.py | 20 ++---- tests/test_wigner_transform.py | 15 +++-- tjpcov/covariance_builder.py | 89 ++++++++++++++------------ tjpcov/covariance_gaussian_fsky.py | 3 +- tjpcov/wigner_transform.py | 3 - 6 files changed, 68 insertions(+), 70 deletions(-) diff --git a/tests/test_covariance_gaussian_fsky.py b/tests/test_covariance_gaussian_fsky.py index 5d3c1f68..131bac75 100644 --- a/tests/test_covariance_gaussian_fsky.py +++ b/tests/test_covariance_gaussian_fsky.py @@ -118,13 +118,14 @@ def test_Fourier_get_covariance_block(cov_fg_fsky, mock_cosmo): tracer_comb1=trs, tracer_comb2=trs, for_real=True ) # 2. Check block - cov2 = cov_fg_fsky.get_covariance_block( + cov2, SN = cov_fg_fsky.get_covariance_block( tracer_comb1=trs, tracer_comb2=trs, for_real=True, lmax=30 ) ell = np.arange(30 + 1) ccltr = ccl_tracers["src0"] cl = ccl.angular_cl(mock_cosmo, ccltr, ccltr, ell) + tracer_noise["src0"] cov = np.diag(2 * cl**2) + cov2 = cov2 + np.diag(np.ones_like(cl)*SN) assert cov2.shape == (ell.size, ell.size) np.testing.assert_allclose(cov2, cov) @@ -150,13 +151,14 @@ def test_Fourier_get_covariance_block(cov_fg_fsky, mock_cosmo): def test_Real_get_fourier_block( cov_rg_fsky, cov_fg_fsky, tracer_comb1, tracer_comb2 ): - cov = cov_rg_fsky._get_fourier_block(tracer_comb1, tracer_comb2) - cov2 = cov_fg_fsky.get_covariance_block( + cov, SN = cov_rg_fsky._get_fourier_block(tracer_comb1, tracer_comb2) + cov2, SN2 = cov_fg_fsky.get_covariance_block( tracer_comb1, tracer_comb2, for_real=True, lmax=cov_rg_fsky.lmax ) norm = np.pi * 4 * cov_rg_fsky.fsky assert np.all(cov == cov2 / norm) + assert np.all(SN == SN2 / norm) def test_smoke_get_covariance(cov_fg_fsky, cov_rg_fsky): diff --git a/tests/test_covariance_real_base.py b/tests/test_covariance_real_base.py index 6022dc8e..16356fdc 100644 --- a/tests/test_covariance_real_base.py +++ b/tests/test_covariance_real_base.py @@ -5,7 +5,7 @@ from functools import partial -from tjpcov.wigner_transform import bin_cov, WignerTransform +from tjpcov.wigner_transform import WignerTransform from tjpcov.covariance_builder import ( CovarianceProjectedReal, CovarianceReal, @@ -115,7 +115,7 @@ def test_get_Wigner_transform(cov_prj_real): wt = cov_prj_real.get_Wigner_transform() assert isinstance(wt, WignerTransform) - assert np.all(wt.ell == np.arange(2, cov_prj_real.lmax + 1)) + assert np.all(wt.ell == np.arange(0, cov_prj_real.lmax + 1)) assert np.all(wt.theta == cov_prj_real.get_binning_info()[0]) assert wt.s1_s2s == [(2, 2), (2, -2), (0, 2), (2, 0), (0, 0)] @@ -158,12 +158,12 @@ def test_build_matrix_from_blocks(cov_prj_real): ) def test_get_covariance_block(cov_prj_real, tracer_comb1, tracer_comb2): lmax = cov_prj_real.lmax - ell = np.arange(2, lmax + 1) + ell = np.arange(0, lmax + 1) fourier_block = np.random.rand(lmax + 1, lmax + 1) # Dynamically override the method on this instance. def override_fourier_block(self, tracer_comb1, tracer_comb2): - return fourier_block + return fourier_block, None cov_prj_real._get_fourier_block = partial( override_fourier_block, cov_prj_real @@ -180,23 +180,15 @@ def override_fourier_block(self, tracer_comb1, tracer_comb2): ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, - cl_cov=fourier_block[2:][:, 2:], + cl_cov=fourier_block, ) gcov_xi_1 = cov_prj_real.get_covariance_block( - tracer_comb1=tracer_comb1, tracer_comb2=tracer_comb2, binned=False + tracer_comb1=tracer_comb1, tracer_comb2=tracer_comb2 ) assert np.max(np.abs((gcov_xi_1 + 1e-100) / (cov + 1e-100) - 1)) < 1e-5 - - gcov_xi_1 = cov_prj_real.get_covariance_block( - tracer_comb1=tracer_comb1, tracer_comb2=tracer_comb2, binned=True - ) - - theta, _, theta_edges = cov_prj_real.get_binning_info(in_radians=False) - _, cov = bin_cov(r=theta, r_bins=theta_edges, cov=cov) assert gcov_xi_1.shape == (20, 20) - assert np.max(np.abs((gcov_xi_1 + 1e-100) / (cov + 1e-100) - 1)) < 1e-5 @pytest.mark.parametrize( diff --git a/tests/test_wigner_transform.py b/tests/test_wigner_transform.py index 4609c10d..da5c9af9 100644 --- a/tests/test_wigner_transform.py +++ b/tests/test_wigner_transform.py @@ -7,11 +7,13 @@ def get_WT_kwargs(): lmax = 96 - ell = np.arange(2, lmax + 1) - theta = np.sort(np.pi / ell)[::2] # Force them to have different sizes + ell = np.arange(0, lmax + 1) + theta = np.sort(np.pi / ell[1:])[::2] # Force them to have different sizes + theta_edges = np.append(theta, [np.max(theta) * 1.05]) WT_kwargs = { "ell": ell, "theta": theta, + "theta_edges": theta_edges, "s1_s2": [(2, 2), (2, -2), (0, 2), (2, 0), (0, 0)], } return WT_kwargs @@ -103,14 +105,13 @@ def test_projected_covariance(s1_s2, s1_s2_cross): wt.projected_covariance(wt.ell[10:], mat, s1_s2, s1_s2_cross) th, matb = wt.projected_covariance(wt.ell, mat, s1_s2, s1_s2_cross) - wd_a = wigner_transform.wigner_d(*s1_s2, wt.theta, wt.ell) - wd_b = wigner_transform.wigner_d(*s1_s2_cross, wt.theta, wt.ell) + wd_a = wt.wig_d[*s1_s2] + wd_b = wt.wig_d[*s1_s2_cross] matb_2 = ( - (wd_a * np.sqrt(wt.norm) * wt.grad_ell) + (wd_a * wt.grad_ell) @ mat - @ (wd_b * np.sqrt(wt.norm)).T + @ wd_b.T ) - assert np.all(th == wt.theta) assert np.max(np.abs(matb / matb_2) - 1) < 1e-5 diff --git a/tjpcov/covariance_builder.py b/tjpcov/covariance_builder.py index e58a7d2d..37c3da09 100644 --- a/tjpcov/covariance_builder.py +++ b/tjpcov/covariance_builder.py @@ -995,6 +995,7 @@ def __init__(self, config): "You need to specify the lmax you want to " "compute the Fourier covariance up to" ) + self.cov_type = None @property def fourier(self): @@ -1111,7 +1112,6 @@ def get_covariance_block( tracer_comb2, xi_plus_minus1="plus", xi_plus_minus2="plus", - binned=True, ): """Compute a single covariance matrix for a given pair of xi. @@ -1139,44 +1139,6 @@ def get_covariance_block( if isinstance(s1_s2_2, dict): s1_s2_2 = s1_s2_2[xi_plus_minus2] - # Load Npair - # see https://github.com/LSSTDESC/TXPipe/blob/ - # a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/twopoint_plots.py#L215 - if tracer_comb1 == tracer_comb2: - sacc_file = self.io.get_sacc_file() - data_type = self.get_tracer_comb_data_types(tracer_comb1)[0] - D = sacc_file.get_data_points( - data_type, (tracer_comb1[0], tracer_comb1[1]) - ) - Npair = np.array([d.get_tag("npair") for d in D]) - - if np.any(Npair is None): - # assuming no survey boundaries. - if np.abs(s1_s2_1[0]) == np.abs(s1_s2_1[1]) == 2: - SN *= 2 - - cov_sn = ( - SN - / np.pi - / (WT.theta_edges[1:] ** 2 - WT.theta_edges[:-1] ** 2) - ) - - else: - # catalog level N_pair from treecorr - T_sn = 1 - if tracer_comb1[0] in self.sigma_e: - T_sn *= self.sigma_e[tracer_comb1[0]] ** 2 - if tracer_comb1[1] in self.sigma_e: - T_sn *= self.sigma_e[tracer_comb1[1]] ** 2 - - if (tracer_comb1[0] in self.sigma_e) and ( - tracer_comb1[1] in self.sigma_e - ): - T_sn *= 2 - - # Eq. 64 of https://arxiv.org/abs/2410.06962 - cov_sn = 2 * T_sn / Npair - # Project sample variance term and mixed term. # Removing ell <= 1 is done in legendre.py ell = np.arange(0, self.lmax + 1) @@ -1184,9 +1146,52 @@ def get_covariance_block( ell_cl=ell, s1_s2=s1_s2_1, s1_s2_cross=s1_s2_2, cl_cov=cov ) - # Add pure shot/shape noise contribution. - if tracer_comb1 == tracer_comb2: - cov += np.diag(cov_sn) + if self.cov_type == "gauss": + # denominator in average + dcost = np.cos(WT.theta_edges[1:]) - np.cos(WT.theta_edges[:-1]) + + if SN is not None: + # Load Npair + # see https://github.com/LSSTDESC/TXPipe/blob/ + # a9dfdb7809ac7ed6c162fd3930c643a67afcd881/txpipe/twopoint_plots.py#L215 + if tracer_comb1 == tracer_comb2: + sacc_file = self.io.get_sacc_file() + data_type = self.get_tracer_comb_data_types(tracer_comb1)[ + 0 + ] + D = sacc_file.get_data_points( + data_type, (tracer_comb1[0], tracer_comb1[1]) + ) + Npair = np.array([d.get_tag("npair") for d in D]) + + if Npair[0] is None: + # assuming no survey boundaries. + if np.abs(s1_s2_1[0]) == np.abs(s1_s2_1[1]) == 2: + SN *= 2 + + cov_sn = SN / (np.pi * dcost) + else: + # catalog level N_pair from treecorr + T_sn = 1 + if tracer_comb1[0] in self.sigma_e: + T_sn *= self.sigma_e[tracer_comb1[0]] ** 2 + if tracer_comb1[1] in self.sigma_e: + T_sn *= self.sigma_e[tracer_comb1[1]] ** 2 + + if (tracer_comb1[0] in self.sigma_e) and ( + tracer_comb1[1] in self.sigma_e + ): + T_sn *= 2 + # Eq. 64 of https://arxiv.org/abs/2410.06962 + cov_sn = 2 * T_sn / Npair + + # Add pure shot/shape noise contribution. + cov += np.diag(cov_sn) + + else: + # TODO: Projection of NaMaster covariance via + # Eq. 67 of https://arxiv.org/abs/2012.08568 + pass return cov diff --git a/tjpcov/covariance_gaussian_fsky.py b/tjpcov/covariance_gaussian_fsky.py index 766d4449..ebd4dbe7 100644 --- a/tjpcov/covariance_gaussian_fsky.py +++ b/tjpcov/covariance_gaussian_fsky.py @@ -198,7 +198,7 @@ def get_covariance_block( class RealGaussianFsky(CovarianceProjectedReal): """Class to compute the Real space Gaussian cov. with the Knox formula. - It projects the the Fourier space Gaussian covariance into the real space. + It projects the Fourier space Gaussian covariance into the real space. """ cov_type = "gauss" @@ -219,6 +219,7 @@ def __init__(self, config): # sacc file. self.fourier = FourierGaussianFsky(config) self.fsky = self.fourier.fsky + self.cov_type = "gauss" def _get_fourier_block(self, tracer_comb1, tracer_comb2): """Return the Fourier covariance block for two pair of tracers. diff --git a/tjpcov/wigner_transform.py b/tjpcov/wigner_transform.py index 704a2c19..1bcda8ad 100755 --- a/tjpcov/wigner_transform.py +++ b/tjpcov/wigner_transform.py @@ -209,7 +209,6 @@ def projected_covariance( self, ell_cl, cl_cov, - SN, s1_s2, s1_s2_cross=None, taper=False, @@ -271,8 +270,6 @@ def projected_covariance( optimize=True, ) - # TODO: add SN term - # FIXME: Check normalization return self.theta, cov From c76fe4b24e7d9cb42ae871fcf9f6552c80a7b7a5 Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Sun, 25 Jan 2026 14:26:44 +0900 Subject: [PATCH 09/10] lint --- tests/test_covariance_gaussian_fsky.py | 2 +- tests/test_wigner_transform.py | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/test_covariance_gaussian_fsky.py b/tests/test_covariance_gaussian_fsky.py index 131bac75..1e9c9f88 100644 --- a/tests/test_covariance_gaussian_fsky.py +++ b/tests/test_covariance_gaussian_fsky.py @@ -125,7 +125,7 @@ def test_Fourier_get_covariance_block(cov_fg_fsky, mock_cosmo): ccltr = ccl_tracers["src0"] cl = ccl.angular_cl(mock_cosmo, ccltr, ccltr, ell) + tracer_noise["src0"] cov = np.diag(2 * cl**2) - cov2 = cov2 + np.diag(np.ones_like(cl)*SN) + cov2 = cov2 + np.diag(np.ones_like(cl) * SN) assert cov2.shape == (ell.size, ell.size) np.testing.assert_allclose(cov2, cov) diff --git a/tests/test_wigner_transform.py b/tests/test_wigner_transform.py index da5c9af9..cdc231ad 100644 --- a/tests/test_wigner_transform.py +++ b/tests/test_wigner_transform.py @@ -107,11 +107,7 @@ def test_projected_covariance(s1_s2, s1_s2_cross): th, matb = wt.projected_covariance(wt.ell, mat, s1_s2, s1_s2_cross) wd_a = wt.wig_d[*s1_s2] wd_b = wt.wig_d[*s1_s2_cross] - matb_2 = ( - (wd_a * wt.grad_ell) - @ mat - @ wd_b.T - ) + matb_2 = (wd_a * wt.grad_ell) @ mat @ wd_b.T assert np.all(th == wt.theta) assert np.max(np.abs(matb / matb_2) - 1) < 1e-5 From c124e2ebc6e701d3ccf53d5b07e1bfdbc65aab31 Mon Sep 17 00:00:00 2001 From: RyoTerasawa Date: Sun, 25 Jan 2026 05:53:02 +0000 Subject: [PATCH 10/10] update legendre.py --- tjpcov/legendre.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tjpcov/legendre.py b/tjpcov/legendre.py index 6d5a70fb..36d90d9d 100644 --- a/tjpcov/legendre.py +++ b/tjpcov/legendre.py @@ -169,8 +169,8 @@ def Gp_plus_minus_Gm_binav_dep1(ells, cost_min, cost_max): # this computes all polynomials of order 0 to ell_max+1 and for all ell's lpns_min = lpn(ell[-1] + 1, cost_min)[0][1:] lpns_max = lpn(ell[-1] + 1, cost_max)[0][1:] - dlpns_min = lpn(ell[-1] + 1, cost_min)[1][1:] - dlpns_max = lpn(ell[-1] + 1, cost_max)[1][1:] + dlpns_min = lpn(ell[-1] + 1, cost_min, diff_n=1)[1][1:] + dlpns_max = lpn(ell[-1] + 1, cost_max, diff_n=1)[1][1:] # denominator in average dcost = cost_max - cost_min @@ -383,8 +383,8 @@ def Gp_plus_minus_Gm_binav(ells, cost_min, cost_max): # this computes all polynomials of order 0 to ell_max+1 and for all ell's lpns_min = lpn(ell[-1] + 1, cost_min)[0][1:] lpns_max = lpn(ell[-1] + 1, cost_max)[0][1:] - dlpns_min = lpn(ell[-1] + 1, cost_min)[1][1:] - dlpns_max = lpn(ell[-1] + 1, cost_max)[1][1:] + dlpns_min = lpn(ell[-1] + 1, cost_min, diff_n=1)[1][1:] + dlpns_max = lpn(ell[-1] + 1, cost_max, diff_n=1)[1][1:] # denominator in average dcost = cost_max - cost_min