Skip to content

Commit 43283c7

Browse files
authored
make plots optional for gp smoothing of cov (#185)
1 parent 0046119 commit 43283c7

1 file changed

Lines changed: 65 additions & 47 deletions

File tree

project/ACT_DR6/python/montecarlo/get_mc_corrected_xarrays_covmat_gp.py

Lines changed: 65 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -23,8 +23,17 @@
2323
help="Correct using the cov that used only use these simulations")
2424
parser.add_argument("--iStop", type=int, default=None,
2525
help="Correct using the cov that used only use these simulations")
26+
parser.add_argument("--plot-all", action='store_true', dest='plot_all',
27+
help="Make plots of the GP fits for both the main and block diagonals")
28+
parser.add_argument("--plot-diag", action='store_true', dest='plot_diag',
29+
help="Make plots of the GP fits for the main diagonals only")
2630
args = parser.parse_args()
2731

32+
plot_all = args.plot_all
33+
plot_diag = args.plot_diag
34+
if plot_all and plot_diag:
35+
raise ValueError('plot-all and plot-diag cannot be both True')
36+
2837
d = so_dict.so_dict()
2938
d.read_from_file(args.paramfile)
3039

@@ -33,8 +42,10 @@
3342
covariances_dir = "covariances"
3443
plot_dir = "plots/x_ar_cov"
3544

36-
pspy_utils.create_directory(plot_dir)
37-
pspy_utils.create_directory(plot_dir + '/off_diags')
45+
if plot_all or plot_diag:
46+
pspy_utils.create_directory(plot_dir)
47+
if not plot_diag:
48+
pspy_utils.create_directory(plot_dir + '/off_diags')
3849

3950
binning_file = d["binning_file"]
4051
lmax = d["lmax"]
@@ -128,53 +139,60 @@
128139
np.save(os.path.join(covariances_dir, "x_ar_final_cov_sim_gp.npy"), corrected_mc_cov)
129140

130141
# make plots
131-
for i in range(len(keys)):
132-
name_i, spec_i = keys[i]
133-
for j in range(i, len(keys)):
134-
name_j, spec_j = keys[j]
142+
if plot_all or plot_diag:
143+
for i in range(len(keys)):
144+
name_i, spec_i = keys[i]
135145

136-
if gprs[i, j] is not None:
137-
kern = gprs[i, j].kernel_
138-
else:
139-
kern = ''
140-
print(name_i, spec_i, name_j, spec_j, kern)
141-
142-
sel = np.s_[i*n_bins:(i+1)*n_bins, j*n_bins:(j+1)*n_bins]
143-
plt.errorbar(bin_mean, np.diag(mc_cov_anaflat[sel]), np.diag(var_mc_cov_anaflat[sel])**.5, zorder=0, label='MC')
144-
plt.plot(bin_mean, np.diag(smoothed_mc_cov_anaflat[sel]), zorder=1, label='GP(MC)')
145-
146-
# get idxs
147-
idxs_i = idx_arrs_by_block[i]
148-
idxs_j = idx_arrs_by_block[j]
149-
if idxs_i is not None and idxs_j is not None:
150-
idxs = np.intersect1d(idxs_i, idxs_j)
151-
elif idxs_i is None:
152-
idxs = idxs_j # idxs_i is "all idxs" so use idxs_j (which might also be "all idxs")
153-
else:
154-
idxs = idxs_i # idxs_j is "all idxs" so use idxs_i
155-
if idxs is None:
156-
idxs = np.arange(n_bins, dtype=int)
157-
158-
if len(idxs) > 0:
159-
plt.axvline(bin_mean[idxs[0]], linestyle='--', color='k')
160-
plt.axvline(bin_mean[idxs[-1]], linestyle='--', color='k')
161-
if i == j:
162-
plt.ylim(.5, 1.5)
146+
if plot_diag:
147+
js = [i]
163148
else:
164-
plt.ylim(-.1, .1)
165-
plt.grid()
166-
plt.xlabel("l")
167-
plt.ylabel("ratio")
168-
plt.legend()
169-
plt.title(f"Cov({name_i}_{spec_i}, {name_j}_{spec_j})\n{kern}")
170-
if args.iStart is not None:
171-
if i == j:
172-
plt.savefig(os.path.join(plot_dir, f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}_{iStart}_{iStop}.png"))
149+
js = range(i, len(keys))
150+
151+
for j in js:
152+
name_j, spec_j = keys[j]
153+
154+
if gprs[i, j] is not None:
155+
kern = gprs[i, j].kernel_
173156
else:
174-
plt.savefig(os.path.join(plot_dir, 'off_diags', f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}_{iStart}_{iStop}.png"))
175-
else:
157+
kern = ''
158+
print(name_i, spec_i, name_j, spec_j, kern)
159+
160+
sel = np.s_[i*n_bins:(i+1)*n_bins, j*n_bins:(j+1)*n_bins]
161+
plt.errorbar(bin_mean, np.diag(mc_cov_anaflat[sel]), np.diag(var_mc_cov_anaflat[sel])**.5, zorder=0, label='MC')
162+
plt.plot(bin_mean, np.diag(smoothed_mc_cov_anaflat[sel]), zorder=1, label='GP(MC)')
163+
164+
# get idxs
165+
idxs_i = idx_arrs_by_block[i]
166+
idxs_j = idx_arrs_by_block[j]
167+
if idxs_i is not None and idxs_j is not None:
168+
idxs = np.intersect1d(idxs_i, idxs_j)
169+
elif idxs_i is None:
170+
idxs = idxs_j # idxs_i is "all idxs" so use idxs_j (which might also be "all idxs")
171+
else:
172+
idxs = idxs_i # idxs_j is "all idxs" so use idxs_i
173+
if idxs is None:
174+
idxs = np.arange(n_bins, dtype=int)
175+
176+
if len(idxs) > 0:
177+
plt.axvline(bin_mean[idxs[0]], linestyle='--', color='k')
178+
plt.axvline(bin_mean[idxs[-1]], linestyle='--', color='k')
176179
if i == j:
177-
plt.savefig(os.path.join(plot_dir, f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}.png"))
180+
plt.ylim(.5, 1.5)
181+
else:
182+
plt.ylim(-.1, .1)
183+
plt.grid()
184+
plt.xlabel("l")
185+
plt.ylabel("ratio")
186+
plt.legend()
187+
plt.title(f"Cov({name_i}_{spec_i}, {name_j}_{spec_j})\n{kern}")
188+
if args.iStart is not None:
189+
if i == j:
190+
plt.savefig(os.path.join(plot_dir, f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}_{iStart}_{iStop}.png"))
191+
else:
192+
plt.savefig(os.path.join(plot_dir, 'off_diags', f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}_{iStart}_{iStop}.png"))
178193
else:
179-
plt.savefig(os.path.join(plot_dir, 'off_diags', f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}.png"))
180-
plt.close()
194+
if i == j:
195+
plt.savefig(os.path.join(plot_dir, f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}.png"))
196+
else:
197+
plt.savefig(os.path.join(plot_dir, 'off_diags', f"GP_MC_covmat_smooth_{name_i}_{spec_i}_{name_j}_{spec_j}.png"))
198+
plt.close()

0 commit comments

Comments
 (0)