Skip to content

Commit bff1647

Browse files
author
Serena Giardiello
committed
Merge branch 'zach_piso' into serena_piso
2 parents bc07c69 + a474fcf commit bff1647

8 files changed

Lines changed: 804 additions & 117 deletions

File tree

project/SO/pISO/python/calibration/calibrator.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -55,15 +55,15 @@
5555
log.info(f'Calibrate spectra using calibs yaml file and {calib_test} test.')
5656
with open(args.calibs, "r") as file:
5757
calibs_dict: dict = yaml.safe_load(file)
58-
calibs_survey_arrays = list(calibs_dict['bestfits'].keys())
58+
calibs_survey_arrays = list(calibs_dict.keys())
5959

6060
poleff_suffix = '_poleff' if args.force_poleff else ''
6161

6262
# Compare all sv_ar with the ones in calibs yaml file and put non existing ones to 1.
6363
for sv_ar in survey_arrays:
6464
if sv_ar not in calibs_survey_arrays:
6565
log.info(f'{sv_ar} not in calibs yaml, setting it to 1.')
66-
calibs_dict['bestfits'][sv_ar] = {test: 1. for test in calib_tests}
66+
calibs_dict[sv_ar] = 1.
6767

6868
# Calibrate the spectra and save with _cal suffix
6969
for (sv1, ar1), (sv2, ar2) in itertools.combinations_with_replacement(survey_arrays_tuple, r=2):
@@ -75,7 +75,7 @@
7575
# Load & calib
7676
spec_filename_load = f'{spec_dir}/Dl_{sv_ar1}x{sv_ar2}_{spec_type}{poleff_suffix}.dat'
7777
ls, Dls = so_spectra.read_ps(spec_filename_load, spectra=spectra)
78-
Dls_cal = {spec: Dls[spec] / calibs_dict['bestfits'][sv_ar1][calib_test] / calibs_dict['bestfits'][sv_ar2][calib_test] for spec in spectra}
78+
Dls_cal = {spec: Dls[spec] * calibs_dict[sv_ar1] * calibs_dict[sv_ar2] for spec in spectra}
7979

8080
# Save with _cal suffix
8181
spec_filename_save = f'{spec_dir}/Dl_{sv_ar1}x{sv_ar2}_{spec_type}_calib{poleff_suffix}.dat'

project/SO/pISO/python/calibration/get_calibs.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
import pickle
99
import sys
1010
import yaml
11+
from scipy.optimize import curve_fit
1112
import scipy.stats as ss
1213

1314

project/SO/pISO/python/calibration/get_calibs_cross.py

Lines changed: 214 additions & 93 deletions
Large diffs are not rendered by default.

project/SO/pISO/python/calibration/get_polar_eff_LCDM.py

Lines changed: 45 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
import matplotlib.pyplot as plt
99
import getdist.plots as gdplt
1010
from cobaya.run import run
11+
import pickle
1112
import numpy as np
1213
import sys
1314
import yaml
@@ -25,15 +26,18 @@
2526
spec_dir = d['spec_dir']
2627
cov_dir = d['cov_dir']
2728
bestfit_dir = d['best_fits_dir']
28-
mcm_dir = f'old_cov/mcm/'
29+
# mcm_dir = f'old_cov/mcm/'
30+
mcm_dir = d['mcm_dir']
2931
planck_corr = False
3032
use_leakage_cov = False
3133

3234
if planck_corr:
3335
spec_dir = "spectra_leak_corr_planck_bias_corr"
3436

3537
output_dir = d['poleff_dir']
38+
chains_dir = f"{output_dir}/chains"
3639
pspy_utils.create_directory(output_dir)
40+
pspy_utils.create_directory(chains_dir)
3741
plots_dir = d['plots_dir'] + '/poleff/'
3842
pspy_utils.create_directory(plots_dir)
3943

