Skip to content

Commit ad69efb

Browse files
Louis ThibautLouis Thibaut
Louis Thibaut
authored and
Louis Thibaut
committed
first version of SPT pipeline
1 parent 0ae9b19 commit ad69efb

10 files changed

Lines changed: 907 additions & 0 deletions

‎project/SPT_dev/SPT.rst‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
**************************
2+
Redoing SPT analysis
3+
**************************
4+
5+
In this pipeline, I am putting a bunch of scripts that can be used to redo the SPT analysis at NERSC.
6+
This is currently in dev, note that you need pspy version 1.8.4 to run it.
7+
8+
9+
Running the main pipeline
10+
-------------------------------------------------------
11+
12+
.. code:: shell
13+
14+
salloc --nodes 1 --qos interactive --time 02:00:00 --constraint cpu
15+
OMP_NUM_THREADS=32 srun -n 6 -c 32 --cpu-bind=cores python spt_get_mcm_and_bbl.py global_spt.dict
16+
OMP_NUM_THREADS=64 srun -n 3 -c 64 --cpu-bind=cores python spt_get_alms.py global_spt.dict
17+
OMP_NUM_THREADS=64 srun -n 3 -c 64 --cpu-bind=cores python spt_get_spectra_from_alms.py global_spt.dict
18+
19+
You can then compare with the spt released spectra by running spt_plot_spectra.py, note that you need to install candl and spt_candl_data first.
20+
The pipeline take roughly 30 minutes, all time is spent doing map2alms since spt has 30 splits per frequency.
21+
An alternative dict file using half-missing is provided that allow for much faster investigation

‎project/SPT_dev/global_spt.dict‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
surveys = ["spt"]
2+
arrays_spt = ["095", "150", "220"]
3+
#release_dir = "/global/cfs/cdirs/sobs/users/tlouis/spt_data/"
4+
release_dir = "/Users/louisthibaut/Desktop/projects/spt_fun/"
5+
6+
alm_conv_spt = "IAU"
7+
8+
9+
binning_file = release_dir + "binning_spt.dat"
10+
apply_kspace_filter = False
11+
12+
deconvolve_pixwin = True
13+
pixwin_spt = {"pix": 'HEALPIX', "nside": 8192} #should be replaced with the more accurate pixel window that is provided
14+
15+
16+
apply_alm_mask = True
17+
lmax_mask = 4000
18+
alm_mask_spt_095 = release_dir + "ancillary_products/specific_to_c25/alm_mask_095ghz.fits"
19+
alm_mask_spt_150 = release_dir + "ancillary_products/specific_to_c25/alm_mask_150ghz.fits"
20+
alm_mask_spt_220 = release_dir + "ancillary_products/specific_to_c25/alm_mask_220ghz.fits"
21+
22+
niter = 0
23+
remove_mean = True
24+
binned_mcm = True
25+
lmax = 3500
26+
type = 'Dl'
27+
write_splits_spectra = True
28+
cov_T_E_only = False
29+
use_toeplitz_mcm = False
30+
use_toeplitz_cov = True
31+
32+
# maps
33+
34+
src_free_maps_spt = False
35+
36+
37+
freq_info_spt_095 = {"freq_tag": 95, "passband": "dummy.dat"}
38+
freq_info_spt_150 = {"freq_tag": 150, "passband": "dummy.dat"}
39+
freq_info_spt_220 = {"freq_tag": 220, "passband": "dummy.dat"}
40+
41+
42+
n_splits_spt = 30
43+
maps_spt_095 = [release_dir + "real_data_maps/one_thirtieth/one_thirtieth_bundle%02d_095ghz_inpainted.fits" % (i) for i in range(n_splits_spt)]
44+
maps_spt_150 = [release_dir + "real_data_maps/one_thirtieth/one_thirtieth_bundle%02d_150ghz_inpainted.fits" % (i) for i in range(n_splits_spt)]
45+
maps_spt_220 = [release_dir + "real_data_maps/one_thirtieth/one_thirtieth_bundle%02d_220ghz_inpainted.fits" % (i) for i in range(n_splits_spt)]
46+
47+
48+
# these number are just the one associated with their beam at l=0
49+
# in act we use the convention that bl(l=0)=1, they don't do this in spt
50+
cal_spt_095 = 1 / 1.09647802
51+
cal_spt_150 = 1 / 1.06827578
52+
cal_spt_220 = 1 / 1.11338355
53+
54+
pol_eff_spt_095 = 1.0
55+
pol_eff_spt_150 = 1.0
56+
pol_eff_spt_220 = 1.0
57+
58+
beam_file_dir = release_dir + "ancillary_products/specific_to_c25/"
59+
60+
beam_T_spt_095 = beam_file_dir + "beam_c25_bl_095ghz.txt"
61+
beam_T_spt_150 = beam_file_dir + "beam_c25_bl_150ghz.txt"
62+
beam_T_spt_220 = beam_file_dir + "beam_c25_bl_220ghz.txt"
63+
64+
beam_pol_spt_095 = beam_file_dir + "beam_c25_bl_095ghz.txt"
65+
beam_pol_spt_150 = beam_file_dir + "beam_c25_bl_150ghz.txt"
66+
beam_pol_spt_220 = beam_file_dir + "beam_c25_bl_220ghz.txt"
67+
68+
window_dir = release_dir + "ancillary_products/generally_applicable/"
69+
70+
window_T_spt_095 = window_dir + "pixel_mask_apodized_borders_only.fits"
71+
window_pol_spt_095 = window_dir + "pixel_mask_apodized_borders_only.fits"
72+
73+
window_T_spt_150 = window_dir + "pixel_mask_apodized_borders_only.fits"
74+
window_pol_spt_150 = window_dir + "pixel_mask_apodized_borders_only.fits"
75+
76+
window_T_spt_220 = window_dir + "pixel_mask_apodized_borders_only.fits"
77+
window_pol_spt_220 = window_dir + "pixel_mask_apodized_borders_only.fits"
78+
79+
80+
cosmo_params = {"cosmomc_theta":0.0104059, "logA": 3.053, "ombh2": 0.02258, "omch2": 0.1242, "ns": 0.9666, "Alens": 1.0, "tau": 0.0567}
81+
82+

