Skip to content

Commit 6d84882

Browse files
author
Merry Duparc
committed
update calib script
1 parent a148e72 commit 6d84882

File tree

1 file changed

+160
-47
lines changed

1 file changed

+160
-47
lines changed

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

Lines changed: 160 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
import pickle
99
import sys
1010
import yaml
11+
import scipy.stats as ss
1112

1213

1314
def get_proj_pattern(test, map_set, ref_map_set):
@@ -150,17 +151,6 @@ def get_proj_indices(test):
150151

151152
lmin, lmax = calib_infos["ell_ranges"][map_set]
152153
id = np.where((lb >= lmin) & (lb <= lmax))[0]
153-
consistency.plot_residual(
154-
lb,
155-
res_spectrum,
156-
{"analytical": res_cov},
157-
"TT",
158-
f"{map_set} {test}",
159-
f"{plot_output_dir}/residual_{name}_before",
160-
lrange=id,
161-
l_pow=1,
162-
ylims=calib_infos["y_lims_plot_TT"],
163-
)
164154

165155
# Calibrate the spectra
166156
cal_mean, cal_std = consistency.get_calibration_amplitudes(
@@ -170,6 +160,7 @@ def get_proj_indices(test):
170160
"TT",
171161
id,
172162
f"{chains_dir}/{name}",
163+
Rminus1_cl_stop=0.1,
173164
)
174165

175166
results_dict[test, map_set] = {
@@ -179,93 +170,215 @@ def get_proj_indices(test):
179170
}
180171

181172
calib_vec = np.array([cal_mean**2, cal_mean, 1])
182-
res_spectrum, res_cov = consistency.project_spectra_vec_and_cov(
173+
res_spectrum_cal, res_cov_cal = consistency.project_spectra_vec_and_cov(
183174
spec_vec, full_cov, proj_pattern, calib_vec=calib_vec
184175
)
185176

186177
np.savetxt(
187178
f"{residual_output_dir}/residual_{name}_after.dat",
188-
np.array([lb, res_spectrum]).T,
189-
)
190-
consistency.plot_residual(
191-
lb,
192-
res_spectrum,
193-
{"analytical": res_cov},
194-
"TT",
195-
f"{map_set} {test}",
196-
f"{plot_output_dir}/residual_{name}_after",
197-
lrange=id,
198-
l_pow=1,
199-
ylims=calib_infos[f"y_lims_plot_TT"],
179+
np.array([lb, res_spectrum_cal]).T,
200180
)
181+
182+
### THIS REPLACES CONSISTENCY.PLOT_RESIDUAL
183+
expected_res = 0.
184+
remove_dof = 0.
185+
186+
res_th = np.ones(len(lb)) * expected_res
187+
chi2 = (res_spectrum[id] - res_th[id]) @ np.linalg.inv(res_cov[np.ix_(id, id)]) @ (res_spectrum[id] - res_th[id])
188+
ndof = len(lb[id]) - remove_dof
189+
pte = 1 - ss.chi2(ndof).cdf(chi2)
190+
191+
chi2_cal = (res_spectrum_cal[id] - res_th[id]) @ np.linalg.inv(res_cov_cal[np.ix_(id, id)]) @ (res_spectrum_cal[id] - res_th[id])
192+
ndof_cal = len(lb[id]) - remove_dof
193+
pte_cal = 1 - ss.chi2(ndof_cal).cdf(chi2_cal)
201194

202195
fig, ax = plt.subplots(
203-
2, gridspec_kw={"hspace": 0, "height_ratios": (2, 1)}, figsize=(8, 6)
196+
2,
197+
gridspec_kw={"hspace": 0, "height_ratios": (2, 1)},
198+
figsize=(8, 6),
199+
sharex=True,
204200
)
205201

206-
proj_indices = get_proj_indices(test)
202+
ax[0].set_yscale("log")
207203

204+
proj_indices = get_proj_indices(test)
205+
ell_pow = 1
208206
ax[0].errorbar(
209-
lb,
207+
lb - 6,
210208
ps_dict[
211209
spectra_for_cal[proj_indices[0]][0],
212210
spectra_for_cal[proj_indices[0]][1],
213211
"TT",
214-
],
212+
] * lb**ell_pow,
213+
np.sqrt(
214+
cov_dict[
215+
(
216+
spectra_for_cal[proj_indices[0]][0],
217+
spectra_for_cal[proj_indices[0]][1],
218+
"TT",
219+
),
220+
(
221+
spectra_for_cal[proj_indices[0]][0],
222+
spectra_for_cal[proj_indices[0]][1],
223+
"TT",
224+
),
225+
].diagonal()
226+
) * lb**ell_pow,
215227
color="grey",
216228
label=test[:3],
217-
ls="--",
229+
ls="",
230+
alpha=.6,
231+
marker='.',
232+
mfc='white',
233+
mec='tab:grey',
234+
linewidth=.5,
235+
elinewidth=1,
218236
)
219237
ax[0].errorbar(
220-
lb,
238+
lb + 2,
221239
ps_dict[
222240
spectra_for_cal[proj_indices[1]][0],
223241
spectra_for_cal[proj_indices[1]][1],
224242
"TT",
225-
],
226-
color="blue",
243+
] * lb**ell_pow,
244+
np.sqrt(
245+
cov_dict[
246+
(
247+
spectra_for_cal[proj_indices[1]][0],
248+
spectra_for_cal[proj_indices[1]][1],
249+
"TT",
250+
),
251+
(
252+
spectra_for_cal[proj_indices[1]][0],
253+
spectra_for_cal[proj_indices[1]][1],
254+
"TT",
255+
),
256+
].diagonal()
257+
) * lb**ell_pow,
258+
color="tab:blue",
227259
label=test[-3:],
228-
ls="--",
260+
ls="",
261+
alpha=.6,
262+
marker='.',
263+
mfc='white',
264+
mec='tab:blue',
265+
linewidth=.5,
266+
elinewidth=1,
229267
)
230268