@@ -135,9 +139,12 @@ def get_model(cmb_th, fg_th, Bbl, dust_amp, pol_eff, mode):
135139
spin_pair = "spin2xspin2" if spectrum == "EE" else ("spin0xspin2" if spectrum == "TE" else None)
136140
if spin_pair is None:
137141
raise ValueError("spectrum must be set to either 'EE' or 'TE'")
138-
Bbl = np.load(f"{mcm_dir}/{spec_name}_Bbl_{spin_pair}.npy")
139-
Bbl = Bbl[:n_bins, :lmax_paramfile]
140-
142+
Bbl = np.load(f"{mcm_dir}/{spec_name}_Bbl.npy")
143+
if d["binned_mcm"]:
144+
Bbl = Bbl[3, :n_bins, :lmax_paramfile] # ? Bbl is 00, 02, 20, ++, -- so we take ++ ?
145+
else:
146+
Bbl = Bbl[:n_bins, :lmax_paramfile] # ? Bbl is 00, 02, 20, ++, -- so we take ++ ?
147+
141148
# Get theory
142149
cmb_th = ps_th[spectrum][:lmax_paramfile]
143150
fg_th = fg_dict[spectrum.lower(), "dust", f"{sv_ar}", f"{sv_ar}"][:lmax_paramfile]
@@ -148,7 +155,7 @@ def loglike(pol_eff, dust_amp):
148155
residual = ps - theory
149156
chi2 = residual @ invcov @ residual
150157
return -0.5 * chi2
151-
158+
152159
loc, scale = dust_priors[spectrum]["loc"], dust_priors[spectrum]["scale"]
153160