‎project/SPT_dev/global_spt_half.dict‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
surveys = ["spt"]
2+
arrays_spt = ["095", "150", "220"]
3+
release_dir = "/global/cfs/cdirs/sobs/users/tlouis/spt_data/"
4+
alm_conv_spt = "IAU"
5+
6+
7+
binning_file = release_dir + "binning_spt.dat"
8+
apply_kspace_filter = False
9+
10+
deconvolve_pixwin = True
11+
pixwin_spt = {"pix": 'HEALPIX', "nside": 8192} #should be replaced with the more accurate pixel window that is provided
12+
13+
14+
apply_alm_mask = True
15+
lmax_mask = 4000
16+
alm_mask_spt_095 = release_dir + "ancillary_products/specific_to_c25/alm_mask_095ghz.fits"
17+
alm_mask_spt_150 = release_dir + "ancillary_products/specific_to_c25/alm_mask_150ghz.fits"
18+
alm_mask_spt_220 = release_dir + "ancillary_products/specific_to_c25/alm_mask_220ghz.fits"
19+
20+
niter = 0
21+
remove_mean = True
22+
binned_mcm = True
23+
lmax = 3500
24+
type = 'Dl'
25+
write_splits_spectra = True
26+
cov_T_E_only = False
27+
use_toeplitz_mcm = False
28+
use_toeplitz_cov = True
29+
30+
# maps
31+
32+
src_free_maps_spt = False
33+
34+
35+
freq_info_spt_095 = {"freq_tag": 95, "passband": "dummy.dat"}
36+
freq_info_spt_150 = {"freq_tag": 150, "passband": "dummy.dat"}
37+
freq_info_spt_220 = {"freq_tag": 220, "passband": "dummy.dat"}
38+
39+
40+
n_splits_spt = 2
41+
maps_spt_095 = [release_dir + "real_data_maps/half/half_bundle%d_095ghz_inpainted.fits" % (i) for i in range(n_splits_spt)]
42+
maps_spt_150 = [release_dir + "real_data_maps/half/half_bundle%d_150ghz_inpainted.fits" % (i) for i in range(n_splits_spt)]
43+
maps_spt_220 = [release_dir + "real_data_maps/half/half_bundle%d_220ghz_inpainted.fits" % (i) for i in range(n_splits_spt)]
44+
45+
46+
# these number are just the one associated with their beam at l=0
47+
# in act we use the convention that bl(l=0)=1, they don't do this in spt
48+
cal_spt_095 = 1 / 1.09647802
49+
cal_spt_150 = 1 / 1.06827578
50+
cal_spt_220 = 1 / 1.11338355
51+
52+
pol_eff_spt_095 = 1.0
53+
pol_eff_spt_150 = 1.0
54+
pol_eff_spt_220 = 1.0
55+
56+
beam_file_dir = release_dir + "ancillary_products/specific_to_c25/"
57+
58+
beam_T_spt_095 = beam_file_dir + "beam_c25_bl_095ghz.txt"
59+
beam_T_spt_150 = beam_file_dir + "beam_c25_bl_150ghz.txt"
60+
beam_T_spt_220 = beam_file_dir + "beam_c25_bl_220ghz.txt"
61+
62+
beam_pol_spt_095 = beam_file_dir + "beam_c25_bl_095ghz.txt"
63+
beam_pol_spt_150 = beam_file_dir + "beam_c25_bl_150ghz.txt"
64+
beam_pol_spt_220 = beam_file_dir + "beam_c25_bl_220ghz.txt"
65+
66+
window_dir = release_dir + "ancillary_products/generally_applicable/"
67+
68+
window_T_spt_095 = window_dir + "pixel_mask_apodized_borders_only.fits"
69+
window_pol_spt_095 = window_dir + "pixel_mask_apodized_borders_only.fits"
70+
71+
window_T_spt_150 = window_dir + "pixel_mask_apodized_borders_only.fits"
72+
window_pol_spt_150 = window_dir + "pixel_mask_apodized_borders_only.fits"
73+
74+
window_T_spt_220 = window_dir + "pixel_mask_apodized_borders_only.fits"
75+
window_pol_spt_220 = window_dir + "pixel_mask_apodized_borders_only.fits"
76+
77+
78+
cosmo_params = {"cosmomc_theta":0.0104059, "logA": 3.053, "ombh2": 0.02258, "omch2": 0.1242, "ns": 0.9666, "Alens": 1.0, "tau": 0.0567}
79+
80+