231269
ax[0].errorbar(
232-
lb,
270+
lb - 2,
233271
ps_dict[
234272
spectra_for_cal[proj_indices[0]][0],
235273
spectra_for_cal[proj_indices[0]][1],
236274
"TT",
237-
],
275+
]
276+
* calib_vec[proj_indices[0]]
277+
* lb**ell_pow,
278+
np.sqrt(
279+
cov_dict[
280+
(
281+
spectra_for_cal[proj_indices[0]][0],
282+
spectra_for_cal[proj_indices[0]][1],
283+
"TT",
284+
),
285+
(
286+
spectra_for_cal[proj_indices[0]][0],
287+
spectra_for_cal[proj_indices[0]][1],
288+
"TT",
289+
),
290+
].diagonal()
291+
) * calib_vec[proj_indices[0]]
292+
* lb**ell_pow,
238293
color="grey",
239-
label=test[:3],
240-
ls="--",
294+
label=test[:3] + ' cal',
295+
ls="-",
296+
marker='.',
297+
mfc='white',
298+
mec='tab:grey',
299+
linewidth=.5,
300+
elinewidth=1,
241301
)
242302
ax[0].errorbar(
243-
lb,
303+
lb + 6,
244304
ps_dict[
245305
spectra_for_cal[proj_indices[1]][0],
246306
spectra_for_cal[proj_indices[1]][1],
247307
"TT",
248-
],
249-
color="blue",
250-
label=test[-3:],
251-
ls="--",
308+
]
309+
* calib_vec[proj_indices[1]]
310+
* lb**ell_pow,
311+
np.sqrt(
312+
cov_dict[
313+
(
314+
spectra_for_cal[proj_indices[1]][0],
315+
spectra_for_cal[proj_indices[1]][1],
316+
"TT",
317+
),
318+
(
319+
spectra_for_cal[proj_indices[1]][0],
320+
spectra_for_cal[proj_indices[1]][1],
321+
"TT",
322+
),
323+
].diagonal()
324+
) * calib_vec[proj_indices[1]]
325+
* lb**ell_pow,
326+
color="tab:blue",
327+
label=test[-3:]+ ' cal',
328+
ls="-",
329+
marker='.',
330+
mfc='white',
331+
mec='tab:blue',
332+
# linewidth=.5,
333+
# elinewidth=1,
252334
)
335+
336+
ax[1].errorbar(lb, res_spectrum / np.sqrt(res_cov.diagonal()),
337+
yerr=np.sqrt(res_cov.diagonal()) / np.sqrt(res_cov.diagonal()),
338+
ls="None", marker = ".",
339+
linewidth=.5,
340+
alpha=.5,
341+
color="tab:blue",
342+
label=f"{test} [$\chi^2 = {{{chi2:.1f}}}/{{{ndof}}}$ (${{{pte:.4f}}}$)]")
343+
ax[1].errorbar(lb, res_spectrum_cal / np.sqrt(res_cov_cal.diagonal()),
344+
yerr=np.sqrt(res_cov_cal.diagonal()) / np.sqrt(res_cov_cal.diagonal()),
345+
ls="None", marker = ".",
346+
color="tab:blue",
347+
label=f"{test} cal [$\chi^2 = {{{chi2_cal:.1f}}}/{{{ndof_cal}}}$ (${{{pte_cal:.4f}}}$)]")
348+
ax[1].axhline(0, color='black', zorder=-10)
349+
350+
ax[1].axvspan(xmin=0, xmax=lmin,
351+
color="gray", alpha=0.7, zorder=-20)
352+
ax[1].axvspan(xmin=lmax, xmax=10000,
353+
color="gray", alpha=0.7, zorder=-20)
354+
355+
ax[0].legend(title=f'A={ref_map_set} B={map_set}')
356+
ax[0].set_xlim(0, lmax + 500)
357+
ax[0].set_ylim(4e5, 4e6)
358+
ax[0].set_ylabel(fr"$\ell^{{{ell_pow}}} D_\ell^\mathrm{{TT}}$", fontsize=18)
359+
ax[0].set_title(f'cal={cal_mean:.3f}+-{cal_std:.3f}')
360+
361+
ax[1].legend(loc='lower right')
362+
ax[1].set_ylim(-8, 7)
363+
ax[1].set_ylabel(fr"$\Delta D_\ell^\mathrm{{TT}} / \sigma(\Delta D_\ell^\mathrm{{TT}})$", fontsize=16)
364+
ax[1].set_xlabel(r"$\ell$", fontsize=20)
253365