154161
# Prepare MCMC sampling
@@ -162,7 +169,7 @@ def loglike(pol_eff, dust_amp):
162169
"min": 0.5,
163170
"max": 1.5
164171
},
165-
"latex": r"\epsilon_\mathrm{pol}^{%s}" % f"{sv_ar}".replace("_", "\_")
172+
"latex": r"\epsilon_\mathrm{pol}^{%s}" % f"{sv_ar}".replace("_",r"\_")
166173
},
167174
"dust_amp": {
168175
"prior": {
@@ -176,17 +183,17 @@ def loglike(pol_eff, dust_amp):
176183
"sampler": {
177184
"mcmc": {
178185
"max_tries": 10**6,
179-
"Rminus1_stop": 0.015,
180-
"Rminus1_cl_stop": 0.015
186+
"Rminus1_stop": 0.005,
187+
# "Rminus1_cl_stop": 0.015
181188
}
182189
},
183-
"output": f"{output_dir}/chain_{spectrum}_{sv_ar}",
190+
"output": f"{chains_dir}/chain_{spectrum}_{sv_ar}",
184191
"force": True,
185192
}
186193

187194
updated_info, sampler = run(info)
188195

189-
samples = loadMCSamples(f"{output_dir}/chain_{spectrum}_{sv_ar}", settings={"ignore_rows": 0.5})
196+
samples = loadMCSamples(f"{chains_dir}/chain_{spectrum}_{sv_ar}", settings={"ignore_rows": 0.5})
190197
pol_eff_mean[sv_ar, spectrum] = samples.mean("pol_eff")
191198
pol_eff_std[sv_ar, spectrum] = np.sqrt(samples.cov(["pol_eff"])[0, 0])
192199
dust_mean[sv_ar, spectrum] = samples.mean("dust_amp")
@@ -298,7 +305,7 @@ def loglike(pol_eff, dust_amp):
298305
plt.close()
299306

300307
# Save results to a yaml file
301-
poleffs_to_save = {
308+
poleffs_dict = {
302309
'bestfits' : {
303310
sv_ar: {
304311
mode: float(pol_eff_mean[sv_ar, mode])
@@ -315,10 +322,34 @@ def loglike(pol_eff, dust_amp):
315322
},
316323
}
317324

318-
file = open(f"{d['calib_dir']}/poleffs_dict.yaml", "w")
319-
yaml.dump(poleffs_to_save, file)
325+
default_calibs = {f'{sv}_{ar}': d[f"pol_eff_{sv}_{ar}"] for sv in d["surveys"] for ar in d[f"arrays_{sv}"]}
326+
results_dict = {
327+
sv_ar: (per_mode_poleff["EE"], per_mode_err["EE"])
328+
for (sv_ar, per_mode_poleff), per_mode_err in zip(poleffs_dict["bestfits"].items(), poleffs_dict["std"].values())
329+
}
330+
331+
poleff_to_save = {
332+
sv_ar: float(poleff)
333+
for sv_ar, (poleff, _) in results_dict.items()
334+
}
335+
stds_to_save = {
336+
sv_ar: float(std)
337+
for sv_ar, (_, std) in results_dict.items()
338+
}
339+
340+
file = open(f"{output_dir}/poleff_dict.yaml", "w")
341+
yaml.dump(poleff_to_save, file)
320342
file.close()
321343

344+
file = open(f"{output_dir}/poleff_errs_dict.yaml", "w")
345+
yaml.dump(stds_to_save, file)
346+
file.close()
347+
348+
with open(f"{output_dir}/poleff_dict.pickle", 'wb') as handle:
349+
pickle.dump(poleff_to_save, handle, protocol=pickle.HIGHEST_PROTOCOL)
350+
with open(f"{output_dir}/poleff_errs_dict.pickle", 'wb') as handle:
351+
pickle.dump(stds_to_save, handle, protocol=pickle.HIGHEST_PROTOCOL)
352+
322353

323354
# Plot and print results
324355
color_list = ["blue", "red", "green"]
@@ -330,7 +361,7 @@ def loglike(pol_eff, dust_amp):
330361
print(f"**************")
331362
p_eff, std = pol_eff_mean[sv_ar, mode], pol_eff_std[sv_ar, mode]
332363
print(f"{sv_ar} {mode} {p_eff} {std}")
333-
cal, std = poleffs_to_save["bestfits"][sv_ar][mode], poleffs_to_save["std"][sv_ar][mode]
364+
cal, std = poleffs_dict["bestfits"][sv_ar][mode], poleffs_dict["std"][sv_ar][mode]
334365
print(f"{mode}, cal: {cal:.5f}, sigma cal: {std:.5f}")
335366

336367
ax.errorbar(
Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
import yaml
2+
import numpy as np
3+
import pylab as plt
4+
import argparse
5+
import os
6+
from pspy import so_dict
7+
8+
9+
parser = argparse.ArgumentParser()
10+
parser.add_argument('odir', help="Where to save figures")
11+
parser.add_argument('--calib-yamls', nargs="+",
12+
help='A list of calib resutls yamls fils to plot')
13+
parser.add_argument('--paramfiles', nargs="+",
14+
help='A list of paramfiles to plot calibs (overrides --calib-yamls)')
15+
parser.add_argument('--arrays', nargs="+",
16+
help='list of arrays to plot (None = all LAT ISO)', default=None)
17+
parser.add_argument('--names', nargs="+", help='optionnally give names to these calibs', default=None)
18+
parser.add_argument('--colors', nargs="+", help='optionnally give colors to these calibs', default=None)
19+
parser.add_argument('--test', help='Which test to use', default='_all_tubes_Pl143')
20+
parser.add_argument('--ylabel', help='Ylabel for the figure', default='Calibration Factor')
21+
parser.add_argument('--plot-line', help='axhline for figure, is also the value to compare to get nsigmas', default=1, type=float)
22+
parser.add_argument('--plot-mean', help='plot mean of all data points (ASSUME UNCORRELATED POINTS)', default=False)
23+
24+
25+
args = parser.parse_args()
26+
27+
test_names = args.names or [os.path.basename(os.path.dirname(yaml_fn)) for yaml_fn in args.calib_yamls]
28+
test_colors = args.colors or [f"C{c}" for c, yaml_fn in enumerate(args.calib_yamls)]
29+
30+
defaults_calibs = {
31+
"i1_f090":1,
32+
"i1_f150":1,
33+
"i3_f090":1,
34+
"i3_f150":1,
35+
"i4_f090":1,
36+
"i4_f150":1,
37+
"i6_f090":1,
38+
"i6_f150":1,
39+
"c1_f220":1,
40+
"c1_f280":1,
41+
"i5_f220":1,
42+
"i5_f280":1,
43+
}
44+
45+
if args.paramfiles is None:
46+
calib_yamls = args.calib_yamls
47+
defaults = None
48+
else:
49+
defaults = {}
50+
calib_yamls = []
51+
for name, prmfl in zip(args.names, args.paramfiles):
52+
d = so_dict.so_dict()
53+
d.read_from_file(prmfl)
54+
calib_yamls.append(f"calib/{d['run_name']}/calibs_dict{args.test}.yaml")
55+
56+
results_dict = {}
57+
errs_dict = {}
58+
for name, yaml_fn in zip(args.names, calib_yamls):
59+
with open(yaml_fn, 'r') as f:
60+
results_dict[name] = yaml.safe_load(f)
61+
with open(yaml_fn.replace('_dict', '_errs_dict'), 'r') as f:
62+
errs_dict[name] = yaml.safe_load(f)
63+
64+
if args.paramfiles is not None:
65+
for name, prmfl in zip(args.names, args.paramfiles):
66+
d = so_dict.so_dict()
67+
d.read_from_file(prmfl)
68+
defaults[name] = {ar[-7:]: d[f"cal_{ar}"] for ar in results_dict[name].keys()}
69+
70+
arrays_set = args.arrays or ['c1_f220', 'c1_f280', 'i1_f090', 'i1_f150', 'i3_f090', 'i3_f150', 'i4_f090', 'i4_f150', 'i5_f220', 'i5_f280', 'i6_f090', 'i6_f150']
71+
72+
fig, ax = plt.subplots(figsize=(12, 6))
73+
ax.axhline(args.plot_line, color='grey', ls='--', lw=1)
74+
75+
x_shift = np.linspace(-.02*len(test_names), .02*len(test_names), len(test_names))
76+
77+
for j, (test, results_subdict) in enumerate(results_dict.items()):
78+
for i, (name, cal) in enumerate(results_subdict.items()):
79+
std = errs_dict[test][name]
80+
ax.errorbar(
81+
i + x_shift[j],
82+
cal,
83+
std,
84+
color=test_colors[j],
85+
marker=".",
86+
ls="",
87+
label=test if i == 0 else None,
88+
markersize=6.5,
89+
# markeredgewidth=2,
90+
fillstyle='none',
91+
)
92+
if defaults is not None:
93+
ax.plot(
94+
i + x_shift[j],
95+
defaults[test][arrays_set[i]],
96+
color=test_colors[j],
97+
marker="_",
98+
ls="",
99+
markersize=6.5,
100+
alpha=1,
101+
markeredgewidth=2,
102+
)
103+
if args.plot_mean:
104+
mean_var = 1 / np.sum([1 / errs_dict[test][name]**2 for name in results_subdict.keys()])
105+
mean_value = np.sum([results_subdict[name] / errs_dict[test][name]**2 for name in results_subdict.keys()]) * mean_var
106+
ax.errorbar(
107+
i+1 + x_shift[j],
108+
mean_value,
109+
np.sqrt(mean_var),
110+
color=test_colors[j],
111+
marker=".",
112+
ls="",
113+
markersize=6.5,
114+
# markeredgewidth=2,
115+
fillstyle='none',
116+
)
117+
print(f"{test}: {mean_value=:.5f}±{np.sqrt(mean_var):.5f}")
118+
# ax.set_ylim(.79, 1.02)
119+
ax.legend(fontsize=15)
120+
121+
xlabels = [ar[-7:] for ar in arrays_set]
122+
if args.plot_mean:
123+
xlabels += ["Mean"]
124+
x = np.arange(0, len(xlabels))
125+
ax.set_xticks(x, xlabels)
126+
ax.set_ylabel(args.ylabel, fontsize=18)
127+
plt.tight_layout()
128+
plt.savefig(f"{args.odir}/{args.ylabel.replace(" ", "_")}_summary_{'_'.join(test_names[:len(results_dict.keys())])}.pdf", bbox_inches="tight")
129+
plt.clf()
130+
plt.close()
131+
132+
# Claude wrote that
133+
def print_table(headers, rows, title=None):
134+
# Column widths: max of header or any cell value
135+
col_widths = [
136+
max(len(str(headers[i])), max((len(str(row[i])) for row in rows), default=0))
137+
for i in range(len(headers))
138+
]
139+
140+
# Box-drawing pieces
141+
top = "┌" + "┬".join("─" * (w + 2) for w in col_widths) + "┐"
142+
mid = "├" + "┼".join("─" * (w + 2) for w in col_widths) + "┤"
143+
bottom = "└" + "┴".join("─" * (w + 2) for w in col_widths) + "┘"
144+
145+
total_width = sum(w + 3 for w in col_widths) + 1
146+
147+
def row_line(cells, widths):
148+
return "│" + "│".join(f" {str(c):<{w}} " for c, w in zip(cells, widths)) + "│"
149+
150+
# Print
151+
if title:
152+
print("┌" + "─" * (total_width - 2) + "┐")
153+
print("│" + title.center(total_width - 2) + "│")
154+
155+
print(top)
156+
print(row_line(headers, col_widths))
157+
print(mid)
158+
for row in rows:
159+
print(row_line(row, col_widths))
160+
print(bottom)
161+
162+
for j, (test, results_subdict) in enumerate(results_dict.items()):
163+
for i, (name, cal) in enumerate(results_subdict.items()):
164+
std = errs_dict[test][name]
165+
166+
headers = ["params"] + [test for test in results_dict.keys()]
167+
rows = [[name] + [f"{results_dict[test][name]:.4f}±{errs_dict[test][name]:.4f} ({(results_dict[test][name] - args.plot_line) / errs_dict[test][name]:.1f}σ)" for test in results_dict.keys()] for name in results_dict[test].keys()]
168+
169+
print_table(headers, rows, title=args.ylabel)

project/SO/pISO/python/get_windows.py

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,6 @@
55
we also produce a window that include the pixel weighting.
66
The different masks are apodized.
77
We also produce a kspace-mask that will later be used for the kspace filtering operation, in order to remove the edges of the survey and avoid nasty pixels.
8-
MODIFIED VERSION, SKIP IVAR AND XLINK
98
"""
109

1110
import sys
@@ -23,6 +22,7 @@
2322
d.read_from_file(sys.argv[1])
2423
log = log.get_logger(**d)
2524

25+
use_weight_mask = True
2626

2727
surveys = d["surveys"]
2828
# the apodisation length of the final survey mask
@@ -137,14 +137,14 @@ def sa(emap):
137137
baseline_mask = so_map.car_template_from_shape_wcs(1, *template_geom, dtype=np.float32)
138138
baseline_mask.data = enmap.read_map(baseline_mask_fn, geometry=template_geom).astype(np.float32, copy=False)
139139
baseline_apod = d[f'baseline_apods_{winname}'][i]
140-
baseline_mask = so_window.create_apodization(
141-
baseline_mask, "C1", baseline_apod, use_rmax=True
142-
)
140+
if baseline_apod is not None:
141+
baseline_mask = so_window.create_apodization(
142+
baseline_mask, "C1", baseline_apod, use_rmax=True
143+
)
143144
my_masks["baseline"].data[:] *= baseline_mask.data
144145

145146
my_masks["baseline"].data = my_masks["baseline"].data.astype(np.float32, copy=False)
146-
my_masks["baseline"].write_map(f"{window_dir}/window_{winname}_{"baseline"}.fits")
147-
147+
my_masks["baseline"].write_map(f"{window_dir}/window_{winname}_baseline.fits")
148148
log.info(f"[{task}] joint baseline mask solid angle: {sa(my_masks["baseline"].data)}")
149149

150150
# Plot baseline and kspace windows

0 commit comments

Comments
 (0)