‎project/SPT_dev/spt_binning_file.py‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
import numpy as np
2+
import scipy.stats as stats
3+
4+
ell = np.arange(4001)
5+
6+
bin_edges = np.arange(0, 5000, 50)
7+
8+
lmax = bin_edges[1:].copy()
9+
lmin = bin_edges[:-1].copy()
10+
lmin[0]= 2
11+
lmin[1:] += 1
12+
13+
lmean = (lmin + lmax)/2
14+
15+
np.savetxt("binning_spt.dat", np.transpose([lmin, lmax, lmean]))
16+

‎project/SPT_dev/spt_get_alms.py‎

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
1+
"""
2+
This script compute all alms write them to disk.
3+
It uses the window function provided in the dictionnary file.
4+
Optionally, it applies a calibration to the maps, a kspace filter and deconvolve the CAR pixel window function.
5+
"""
6+
7+
import sys
8+
import time
9+
10+
import numpy as np
11+
import healpy as hp
12+
from pixell import enmap
13+
from pspipe_utils import kspace, log, misc, pspipe_list
14+
from pspy import pspy_utils, so_dict, so_mpi, sph_tools, so_map
15+
16+
17+
d = so_dict.so_dict()
18+
d.read_from_file(sys.argv[1])
19+
log = log.get_logger(**d)
20+
21+
surveys = d["surveys"]
22+
lmax = d["lmax"]
23+
deconvolve_pixwin = d["deconvolve_pixwin"]
24+
niter = d["niter"]
25+
apply_kspace_filter = d["apply_kspace_filter"]
26+
apply_alm_mask = d["apply_alm_mask"]
27+
28+
alms_dir = "alms"
29+
pspy_utils.create_directory(alms_dir)
30+
31+
n_ar, sv_list, ar_list = pspipe_list.get_arrays_list(d)
32+
33+
log.info(f"number of arrays for the mpi loop : {n_ar}")
34+
so_mpi.init(True)
35+
subtasks = so_mpi.taskrange(imin=0, imax=n_ar-1)
36+
37+
for task in subtasks:
38+
task = int(task)
39+
sv, ar = sv_list[task], ar_list[task]
40+
alm_conv = d[f"alm_conv_{sv}"]
41+
42+
log.info(f"[{task}] Computing alm for '{sv}' survey and '{ar}' array")
43+
44+
win_T = so_map.read_map(d[f"window_T_{sv}_{ar}"])
45+
win_pol = so_map.read_map(d[f"window_pol_{sv}_{ar}"])
46+
47+
window_tuple = (win_T, win_pol)
48+
49+
50+
if win_T.pixel == "CAR":
51+
win_kspace = so_map.read_map(d[f"window_kspace_{sv}_{ar}"])
52+
53+
if apply_kspace_filter:
54+
ks_f = d[f"k_filter_{sv}"]
55+
filter = kspace.get_kspace_filter(win_T, ks_f)
56+
57+
inv_pixwin_lxly = None
58+
if deconvolve_pixwin:
59+
if d[f"pixwin_{sv}"]["pix"] == "CAR":
60+
# compute the CAR pixel function in fourier space
61+
wy, wx = enmap.calc_window(win_T.data.shape, order=d[f"pixwin_{sv}"]["order"])
62+
inv_pixwin_lxly = (wy[:,None] * wx[None,:]) ** (-1)
63+
64+
65+
maps = d[f"maps_{sv}_{ar}"]
66+
cal, pol_eff = d[f"cal_{sv}_{ar}"], d[f"pol_eff_{sv}_{ar}"]
67+
68+
t0 = time.time()
69+
for k, map in enumerate(maps):
70+
71+
if win_T.pixel == "CAR":
72+
split = so_map.read_map(map, geometry=win_T.data.geometry)
73+
74+
if d[f"src_free_maps_{sv}"] == True:
75+
ps_map_name = map.replace("_srcfree.fits", ".fits")
76+
if ps_map_name == map:
77+
raise ValueError("name should contain srcfree, check map names!")
78+
ps_map = so_map.read_map(ps_map_name, geometry=win_T.data.geometry)
79+
ps_map.data -= split.data
80+
81+
82+
ps_mask = so_map.read_map(d[f"ps_mask_{sv}_{ar}"])
83+
ps_map.data *= ps_mask.data
84+
split.data += ps_map.data
85+
86+
if apply_kspace_filter:
87+
log.info(f"[{task}] apply kspace filter on {map}")
88+
split = kspace.filter_map(split,
89+
filter,
90+
win_kspace,
91+
inv_pixwin=inv_pixwin_lxly,
92+
weighted_filter=ks_f["weighted"],
93+
use_ducc_rfft=True)
94+
95+
else:
96+
log.info(f"[{task}] WARNING: no kspace filter is applied")
97+
if (deconvolve_pixwin) & (inv_pixwin_lxly is not None):
98+
split = so_map.fourier_convolution(split,
99+
inv_pixwin_lxly,
100+
window=win_kspace,
101+
use_ducc_rfft=True)
102+
103+
elif win_T.pixel == "HEALPIX":
104+
split = so_map.read_map(map)
105+
106+
split = split.calibrate(cal=cal, pol_eff=pol_eff)
107+
108+
if d["remove_mean"] == True:
109+
split = split.subtract_mean(window_tuple)
110+
111+
master_alms = sph_tools.get_alms(split, window_tuple, niter, lmax, alm_conv=alm_conv)
112+
113+
if apply_alm_mask == True:
114+
alm_mask = hp.read_alm(d[f"alm_mask_{sv}_{ar}"], hdu=1)
115+
alm_mask = hp.sphtfunc.resize_alm(alm_mask, d["lmax_mask"], d["lmax_mask"], lmax, lmax)
116+
master_alms *= alm_mask
117+
118+
np.save(f"{alms_dir}/alms_{sv}_{ar}_{k}.npy", master_alms)
119+
120+
log.info(f"[{task}] execution time {time.time() - t0} seconds")