254366
plt.savefig(f"{plot_output_dir}/calib_full_{name}.png")
367+
plt.close()
255368

256369

257370
# plot the cal factors
258371
color_list = ["blue", "red", "green"]
259-
372+
fig, ax = plt.subplots(figsize=(8, 6))
260373
for i, (map_set, ref_map_set) in enumerate(calib_infos[f"ref_map_sets"].items()):
261374
print(f"**************")
262375
print(f"calibration {map_set} with {ref_map_set}")
263376

264377
for j, test in enumerate(tests):
265378
cal, std = results_dict[test, map_set]["calibs"]
266-
print(f"{test}, cal: {cal}, sigma cal: {std}")
379+
print(f"{test}, cal: {cal:.4f}, sigma cal: {std:.5f}")
267380

268-
plt.errorbar(
381+
ax.errorbar(
269382
i + 0.9 + j * 0.1,
270383
cal,
271384
std,
@@ -278,10 +391,10 @@ def get_proj_indices(test):
278391
)
279392

280393
if i == 0:
281-
plt.legend(fontsize=15)
394+
ax.legend(fontsize=15)
282395

283396
x = np.arange(1, len(calib_infos[f"ref_map_sets"]) + 1)
284-
plt.xticks(x, calib_infos[f"ref_map_sets"].keys())
397+
ax.set_xticks(x, calib_infos[f"ref_map_sets"].keys())
285398
# plt.ylim(0.967, 1.06)
286399
plt.tight_layout()
287400
plt.savefig(f"{plot_output_dir}/calibs_summary.pdf", bbox_inches="tight")

0 commit comments

Comments
 (0)