From 5c7a6403fa20b9cceae4467a1e432907ae990b95 Mon Sep 17 00:00:00 2001 From: xzackli Date: Sun, 9 Jul 2023 01:37:05 -0400 Subject: [PATCH 01/23] generalize to EEEE --- pspy/so_cov.py | 172 ++++++++++++++++++-------------------- pspy/tests/test_so_cov.py | 14 +++- 2 files changed, 90 insertions(+), 96 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index 257efc4..42b7d40 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1083,96 +1083,83 @@ def measure_white_noise_level(var, mask, dtype=np.float64): noise_tot = np.sum(buffer) return noise_tot / mask_tot -# only works for TTTT for now -def generate_aniso_couplings(survey_name, win, var, lmax, niter=0): - """ - survey_name : list of names of the four splits - win : dict with windows with keys 'Ta', "Tb', ... - var : dict with variance maps, with keys 'Ta', 'Tb' - """ - survey_id = ["Ta", "Tb", "Tc", "Td"] - # set up some tools to work with the survey names - id2name = dict(zip(survey_id, survey_name)) +def masked_variances(survey_name, win, var): var_a, var_b, var_c, var_d = ( win['Ta'].copy(), win['Tb'].copy(), win['Tc'].copy(), win['Td'].copy()) - pixsizes = win['Ta'].data.pixsizemap() # assuming they're the same footprint var_a.data *= np.sqrt(pixsizes * var['Ta'].data) var_b.data *= np.sqrt(pixsizes * var['Tb'].data) var_c.data *= np.sqrt(pixsizes * var['Tc'].data) var_d.data *= np.sqrt(pixsizes * var['Td'].data) - - zero_var = win['Ta'].copy() - zero_var.data *= 0 - - # allows indexing with the delta function - v_a = (zero_var, var_a) - v_b = (zero_var, var_b) - v_c = (zero_var, var_c) - v_d = (zero_var, var_d) - - coupling_1 = cov_coupling_spin0( - {'Ta': win['Ta'], - 'Tb': win['Tb'], - 'Tc': win['Tc'], - 'Td': win['Td']}, - lmax, niter)['TaTcTbTd'] - - coupling_2 = cov_coupling_spin0( - {'Ta': win['Ta'], - 'Tb': win['Tb'], - 'Tc': win['Tc'], - 'Td': win['Td']}, - lmax, niter)['TaTdTbTc'] - - coupling_3 = cov_coupling_spin0( - {'Ta': win['Ta'], - 'Tb': v_b[delta2(id2name['Tb'], id2name['Td'])], - 'Tc': win['Tc'], - 'Td': v_d[delta2(id2name['Tb'], id2name['Td'])]}, - lmax, niter)['TaTcTbTd'] - - coupling_4 = cov_coupling_spin0( - {'Ta': v_a[delta2(id2name['Ta'], id2name['Tc'])], - 'Tb': win['Tb'], - 'Tc': v_c[delta2(id2name['Ta'], id2name['Tc'])], - 'Td': win['Td']}, - lmax, niter)['TaTcTbTd'] - - coupling_5 = cov_coupling_spin0( - {'Ta': win['Ta'], - 'Tb': v_a[delta2(id2name['Tb'], id2name['Tc'])], - 'Tc': v_c[delta2(id2name['Tb'], id2name['Tc'])], - 'Td': win['Td']}, - lmax, niter)['TaTdTbTc'] - - coupling_6 = cov_coupling_spin0( - {'Ta': v_a[delta2(id2name['Ta'], id2name['Td'])], - 'Tb': win['Tb'], - 'Tc': win['Tc'], - 'Td': v_d[delta2(id2name['Ta'], id2name['Td'])]}, - lmax, niter)['TaTdTbTc'] - - coupling_7 = cov_coupling_spin0( - {'Ta': v_a[delta2(id2name['Ta'], id2name['Tc'])], - 'Tb': v_b[delta2(id2name['Tb'], id2name['Td'])], - 'Tc': v_c[delta2(id2name['Ta'], id2name['Tc'])], - 'Td': v_d[delta2(id2name['Tb'], id2name['Td'])]}, - lmax, niter)['TaTcTbTd'] - - coupling_8 = cov_coupling_spin0( - {'Ta': v_a[delta2(id2name['Ta'], id2name['Td'])], - 'Tb': v_b[delta2(id2name['Tb'], id2name['Tc'])], - 'Tc': v_c[delta2(id2name['Tb'], id2name['Tc'])], - 'Td': v_d[delta2(id2name['Ta'], id2name['Td'])]}, - lmax, niter)['TaTcTbTd'] - - return [coupling_1, coupling_2, coupling_3, coupling_4, - coupling_5, coupling_6, coupling_7, coupling_8] - - -def cov_spin0_aniso1(Clth, Rl, couplings, binning_file, lmax, + return dict(zip(survey_name, (var_a, var_b, var_c, var_d))) + + + +def _same_pol(a, b): # i.e. 'Ta' and 'Tb' + return a[0] == b[0] + + +def get_zero_map(template): + c = template.copy() + c.data *= 0 + return c + +def _product_of_so_map(a, b): + c = a.copy() + c.data *= b.data + return c + +def make_weighted_variance_map(window, variance): + var = window.copy() # window^2 * variance * pixsize + var.data *= window.data + var.data *= window.data.pixsizemap() * variance.data + return var + + +def organize_covmat_products(field_names, survey_names, weighted_variances, windows): + i, j, p, q = field_names + survey = dict(zip(field_names, survey_names)) + zero_map = get_zero_map(weighted_variances[i]) + win = lambda a, b : _product_of_so_map(windows[a], windows[b]) + var = lambda a, b : weighted_variances[a] if ( + _same_pol(a, b) and survey[a] == survey[b]) else zero_map + return win, var + +def generate_aniso_couplings_TTTT(fields, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = fields + win, var = organize_covmat_products(fields, surveys, weighted_var, windows) + coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( + prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + return [ + coupling(win(i, p), win(j, q)), + coupling(win(i, q), win(j, p)), + coupling(win(i, p), var(j, q)), + coupling(win(j, q), var(i, p)), + coupling(win(i, q), var(j, p)), + coupling(win(j, p), var(i, q)), + coupling(var(i, p), var(j, q)), + coupling(var(i, q), var(j, p))] + + +def generate_aniso_couplings_EEEE(fields, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = fields + win, var = organize_covmat_products(fields, surveys, weighted_var, windows) + coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( + prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[3,:,:] / (4 * np.pi) + return [ + coupling(win(i, p), win(j, q)), + coupling(win(i, q), win(j, p)), + coupling(win(i, p), var(j, q)), + coupling(win(j, q), var(i, p)), + coupling(win(i, q), var(j, p)), + coupling(win(j, p), var(i, q)), + coupling(var(i, p), var(j, q)), + coupling(var(i, q), var(j, p))] + + +# for TTTT and EEEE +def cov_spin0_aniso_same_pol(fields, Clth, Rl, couplings, binning_file, lmax, mbb_inv_ab, mbb_inv_cd, binned_mcm=True): """Estimate the analytic covariance in the case of anisotropic noise pixel variance, using the theoretical Cls, the ratios of noise to white noise power spectra, the @@ -1181,7 +1168,8 @@ def cov_spin0_aniso1(Clth, Rl, couplings, binning_file, lmax, Parameters ---------- - + fields: list + List of fields / survey IDs, i.e. ['Ta', 'Tb', 'Tc', 'Td'] Clth: dictionary A dictionary of theoretical power spectrum (auto and cross) for the different split combinations ('TaTb' etc) @@ -1201,17 +1189,17 @@ def cov_spin0_aniso1(Clth, Rl, couplings, binning_file, lmax, specify if the mode coupling matrices are binned or not """ - + i, j, p, q = fields geom = lambda cl : symmetrize(cl, mode="geo") - cov = geom(Clth["TaTc"]) * geom(Clth["TbTd"]) * couplings[0], - cov += geom(Clth["TaTd"]) * geom(Clth["TbTc"]) * couplings[1] - cov += geom(Rl['Tb']) * geom(Rl['Td']) * geom(Clth["TaTc"]) * couplings[2] - cov += geom(Rl['Ta']) * geom(Rl['Tc']) * geom(Clth["TbTd"]) * couplings[3] - cov += geom(Rl['Tb']) * geom(Rl['Tc']) * geom(Clth["TaTd"]) * couplings[4] - cov += geom(Rl['Ta']) * geom(Rl['Td']) * geom(Clth["TbTc"]) * couplings[5] - cov += geom(Rl['Ta']) * geom(Rl['Tb']) * geom(Rl['Tc']) * geom(Rl['Td']) * couplings[6] - cov += geom(Rl['Ta']) * geom(Rl['Tb']) * geom(Rl['Tc']) * geom(Rl['Td']) * couplings[7] + cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0], + cov += geom(Clth[i+q]) * geom(Clth[j+p]) * couplings[1] + cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i +p]) * couplings[2] + cov += geom(Rl[i]) * geom(Rl[p]) * geom(Clth[j+q]) * couplings[3] + cov += geom(Rl[j]) * geom(Rl[p]) * geom(Clth[i+q]) * couplings[4] + cov += geom(Rl[i]) * geom(Rl[q]) * geom(Clth[j+p]) * couplings[5] + cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[6] + cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[7] if binned_mcm == True: analytic_cov = bin_mat(cov, binning_file, lmax) diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 5146983..e750f49 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -38,13 +38,19 @@ def test_covmat(self, verbose=False): survey_id = ["Ta", "Tb", "Tc", "Td"] survey_name = ["split_0", "split_1", "split_0", "split_1"] win = {'Ta': window, 'Tb': window, 'Tc': window, 'Td': window} - var = {'Ta': var0, 'Tb': var1, 'Tc': var0, 'Td': var1} + var = { + 'Ta': so_cov.make_weighted_variance_map(window, var0), + 'Tb': so_cov.make_weighted_variance_map(window, var1), + 'Tc': so_cov.make_weighted_variance_map(window, var0), + 'Td': so_cov.make_weighted_variance_map(window, var1) + } # generate mcm mbb_inv, Bbl = so_mcm.mcm_and_bbl_spin0(window, binning_file, lmax=lmax, type="Dl", niter=0, binned_mcm = False) - couplings = so_cov.generate_aniso_couplings(survey_name, win, var, lmax) + couplings = so_cov.generate_aniso_couplings_TTTT( + survey_id, survey_name, win, var, lmax) id2spec = {'Ta': 'dr6&pa6_f150_s1', 'Tb': 'dr6&pa6_f150_s4', 'Tc': 'dr6&pa6_f150_s1', 'Td': 'dr6&pa6_f150_s4'} @@ -71,8 +77,8 @@ def test_covmat(self, verbose=False): Clth_dict[id1 + id2] = cl_dict[key12][:lmax] Clth_dict[id1 + id2] = Clth_dict[id1 + id2][2:] - cov_e = so_cov.cov_spin0_aniso1( - Clth_dict, Rl_dict, couplings, + cov_e = so_cov.cov_spin0_aniso_same_pol( + survey_id, Clth_dict, Rl_dict, couplings, binning_file, lmax, mbb_inv, mbb_inv, binned_mcm=False) num_diag = 30 # check first 30 off-diagonals From 25b6742fd13d06428924da5e78c8e881a6842560 Mon Sep 17 00:00:00 2001 From: xzackli Date: Sun, 9 Jul 2023 14:52:50 -0400 Subject: [PATCH 02/23] EEEE test --- pspy/so_cov.py | 105 +++++++++++++++++++------------------- pspy/tests/test_so_cov.py | 70 +++++++++++++++++++++---- 2 files changed, 113 insertions(+), 62 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index 8cd333e..b1cc128 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1119,18 +1119,17 @@ def make_weighted_variance_map(window, variance): return var -def organize_covmat_products(field_names, survey_names, weighted_variances, windows): - i, j, p, q = field_names - survey = dict(zip(field_names, survey_names)) - zero_map = get_zero_map(weighted_variances[i]) +def organize_covmat_products(survey_id, survey_names, weighted_variances, windows): + survey = dict(zip(survey_id, survey_names)) + zero_map = get_zero_map(weighted_variances[survey_id[0]]) win = lambda a, b : _product_of_so_map(windows[a], windows[b]) var = lambda a, b : weighted_variances[a] if ( _same_pol(a, b) and survey[a] == survey[b]) else zero_map return win, var -def generate_aniso_couplings_TTTT(fields, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = fields - win, var = organize_covmat_products(fields, surveys, weighted_var, windows) +def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = survey_id + win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) return [ @@ -1144,11 +1143,12 @@ def generate_aniso_couplings_TTTT(fields, surveys, windows, weighted_var, lmax, coupling(var(i, q), var(j, p))] -def generate_aniso_couplings_EEEE(fields, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = fields - win, var = organize_covmat_products(fields, surveys, weighted_var, windows) +def generate_aniso_couplings_EEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = survey_id + win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( - prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[3,:,:] / (4 * np.pi) + (prod1, prod1), "", lmax, niter, "Cl", (prod2, prod2), # EXTREMELY INEFFICIENT + return_coupling_only=True)[3,:,:] / (4 * np.pi) return [ coupling(win(i, p), win(j, q)), coupling(win(i, q), win(j, p)), @@ -1160,40 +1160,10 @@ def generate_aniso_couplings_EEEE(fields, surveys, windows, weighted_var, lmax, coupling(var(i, q), var(j, p))] -# for TTTT and EEEE -def cov_spin0_aniso_same_pol(fields, Clth, Rl, couplings, binning_file, lmax, - mbb_inv_ab, mbb_inv_cd, binned_mcm=True): - """Estimate the analytic covariance in the case of anisotropic noise pixel variance, - using the theoretical Cls, the ratios of noise to white noise power spectra, the - coupling matrices estimated from the masks and variance maps, and the mode-coupling - matrices. - - Parameters - ---------- - fields: list - List of fields / survey IDs, i.e. ['Ta', 'Tb', 'Tc', 'Td'] - Clth: dictionary - A dictionary of theoretical power spectrum (auto and cross) for the different split - combinations ('TaTb' etc) - Rl: dictionary - a dictionary containing the ratios of noise power spectra to white noise - couplings: list of arrays - a list containing the eight coupling matrices involved in this covariance - binning_file: data file - a binning file with format bin low, bin high, bin mean - lmax: int - the maximum multipole to consider - mbb_inv_ab: 2d array - the inverse mode coupling matrix for the 'TaTb' power spectrum - mbb_inv_cd: 2d array - the inverse mode coupling matrix for the 'TcTd' power spectrum - binned_mcm: boolean - specify if the mode coupling matrices are binned or not - - """ - i, j, p, q = fields +# for TTTT, EEEE +def coupled_cov_aniso_same_pol(survey_id, Clth, Rl, couplings): + i, j, p, q = survey_id geom = lambda cl : symmetrize(cl, mode="geo") - cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0], cov += geom(Clth[i+q]) * geom(Clth[j+p]) * couplings[1] cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i +p]) * couplings[2] @@ -1202,12 +1172,43 @@ def cov_spin0_aniso_same_pol(fields, Clth, Rl, couplings, binning_file, lmax, cov += geom(Rl[i]) * geom(Rl[q]) * geom(Clth[j+p]) * couplings[5] cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[6] cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[7] - - if binned_mcm == True: - analytic_cov = bin_mat(cov, binning_file, lmax) - analytic_cov = mbb_inv_ab @ analytic_cov @ mbb_inv_cd.T - else: - full_analytic_cov = mbb_inv_ab @ cov @ mbb_inv_cd.T - analytic_cov = bin_mat(full_analytic_cov, binning_file, lmax) - return cov + + +def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = survey_id + win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) + coupling_TT = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( + prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + coupling_TE = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( + prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[1,:,:] / (4 * np.pi) + return [ + coupling_TE(win(i, p), win(j, q)), + coupling_TT(win(i, q), win(j, p)), + coupling_TE(win(i, p), var(j, q)), + coupling_TE(win(j, q), var(i, p)), + coupling_TE(var(i, p), var(j, q))] + +# todo: need to check arithmetic mean routine + + +# def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): +# i, j, p, q = survey_id +# win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) +# coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( +# prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[3,:,:] / (4 * np.pi) +# return [ +# coupling(win(i, p), win(j, q)), +# coupling(win(i, q), win(j, p)), +# coupling(win(j, q), var(i, p)), +# coupling(win(i, q), var(j, p))] + +# def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): +# i, j, p, q = survey_id +# win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) +# coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( +# prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[3,:,:] / (4 * np.pi) +# return [ +# coupling(win(i, p), win(j, q)), +# coupling(win(i, q), win(j, p))] + diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index e750f49..778b4b5 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -16,7 +16,8 @@ def setUp(self, verbose=False): var1 = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "ivar2_adjusted.fits")) self.cl_dict = np.load(os.path.join(TEST_DATA_PREFIX, "cl_dict.npy")) self.nl_dict = np.load(os.path.join(TEST_DATA_PREFIX, "nl_dict.npy")) - self.cov_ref = np.load(os.path.join(TEST_DATA_PREFIX, "analytic_cov_ll2.npy")) + self.cov_ref_TTTT = np.load(os.path.join(TEST_DATA_PREFIX, "cov_TTTT.npy")) + self.cov_ref_EEEE = np.load(os.path.join(TEST_DATA_PREFIX, "cov_EEEE.npy")) var0.data = 1 / var0.data var1.data = 1 / var1.data @@ -29,7 +30,7 @@ def setUp(self, verbose=False): - def test_covmat(self, verbose=False): + def test_covmat_TTTT(self, verbose=False): lmax, binning_file = self.lmax, self.binning_file window, var0, var1 = self.window, self.var0, self.var1 cl_dict, nl_dict = self.cl_dict, self.nl_dict @@ -45,10 +46,6 @@ def test_covmat(self, verbose=False): 'Td': so_cov.make_weighted_variance_map(window, var1) } - # generate mcm - mbb_inv, Bbl = so_mcm.mcm_and_bbl_spin0(window, binning_file, lmax=lmax, - type="Dl", niter=0, binned_mcm = False) - couplings = so_cov.generate_aniso_couplings_TTTT( survey_id, survey_name, win, var, lmax) @@ -77,13 +74,66 @@ def test_covmat(self, verbose=False): Clth_dict[id1 + id2] = cl_dict[key12][:lmax] Clth_dict[id1 + id2] = Clth_dict[id1 + id2][2:] - cov_e = so_cov.cov_spin0_aniso_same_pol( - survey_id, Clth_dict, Rl_dict, couplings, - binning_file, lmax, mbb_inv, mbb_inv, binned_mcm=False) + cov_e = so_cov.coupled_cov_aniso_same_pol(survey_id, Clth_dict, Rl_dict, couplings) + + num_diag = 30 # check first 30 off-diagonals + diag_max_errors = [ + np.max(np.abs(np.diag(self.cov_ref_TTTT[2:,2:] / cov_e[0], k) - 1)) + for k in range(0,num_diag) + ] + np.testing.assert_almost_equal(np.zeros(num_diag), + diag_max_errors, 7, + err_msg="aniso covmats don't match") + + def test_covmat_EEEE(self, verbose=False): + lmax, binning_file = self.lmax, self.binning_file + window, var0, var1 = self.window, self.var0, self.var1 + cl_dict, nl_dict = self.cl_dict, self.nl_dict + + # set up scenario + survey_id = ["Ea", "Eb", "Ec", "Ed"] + survey_name = ["split_0", "split_1", "split_0", "split_1"] + win = {'Ea': window, 'Eb': window, 'Ec': window, 'Ed': window} + var = { + 'Ea': so_cov.make_weighted_variance_map(window, var0), + 'Eb': so_cov.make_weighted_variance_map(window, var1), + 'Ec': so_cov.make_weighted_variance_map(window, var0), + 'Ed': so_cov.make_weighted_variance_map(window, var1) + } + + couplings = so_cov.generate_aniso_couplings_EEEE( + survey_id, survey_name, win, var, lmax) + + id2spec = {'Ea': 'dr6&pa6_f150_s1', 'Eb': 'dr6&pa6_f150_s4', + 'Ec': 'dr6&pa6_f150_s1', 'Ed': 'dr6&pa6_f150_s4'} + + white_noise = { + 'Ea': so_cov.measure_white_noise_level(var0.data, window.data), + 'Eb': so_cov.measure_white_noise_level(var1.data, window.data), + 'Ec': so_cov.measure_white_noise_level(var0.data, window.data), + 'Ed': so_cov.measure_white_noise_level(var1.data, window.data) + } + + # compute Rl, could put this in a function maybe + Clth_dict = {} + Rl_dict = {} + for name1, id1 in zip(survey_name, survey_id): + Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:lmax] + / white_noise[id1])[2:] + for name2, id2 in zip(survey_name, survey_id): + spec = id1[0] + id2[0] + key12 = spec + id2spec[id1] + id2spec[id2] + if name1 == name2: + Clth_dict[id1 + id2] = cl_dict[key12][:lmax] + else: + Clth_dict[id1 + id2] = cl_dict[key12][:lmax] + Clth_dict[id1 + id2] = Clth_dict[id1 + id2][2:] + + cov_e = so_cov.coupled_cov_aniso_same_pol(survey_id, Clth_dict, Rl_dict, couplings) num_diag = 30 # check first 30 off-diagonals diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref[2:,2:] / cov_e[0], k) - 1)) + np.max(np.abs(np.diag(self.cov_ref_EEEE[2:,2:] / cov_e[0], k) - 1)) for k in range(0,num_diag) ] np.testing.assert_almost_equal(np.zeros(num_diag), From d92d2cbc435a96b12df6200aa15c9dd06c63f950 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 10 Jul 2023 08:12:47 -0400 Subject: [PATCH 03/23] TETE support --- pspy/so_cov.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index b1cc128..65ee102 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1189,6 +1189,24 @@ def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lma coupling_TE(win(j, q), var(i, p)), coupling_TE(var(i, p), var(j, q))] + + +def symmetrized_arithmetic_mean(cl_a, cl_b): + A = (np.repeat(cl_a[:,np.newaxis], len(cl_a), 1) * + np.repeat(cl_b[np.newaxis,:], len(cl_b), 0)) + return (A + A.T) / 2 + +# for TETE +def coupled_cov_aniso_TETE(survey_id, Clth, Rl, couplings): + i, j, p, q = survey_id + geom = lambda cl : symmetrize(cl, mode="geo") + cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0] + cov += symmetrized_arithmetic_mean(Clth[i+q], Clth[j+p]) * couplings[1] + cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i+p]) * couplings[2] + cov += geom(Rl[i]) * geom(Rl[p]) * geom(Clth[j+q]) * couplings[3] + cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[4] + return cov + # todo: need to check arithmetic mean routine From f804afd534e2c90921e968d37733cf4f837f6eed Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 13 Jul 2023 13:19:21 -0400 Subject: [PATCH 04/23] TEEE and TTTE --- pspy/so_cov.py | 79 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 73 insertions(+), 6 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index 65ee102..0317981 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1164,9 +1164,9 @@ def generate_aniso_couplings_EEEE(survey_id, surveys, windows, weighted_var, lma def coupled_cov_aniso_same_pol(survey_id, Clth, Rl, couplings): i, j, p, q = survey_id geom = lambda cl : symmetrize(cl, mode="geo") - cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0], + cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0] cov += geom(Clth[i+q]) * geom(Clth[j+p]) * couplings[1] - cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i +p]) * couplings[2] + cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i+p]) * couplings[2] cov += geom(Rl[i]) * geom(Rl[p]) * geom(Clth[j+q]) * couplings[3] cov += geom(Rl[j]) * geom(Rl[p]) * geom(Clth[i+q]) * couplings[4] cov += geom(Rl[i]) * geom(Rl[q]) * geom(Clth[j+p]) * couplings[5] @@ -1190,23 +1190,90 @@ def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lma coupling_TE(var(i, p), var(j, q))] - -def symmetrized_arithmetic_mean(cl_a, cl_b): +def arith(cl_a, cl_b): A = (np.repeat(cl_a[:,np.newaxis], len(cl_a), 1) * np.repeat(cl_b[np.newaxis,:], len(cl_b), 0)) - return (A + A.T) / 2 + return A # for TETE def coupled_cov_aniso_TETE(survey_id, Clth, Rl, couplings): i, j, p, q = survey_id geom = lambda cl : symmetrize(cl, mode="geo") cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0] - cov += symmetrized_arithmetic_mean(Clth[i+q], Clth[j+p]) * couplings[1] + cov += 0.5 * (np.outer(Clth[i+q], Clth[j+p]) + np.outer(Clth[j+p], Clth[i+q])) * couplings[1] cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i+p]) * couplings[2] cov += geom(Rl[i]) * geom(Rl[p]) * geom(Clth[j+q]) * couplings[3] cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[4] return cov + +def coupled_cov_aniso_TTTE(survey_id, Clth, Rl, couplings): + i, j, p, q = survey_id + geom = lambda cl : symmetrize(cl, mode="geo") + arit = lambda cl : symmetrize(cl, mode="arithm") + cov = geom(Clth[i+p]) * arit(Clth[j+q]) * couplings[0] + cov += geom(Clth[j+p]) * arit(Clth[i+q]) * couplings[1] + cov += geom(Rl[i]) * geom(Rl[p]) * arit(Clth[j+q]) * couplings[2] + cov += geom(Rl[j]) * geom(Rl[p]) * arit(Clth[i+q]) * couplings[3] + return cov + + +def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = survey_id + win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) + coupling_TT = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( + prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + return [ + coupling_TT(win(i, p), win(j, q)), + coupling_TT(win(i, q), win(j, p)), + coupling_TT(win(j, q), var(i, p)), + coupling_TT(win(i, q), var(j, p)) + ] + + +def coupled_cov_aniso_TTEE(survey_id, Clth, Rl, couplings): + i, j, p, q = survey_id + cov = 0.5 * (np.outer(Clth[i+p], Clth[j+q]) + np.outer(Clth[j+q], Clth[i+p])) * couplings[0] + cov += 0.5 * (np.outer(Clth[i+q], Clth[j+p]) + np.outer(Clth[j+p], Clth[i+q])) * couplings[1] + return cov + + +def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = survey_id + win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) + coupling_TT = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( + prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + return [ + coupling_TT(win(i, p), win(j, q)), + coupling_TT(win(i, q), win(j, p)) + ] + + +def coupled_cov_aniso_TEEE(survey_id, Clth, Rl, couplings): + i, j, p, q = survey_id + geom = lambda cl : symmetrize(cl, mode="geo") + arit = lambda cl : symmetrize(cl, mode="arithm") + cov = geom(Clth[j+q]) * arit(Clth[i+p]) * couplings[0] + cov += geom(Clth[j+p]) * arit(Clth[i+q]) * couplings[1] + cov += geom(Rl[j]) * geom(Rl[q]) * arit(Clth[i+p]) * couplings[2] + cov += geom(Rl[j]) * geom(Rl[p]) * arit(Clth[i+q]) * couplings[3] + return cov + + +def generate_aniso_couplings_TEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): + i, j, p, q = survey_id + win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) + coupling_EE = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( + (prod1, prod1), "", lmax, niter, "Cl", (prod2, prod2), # EXTREMELY INEFFICIENT + return_coupling_only=True)[3,:,:] / (4 * np.pi) + return [ + coupling_EE(win(i, p), win(j, q)), + coupling_EE(win(i, q), win(j, p)), + coupling_EE(win(i, p), var(j, q)), + coupling_EE(win(i, q), var(j, p)) + ] + + # todo: need to check arithmetic mean routine From 46c85f1355001a12688d5c36dcb50a0bfbb085d3 Mon Sep 17 00:00:00 2001 From: xzackli Date: Tue, 18 Jul 2023 14:12:06 -0400 Subject: [PATCH 05/23] fix test index --- pspy/tests/test_so_cov.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 778b4b5..be532f5 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -78,7 +78,7 @@ def test_covmat_TTTT(self, verbose=False): num_diag = 30 # check first 30 off-diagonals diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref_TTTT[2:,2:] / cov_e[0], k) - 1)) + np.max(np.abs(np.diag(self.cov_ref_TTTT[2:,2:] / cov_e, k) - 1)) for k in range(0,num_diag) ] np.testing.assert_almost_equal(np.zeros(num_diag), @@ -133,7 +133,7 @@ def test_covmat_EEEE(self, verbose=False): num_diag = 30 # check first 30 off-diagonals diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref_EEEE[2:,2:] / cov_e[0], k) - 1)) + np.max(np.abs(np.diag(self.cov_ref_EEEE[2:,2:] / cov_e, k) - 1)) for k in range(0,num_diag) ] np.testing.assert_almost_equal(np.zeros(num_diag), From a26b7cb2d1c68a2f49e4bb746597f7736a0ad884 Mon Sep 17 00:00:00 2001 From: xzackli Date: Fri, 28 Jul 2023 07:28:57 -0400 Subject: [PATCH 06/23] fix survey names in test --- pspy/tests/test_so_cov.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index be532f5..2b2e5a9 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -37,7 +37,7 @@ def test_covmat_TTTT(self, verbose=False): # set up scenario survey_id = ["Ta", "Tb", "Tc", "Td"] - survey_name = ["split_0", "split_1", "split_0", "split_1"] + survey_name = ["s1", "s4", "s1", "s4"] win = {'Ta': window, 'Tb': window, 'Tc': window, 'Td': window} var = { 'Ta': so_cov.make_weighted_variance_map(window, var0), @@ -92,7 +92,7 @@ def test_covmat_EEEE(self, verbose=False): # set up scenario survey_id = ["Ea", "Eb", "Ec", "Ed"] - survey_name = ["split_0", "split_1", "split_0", "split_1"] + survey_name = ["s1", "s4", "s1", "s4"] win = {'Ea': window, 'Eb': window, 'Ec': window, 'Ed': window} var = { 'Ea': so_cov.make_weighted_variance_map(window, var0), From e55ff4962b70e7c72ef9d0b76a73a3ba86cf19fc Mon Sep 17 00:00:00 2001 From: xzackli Date: Fri, 28 Jul 2023 08:47:25 -0400 Subject: [PATCH 07/23] always use var first, then window in args --- pspy/so_cov.py | 2 +- pspy/tests/test_so_cov.py | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index 0317981..4556dc1 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1112,7 +1112,7 @@ def _product_of_so_map(a, b): c.data *= b.data return c -def make_weighted_variance_map(window, variance): +def make_weighted_variance_map(variance, window): var = window.copy() # window^2 * variance * pixsize var.data *= window.data var.data *= window.data.pixsizemap() * variance.data diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 2b2e5a9..6ee9794 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -40,10 +40,10 @@ def test_covmat_TTTT(self, verbose=False): survey_name = ["s1", "s4", "s1", "s4"] win = {'Ta': window, 'Tb': window, 'Tc': window, 'Td': window} var = { - 'Ta': so_cov.make_weighted_variance_map(window, var0), - 'Tb': so_cov.make_weighted_variance_map(window, var1), - 'Tc': so_cov.make_weighted_variance_map(window, var0), - 'Td': so_cov.make_weighted_variance_map(window, var1) + 'Ta': so_cov.make_weighted_variance_map(var0, window), + 'Tb': so_cov.make_weighted_variance_map(var1, window), + 'Tc': so_cov.make_weighted_variance_map(var0, window), + 'Td': so_cov.make_weighted_variance_map(var1, window) } couplings = so_cov.generate_aniso_couplings_TTTT( @@ -95,10 +95,10 @@ def test_covmat_EEEE(self, verbose=False): survey_name = ["s1", "s4", "s1", "s4"] win = {'Ea': window, 'Eb': window, 'Ec': window, 'Ed': window} var = { - 'Ea': so_cov.make_weighted_variance_map(window, var0), - 'Eb': so_cov.make_weighted_variance_map(window, var1), - 'Ec': so_cov.make_weighted_variance_map(window, var0), - 'Ed': so_cov.make_weighted_variance_map(window, var1) + 'Ea': so_cov.make_weighted_variance_map(var0, window), + 'Eb': so_cov.make_weighted_variance_map(var1, window), + 'Ec': so_cov.make_weighted_variance_map(var0, window), + 'Ed': so_cov.make_weighted_variance_map(var1, window) } couplings = so_cov.generate_aniso_couplings_EEEE( From 900431ea27093991505d73f113a636e9317412f5 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 12:37:32 -0700 Subject: [PATCH 08/23] utility coupling functions --- pspy/so_cov.py | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index 4556dc1..a725b94 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1127,11 +1127,26 @@ def organize_covmat_products(survey_id, survey_names, weighted_variances, window _same_pol(a, b) and survey[a] == survey[b]) else zero_map return win, var + +def coupling_00(m1, m2, lmax, niter=0): + return so_mcm.mcm_and_bbl_spin0( + m1, "", lmax, niter, "Cl", m2, return_coupling_only=True) / (4 * np.pi) + +def coupling_02(m1, m2, lmax, niter=0): + return so_mcm.mcm_and_bbl_spin0and2( + m1, "", lmax, niter, "Cl", m2, return_coupling_only=True)[1,:,:] / (4 * np.pi) + +# ++ coupling +def coupling_PP(m1, m2, lmax, niter=0): + return so_mcm.mcm_and_bbl_spin0and2( + (m1, m1), "", lmax, niter, "Cl", (m2, m2), # EXTREMELY INEFFICIENT + return_coupling_only=True)[3,:,:] / (4 * np.pi) + + def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( - prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + coupling = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) return [ coupling(win(i, p), win(j, q)), coupling(win(i, q), win(j, p)), @@ -1146,9 +1161,7 @@ def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lma def generate_aniso_couplings_EEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( - (prod1, prod1), "", lmax, niter, "Cl", (prod2, prod2), # EXTREMELY INEFFICIENT - return_coupling_only=True)[3,:,:] / (4 * np.pi) + coupling = lambda m1, m2 : coupling_PP(m1, m2, lmax, niter=niter) return [ coupling(win(i, p), win(j, q)), coupling(win(i, q), win(j, p)), @@ -1178,10 +1191,8 @@ def coupled_cov_aniso_same_pol(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( - prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) - coupling_TE = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( - prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[1,:,:] / (4 * np.pi) + coupling_TT = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) + coupling_TE = lambda m1, m2 : coupling_02(m1, m2, lmax, niter=niter) return [ coupling_TE(win(i, p), win(j, q)), coupling_TT(win(i, q), win(j, p)), @@ -1221,8 +1232,7 @@ def coupled_cov_aniso_TTTE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( - prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + coupling_TT = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) return [ coupling_TT(win(i, p), win(j, q)), coupling_TT(win(i, q), win(j, p)), @@ -1241,8 +1251,7 @@ def coupled_cov_aniso_TTEE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0( - prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True) / (4 * np.pi) + coupling_TT = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) return [ coupling_TT(win(i, p), win(j, q)), coupling_TT(win(i, q), win(j, p)) @@ -1263,9 +1272,7 @@ def coupled_cov_aniso_TEEE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_EE = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( - (prod1, prod1), "", lmax, niter, "Cl", (prod2, prod2), # EXTREMELY INEFFICIENT - return_coupling_only=True)[3,:,:] / (4 * np.pi) + coupling_EE = lambda m1, m2 : coupling_PP(m1, m2, lmax, niter=niter) return [ coupling_EE(win(i, p), win(j, q)), coupling_EE(win(i, q), win(j, p)), From af11108c617b68747a0eab85ff7c732687ef925c Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 13:55:14 -0700 Subject: [PATCH 09/23] dedicated coupling matrix utility functions --- pspy/mcm_fortran/mcm_fortran.f90 | 189 +++++++++++++++++++++++++++++++ pspy/so_cov.py | 36 +++--- pspy/so_mcm.py | 69 +++++++++++ pspy/tests/test_so_cov.py | 4 +- 4 files changed, 274 insertions(+), 24 deletions(-) diff --git a/pspy/mcm_fortran/mcm_fortran.f90 b/pspy/mcm_fortran/mcm_fortran.f90 index 92b4d61..686d4d8 100644 --- a/pspy/mcm_fortran/mcm_fortran.f90 +++ b/pspy/mcm_fortran/mcm_fortran.f90 @@ -337,5 +337,194 @@ subroutine binning_matrix(mcm, binLo, binHi, binsize, bbl, doDl) end subroutine +! sub-blocks of the coupling matrix + +subroutine calc_coupling_elem_spin2_pp(wcl, l1, l2, elem) + implicit none + integer, intent(in) :: l1, l2 + real(8), intent(in) :: wcl(:) + real(8), intent(inout):: elem + real(8) :: thrcof1(2*size(wcl)), l1f, l2f + integer :: info, l3, wlmin, wlmax, i, lmax + lmax = size(wcl)-1 ! wcl starts at 0 + + call drc3jj(dble(l1), dble(l2), -2d0, 2d0, l1f, l2f, thrcof1, size(thrcof1), info) + + wlmin = int(l1f) + wlmax = min(lmax, int(l2f)) + elem = 0 + do l3 = wlmin, wlmax + i = l3 - wlmin + 1 + elem = elem + wcl(l3 + 1) * thrcof1(i)**2 * (1 + (-1)**(l1 + l2 + l3)) / 2 + end do +end subroutine + + +subroutine calc_coupling_elem_spin02(wcl, l1, l2, elem) + implicit none + integer, intent(in) :: l1, l2 + real(8), intent(in) :: wcl(:) + real(8), intent(inout):: elem + real(8) :: thrcof0(2*size(wcl)), thrcof1(2*size(wcl)), l1f, l2f + integer :: info, l3, wlmin, wlmax, i, lmax + lmax = size(wcl)-1 ! wcl starts at 0 + + call drc3jj(dble(l1), dble(l2), 0d0, 0d0, l1f, l2f, thrcof0, size(thrcof0), info) + call drc3jj(dble(l1), dble(l2), -2d0, 2d0, l1f, l2f, thrcof1, size(thrcof1), info) + + wlmin = int(l1f) + wlmax = min(lmax, int(l2f)) + elem = 0 + do l3 = wlmin, wlmax + i = l3 - wlmin + 1 + elem = elem + wcl(l3 + 1) * thrcof0(i) * thrcof1(i) + end do +end subroutine + + +subroutine calc_coupling_elem_spin2_mm(wcl, l1, l2, elem) + implicit none + integer, intent(in) :: l1, l2 + real(8), intent(in) :: wcl(:) + real(8), intent(inout):: elem + real(8) :: thrcof1(2*size(wcl)), l1f, l2f + integer :: info, l3, wlmin, wlmax, i, lmax + lmax = size(wcl)-1 ! wcl starts at 0 + + call drc3jj(dble(l1), dble(l2), -2d0, 2d0, l1f, l2f, thrcof1, size(thrcof1), info) + + wlmin = int(l1f) + wlmax = min(lmax, int(l2f)) + elem = 0 + do l3 = wlmin, wlmax + i = l3 - wlmin + 1 + elem = elem + wcl(l3 + 1) * thrcof1(i)**2 * (1 - (-1)**(l1 + l2 + l3)) / 2 + end do +end subroutine + + +subroutine calc_coupling_spin02(wcl, l_exact, l_band, l_toeplitz, coupling) + implicit none + real(8), intent(in) :: wcl(:) + integer, intent(in) :: l_exact , l_band, l_toeplitz + real(8), intent(inout) :: coupling(:,:) + integer :: l1, l2, nlmax, lmax_band + + nlmax = size(coupling,1) - 1 + + !$omp parallel do private(l2, l1) schedule(dynamic) + do l1 = 2, min(nlmax, l_exact) + do l2 = l1, nlmax + call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_exact .lt. nlmax) then + !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) + do l1 = l_exact + 1, l_toeplitz + if (l1 .lt. l_toeplitz) then + lmax_band = min(l1 + l_band, nlmax) + else + lmax_band = nlmax + end if + + do l2 = l1, lmax_band + call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_toeplitz .lt. nlmax) then + !$omp parallel do + do l1 = l_toeplitz + 1, nlmax + call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1-1, l1-1)) + end do + end if + end if + +end subroutine + + +subroutine calc_coupling_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupling) + implicit none + real(8), intent(in) :: wcl(:) + integer, intent(in) :: l_exact , l_band, l_toeplitz + real(8), intent(inout) :: coupling(:,:) + integer :: l1, l2, nlmax, lmax_band + + nlmax = size(coupling,1) - 1 + + !$omp parallel do private(l2, l1) schedule(dynamic) + do l1 = 2, min(nlmax, l_exact) + do l2 = l1, nlmax + call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_exact .lt. nlmax) then + !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) + do l1 = l_exact + 1, l_toeplitz + if (l1 .lt. l_toeplitz) then + lmax_band = min(l1 + l_band, nlmax) + else + lmax_band = nlmax + end if + + do l2 = l1, lmax_band + call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_toeplitz .lt. nlmax) then + !$omp parallel do + do l1 = l_toeplitz + 1, nlmax + call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1-1, l1-1)) + end do + end if + end if + +end subroutine + + +subroutine calc_coupling_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupling) + implicit none + real(8), intent(in) :: wcl(:) + integer, intent(in) :: l_exact , l_band, l_toeplitz + real(8), intent(inout) :: coupling(:,:) + integer :: l1, l2, nlmax, lmax_band + + nlmax = size(coupling,1) - 1 + + !$omp parallel do private(l2, l1) schedule(dynamic) + do l1 = 2, min(nlmax, l_exact) + do l2 = l1, nlmax + call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_exact .lt. nlmax) then + !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) + do l1 = l_exact + 1, l_toeplitz + if (l1 .lt. l_toeplitz) then + lmax_band = min(l1 + l_band, nlmax) + else + lmax_band = nlmax + end if + + do l2 = l1, lmax_band + call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_toeplitz .lt. nlmax) then + !$omp parallel do + do l1 = l_toeplitz + 1, nlmax + call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1-1, l1-1)) + end do + end if + end if + +end subroutine + + end module diff --git a/pspy/so_cov.py b/pspy/so_cov.py index a725b94..de52e1a 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1128,25 +1128,11 @@ def organize_covmat_products(survey_id, survey_names, weighted_variances, window return win, var -def coupling_00(m1, m2, lmax, niter=0): - return so_mcm.mcm_and_bbl_spin0( - m1, "", lmax, niter, "Cl", m2, return_coupling_only=True) / (4 * np.pi) - -def coupling_02(m1, m2, lmax, niter=0): - return so_mcm.mcm_and_bbl_spin0and2( - m1, "", lmax, niter, "Cl", m2, return_coupling_only=True)[1,:,:] / (4 * np.pi) - -# ++ coupling -def coupling_PP(m1, m2, lmax, niter=0): - return so_mcm.mcm_and_bbl_spin0and2( - (m1, m1), "", lmax, niter, "Cl", (m2, m2), # EXTREMELY INEFFICIENT - return_coupling_only=True)[3,:,:] / (4 * np.pi) - - def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) + coupling = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling(win(i, p), win(j, q)), coupling(win(i, q), win(j, p)), @@ -1161,7 +1147,8 @@ def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lma def generate_aniso_couplings_EEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda m1, m2 : coupling_PP(m1, m2, lmax, niter=niter) + coupling = lambda m1, m2 : so_mcm.coupling_matrix("++", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling(win(i, p), win(j, q)), coupling(win(i, q), win(j, p)), @@ -1191,8 +1178,10 @@ def coupled_cov_aniso_same_pol(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) - coupling_TE = lambda m1, m2 : coupling_02(m1, m2, lmax, niter=niter) + coupling_TT = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) + coupling_TE = lambda m1, m2 : so_mcm.coupling_matrix("02", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_TE(win(i, p), win(j, q)), coupling_TT(win(i, q), win(j, p)), @@ -1232,7 +1221,8 @@ def coupled_cov_aniso_TTTE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) + coupling_TT = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_TT(win(i, p), win(j, q)), coupling_TT(win(i, q), win(j, p)), @@ -1251,7 +1241,8 @@ def coupled_cov_aniso_TTEE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : coupling_00(m1, m2, lmax, niter=niter) + coupling_TT = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_TT(win(i, p), win(j, q)), coupling_TT(win(i, q), win(j, p)) @@ -1272,7 +1263,8 @@ def coupled_cov_aniso_TEEE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_EE = lambda m1, m2 : coupling_PP(m1, m2, lmax, niter=niter) + coupling_EE = lambda m1, m2 : so_mcm.coupling_matrix("++", win1=m1, win2=m2, + lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_EE(win(i, p), win(j, q)), coupling_EE(win(i, q), win(j, p)), diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index 2163860..f29d14b 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -470,3 +470,72 @@ def read_coupling(prefix, spin_pairs=None): Bbl = np.load(prefix + "_Bbl.npy") return mode_coupling_inv, Bbl + + +def coupling_matrix( + mode, + win1, + lmax, + niter, + win2=None, + input_alm=False, + l_exact=None, + l_toep=None, + l_band=None, + l3_pad=2000): + + """Get the coupling matrix for spin0 fields + Parameters + ---------- + win1: so_map (or alm) + the window function of survey 1, if input_alm=True, expect wlm1 + lmax: integer + the maximum multipole to consider for the spectra computation + win2: so_map (or alm) + the window function of survey 2, if input_alm=True, expect wlm2 + niter: int + specify the number of iteration in map2alm + l_toep: int + l_band: int + l_exact: int + """ + + if input_alm == False: + l_max_limit = win1.get_lmax_limit() + if lmax > l_max_limit: raise ValueError("the requested lmax is too high with respect to the map pixellisation") + maxl = np.minimum(lmax + l3_pad, l_max_limit) + + win1 = sph_tools.map2alm(win1, niter=niter, lmax=maxl) + if win2 is not None: + win2 = sph_tools.map2alm(win2, niter=niter, lmax=maxl) + + if win2 is None: + wcl = hp.alm2cl(win1) + else: + wcl = hp.alm2cl(win1, win2) + + l = np.arange(len(wcl)) + wcl *= (2 * l + 1) + + + mcm = np.zeros((lmax, lmax)) + + if l_toep is None: l_toep = lmax + if l_band is None: l_band = lmax + if l_exact is None: l_exact = lmax + + if mode == "00": + mcm_fortran.calc_coupling_spin0(wcl, l_exact, l_band, l_toep, mcm.T) + elif mode == "02": + mcm_fortran.calc_coupling_spin02(wcl, l_exact, l_band, l_toep, mcm.T) + elif mode == "++": + mcm_fortran.calc_coupling_spin2_pp(wcl, l_exact, l_band, l_toep, mcm.T) + elif mode == "--": + mcm_fortran.calc_coupling_spin2_mm(wcl, l_exact, l_band, l_toep, mcm.T) + + if l_toep < lmax: + mcm = format_toepliz_fortran2(mcm, l_toep, l_exact, lmax) + + mcm_fortran.fill_upper(mcm.T) + + return mcm[:lmax - 2, :lmax - 2] diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 6ee9794..0812723 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -19,8 +19,8 @@ def setUp(self, verbose=False): self.cov_ref_TTTT = np.load(os.path.join(TEST_DATA_PREFIX, "cov_TTTT.npy")) self.cov_ref_EEEE = np.load(os.path.join(TEST_DATA_PREFIX, "cov_EEEE.npy")) - var0.data = 1 / var0.data - var1.data = 1 / var1.data + var0.data = np.reciprocal(var0.data,where=(var0.data!=0)) + var1.data = np.reciprocal(var1.data,where=(var1.data!=0)) var0.data[np.isinf(var0.data)] = 0. var1.data[np.isinf(var1.data)] = 0. self.var0, self.var1 = var0, var1 From 6846d310136383bf101fc01a2a1651990f6f02d5 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 14:10:51 -0700 Subject: [PATCH 10/23] clean up coupling matrix block, restore test inverse --- pspy/so_cov.py | 14 +++++++------- pspy/so_mcm.py | 6 ++++-- pspy/tests/test_so_cov.py | 4 ++-- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index de52e1a..e1e9e24 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1131,7 +1131,7 @@ def organize_covmat_products(survey_id, survey_names, weighted_variances, window def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + coupling = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling(win(i, p), win(j, q)), @@ -1147,7 +1147,7 @@ def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lma def generate_aniso_couplings_EEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda m1, m2 : so_mcm.coupling_matrix("++", win1=m1, win2=m2, + coupling = lambda m1, m2 : so_mcm.coupling_block("++", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling(win(i, p), win(j, q)), @@ -1178,9 +1178,9 @@ def coupled_cov_aniso_same_pol(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + coupling_TT = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) - coupling_TE = lambda m1, m2 : so_mcm.coupling_matrix("02", win1=m1, win2=m2, + coupling_TE = lambda m1, m2 : so_mcm.coupling_block("02", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_TE(win(i, p), win(j, q)), @@ -1221,7 +1221,7 @@ def coupled_cov_aniso_TTTE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + coupling_TT = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_TT(win(i, p), win(j, q)), @@ -1241,7 +1241,7 @@ def coupled_cov_aniso_TTEE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : so_mcm.coupling_matrix("00", win1=m1, win2=m2, + coupling_TT = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_TT(win(i, p), win(j, q)), @@ -1263,7 +1263,7 @@ def coupled_cov_aniso_TEEE(survey_id, Clth, Rl, couplings): def generate_aniso_couplings_TEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): i, j, p, q = survey_id win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_EE = lambda m1, m2 : so_mcm.coupling_matrix("++", win1=m1, win2=m2, + coupling_EE = lambda m1, m2 : so_mcm.coupling_block("++", win1=m1, win2=m2, lmax=lmax, niter=niter) / (4 * np.pi) return [ coupling_EE(win(i, p), win(j, q)), diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index f29d14b..2a53541 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -472,7 +472,7 @@ def read_coupling(prefix, spin_pairs=None): return mode_coupling_inv, Bbl -def coupling_matrix( +def coupling_block( mode, win1, lmax, @@ -484,9 +484,11 @@ def coupling_matrix( l_band=None, l3_pad=2000): - """Get the coupling matrix for spin0 fields + """Get the coupling matrix for two fields Parameters ---------- + mode: string + one of "00", "02", "++" or "--" win1: so_map (or alm) the window function of survey 1, if input_alm=True, expect wlm1 lmax: integer diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 0812723..6ee9794 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -19,8 +19,8 @@ def setUp(self, verbose=False): self.cov_ref_TTTT = np.load(os.path.join(TEST_DATA_PREFIX, "cov_TTTT.npy")) self.cov_ref_EEEE = np.load(os.path.join(TEST_DATA_PREFIX, "cov_EEEE.npy")) - var0.data = np.reciprocal(var0.data,where=(var0.data!=0)) - var1.data = np.reciprocal(var1.data,where=(var1.data!=0)) + var0.data = 1 / var0.data + var1.data = 1 / var1.data var0.data[np.isinf(var0.data)] = 0. var1.data[np.isinf(var1.data)] = 0. self.var0, self.var1 = var0, var1 From d36b5c0cced3ad4665919f866a034bbb0a9edaea Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 14:25:42 -0700 Subject: [PATCH 11/23] further move away from pspy default block in 00 --- pspy/mcm_fortran/mcm_fortran.f90 | 48 ++++++++++++++++++++++++++++++-- pspy/so_mcm.py | 8 +++--- 2 files changed, 49 insertions(+), 7 deletions(-) diff --git a/pspy/mcm_fortran/mcm_fortran.f90 b/pspy/mcm_fortran/mcm_fortran.f90 index 686d4d8..73dca15 100644 --- a/pspy/mcm_fortran/mcm_fortran.f90 +++ b/pspy/mcm_fortran/mcm_fortran.f90 @@ -403,7 +403,49 @@ subroutine calc_coupling_elem_spin2_mm(wcl, l1, l2, elem) end subroutine -subroutine calc_coupling_spin02(wcl, l_exact, l_band, l_toeplitz, coupling) +! unlike calc_coupling_spin0, this keeps ell=0,1,lmax +subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling) + implicit none + real(8), intent(in) :: wcl(:) + integer, intent(in) :: l_exact , l_band, l_toeplitz + real(8), intent(inout) :: coupling(:,:) + integer :: l1, l2, nlmax, lmax_band + + nlmax = size(coupling,1) - 1 + + !$omp parallel do private(l2, l1) schedule(dynamic) + do l1 = 2, min(nlmax, l_exact) + do l2 = l1, nlmax + call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_exact .lt. nlmax) then + !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) + do l1 = l_exact + 1, l_toeplitz + if (l1 .lt. l_toeplitz) then + lmax_band = min(l1 + l_band, nlmax) + else + lmax_band = nlmax + end if + + do l2 = l1, lmax_band + call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1)) + end do + end do + + if (l_toeplitz .lt. nlmax) then + !$omp parallel do + do l1 = l_toeplitz + 1, nlmax + call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1-1, l1-1)) + end do + end if + end if + +end subroutine + + +subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling) implicit none real(8), intent(in) :: wcl(:) integer, intent(in) :: l_exact , l_band, l_toeplitz @@ -444,7 +486,7 @@ subroutine calc_coupling_spin02(wcl, l_exact, l_band, l_toeplitz, coupling) end subroutine -subroutine calc_coupling_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupling) +subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupling) implicit none real(8), intent(in) :: wcl(:) integer, intent(in) :: l_exact , l_band, l_toeplitz @@ -485,7 +527,7 @@ subroutine calc_coupling_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupling) end subroutine -subroutine calc_coupling_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupling) +subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupling) implicit none real(8), intent(in) :: wcl(:) integer, intent(in) :: l_exact , l_band, l_toeplitz diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index 2a53541..c294d98 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -527,13 +527,13 @@ def coupling_block( if l_exact is None: l_exact = lmax if mode == "00": - mcm_fortran.calc_coupling_spin0(wcl, l_exact, l_band, l_toep, mcm.T) + mcm_fortran.calc_coupling_block_spin0(wcl, l_exact, l_band, l_toep, mcm.T) elif mode == "02": - mcm_fortran.calc_coupling_spin02(wcl, l_exact, l_band, l_toep, mcm.T) + mcm_fortran.calc_coupling_block_spin02(wcl, l_exact, l_band, l_toep, mcm.T) elif mode == "++": - mcm_fortran.calc_coupling_spin2_pp(wcl, l_exact, l_band, l_toep, mcm.T) + mcm_fortran.calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toep, mcm.T) elif mode == "--": - mcm_fortran.calc_coupling_spin2_mm(wcl, l_exact, l_band, l_toep, mcm.T) + mcm_fortran.calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toep, mcm.T) if l_toep < lmax: mcm = format_toepliz_fortran2(mcm, l_toep, l_exact, lmax) From e5caf539419ce15b7844ef2f1316823c69861170 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 15:59:19 -0700 Subject: [PATCH 12/23] change ell limits to [0, lmax] --- pspy/mcm_fortran/mcm_fortran.f90 | 33 ++++++++++++++++++-------------- pspy/so_mcm.py | 10 +++++++--- 2 files changed, 26 insertions(+), 17 deletions(-) diff --git a/pspy/mcm_fortran/mcm_fortran.f90 b/pspy/mcm_fortran/mcm_fortran.f90 index 73dca15..6d8e9c1 100644 --- a/pspy/mcm_fortran/mcm_fortran.f90 +++ b/pspy/mcm_fortran/mcm_fortran.f90 @@ -414,9 +414,9 @@ subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling) nlmax = size(coupling,1) - 1 !$omp parallel do private(l2, l1) schedule(dynamic) - do l1 = 2, min(nlmax, l_exact) + do l1 = 0, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do @@ -430,14 +430,14 @@ subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling) end if do l2 = l1, lmax_band - call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1-1, l1-1)) + call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1+1, l1+1)) end do end if end if @@ -457,9 +457,11 @@ subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling !$omp parallel do private(l2, l1) schedule(dynamic) do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do + coupling(1,1) = 1 + coupling(2,2) = 1 if (l_exact .lt. nlmax) then !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) @@ -471,14 +473,14 @@ subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling end if do l2 = l1, lmax_band - call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1-1, l1-1)) + call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1+1, l1+1)) end do end if end if @@ -498,9 +500,11 @@ subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupli !$omp parallel do private(l2, l1) schedule(dynamic) do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do + coupling(1,1) = 1 + coupling(2,2) = 1 if (l_exact .lt. nlmax) then !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) @@ -512,14 +516,14 @@ subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupli end if do l2 = l1, lmax_band - call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1-1, l1-1)) + call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1+1, l1+1)) end do end if end if @@ -539,9 +543,11 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli !$omp parallel do private(l2, l1) schedule(dynamic) do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do + coupling(1,1) = 1 + coupling(2,2) = 1 if (l_exact .lt. nlmax) then !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) @@ -553,14 +559,14 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli end if do l2 = l1, lmax_band - call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1)) + call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1+1, l2+1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1-1, l1-1)) + call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1+1, l1+1)) end do end if end if @@ -569,4 +575,3 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli end module - diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index c294d98..b46bfcd 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -482,7 +482,8 @@ def coupling_block( l_exact=None, l_toep=None, l_band=None, - l3_pad=2000): + l3_pad=2000, + legacy_thibaut=True): """Get the coupling matrix for two fields Parameters @@ -520,7 +521,7 @@ def coupling_block( wcl *= (2 * l + 1) - mcm = np.zeros((lmax, lmax)) + mcm = np.zeros((lmax+1, lmax+1)) if l_toep is None: l_toep = lmax if l_band is None: l_band = lmax @@ -540,4 +541,7 @@ def coupling_block( mcm_fortran.fill_upper(mcm.T) - return mcm[:lmax - 2, :lmax - 2] + if legacy_thibaut: + return mcm[2:-1, 2:-1] + else: + return mcm From 685116a453beb0f599710b2891ac7fd9145f1911 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 17:37:01 -0700 Subject: [PATCH 13/23] extend ells to [0,lmax] --- pspy/so_mcm.py | 2 +- pspy/tests/test_so_cov.py | 24 +++++++++++------------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index b46bfcd..4f6701a 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -483,7 +483,7 @@ def coupling_block( l_toep=None, l_band=None, l3_pad=2000, - legacy_thibaut=True): + legacy_thibaut=False): """Get the coupling matrix for two fields Parameters diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 6ee9794..06c5dc9 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -10,7 +10,7 @@ class SOCovmatTests(unittest.TestCase): def setUp(self, verbose=False): # read in downgraded windows and ivar - self.lmax = 1350 + self.lmax = 1349 self.window = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "window_adjusted.fits")) var0 = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "ivar1_adjusted.fits")) var1 = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "ivar2_adjusted.fits")) @@ -63,22 +63,21 @@ def test_covmat_TTTT(self, verbose=False): Clth_dict = {} Rl_dict = {} for name1, id1 in zip(survey_name, survey_id): - Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:lmax] - / white_noise[id1])[2:] + Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:(lmax+1)] + / white_noise[id1]) for name2, id2 in zip(survey_name, survey_id): spec = id1[0] + id2[0] key12 = spec + id2spec[id1] + id2spec[id2] if name1 == name2: - Clth_dict[id1 + id2] = cl_dict[key12][:lmax] + Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] else: - Clth_dict[id1 + id2] = cl_dict[key12][:lmax] - Clth_dict[id1 + id2] = Clth_dict[id1 + id2][2:] + Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] cov_e = so_cov.coupled_cov_aniso_same_pol(survey_id, Clth_dict, Rl_dict, couplings) num_diag = 30 # check first 30 off-diagonals diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref_TTTT[2:,2:] / cov_e, k) - 1)) + np.max(np.abs(np.diag(self.cov_ref_TTTT / cov_e, k) - 1)) for k in range(0,num_diag) ] np.testing.assert_almost_equal(np.zeros(num_diag), @@ -118,22 +117,21 @@ def test_covmat_EEEE(self, verbose=False): Clth_dict = {} Rl_dict = {} for name1, id1 in zip(survey_name, survey_id): - Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:lmax] - / white_noise[id1])[2:] + Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:(lmax+1)] + / white_noise[id1]) for name2, id2 in zip(survey_name, survey_id): spec = id1[0] + id2[0] key12 = spec + id2spec[id1] + id2spec[id2] if name1 == name2: - Clth_dict[id1 + id2] = cl_dict[key12][:lmax] + Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] else: - Clth_dict[id1 + id2] = cl_dict[key12][:lmax] - Clth_dict[id1 + id2] = Clth_dict[id1 + id2][2:] + Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] cov_e = so_cov.coupled_cov_aniso_same_pol(survey_id, Clth_dict, Rl_dict, couplings) num_diag = 30 # check first 30 off-diagonals diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref_EEEE[2:,2:] / cov_e, k) - 1)) + np.max(np.abs(np.diag(self.cov_ref_EEEE[2:,2:] / cov_e[2:,2:], k) - 1)) for k in range(0,num_diag) ] np.testing.assert_almost_equal(np.zeros(num_diag), From 85fa14f51c89059c8860f268c7e8e575c518dfb3 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 14 Sep 2023 17:47:37 -0700 Subject: [PATCH 14/23] rename legacy parameter --- pspy/so_mcm.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index 4f6701a..ca8ce5a 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -483,7 +483,7 @@ def coupling_block( l_toep=None, l_band=None, l3_pad=2000, - legacy_thibaut=False): + legacy_ell_range=False): """Get the coupling matrix for two fields Parameters @@ -541,7 +541,7 @@ def coupling_block( mcm_fortran.fill_upper(mcm.T) - if legacy_thibaut: + if legacy_ell_range: return mcm[2:-1, 2:-1] else: return mcm From ea4483a93582a99f62e31eec81f1fbfddb845984 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 11 Dec 2023 21:16:59 -0800 Subject: [PATCH 15/23] make niter optional if alms are passed --- pspy/so_mcm.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index ca8ce5a..35c8903 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -476,7 +476,7 @@ def coupling_block( mode, win1, lmax, - niter, + niter=None, win2=None, input_alm=False, l_exact=None, @@ -504,6 +504,8 @@ def coupling_block( """ if input_alm == False: + if niter is None: + raise ValueError("niter required") l_max_limit = win1.get_lmax_limit() if lmax > l_max_limit: raise ValueError("the requested lmax is too high with respect to the map pixellisation") maxl = np.minimum(lmax + l3_pad, l_max_limit) From c4ff4f9389ebd5e7b52874c86a0056041552cb31 Mon Sep 17 00:00:00 2001 From: Zachary Atkins Date: Tue, 19 Dec 2023 19:25:07 -0500 Subject: [PATCH 16/23] new analytical tf methods --- pspy/so_map_preprocessing.py | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/pspy/so_map_preprocessing.py b/pspy/so_map_preprocessing.py index 0fe34b5..1e6e623 100644 --- a/pspy/so_map_preprocessing.py +++ b/pspy/so_map_preprocessing.py @@ -1,6 +1,8 @@ import numpy as np, pylab as plt +from scipy.integrate import dblquad from pspy import flat_tools, pspy_utils from pixell import enmap, utils +from itertools import product def build_std_filter(shape, wcs, vk_mask, hk_mask, dtype=np.float64): ly, lx = enmap.laxes(shape, wcs, method="auto") @@ -68,7 +70,35 @@ def analytical_tf(map, filter, binning_file, lmax): return bin_c, tf - +def analytical_tf_vkhk(vk_mask, hk_mask, ells, dtype=np.float64): + out = np.zeros(ells.shape, dtype=dtype) + for vk, hk in product(vk_mask, hk_mask): + rk = np.sqrt(vk**2 + hk**2) + out[ells > rk] += .25 - (np.abs(np.arcsin(vk/ells[ells > rk])) + np.abs(np.arcsin(hk/ells[ells > rk]))) / (2*np.pi) + return out + +def analytical_tf_yxint(yxfilter_func, ells, dtype=np.float64): + def x(r, phi): + return r*np.cos(phi) + def y(r, phi): + return r*np.sin(phi) + def rphifilter_func(r, phi): + return yxfilter_func(y(r, phi), x(r, phi)) * r + + assert ells.ndim == 1, 'expect 1d ells' + assert np.all(ells >= 0), 'expect non-negative ells' + out = np.zeros(ells.shape, dtype=dtype) + for i, ell in enumerate(ells): + if i % 10 == 0: print(ell) + if ell == 0: + rlow = ell + else: + rlow = ell - 0.5 + rhigh = ell + 0.5 + + out[i] = dblquad(rphifilter_func, 0, 2*np.pi, rlow, rhigh)[0] / (np.pi * (rhigh**2 - rlow**2)) + + return out #def analytical_tf_old(map, binning_file, lmax, vk_mask=None, hk_mask=None): # import time From 983ba20c359ecaec0393429a3ec972a5a098b850 Mon Sep 17 00:00:00 2001 From: Zachary Atkins Date: Thu, 11 Jan 2024 10:55:51 -0500 Subject: [PATCH 17/23] more analytical tf updates --- pspy/so_map_preprocessing.py | 46 ++++++++++++++++++++++++------------ 1 file changed, 31 insertions(+), 15 deletions(-) diff --git a/pspy/so_map_preprocessing.py b/pspy/so_map_preprocessing.py index 1e6e623..3b64ad3 100644 --- a/pspy/so_map_preprocessing.py +++ b/pspy/so_map_preprocessing.py @@ -1,5 +1,5 @@ import numpy as np, pylab as plt -from scipy.integrate import dblquad +from scipy.integrate import quad from pspy import flat_tools, pspy_utils from pixell import enmap, utils from itertools import product @@ -71,32 +71,48 @@ def analytical_tf(map, filter, binning_file, lmax): return bin_c, tf def analytical_tf_vkhk(vk_mask, hk_mask, ells, dtype=np.float64): + if vk_mask is None: + vk_mask = [0, 0] + if hk_mask is None: + hk_mask = [0, 0] + assert len(vk_mask) == 2 and len(hk_mask) == 2, \ + f'{len(vk_mask)=} and {len(hk_mask)=} should both be 2' + out = np.zeros(ells.shape, dtype=dtype) for vk, hk in product(vk_mask, hk_mask): rk = np.sqrt(vk**2 + hk**2) out[ells > rk] += .25 - (np.abs(np.arcsin(vk/ells[ells > rk])) + np.abs(np.arcsin(hk/ells[ells > rk]))) / (2*np.pi) return out -def analytical_tf_yxint(yxfilter_func, ells, dtype=np.float64): - def x(r, phi): - return r*np.cos(phi) - def y(r, phi): - return r*np.sin(phi) - def rphifilter_func(r, phi): - return yxfilter_func(y(r, phi), x(r, phi)) * r +def analytical_tf_yxint(yxfilter_func, ells, x_fac=1, y_fac=1, dtype=np.float64): + def x(r, t): + return r*x_fac*np.cos(t) + def y(r, t): + return r*y_fac*np.sin(t) + + if x_fac == y_fac: + def v(r, t): + return r*x_fac + else: + def v(r, t): + return r*np.sqrt(x_fac**2 * np.sin(t)**2 + y_fac**2 * np.cos(t)**2) assert ells.ndim == 1, 'expect 1d ells' assert np.all(ells >= 0), 'expect non-negative ells' out = np.zeros(ells.shape, dtype=dtype) for i, ell in enumerate(ells): - if i % 10 == 0: print(ell) - if ell == 0: - rlow = ell + if i % 100 == 0: + print(ell) + + def vell(t): + return v(ell, t) + def rphifilter_func(t): + return yxfilter_func(y(ell, t), x(ell, t)) * vell(t) + + if x_fac == y_fac: + out[i] = quad(rphifilter_func, 0, 2*np.pi)[0] / (2*np.pi*ell*x_fac) else: - rlow = ell - 0.5 - rhigh = ell + 0.5 - - out[i] = dblquad(rphifilter_func, 0, 2*np.pi, rlow, rhigh)[0] / (np.pi * (rhigh**2 - rlow**2)) + out[i] = quad(rphifilter_func, 0, 2*np.pi)[0] / quad(vell, 0, 2*np.pi)[0] return out From c776d2fc9f51e24034f93611b416ca8f5cb43eb7 Mon Sep 17 00:00:00 2001 From: Zachary Atkins Date: Tue, 27 Feb 2024 17:28:46 -0500 Subject: [PATCH 18/23] allow kspace filter constructors to upgrade the y-direction --- pspy/so_map_preprocessing.py | 40 ++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/pspy/so_map_preprocessing.py b/pspy/so_map_preprocessing.py index 3b64ad3..739161b 100644 --- a/pspy/so_map_preprocessing.py +++ b/pspy/so_map_preprocessing.py @@ -4,11 +4,40 @@ from pixell import enmap, utils from itertools import product -def build_std_filter(shape, wcs, vk_mask, hk_mask, dtype=np.float64): +def build_std_filter(shape, wcs, vk_mask, hk_mask, dtype=np.float64, + shape_y=None, wcs_y=None): + """Construct a "standard," binary-cross style filter. + + Parameters + ---------- + shape : (ny, nx) tuple + Shape of map and therefore shape of DFT. + wcs : astropy.wcs.WCS + wcs describing the map, used in enmap.laxes + vk_mask : (kx0, kx1) tuple + Modes between kx0 and kx1 (exclusive) are set to 0, otherwise 1. + hk_mask : (ky0, ky1) tuple + Modes between ky0 and ky1 (exclusive) are set to 0, otherwise 1. + dtype : np.dtype, optional + Type of output filter array, by default np.float64. + shape_y : (ny, nx) tuple, optional + If provided, together with wcs_y, used to override the y-spacing and + shape in Fourier space in enmap.laxes, by default None. + wcs_y : astropy.wcs.WCS, optional + If provided, together with shape_y, used to override the y-spacing and + shape in Fourier space in enmap.laxes, by default None. + + Returns + ------- + (ny, nx) enmap.ndmap + The binary-cross filter. + """ ly, lx = enmap.laxes(shape, wcs, method="auto") + if shape_y is not None and wcs_y is not None: + ly, _ = enmap.laxes(shape_y, wcs_y, method="auto") ly = ly.astype(dtype) lx = lx.astype(dtype) - filter = enmap.ones(shape[-2:], wcs, dtype) + filter = enmap.ones((ly.size, lx.size), wcs, dtype) if vk_mask is not None: id_vk = np.where((lx > vk_mask[0]) & (lx < vk_mask[1])) filter[:, id_vk] = 0 @@ -27,8 +56,11 @@ def apply_std_filter(imap, filter): return filtered_map -def build_sigurd_filter(shape, wcs, lbounds, dtype=np.float64): +def build_sigurd_filter(shape, wcs, lbounds, dtype=np.float64, + shape_y=None, wcs_y=None): ly, lx = enmap.laxes(shape, wcs, method="auto") + if shape_y is not None and wcs_y is not None: + ly, _ = enmap.laxes(shape_y, wcs_y, method="auto") ly = ly.astype(dtype) lx = lx.astype(dtype) lbounds = np.asarray(lbounds) @@ -36,7 +68,7 @@ def build_sigurd_filter(shape, wcs, lbounds, dtype=np.float64): lbounds = np.broadcast_to(lbounds, (1,2)) if lbounds.ndim > 2 or lbounds.shape[-1] != 2: raise ValueError("lbounds must be [:,{ly,lx}]") - filter = enmap.ones(shape[-2:], wcs, dtype) + filter = enmap.ones((ly.size, lx.size), wcs, dtype) # Apply the filters for i , (ycut, xcut) in enumerate(lbounds): filter *= 1-(np.exp(-0.5*(ly/ycut)**2)[:,None]*np.exp(-0.5*(lx/xcut)**2)[None,:]) From 5be212dcb109ea515a2834ae9ef82fa70f4c70f1 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 1 Apr 2024 14:00:01 -0700 Subject: [PATCH 19/23] fix toeplitz issues with custom couplings --- pspy/mcm_fortran/mcm_fortran.f90 | 32 +++++++++++++------------------- pspy/so_mcm.py | 12 ++++++++---- 2 files changed, 21 insertions(+), 23 deletions(-) diff --git a/pspy/mcm_fortran/mcm_fortran.f90 b/pspy/mcm_fortran/mcm_fortran.f90 index 6d8e9c1..bfb5a62 100644 --- a/pspy/mcm_fortran/mcm_fortran.f90 +++ b/pspy/mcm_fortran/mcm_fortran.f90 @@ -414,9 +414,9 @@ subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling) nlmax = size(coupling,1) - 1 !$omp parallel do private(l2, l1) schedule(dynamic) - do l1 = 0, min(nlmax, l_exact) + do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do @@ -430,14 +430,14 @@ subroutine calc_coupling_block_spin0(wcl, l_exact, l_band, l_toeplitz, coupling) end if do l2 = l1, lmax_band - call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin0(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1+1, l1+1)) + call calc_coupling_elem_spin0(wcl, l1, l1, coupling(l1-1, l1-1)) end do end if end if @@ -457,11 +457,9 @@ subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling !$omp parallel do private(l2, l1) schedule(dynamic) do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do - coupling(1,1) = 1 - coupling(2,2) = 1 if (l_exact .lt. nlmax) then !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) @@ -473,14 +471,14 @@ subroutine calc_coupling_block_spin02(wcl, l_exact, l_band, l_toeplitz, coupling end if do l2 = l1, lmax_band - call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin02(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1+1, l1+1)) + call calc_coupling_elem_spin02(wcl, l1, l1, coupling(l1-1, l1-1)) end do end if end if @@ -500,11 +498,9 @@ subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupli !$omp parallel do private(l2, l1) schedule(dynamic) do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do - coupling(1,1) = 1 - coupling(2,2) = 1 if (l_exact .lt. nlmax) then !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) @@ -516,14 +512,14 @@ subroutine calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toeplitz, coupli end if do l2 = l1, lmax_band - call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin2_pp(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1+1, l1+1)) + call calc_coupling_elem_spin2_pp(wcl, l1, l1, coupling(l1-1, l1-1)) end do end if end if @@ -543,11 +539,9 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli !$omp parallel do private(l2, l1) schedule(dynamic) do l1 = 2, min(nlmax, l_exact) do l2 = l1, nlmax - call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do - coupling(1,1) = 1 - coupling(2,2) = 1 if (l_exact .lt. nlmax) then !$omp parallel do private(l2, l1, lmax_band) schedule(dynamic) @@ -559,14 +553,14 @@ subroutine calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toeplitz, coupli end if do l2 = l1, lmax_band - call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1+1, l2+1)) + call calc_coupling_elem_spin2_mm(wcl, l1, l2, coupling(l1-1, l2-1)) end do end do if (l_toeplitz .lt. nlmax) then !$omp parallel do do l1 = l_toeplitz + 1, nlmax - call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1+1, l1+1)) + call calc_coupling_elem_spin2_mm(wcl, l1, l1, coupling(l1-1, l1-1)) end do end if end if diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index 35c8903..cd0a864 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -522,8 +522,7 @@ def coupling_block( l = np.arange(len(wcl)) wcl *= (2 * l + 1) - - mcm = np.zeros((lmax+1, lmax+1)) + mcm = np.zeros((lmax, lmax)) if l_toep is None: l_toep = lmax if l_band is None: l_band = lmax @@ -542,8 +541,13 @@ def coupling_block( mcm = format_toepliz_fortran2(mcm, l_toep, l_exact, lmax) mcm_fortran.fill_upper(mcm.T) + + mcm_full = np.zeros((lmax+1, lmax+1)) + mcm_full[2:, 2:] = mcm[:-1,:-1] + mcm_full[0,0] = 1.0 + mcm_full[1,1] = 1.0 if legacy_ell_range: - return mcm[2:-1, 2:-1] + return mcm_full[2:-1, 2:-1] else: - return mcm + return mcm_full From 0509f451b059f4295f45435c9953ee9f29a42a43 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 1 Apr 2024 14:12:35 -0700 Subject: [PATCH 20/23] add guardrails in custom block against unsupported modes --- pspy/so_mcm.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index cd0a864..638b112 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -536,7 +536,9 @@ def coupling_block( mcm_fortran.calc_coupling_block_spin2_pp(wcl, l_exact, l_band, l_toep, mcm.T) elif mode == "--": mcm_fortran.calc_coupling_block_spin2_mm(wcl, l_exact, l_band, l_toep, mcm.T) - + else: + raise ValueError(f"mode \"{mode}\" not one of (\"00\", \"02\", \"++\", \"--\")") + if l_toep < lmax: mcm = format_toepliz_fortran2(mcm, l_toep, l_exact, lmax) From a57cf72664c8b463ddc25b7cf73d62cbe5230395 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 1 Apr 2024 14:57:51 -0700 Subject: [PATCH 21/23] add unit test for custom couplingsg --- pspy/tests/test_so_cov.py | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index 06c5dc9..ec48116 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -2,7 +2,7 @@ import tempfile, unittest, os import numpy as np -from pspy import so_cov, so_map, so_mcm, pspy_utils +from pspy import so_cov, so_map, so_mcm, pspy_utils, so_window TEST_DATA_PREFIX = os.path.join(os.path.dirname(__file__), "data", "example_newcov") @@ -138,6 +138,36 @@ def test_covmat_EEEE(self, verbose=False): diag_max_errors, 7, err_msg="aniso covmats don't match") + def test_custom_couplings(self, verbose=False): + ra0, ra1, dec0, dec1 = -25, 25, -15, 15 + res = 2 + lmax = 500 + l_exact, l_band, l_toep = 50, 100, 200 + + binary_car = so_map.car_template(1, ra0, ra1, dec0, dec1, res) + binary_car.data[:] = 0 + binary_car.data[1:-1, 1:-1] = 1 + window = so_window.create_apodization(binary_car, apo_type="C1", apo_radius_degree=2) + window2 = so_window.create_apodization(binary_car, apo_type="C2", apo_radius_degree=4) + + default_coup_02 = so_mcm.mcm_and_bbl_spin0and2( + (window, window2), "", lmax=lmax, niter=0, l_exact=l_exact, l_band=l_band, l_toep=l_toep, return_coupling_only=True) + + xi = so_mcm.coupling_block("00", window, lmax, niter=0, l_exact=l_exact, l_band=l_band, l_toep=l_toep) + np.testing.assert_almost_equal( default_coup_02[0,:,:], xi[2:-1,2:-1]) + + xi = so_mcm.coupling_block("02", window, lmax, niter=0, win2=window2, l_exact=l_exact, l_band=l_band, l_toep=l_toep) + np.testing.assert_almost_equal( default_coup_02[1,:,:], xi[2:-1,2:-1]) + + xi = so_mcm.coupling_block("02", window2, lmax, niter=0, win2=window, l_exact=l_exact, l_band=l_band, l_toep=l_toep) + np.testing.assert_almost_equal( default_coup_02[2,:,:], xi[2:-1,2:-1]) + + xi = so_mcm.coupling_block("++", window2, lmax, niter=0, win2=window2, l_exact=l_exact, l_band=l_band, l_toep=l_toep) + np.testing.assert_almost_equal( default_coup_02[3,:,:], xi[2:-1,2:-1]) + + xi = so_mcm.coupling_block("--", window2, lmax, niter=0, win2=window2, l_exact=l_exact, l_band=l_band, l_toep=l_toep) + np.testing.assert_almost_equal( default_coup_02[4,:,:], xi[2:-1,2:-1]) + if __name__ == "__main__": unittest.main() From 59f5c205b417ea35a36d7f5c66f08590f469fd1a Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 1 Apr 2024 15:17:51 -0700 Subject: [PATCH 22/23] remove aniso covmat stuff, functionality is now in PSpipe --- pspy/so_cov.py | 212 -------------------------------------- pspy/tests/test_so_cov.py | 130 ----------------------- 2 files changed, 342 deletions(-) diff --git a/pspy/so_cov.py b/pspy/so_cov.py index c9b1ddd..73e377e 100644 --- a/pspy/so_cov.py +++ b/pspy/so_cov.py @@ -1084,215 +1084,3 @@ def measure_white_noise_level(var, mask, dtype=np.float64): buffer *= var.pixsizemap().flatten() noise_tot = np.sum(buffer) return noise_tot / mask_tot - - -def masked_variances(survey_name, win, var): - var_a, var_b, var_c, var_d = ( - win['Ta'].copy(), win['Tb'].copy(), win['Tc'].copy(), win['Td'].copy()) - pixsizes = win['Ta'].data.pixsizemap() # assuming they're the same footprint - var_a.data *= np.sqrt(pixsizes * var['Ta'].data) - var_b.data *= np.sqrt(pixsizes * var['Tb'].data) - var_c.data *= np.sqrt(pixsizes * var['Tc'].data) - var_d.data *= np.sqrt(pixsizes * var['Td'].data) - return dict(zip(survey_name, (var_a, var_b, var_c, var_d))) - - - -def _same_pol(a, b): # i.e. 'Ta' and 'Tb' - return a[0] == b[0] - - -def get_zero_map(template): - c = template.copy() - c.data *= 0 - return c - -def _product_of_so_map(a, b): - c = a.copy() - c.data *= b.data - return c - -def make_weighted_variance_map(variance, window): - var = window.copy() # window^2 * variance * pixsize - var.data *= window.data - var.data *= window.data.pixsizemap() * variance.data - return var - - -def organize_covmat_products(survey_id, survey_names, weighted_variances, windows): - survey = dict(zip(survey_id, survey_names)) - zero_map = get_zero_map(weighted_variances[survey_id[0]]) - win = lambda a, b : _product_of_so_map(windows[a], windows[b]) - var = lambda a, b : weighted_variances[a] if ( - _same_pol(a, b) and survey[a] == survey[b]) else zero_map - return win, var - - -def generate_aniso_couplings_TTTT(survey_id, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = survey_id - win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - return [ - coupling(win(i, p), win(j, q)), - coupling(win(i, q), win(j, p)), - coupling(win(i, p), var(j, q)), - coupling(win(j, q), var(i, p)), - coupling(win(i, q), var(j, p)), - coupling(win(j, p), var(i, q)), - coupling(var(i, p), var(j, q)), - coupling(var(i, q), var(j, p))] - - -def generate_aniso_couplings_EEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = survey_id - win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling = lambda m1, m2 : so_mcm.coupling_block("++", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - return [ - coupling(win(i, p), win(j, q)), - coupling(win(i, q), win(j, p)), - coupling(win(i, p), var(j, q)), - coupling(win(j, q), var(i, p)), - coupling(win(i, q), var(j, p)), - coupling(win(j, p), var(i, q)), - coupling(var(i, p), var(j, q)), - coupling(var(i, q), var(j, p))] - - -# for TTTT, EEEE -def coupled_cov_aniso_same_pol(survey_id, Clth, Rl, couplings): - i, j, p, q = survey_id - geom = lambda cl : symmetrize(cl, mode="geo") - cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0] - cov += geom(Clth[i+q]) * geom(Clth[j+p]) * couplings[1] - cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i+p]) * couplings[2] - cov += geom(Rl[i]) * geom(Rl[p]) * geom(Clth[j+q]) * couplings[3] - cov += geom(Rl[j]) * geom(Rl[p]) * geom(Clth[i+q]) * couplings[4] - cov += geom(Rl[i]) * geom(Rl[q]) * geom(Clth[j+p]) * couplings[5] - cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[6] - cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[7] - return cov - - -def generate_aniso_couplings_TETE(survey_id, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = survey_id - win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - coupling_TE = lambda m1, m2 : so_mcm.coupling_block("02", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - return [ - coupling_TE(win(i, p), win(j, q)), - coupling_TT(win(i, q), win(j, p)), - coupling_TE(win(i, p), var(j, q)), - coupling_TE(win(j, q), var(i, p)), - coupling_TE(var(i, p), var(j, q))] - - -def arith(cl_a, cl_b): - A = (np.repeat(cl_a[:,np.newaxis], len(cl_a), 1) * - np.repeat(cl_b[np.newaxis,:], len(cl_b), 0)) - return A - -# for TETE -def coupled_cov_aniso_TETE(survey_id, Clth, Rl, couplings): - i, j, p, q = survey_id - geom = lambda cl : symmetrize(cl, mode="geo") - cov = geom(Clth[i+p]) * geom(Clth[j+q]) * couplings[0] - cov += 0.5 * (np.outer(Clth[i+q], Clth[j+p]) + np.outer(Clth[j+p], Clth[i+q])) * couplings[1] - cov += geom(Rl[j]) * geom(Rl[q]) * geom(Clth[i+p]) * couplings[2] - cov += geom(Rl[i]) * geom(Rl[p]) * geom(Clth[j+q]) * couplings[3] - cov += geom(Rl[i]) * geom(Rl[j]) * geom(Rl[p]) * geom(Rl[q]) * couplings[4] - return cov - - -def coupled_cov_aniso_TTTE(survey_id, Clth, Rl, couplings): - i, j, p, q = survey_id - geom = lambda cl : symmetrize(cl, mode="geo") - arit = lambda cl : symmetrize(cl, mode="arithm") - cov = geom(Clth[i+p]) * arit(Clth[j+q]) * couplings[0] - cov += geom(Clth[j+p]) * arit(Clth[i+q]) * couplings[1] - cov += geom(Rl[i]) * geom(Rl[p]) * arit(Clth[j+q]) * couplings[2] - cov += geom(Rl[j]) * geom(Rl[p]) * arit(Clth[i+q]) * couplings[3] - return cov - - -def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = survey_id - win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - return [ - coupling_TT(win(i, p), win(j, q)), - coupling_TT(win(i, q), win(j, p)), - coupling_TT(win(j, q), var(i, p)), - coupling_TT(win(i, q), var(j, p)) - ] - - -def coupled_cov_aniso_TTEE(survey_id, Clth, Rl, couplings): - i, j, p, q = survey_id - cov = 0.5 * (np.outer(Clth[i+p], Clth[j+q]) + np.outer(Clth[j+q], Clth[i+p])) * couplings[0] - cov += 0.5 * (np.outer(Clth[i+q], Clth[j+p]) + np.outer(Clth[j+p], Clth[i+q])) * couplings[1] - return cov - - -def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = survey_id - win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_TT = lambda m1, m2 : so_mcm.coupling_block("00", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - return [ - coupling_TT(win(i, p), win(j, q)), - coupling_TT(win(i, q), win(j, p)) - ] - - -def coupled_cov_aniso_TEEE(survey_id, Clth, Rl, couplings): - i, j, p, q = survey_id - geom = lambda cl : symmetrize(cl, mode="geo") - arit = lambda cl : symmetrize(cl, mode="arithm") - cov = geom(Clth[j+q]) * arit(Clth[i+p]) * couplings[0] - cov += geom(Clth[j+p]) * arit(Clth[i+q]) * couplings[1] - cov += geom(Rl[j]) * geom(Rl[q]) * arit(Clth[i+p]) * couplings[2] - cov += geom(Rl[j]) * geom(Rl[p]) * arit(Clth[i+q]) * couplings[3] - return cov - - -def generate_aniso_couplings_TEEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): - i, j, p, q = survey_id - win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) - coupling_EE = lambda m1, m2 : so_mcm.coupling_block("++", win1=m1, win2=m2, - lmax=lmax, niter=niter) / (4 * np.pi) - return [ - coupling_EE(win(i, p), win(j, q)), - coupling_EE(win(i, q), win(j, p)), - coupling_EE(win(i, p), var(j, q)), - coupling_EE(win(i, q), var(j, p)) - ] - - -# todo: need to check arithmetic mean routine - - -# def generate_aniso_couplings_TTTE(survey_id, surveys, windows, weighted_var, lmax, niter=0): -# i, j, p, q = survey_id -# win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) -# coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( -# prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[3,:,:] / (4 * np.pi) -# return [ -# coupling(win(i, p), win(j, q)), -# coupling(win(i, q), win(j, p)), -# coupling(win(j, q), var(i, p)), -# coupling(win(i, q), var(j, p))] - -# def generate_aniso_couplings_TTEE(survey_id, surveys, windows, weighted_var, lmax, niter=0): -# i, j, p, q = survey_id -# win, var = organize_covmat_products(survey_id, surveys, weighted_var, windows) -# coupling = lambda prod1, prod2 : so_mcm.mcm_and_bbl_spin0and2( -# prod1, "", lmax, niter, "Cl", prod2, return_coupling_only=True)[3,:,:] / (4 * np.pi) -# return [ -# coupling(win(i, p), win(j, q)), -# coupling(win(i, q), win(j, p))] - diff --git a/pspy/tests/test_so_cov.py b/pspy/tests/test_so_cov.py index ec48116..a7da0e9 100644 --- a/pspy/tests/test_so_cov.py +++ b/pspy/tests/test_so_cov.py @@ -7,136 +7,6 @@ TEST_DATA_PREFIX = os.path.join(os.path.dirname(__file__), "data", "example_newcov") class SOCovmatTests(unittest.TestCase): - def setUp(self, verbose=False): - - # read in downgraded windows and ivar - self.lmax = 1349 - self.window = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "window_adjusted.fits")) - var0 = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "ivar1_adjusted.fits")) - var1 = so_map.read_map(os.path.join(TEST_DATA_PREFIX, "ivar2_adjusted.fits")) - self.cl_dict = np.load(os.path.join(TEST_DATA_PREFIX, "cl_dict.npy")) - self.nl_dict = np.load(os.path.join(TEST_DATA_PREFIX, "nl_dict.npy")) - self.cov_ref_TTTT = np.load(os.path.join(TEST_DATA_PREFIX, "cov_TTTT.npy")) - self.cov_ref_EEEE = np.load(os.path.join(TEST_DATA_PREFIX, "cov_EEEE.npy")) - - var0.data = 1 / var0.data - var1.data = 1 / var1.data - var0.data[np.isinf(var0.data)] = 0. - var1.data[np.isinf(var1.data)] = 0. - self.var0, self.var1 = var0, var1 - - self.binning_file = os.path.join(TEST_DATA_PREFIX, "binning.dat") # ignored - # pspy_utils.create_binning_file(bin_size=40, n_bins=100, file_name=self.binning_file) - - - - def test_covmat_TTTT(self, verbose=False): - lmax, binning_file = self.lmax, self.binning_file - window, var0, var1 = self.window, self.var0, self.var1 - cl_dict, nl_dict = self.cl_dict, self.nl_dict - - # set up scenario - survey_id = ["Ta", "Tb", "Tc", "Td"] - survey_name = ["s1", "s4", "s1", "s4"] - win = {'Ta': window, 'Tb': window, 'Tc': window, 'Td': window} - var = { - 'Ta': so_cov.make_weighted_variance_map(var0, window), - 'Tb': so_cov.make_weighted_variance_map(var1, window), - 'Tc': so_cov.make_weighted_variance_map(var0, window), - 'Td': so_cov.make_weighted_variance_map(var1, window) - } - - couplings = so_cov.generate_aniso_couplings_TTTT( - survey_id, survey_name, win, var, lmax) - - id2spec = {'Ta': 'dr6&pa6_f150_s1', 'Tb': 'dr6&pa6_f150_s4', - 'Tc': 'dr6&pa6_f150_s1', 'Td': 'dr6&pa6_f150_s4'} - - white_noise = { - 'Ta': so_cov.measure_white_noise_level(var0.data, window.data), - 'Tb': so_cov.measure_white_noise_level(var1.data, window.data), - 'Tc': so_cov.measure_white_noise_level(var0.data, window.data), - 'Td': so_cov.measure_white_noise_level(var1.data, window.data) - } - - # compute Rl, could put this in a function maybe - Clth_dict = {} - Rl_dict = {} - for name1, id1 in zip(survey_name, survey_id): - Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:(lmax+1)] - / white_noise[id1]) - for name2, id2 in zip(survey_name, survey_id): - spec = id1[0] + id2[0] - key12 = spec + id2spec[id1] + id2spec[id2] - if name1 == name2: - Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] - else: - Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] - - cov_e = so_cov.coupled_cov_aniso_same_pol(survey_id, Clth_dict, Rl_dict, couplings) - - num_diag = 30 # check first 30 off-diagonals - diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref_TTTT / cov_e, k) - 1)) - for k in range(0,num_diag) - ] - np.testing.assert_almost_equal(np.zeros(num_diag), - diag_max_errors, 7, - err_msg="aniso covmats don't match") - - def test_covmat_EEEE(self, verbose=False): - lmax, binning_file = self.lmax, self.binning_file - window, var0, var1 = self.window, self.var0, self.var1 - cl_dict, nl_dict = self.cl_dict, self.nl_dict - - # set up scenario - survey_id = ["Ea", "Eb", "Ec", "Ed"] - survey_name = ["s1", "s4", "s1", "s4"] - win = {'Ea': window, 'Eb': window, 'Ec': window, 'Ed': window} - var = { - 'Ea': so_cov.make_weighted_variance_map(var0, window), - 'Eb': so_cov.make_weighted_variance_map(var1, window), - 'Ec': so_cov.make_weighted_variance_map(var0, window), - 'Ed': so_cov.make_weighted_variance_map(var1, window) - } - - couplings = so_cov.generate_aniso_couplings_EEEE( - survey_id, survey_name, win, var, lmax) - - id2spec = {'Ea': 'dr6&pa6_f150_s1', 'Eb': 'dr6&pa6_f150_s4', - 'Ec': 'dr6&pa6_f150_s1', 'Ed': 'dr6&pa6_f150_s4'} - - white_noise = { - 'Ea': so_cov.measure_white_noise_level(var0.data, window.data), - 'Eb': so_cov.measure_white_noise_level(var1.data, window.data), - 'Ec': so_cov.measure_white_noise_level(var0.data, window.data), - 'Ed': so_cov.measure_white_noise_level(var1.data, window.data) - } - - # compute Rl, could put this in a function maybe - Clth_dict = {} - Rl_dict = {} - for name1, id1 in zip(survey_name, survey_id): - Rl_dict[id1] = np.sqrt(nl_dict[id1[0] + id1[0] + id2spec[id1] + id2spec[id1]][:(lmax+1)] - / white_noise[id1]) - for name2, id2 in zip(survey_name, survey_id): - spec = id1[0] + id2[0] - key12 = spec + id2spec[id1] + id2spec[id2] - if name1 == name2: - Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] - else: - Clth_dict[id1 + id2] = cl_dict[key12][:(lmax+1)] - - cov_e = so_cov.coupled_cov_aniso_same_pol(survey_id, Clth_dict, Rl_dict, couplings) - - num_diag = 30 # check first 30 off-diagonals - diag_max_errors = [ - np.max(np.abs(np.diag(self.cov_ref_EEEE[2:,2:] / cov_e[2:,2:], k) - 1)) - for k in range(0,num_diag) - ] - np.testing.assert_almost_equal(np.zeros(num_diag), - diag_max_errors, 7, - err_msg="aniso covmats don't match") def test_custom_couplings(self, verbose=False): ra0, ra1, dec0, dec1 = -25, 25, -15, 15 From 0a533741bc2d162bc963349dce12b257e3091c0d Mon Sep 17 00:00:00 2001 From: Zachary Atkins Date: Sat, 20 Apr 2024 17:29:24 -0400 Subject: [PATCH 23/23] calculate the lmax'th row and column in so_mcm.coupling_block --- pspy/so_mcm.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/pspy/so_mcm.py b/pspy/so_mcm.py index 590611b..4e9f494 100644 --- a/pspy/so_mcm.py +++ b/pspy/so_mcm.py @@ -525,11 +525,11 @@ def coupling_block( l = np.arange(len(wcl)) wcl *= (2 * l + 1) - mcm = np.zeros((lmax, lmax)) + mcm = np.zeros((lmax + 1, lmax + 1)) - if l_toep is None: l_toep = lmax - if l_band is None: l_band = lmax - if l_exact is None: l_exact = lmax + if l_toep is None: l_toep = lmax + 1 + if l_band is None: l_band = lmax + 1 + if l_exact is None: l_exact = lmax + 1 if mode == "00": mcm_fortran.calc_coupling_block_spin0(wcl, l_exact, l_band, l_toep, mcm.T) @@ -547,10 +547,10 @@ def coupling_block( mcm_fortran.fill_upper(mcm.T) - mcm_full = np.zeros((lmax+1, lmax+1)) - mcm_full[2:, 2:] = mcm[:-1,:-1] - mcm_full[0,0] = 1.0 - mcm_full[1,1] = 1.0 + mcm_full = np.zeros((lmax + 1, lmax + 1)) + mcm_full[2:, 2:] = mcm[:-2, :-2] + mcm_full[0, 0] = 1.0 + mcm_full[1, 1] = 1.0 if legacy_ell_range: return mcm_full[2:-1, 2:-1]