‎project/SPT_dev/spt_get_mcm_and_bbl.py‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
"""
2+
This script computes the mode coupling matrices and the binning matrices Bbl
3+
for the different surveys and arrays.
4+
"""
5+
6+
import sys
7+
8+
from pspipe_utils import log, pspipe_list, misc
9+
from pspy import pspy_utils, so_dict, so_mcm, so_mpi, so_map
10+
11+
d = so_dict.so_dict()
12+
d.read_from_file(sys.argv[1])
13+
log = log.get_logger(**d)
14+
15+
mcm_dir = "mcms"
16+
pspy_utils.create_directory(mcm_dir)
17+
18+
surveys = d["surveys"]
19+
lmax = d["lmax"]
20+
binned_mcm = d["binned_mcm"]
21+
22+
if d["use_toeplitz_mcm"] == True:
23+
log.info("we will use the toeplitz approximation")
24+
l_exact, l_band, l_toep = 800, 2000, 2750
25+
else:
26+
l_exact, l_band, l_toep = None, None, None
27+
28+
n_mcms, sv1_list, ar1_list, sv2_list, ar2_list = pspipe_list.get_spectra_list(d)
29+
30+
log.info(f"number of mcm matrices to compute : {n_mcms}")
31+
so_mpi.init(True)
32+
subtasks = so_mpi.taskrange(imin=0, imax=n_mcms - 1)
33+
for task in subtasks:
34+
task = int(task)
35+
sv1, ar1, sv2, ar2 = sv1_list[task], ar1_list[task], sv2_list[task], ar2_list[task]
36+
37+
log.info(f"[{task:02d}] mcm matrix for {sv1}_{ar1} x {sv2}_{ar2}")
38+
39+
l, bl1 = misc.read_beams(d[f"beam_T_{sv1}_{ar1}"], d[f"beam_pol_{sv1}_{ar1}"])
40+
41+
win1_T = so_map.read_map(d[f"window_T_{sv1}_{ar1}"])
42+
win1_pol = so_map.read_map(d[f"window_pol_{sv1}_{ar1}"])
43+
44+
l, bl2 = misc.read_beams(d[f"beam_T_{sv2}_{ar2}"], d[f"beam_pol_{sv2}_{ar2}"])
45+
46+
win2_T = so_map.read_map(d[f"window_T_{sv2}_{ar2}"])
47+
win2_pol = so_map.read_map(d[f"window_pol_{sv2}_{ar2}"])
48+
49+
mbb_inv, Bbl = so_mcm.mcm_and_bbl_spin0and2(win1=(win1_T, win1_pol),
50+
win2=(win2_T, win2_pol),
51+
bl1=(bl1["T"], bl1["E"]),
52+
bl2=(bl2["T"], bl2["E"]),
53+
binning_file=d["binning_file"],
54+
niter=d["niter"],
55+
lmax=d["lmax"],
56+
type=d["type"],
57+
l_exact=l_exact,
58+
l_band=l_band,
59+
l_toep=l_toep,
60+
binned_mcm=binned_mcm,
61+
save_file=f"{mcm_dir}/{sv1}_{ar1}x{sv2}_{ar2}")

0 commit comments

Comments
 (0)