From 86f96166fbb028778c3bb7ce4724ea805b884557 Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 12:50:41 +0100 Subject: [PATCH 01/10] garrido black magic --- .github/workflows/wheels.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index e22cff8..d1a34df 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -17,7 +17,7 @@ jobs: - name: Build sdist and wheel run: | - python -m pip install wheel + python -m pip install wheel setuptools python setup.py sdist bdist_wheel - uses: actions/upload-artifact@v3 From 18f74a3a5d49dc36ef95db175f1f0457636b4418 Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 12:53:38 +0100 Subject: [PATCH 02/10] now assume T and P beam can differ, also some cosmetix --- pspipe_utils/best_fits.py | 84 +++++++++++++++++++++------------- pspipe_utils/misc.py | 48 +++++++++++++++++++ tutorials/tuto_covariances.py | 8 ++-- tutorials/tuto_kspace.py | 81 ++++++++++++++++---------------- tutorials/tuto_prepare_data.py | 15 ++++-- tutorials/tuto_simulation.py | 17 +++---- tutorials/tuto_split_nulls.py | 5 +- 7 files changed, 167 insertions(+), 91 deletions(-) diff --git a/pspipe_utils/best_fits.py b/pspipe_utils/best_fits.py index 6803124..d71063f 100644 --- a/pspipe_utils/best_fits.py +++ b/pspipe_utils/best_fits.py @@ -4,6 +4,7 @@ import numpy as np from mflike import theoryforge as th_mflike from pspy import pspy_utils, so_spectra +from pspipe_utils import misc def cmb_dict_from_file(f_name_cmb, lmax, spectra, lmin=2): """ @@ -31,7 +32,7 @@ def cmb_dict_from_file(f_name_cmb, lmax, spectra, lmin=2): return l_cmb, cmb_dict -def fg_dict_from_files(f_name_fg, array_list, lmax, spectra, lmin=2, f_name_cmb=None): +def fg_dict_from_files(f_name_fg, map_set_list, lmax, spectra, lmin=2, f_name_cmb=None): """ create a fg power spectrum dict from files @@ -39,10 +40,8 @@ def fg_dict_from_files(f_name_fg, array_list, lmax, spectra, lmin=2, f_name_cmb= __________ f_name_fg: string a template for the name of the fg power spectra files - array_list: list - list of all arrays - spectra: list - the list of spectra ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"] + map_set_list: list + list of all map set we will consider, format is {survey}_{array} lmax: integer the maximum multipole to consider (not inclusive) spectra: list @@ -57,20 +56,20 @@ def fg_dict_from_files(f_name_fg, array_list, lmax, spectra, lmin=2, f_name_cmb= l_cmb, cmb_dict = cmb_dict_from_file(f_name_cmb, lmax, spectra, lmin) fg_dict = {} - for i, ar1 in enumerate(array_list): - for j, ar2 in enumerate(array_list): - if i > j: array_tuple = (ar2, ar1) - else: array_tuple = (ar1, ar2) + for i, ms_1 in enumerate(map_set_list): + for j, ms_2 in enumerate(map_set_list): + if i > j: ms_tuple = (ms_2, ms_1) + else: ms_tuple = (ms_1, ms_2) - l_fg, fg = so_spectra.read_ps(f_name_fg.format(*array_tuple), spectra=spectra) + l_fg, fg = so_spectra.read_ps(f_name_fg.format(*ms_tuple), spectra=spectra) id_fg = np.where((l_fg >= lmin) & (l_fg < lmax)) - fg_dict[ar1, ar2] = {} + fg_dict[ms_1, ms_2] = {} for spec in spectra: if i > j: spec = spec[::-1] - fg_dict[ar1, ar2][spec] = fg[spec][id_fg] + fg_dict[ms_1, ms_2][spec] = fg[spec][id_fg] if f_name_cmb is not None: - fg_dict[ar1, ar2][spec] += cmb_dict[spec] + fg_dict[ms_1, ms_2][spec] += cmb_dict[spec] l_fg = l_fg[id_fg] @@ -120,14 +119,16 @@ def noise_dict_from_files(f_name_noise, sv_list, arrays, lmax, spectra, n_splits return l_noise, nl_dict -def beam_dict_from_files(f_name_beam, sv_list, arrays, lmax, lmin=2): +def beam_dict_from_files(f_name_beam_T, f_name_beam_pol, sv_list, arrays, lmax, lmin=2): """ create a beam dict from files Parameters __________ - f_name_beam: string - a template for the name of the beam files + f_name_beam_T: string + a template for the name of the temperature beam files + f_name_beam_pol: string + a template for the name of the polarisation beam files sv_list: list list of the surveys to consider arrays: dict @@ -141,16 +142,22 @@ def beam_dict_from_files(f_name_beam, sv_list, arrays, lmax, lmin=2): bl_dict = {} for sv in sv_list: for ar in arrays[sv]: - l_beam, bl = pspy_utils.read_beam_file(f_name_beam.format(sv, ar)) + + l_beam, bl = misc.read_beams(f_name_beam_T.format(sv, ar), + f_name_beam_pol.format(sv, ar)) + id_beam = np.where((l_beam >= lmin) & (l_beam < lmax)) - bl_dict[sv, ar] = bl[id_beam] + + bl_dict[sv, ar] = {} + for field in ["T", "E", "B"]: + bl_dict[sv, ar][field] = bl[field][id_beam] l_beam = l_beam[id_beam] return l_beam, bl_dict -def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, nl_dict=None, bl_dict=None): +def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", nl_dict=None, bl_dict=None): """ This function prepare all best fit corresponding to the spec_name_list. the ps_all_th and nl_all_th are in particular useful for the analytical covariance computation @@ -166,6 +173,13 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, nl_dict=N the cmb ps (format is [spec] fg_dict: dict of dict the fg ps (format is [sv1_ar1,sv2_ar2][spec]) + spectra: list + the list of spectra ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"] + char: string + serve to split the map_set into survey and array + a bit annoying to have to keep this, in principle we could work only + with map_set, but noise for different survey should be set to zero since + this entry doesn't exist nl_dict: dict of dict the noise ps (format is [sv, ar1, ar2][spec]) bl_dict: dict @@ -175,33 +189,37 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, nl_dict=N ps_all_th, nl_all_th = {}, {} for spec_name in spec_name_list: - na, nb = spec_name.split("x") - if len(na.split("&")) == 2: - sv_a, ar_a = na.split("&") - sv_b, ar_b = nb.split("&") + ms_a, ms_b = spec_name.split("x") + if len(msa.split(char)) == 2: + sv_a, ar_a = ms_a.split(char) + sv_b, ar_b = ms_b.split(char) noise_key_a = ar_a noise_key_b = ar_b - elif len(na.split("&")) == 3: - sv_a, ar_a, split_a = na.split("&") - sv_b, ar_b, split_b = nb.split("&") + elif len(msa.split("&")) == 3: + sv_a, ar_a, split_a = ms_a.split(char) + sv_b, ar_b, split_b = ms_b.split(char) noise_key_a = f"{ar_a}_{split_a}" noise_key_b = f"{ar_b}_{split_b}" for spec in spectra: - ps_all_th[na, nb, spec] = cmb_dict[spec] + fg_dict[f"{sv_a}_{ar_a}", f"{sv_b}_{ar_b}"][spec] + + ps_all_th[ms_a, ms_b, spec] = cmb_dict[spec] + fg_dict[f"{sv_a}_{ar_a}", f"{sv_b}_{ar_b}"][spec] + ps_all_th[ms_b, ms_a, spec] = ps_all_th[ms_a, ms_b, spec].copy() if bl_dict is not None: - ps_all_th[na, nb, spec] *= bl_dict[sv_a, ar_a] * bl_dict[sv_b, ar_b] - - ps_all_th[nb, na, spec] = ps_all_th[na, nb, spec] + X, Y = spec + ps_all_th[ms_a, ms_b, spec] *= bl_dict[sv_a, ar_a][X] * bl_dict[sv_b, ar_b][Y] + if ms_a != ms_b: + # the if avoid a repetition in the case ms_a == ms_b + ps_all_th[ms_b, ms_a, spec] *= bl_dict[sv_b, ar_b][X] * bl_dict[sv_a, ar_a][Y] if nl_dict is not None: if sv_a == sv_b: - nl_all_th[na, nb, spec] = nl_dict[sv_a, noise_key_a, noise_key_b][spec] + nl_all_th[ms_a, ms_b, spec] = nl_dict[sv_a, noise_key_a, noise_key_b][spec] else: - nl_all_th[na, nb, spec] = 0 + nl_all_th[ms_a, ms_b, spec] = 0 - nl_all_th[nb, na, spec] = nl_all_th[na, nb, spec] + nl_all_th[ms_b, ms_a, spec] = nl_all_th[ms_a, ms_b, spec] if nl_dict is not None: return l_th, ps_all_th, nl_all_th diff --git a/pspipe_utils/misc.py b/pspipe_utils/misc.py index 968dd61..4fe5596 100644 --- a/pspipe_utils/misc.py +++ b/pspipe_utils/misc.py @@ -2,11 +2,59 @@ Some general utility functions that do not belong to other files. """ +from pspy import pspy_utils +from pixell import curvedsky def str_replace(my_str, old, new): + """ + just like replace but check that the replacement actually happened + + Parameters + __________ + my_str: string + the string in which the replacment will happen + old: string + old part of the string to be replaced + new: string + what will replace old + """ + my_new_str = my_str my_new_str = my_str.replace(old, new) if my_new_str == my_str: error = f" the name '{my_str}' does not contain '{old}' so I can't replace '{old}' by '{new}'" raise NameError(error) return my_new_str + +def read_beams(f_name_beam_T, f_name_beam_pol): + """ + read T and pol beams and return a beam dictionnary with entry T, E, B + + Parameters + __________ + f_name_beam_T: string + the filename of the temperature beam file + f_name_beam_pol: string + the filename of the polarisation beam file + """ + + bl = {} + l, bl["T"] = pspy_utils.read_beam_file(f_name_beam_T) + l, bl["E"] = pspy_utils.read_beam_file(f_name_beam_pol) + bl["B"] = bl["E"] + return l, bl + +def apply_beams(alms, bl): + """ + apply T and pol beams to alms + + Parameters + __________ + alms: 2d array + array of alms, alms[0]=alm_T, alms[1]=alm_E, alms[2]=alm_B + bl: dict + dictionnary containing T and E,B beams + """ + for i, f in enumerate(["T", "E", "B"]): + alms[i] = curvedsky.almxfl(alms[i], bl[f]) + return alms diff --git a/tutorials/tuto_covariances.py b/tutorials/tuto_covariances.py index 90e7ae7..26a966e 100644 --- a/tutorials/tuto_covariances.py +++ b/tutorials/tuto_covariances.py @@ -2,10 +2,11 @@ This script generate covariance matrices """ import numpy as np -from pspipe_utils import best_fits, pspipe_list +from pspipe_utils import pspipe_list, best_fits from pspy import pspy_utils, so_cov, so_mcm, so_dict import time, sys + d = so_dict.so_dict() d.read_from_file(sys.argv[1]) @@ -37,12 +38,13 @@ f_name_cmb = tuto_data_dir + "/cmb.dat" f_name_noise = tuto_data_dir + "/mean_{}x{}_{}_noise.dat" f_name_fg = tuto_data_dir + "/fg_{}x{}.dat" -f_name_beam = tuto_data_dir + "/beam_{}_{}.dat" +f_name_beam_T = tuto_data_dir + "/beam_{}_{}_T.dat" +f_name_beam_pol = tuto_data_dir + "/beam_{}_{}_pol.dat" l_cmb, cmb_dict = best_fits.cmb_dict_from_file(f_name_cmb, lmax, spectra) l_fg, fg_dict = best_fits.fg_dict_from_files(f_name_fg, arrays_list, lmax, spectra) l_noise, nl_dict = best_fits.noise_dict_from_files(f_name_noise, surveys, arrays, lmax, spectra, n_splits=n_splits) -l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam, surveys, arrays, lmax) +l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam_T, f_name_beam_pol, surveys, arrays, lmax) l_cmb, ps_all_th, nl_all_th = best_fits.get_all_best_fit(spec_name_list, l_cmb, diff --git a/tutorials/tuto_kspace.py b/tutorials/tuto_kspace.py index 44ce2c8..3413bf4 100644 --- a/tutorials/tuto_kspace.py +++ b/tutorials/tuto_kspace.py @@ -6,7 +6,7 @@ import pylab as plt import numpy as np import pspipe_utils -from pspipe_utils import simulation, best_fits, kspace +from pspipe_utils import simulation, best_fits, kspace, misc from pixell import curvedsky, enmap import os, time @@ -20,7 +20,7 @@ ncomp = 3 niter = 0 -n_sims = 40 +n_sims = 4 spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"] binning_file = f"{data_path}/binning_files/BIN_ACTPOL_50_4_SC_large_bin_at_low_ell" vk_mask = [-90, 90] @@ -52,23 +52,26 @@ survey = "dr6" arrays = ["pa4_f150", "pa4_f220", "pa5_f090", "pa5_f150", "pa6_f090", "pa6_f150"] nu_effs = {} -nu_effs["dr6", "pa4_f150"] = 150 -nu_effs["dr6", "pa4_f220"] = 220 -nu_effs["dr6", "pa5_f090"] = 90 -nu_effs["dr6", "pa5_f150"] = 150 -nu_effs["dr6", "pa6_f090"] = 90 -nu_effs["dr6", "pa6_f150"] = 150 +nu_effs["dr6_pa4_f150"] = 150 +nu_effs["dr6_pa4_f220"] = 220 +nu_effs["dr6_pa5_f090"] = 90 +nu_effs["dr6_pa5_f150"] = 150 +nu_effs["dr6_pa6_f090"] = 90 +nu_effs["dr6_pa6_f150"] = 150 ########################################################################################################## # let's generate the filter, the spactial window function and compute the associated mode coupling matrices ########################################################################################################### -window_tuple, mcm_inv, Bbl, filter_std, binary, bl = {}, {}, {}, {}, {}, {} +window_tuple, mcm_inv, Bbl, filter_std, binary, bl, passbands = {}, {}, {}, {}, {}, {}, {} spec_name_list = [] +map_set_list = [] for ar in arrays: - binary[ar] = so_map.read_map(f"{data_path}/binaries/binary_dr6_{ar}_downgraded.fits") + map_set = f"{survey}_{ar}" + + binary[ar] = so_map.read_map(f"{data_path}/binaries/binary_{map_set}_downgraded.fits") lmax = int(binary[ar].get_lmax_limit()) template = binary[ar].copy() @@ -86,30 +89,34 @@ mask = so_window.create_apodization(mask, apo_type="C1", apo_radius_degree=0.3) window.data *= mask.data - window.plot(file_name=f"{test_dir}/window_dr6_{ar}") + window.plot(file_name=f"{test_dir}/window_{map_set}") - binary[ar].plot(file_name=f"{test_dir}/binary_{ar}") + binary[ar].plot(file_name=f"{test_dir}/binary_{map_set}") window_tuple[ar] = (window, window) - l, bl[ar] = pspy_utils.read_beam_file(f"{data_path}/beams/coadd_{ar}_night_beam_tform_jitter_cmb.txt") + l, bl[ar] = misc.read_beams(f"{data_path}/beams/coadd_{ar}_night_beam_tform_jitter_cmb.txt", + f"{data_path}/beams/coadd_{ar}_night_beam_tform_jitter_cmb.txt") mcm_inv[ar], Bbl[ar] = so_mcm.mcm_and_bbl_spin0and2(window_tuple[ar], - bl1 = (bl[ar], bl[ar]), - binning_file = binning_file, - lmax=lmax, - type=type, - niter=niter, - binned_mcm=binned_mcm) + bl1 = (bl[ar]["T"], bl[ar]["E"]), + binning_file = binning_file, + lmax=lmax, + type=type, + niter=niter, + binned_mcm=binned_mcm) filter_std[ar] = so_map_preprocessing.build_std_filter(template.data.shape, template.data.wcs, vk_mask, hk_mask, dtype=np.float64) - - spec_name = f"{survey}&{ar}x{survey}&{ar}" - spec_name_list += [spec_name] + + passbands[f"{map_set}"] = np.array([nu_effs[f"{map_set}"]]), np.array([1]) + spec_name_list += [f"{survey}&{ar}x{survey}&{ar}"] + map_set_list += [map_set] + + ############################################################################ # let's prepare the theory and foreground matrix used to generate simulation ############################################################################ @@ -120,22 +127,17 @@ so_spectra.write_ps(f_name_cmb, l_th, ps_dict, type, spectra=spectra) ps_mat = simulation.cmb_matrix_from_file(f_name_cmb, lmax, spectra) -freq_list = [] -for ar in arrays: - freq_list += [nu_effs[survey, ar]] -freq_list = list(dict.fromkeys(freq_list)) - -fg_dict = best_fits.get_foreground_dict(l_th, freq_list, fg_components, fg_params, fg_norm=None) -fg= {} -for freq1 in freq_list: - for freq2 in freq_list: - fg[freq1, freq2] = {} + +fg_dict = best_fits.get_foreground_dict(l_th, passbands, fg_components, fg_params, fg_norm=None) +for ms1 in map_set_list: + for ms2 in map_set_list: + fg = {} for spec in spectra: - fg[freq1,freq2][spec] = fg_dict[spec.lower(), "all", freq1, freq2] - so_spectra.write_ps(f"{test_dir}/fg_{freq1}x{freq2}.dat", l_th, fg[freq1,freq2], type, spectra=spectra) + fg[spec] = fg_dict[spec.lower(), "all", ms1, ms2] + so_spectra.write_ps(f"{test_dir}/fg_{ms1}x{ms2}.dat", l_th, fg, type, spectra=spectra) f_name_fg = test_dir + "/fg_{}x{}.dat" -l, fg_mat = simulation.foreground_matrix_from_files(f_name_fg, freq_list, lmax, spectra) +l, fg_mat = simulation.foreground_matrix_from_files(f_name_fg, map_set_list, lmax, spectra) ################################# # let's generate the simulations @@ -151,16 +153,16 @@ t = time.time() alms_cmb = curvedsky.rand_alm(ps_mat, lmax=lmax, dtype="complex64") - fglms = simulation.generate_fg_alms(fg_mat, freq_list, lmax) + fglms = simulation.generate_fg_alms(fg_mat, map_set_list, lmax) for scenario in scenarios: for ar in arrays: alms_beamed = alms_cmb.copy() - alms_beamed += fglms[nu_effs[survey, ar]] + alms_beamed += fglms[f"{survey}_{ar}"] - alms_beamed = curvedsky.almxfl(alms_beamed, bl[ar]) + alms_beamed = misc.apply_beams(alms_beamed, bl[ar]) if scenario == "noE": alms_beamed[1] *= 0 if scenario == "noB": alms_beamed[2] *= 0 @@ -240,12 +242,11 @@ plt.close() l_cmb, cmb_dict = best_fits.cmb_dict_from_file(f_name_cmb, lmax + 2, spectra) -l_fg, fg_dict = best_fits.fg_dict_from_files(f_name_fg, freq_list, lmax + 2, spectra) +l_fg, fg_dict = best_fits.fg_dict_from_files(f_name_fg, map_set_list, lmax + 2, spectra) l_cmb, ps_all_th = best_fits.get_all_best_fit(spec_name_list, l_cmb, cmb_dict, fg_dict, - nu_effs, spectra) for ar in arrays: diff --git a/tutorials/tuto_prepare_data.py b/tutorials/tuto_prepare_data.py index 0b83874..89ed5d7 100644 --- a/tutorials/tuto_prepare_data.py +++ b/tutorials/tuto_prepare_data.py @@ -37,9 +37,14 @@ rms_uKarcmin_T[sv, f"{ar}x{ar}"] = d[f"rms_uKarcmin_T_{sv}_{ar}"] fwhm = d[f"beam_fwhm_{sv}_{ar}"] - l_beam, bl[sv, ar] = pspy_utils.beam_from_fwhm(fwhm, lmax_sim) - f_name = f"{tuto_data_dir}/beam_{sv}_{ar}.dat" - np.savetxt(f_name, np.transpose([l_beam, bl[sv, ar]])) + l_beam, bl[sv, ar, "T"] = pspy_utils.beam_from_fwhm(fwhm, lmax_sim) + l_beam, bl[sv, ar, "E"] = pspy_utils.beam_from_fwhm(fwhm * 1.3, lmax_sim) + + f_name = f"{tuto_data_dir}/beam_{sv}_{ar}_T.dat" + np.savetxt(f_name, np.transpose([l_beam, bl[sv, ar, "T"]])) + + f_name = f"{tuto_data_dir}/beam_{sv}_{ar}_pol.dat" + np.savetxt(f_name, np.transpose([l_beam, bl[sv, ar, "E"]])) for sv in surveys: @@ -124,7 +129,7 @@ # we need both the window for T and pol, here we assume they are the same window_tuple_a = (window[sv_a, ar_a], window[sv_a, ar_a]) - bl_a = (bl[sv_a, ar_a], bl[sv_a, ar_a]) + bl_a = (bl[sv_a, ar_a, "T"], bl[sv_a, ar_a, "E"]) for id_sv_b, sv_b in enumerate(surveys): for id_ar_b, ar_b in enumerate(arrays[sv_b]): @@ -138,7 +143,7 @@ # the mode coupling matrices for sv1 x sv2 are the same as the one for sv2 x sv1 # identically, within a survey, the mode coupling matrix for ar1 x ar2 = ar2 x ar1 window_tuple_b = (window[sv_b, ar_b], window[sv_b, ar_b]) - bl_b = (bl[sv_b, ar_b], bl[sv_b, ar_b]) + bl_b = (bl[sv_b, ar_b, "T"], bl[sv_b, ar_b, "E"]) mbb_inv, Bbl = so_mcm.mcm_and_bbl_spin0and2(window_tuple_a, win2 = window_tuple_b, diff --git a/tutorials/tuto_simulation.py b/tutorials/tuto_simulation.py index e5187ae..06e6d20 100644 --- a/tutorials/tuto_simulation.py +++ b/tutorials/tuto_simulation.py @@ -2,7 +2,7 @@ """ import numpy as np import pylab as plt -from pspipe_utils import simulation, best_fits, pspipe_list +from pspipe_utils import simulation, pspipe_list, misc, best_fits from pspy import pspy_utils, so_spectra, sph_tools, so_map, so_cov, so_window, so_mcm, so_dict from pixell import curvedsky import time, sys @@ -33,14 +33,14 @@ f_name_cmb = tuto_data_dir + "/cmb.dat" f_name_noise = tuto_data_dir + "/mean_{}x{}_{}_noise.dat" f_name_fg = tuto_data_dir + "/fg_{}x{}.dat" -f_name_beam = tuto_data_dir + "/beam_{}_{}.dat" +f_name_beam_T = tuto_data_dir + "/beam_{}_{}_T.dat" +f_name_beam_pol = tuto_data_dir + "/beam_{}_{}_pol.dat" arrays, bl, window = {}, {}, {} for sv in surveys: arrays[sv] = d[f"arrays_{sv}"] for ar in arrays[sv]: - l, bl[sv, ar] = np.loadtxt(f"{f_name_beam.format(sv, ar)}", unpack=True) - #nu_eff[sv, ar] = d[f"nu_eff_{sv}_{ar}"] + l, bl[sv, ar] = misc.read_beams(f_name_beam_T.format(sv, ar), f_name_beam_pol.format(sv, ar)) window[sv, ar] = so_map.read_map(f"{tuto_data_dir}/window_{sv}_{ar}.fits") ps_mat = simulation.cmb_matrix_from_file(f_name_cmb, lmax_sim, spectra) @@ -68,8 +68,7 @@ signal_alms = {} for ar in arrays[sv]: signal_alms[ar] = alms_cmb + fglms[f"{sv}_{ar}"] - for i in range(3): - signal_alms[ar][i] = curvedsky.almxfl(signal_alms[ar][i], bl[sv, ar]) + signal_alms[ar] = misc.apply_beams(signal_alms[ar], bl[sv, ar]) for k in range(n_splits[sv]): noise_alms = simulation.generate_noise_alms(noise_mat[sv], arrays[sv], lmax_sim) @@ -150,7 +149,7 @@ l_cmb, cmb_dict = best_fits.cmb_dict_from_file(f_name_cmb, lmax, spectra) l_fg, fg_dict = best_fits.fg_dict_from_files(f_name_fg, array_list, lmax, spectra) l_noise, nl_dict = best_fits.noise_dict_from_files(f_name_noise, surveys, arrays, lmax, spectra) -l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam, surveys, arrays, lmax) +l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam_T, f_name_beam_pol, surveys, arrays, lmax) l_cmb, ps_all_th, nl_all_th = best_fits.get_all_best_fit(spec_name_list, l_cmb, @@ -173,6 +172,8 @@ plt.figure(figsize=(18, 14)) for count, spec in enumerate(spectra): + X, Y = spec + mc_cov_select = so_cov.selectblock(mc_cov, spectra, n_bins, block=spec+spec) std = np.sqrt(mc_cov_select.diagonal()) plt.subplot(3, 3, count + 1) @@ -188,7 +189,7 @@ plt.title(spec_name) plt.legend(fontsize=14) else: - plt.plot(l_cmb, nl_all_th[f"{sv_a}&{ar_a}", f"{sv_b}&{ar_b}", spec] / (bl_dict[sv_a, ar_a] * bl_dict[sv_b, ar_b])) + plt.plot(l_cmb, nl_all_th[f"{sv_a}&{ar_a}", f"{sv_b}&{ar_b}", spec] / (bl_dict[sv_a, ar_a][X] * bl_dict[sv_b, ar_b][Y])) plt.xlabel("$\ell$", fontsize=14) plt.ylabel("$DN^{%s}_\ell$" % spec, fontsize=16) diff --git a/tutorials/tuto_split_nulls.py b/tutorials/tuto_split_nulls.py index f88710f..8ed8e6b 100644 --- a/tutorials/tuto_split_nulls.py +++ b/tutorials/tuto_split_nulls.py @@ -47,12 +47,13 @@ f_name_cmb = tuto_data_dir + "/cmb.dat" f_name_noise = tuto_data_dir + "/mean_{}x{}_{}_noise.dat" f_name_fg = tuto_data_dir + "/fg_{}x{}.dat" -f_name_beam = tuto_data_dir + "/beam_{}_{}.dat" +f_name_beam_T = tuto_data_dir + "/beam_{}_{}_T.dat" +f_name_beam_pol = tuto_data_dir + "/beam_{}_{}_pol.dat" l_cmb, cmb_dict = best_fits.cmb_dict_from_file(f_name_cmb, lmax, spectra) l_fg, fg_dict = best_fits.fg_dict_from_files(f_name_fg, arrays_list, lmax, spectra) l_noise, nl_dict = best_fits.noise_dict_from_files(f_name_noise, surveys, arrays, lmax, spectra, n_splits=n_splits) -l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam, surveys, arrays, lmax) +l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam_T, f_name_beam_pol, surveys, arrays, lmax) l_cmb, ps_all_th, nl_all_th = best_fits.get_all_best_fit(spec_name_list, l_cmb, From f3b38f3ccc0a14aafbd39cc22568d54075c2f526 Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 13:25:57 +0100 Subject: [PATCH 03/10] fix typo --- pspipe_utils/best_fits.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pspipe_utils/best_fits.py b/pspipe_utils/best_fits.py index d71063f..b7cef5c 100644 --- a/pspipe_utils/best_fits.py +++ b/pspipe_utils/best_fits.py @@ -190,12 +190,12 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", for spec_name in spec_name_list: ms_a, ms_b = spec_name.split("x") - if len(msa.split(char)) == 2: + if len(ms_a.split(char)) == 2: sv_a, ar_a = ms_a.split(char) sv_b, ar_b = ms_b.split(char) noise_key_a = ar_a noise_key_b = ar_b - elif len(msa.split("&")) == 3: + elif len(ms_a.split("&")) == 3: sv_a, ar_a, split_a = ms_a.split(char) sv_b, ar_b, split_b = ms_b.split(char) noise_key_a = f"{ar_a}_{split_a}" From ad854aad681bfed7059d8340c41ac2bfb0f7f9fd Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 19:21:47 +0100 Subject: [PATCH 04/10] useless following Xavier comment --- pspipe_utils/misc.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pspipe_utils/misc.py b/pspipe_utils/misc.py index 4fe5596..8f49570 100644 --- a/pspipe_utils/misc.py +++ b/pspipe_utils/misc.py @@ -19,7 +19,6 @@ def str_replace(my_str, old, new): what will replace old """ - my_new_str = my_str my_new_str = my_str.replace(old, new) if my_new_str == my_str: error = f" the name '{my_str}' does not contain '{old}' so I can't replace '{old}' by '{new}'" From 26665141cc03cb0565f512e70be3c48a2f8c03ad Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 19:24:37 +0100 Subject: [PATCH 05/10] fix ug --- pspipe_utils/best_fits.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pspipe_utils/best_fits.py b/pspipe_utils/best_fits.py index b7cef5c..e5f9dda 100644 --- a/pspipe_utils/best_fits.py +++ b/pspipe_utils/best_fits.py @@ -195,7 +195,7 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", sv_b, ar_b = ms_b.split(char) noise_key_a = ar_a noise_key_b = ar_b - elif len(ms_a.split("&")) == 3: + elif len(ms_a.split(char)) == 3: sv_a, ar_a, split_a = ms_a.split(char) sv_b, ar_b, split_b = ms_b.split(char) noise_key_a = f"{ar_a}_{split_a}" @@ -217,7 +217,7 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", if sv_a == sv_b: nl_all_th[ms_a, ms_b, spec] = nl_dict[sv_a, noise_key_a, noise_key_b][spec] else: - nl_all_th[ms_a, ms_b, spec] = 0 + nl_all_th[ms_a, ms_b, spec] = 0.0 nl_all_th[ms_b, ms_a, spec] = nl_all_th[ms_a, ms_b, spec] From ebe11f9a031fcc818be4c5bd2a86f2a00445897c Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 19:34:49 +0100 Subject: [PATCH 06/10] better name --- pspipe_utils/best_fits.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/pspipe_utils/best_fits.py b/pspipe_utils/best_fits.py index e5f9dda..82d7ea8 100644 --- a/pspipe_utils/best_fits.py +++ b/pspipe_utils/best_fits.py @@ -157,7 +157,7 @@ def beam_dict_from_files(f_name_beam_T, f_name_beam_pol, sv_list, arrays, lmax, return l_beam, bl_dict -def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", nl_dict=None, bl_dict=None): +def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, delimiter="&", nl_dict=None, bl_dict=None): """ This function prepare all best fit corresponding to the spec_name_list. the ps_all_th and nl_all_th are in particular useful for the analytical covariance computation @@ -175,7 +175,7 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", the fg ps (format is [sv1_ar1,sv2_ar2][spec]) spectra: list the list of spectra ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"] - char: string + delimiter: string serve to split the map_set into survey and array a bit annoying to have to keep this, in principle we could work only with map_set, but noise for different survey should be set to zero since @@ -190,14 +190,14 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, char="&", for spec_name in spec_name_list: ms_a, ms_b = spec_name.split("x") - if len(ms_a.split(char)) == 2: - sv_a, ar_a = ms_a.split(char) - sv_b, ar_b = ms_b.split(char) + if len(ms_a.split(delimiter)) == 2: + sv_a, ar_a = ms_a.split(delimiter) + sv_b, ar_b = ms_b.split(delimiter) noise_key_a = ar_a noise_key_b = ar_b - elif len(ms_a.split(char)) == 3: - sv_a, ar_a, split_a = ms_a.split(char) - sv_b, ar_b, split_b = ms_b.split(char) + elif len(ms_a.split(delimiter)) == 3: + sv_a, ar_a, split_a = ms_a.split(delimiter) + sv_b, ar_b, split_b = ms_b.split(delimiter) noise_key_a = f"{ar_a}_{split_a}" noise_key_b = f"{ar_b}_{split_b}" From 3b5eb21b8b9eb35b4c09fc589383e858124c629e Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 19:36:07 +0100 Subject: [PATCH 07/10] better name --- pspipe_utils/pspipe_list.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pspipe_utils/pspipe_list.py b/pspipe_utils/pspipe_list.py index d123b41..002ac04 100644 --- a/pspipe_utils/pspipe_list.py +++ b/pspipe_utils/pspipe_list.py @@ -86,14 +86,14 @@ def get_covariances_list(dict): return ncovs, na_list, nb_list, nc_list, nd_list -def get_spec_name_list(dict, char="&", kind=None, freq_pair=None, remove_same_ar_and_sv=False, return_nu_tag=False): +def get_spec_name_list(dict, delimiter="&", kind=None, freq_pair=None, remove_same_ar_and_sv=False, return_nu_tag=False): """This function creates a list with the name of all spectra we consider Parameters ---------- dict : dict the global dictionnary file used in pspipe - char: str + delimiter: str a character that separate the suvey and array name kind : str if "noise" or "auto" won't return @@ -136,7 +136,7 @@ def get_spec_name_list(dict, char="&", kind=None, freq_pair=None, remove_same_ar if remove_same_ar_and_sv == True: if (sv1 == sv2) & (ar1 == ar2): continue - spec_name_list += [f"{sv1}{char}{ar1}x{sv2}{char}{ar2}"] + spec_name_list += [f"{sv1}{delimiter}{ar1}x{sv2}{delimiter}{ar2}"] nu_tag_list += [(nu_tag1, nu_tag2)] if return_nu_tag == False: From c1bb890ca8580e569cd72a80b4268e475908e95c Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Sun, 29 Oct 2023 19:38:13 +0100 Subject: [PATCH 08/10] bump version --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 854df6f..b10bb5f 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.1.3 +current_version = 0.1.4 commit = True tag = True From 752ce416a179f8f9504f860761dee74f85e4ef7e Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Mon, 30 Oct 2023 15:57:24 +0100 Subject: [PATCH 09/10] add option to choose lmax --- pspipe_utils/misc.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pspipe_utils/misc.py b/pspipe_utils/misc.py index 8f49570..2466455 100644 --- a/pspipe_utils/misc.py +++ b/pspipe_utils/misc.py @@ -25,7 +25,7 @@ def str_replace(my_str, old, new): raise NameError(error) return my_new_str -def read_beams(f_name_beam_T, f_name_beam_pol): +def read_beams(f_name_beam_T, f_name_beam_pol, lmax=None): """ read T and pol beams and return a beam dictionnary with entry T, E, B @@ -38,8 +38,8 @@ def read_beams(f_name_beam_T, f_name_beam_pol): """ bl = {} - l, bl["T"] = pspy_utils.read_beam_file(f_name_beam_T) - l, bl["E"] = pspy_utils.read_beam_file(f_name_beam_pol) + l, bl["T"] = pspy_utils.read_beam_file(f_name_beam_T, lmax=lmax) + l, bl["E"] = pspy_utils.read_beam_file(f_name_beam_pol, lmax=lmax) bl["B"] = bl["E"] return l, bl From ca54b9d072ca658bf3e78046b4fdbc272f651594 Mon Sep 17 00:00:00 2001 From: Thibaut Louis Date: Mon, 30 Oct 2023 15:58:11 +0100 Subject: [PATCH 10/10] add docstring --- pspipe_utils/misc.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pspipe_utils/misc.py b/pspipe_utils/misc.py index 2466455..8549733 100644 --- a/pspipe_utils/misc.py +++ b/pspipe_utils/misc.py @@ -35,6 +35,8 @@ def read_beams(f_name_beam_T, f_name_beam_pol, lmax=None): the filename of the temperature beam file f_name_beam_pol: string the filename of the polarisation beam file + lmax : integer + the maximum multipole to consider (note that usually beam file start at l=0) """ bl = {}