import numpy as np
import copy
import yaml
import pdb
W = "galaxy_density_xi"
GAMMA = "galaxy_shearDensity_xi_t"
GAMMAX = "galaxy_shearDensity_xi_x"
XIP = "galaxy_shear_xi_plus"
XIM = "galaxy_shear_xi_minus"
EE = "galaxy_shear_cl_ee"
DD = "galaxy_density_cl"
ED = "galaxy_shearDensity_cl_e"
types = {
W: ("theta", "lens", "lens"),
GAMMA: ("theta", "source", "lens"),
GAMMAX: ("theta", "source", "lens"),
XIP: ("theta", "source", "source"),
XIM: ("theta", "source", "source"),
EE: ("ell", "source", "source"),
DD: ("ell", "lens", "lens"),
ED: ("ell", "source", "lens"),
}
[docs]
def make_axis(i, j, nx, ny, axes):
import matplotlib.pyplot as plt
if i == 0 and j == 0:
shares = {}
elif j == 0:
shares = {"sharex": axes[0, 0]}
elif i == j:
shares = {"sharey": axes[i, 0]}
else:
shares = {"sharey": axes[i, 0], "sharey": axes[j, j]}
a = plt.subplot(ny, nx, i * ny + j + 1, **shares)
axes[i, j] = a
return a
[docs]
def apply_galaxy_bias_ggl(obs, theory, xi):
theory = copy.deepcopy(theory)
nbin_source = obs["nbin_source"]
nbin_lens = obs["nbin_lens"]
bias = np.zeros(nbin_lens)
# Best fit biases as a function of scale
for i in range(nbin_lens):
ell_obs, cl_obs = obs[DD, i, i]
ell_theory, cl_theory = theory[DD, i, i]
theory_at_obs = np.interp(ell_obs, ell_theory, cl_theory)
b2 = cl_obs / theory_at_obs
b = np.sqrt(b2)
mean_b = np.mean(b**0.5)
cl_theory *= mean_b**2
bias[i] = mean_b
print(f"Bias {i} = {mean_b:.2f}")
# Now apply to GGL
for i in range(nbin_source):
for j in range(nbin_lens):
_, theory_cl = theory[ED, i, j]
theory_cl *= bias[j]
return theory
[docs]
def full_3x2pt_plots(
sacc_files,
labels,
theory_sacc_files=None,
theory_labels=None,
xi=None,
fit_bias=False,
figures=None,
xlogscale=True,
ratios=False,
):
import sacc
sacc_data = []
for f in sacc_files:
if isinstance(f, sacc.Sacc):
sacc_data.append(f)
else:
sacc_data.append(sacc.Sacc.load_fits(f))
obs_data = [extract_observables_plot_data(s, label) for s, label in zip(sacc_data, labels)]
if theory_sacc_files is not None:
sacc_theory = []
for f in theory_sacc_files:
if isinstance(f, sacc.Sacc):
sacc_theory.append(f)
else:
sacc_theory.append(sacc.Sacc.load_fits(f))
obs_theory = [extract_observables_plot_data(s, label) for s, label in zip(sacc_theory, theory_labels)]
else:
obs_theory = []
if fit_bias:
if len(obs_theory) > 1:
print("warning - fitting to just the first set of theory spectra")
fitted_obs_theory = []
for obs, label in zip(obs_data, labels):
theory_fit = apply_galaxy_bias_ggl(obs, obs_theory[0], xi)
theory_fit["name"] = f"Fit to {label}"
fitted_obs_theory.append(theory_fit)
obs_theory = fitted_obs_theory
if figures is None:
figures = {}
for t in types:
f = figures.get(t)
output_figures = {}
for t in types:
if any(obs[t] for obs in obs_data):
print(f"Making Plot {t}")
f = figures.get(t)
output_figures[t] = make_plot(t, obs_data, obs_theory, fig=f, xlogscale=xlogscale, ratios=ratios)
return output_figures
[docs]
def axis_setup(a, i, j, ny, ymin, ymax, name):
import matplotlib.pyplot as plt
if j > 0:
plt.setp(a.get_yticklabels(), visible=False)
else:
a.set_ylabel(f"${name}$")
if i < ny - 1:
plt.setp(a.get_xticklabels(), visible=False)
if name.startswith(r"C_\ell"):
a.set_xlabel(r"$C_\ell$")
else:
a.set_xlabel(r"$\theta$ [arcmin]")
a.tick_params(axis="both", which="major", length=10, direction="in")
a.tick_params(axis="both", which="minor", length=5, direction="in")
# Fix
a.text(0.1, 0.1, f"{j} - {i}", transform=a.transAxes)
if i == j == 0:
a.legend()
a.set_ylim(ymin, ymax)
[docs]
def make_plot(corr, obs_data, obs_theory, fig=None, xlogscale=True, ratios=False):
import matplotlib.pyplot as plt
nbin_source = obs_data[0]["nbin_source"]
nbin_lens = obs_data[0]["nbin_lens"]
ny = nbin_source if types[corr][1] == "source" else nbin_lens
nx = nbin_source if types[corr][2] == "source" else nbin_lens
if corr == XIP:
name = r"\xi_+(\theta)"
ymin = 5e-7
ymax = 9e-5
auto_only = False
half_only = True
elif corr == XIM:
name = r"\xi_-(\theta)"
ymin = 5e-7
ymax = 9e-5
auto_only = False
half_only = True
elif corr == GAMMA:
ymin = 5e-7
ymax = 2e-2
name = r"\gamma_T(\theta)"
auto_only = False
half_only = False
elif corr == W:
ymin = 2e-4
ymax = 1e-0
name = r"w(\theta)"
auto_only = True
half_only = False
elif corr == EE:
name = r"C_\ell^{EE}"
ymin = 2e-12
ymax = 9e-8
auto_only = False
half_only = True
elif corr == ED:
ymin = 2e-10
ymax = 2e-6
name = r"C_\ell^{ED}"
auto_only = False
half_only = False
elif corr == DD:
ymin = 2e-8
ymax = 1e-4
name = r"C_\ell^{DD}"
auto_only = True
half_only = False
elif corr == GAMMAX:
ymin = 5e-7
ymax = 2e-2
name = r"\gamma_X(\theta)"
auto_only = False
half_only = False
if ratios:
name += "\ \mathrm{ ratios}"
plt.rcParams["font.size"] = 14
f = fig if fig is not None else plt.figure(figsize=(nx * 3.5, ny * 3))
ax = {}
axes = f.subplots(ny, nx, sharex="col", sharey="row", squeeze=False)
for i in range(ny):
if auto_only:
J = [i]
elif half_only:
J = range(i + 1)
else:
J = range(nx)
for j in range(nx):
a = axes[i, j]
if j not in J:
f.delaxes(a)
continue
for index, obs in enumerate(obs_data):
res = obs[(corr, i, j)]
if ratios:
res_theory = obs_theory[index][(corr, i, j)]
if len(res) == 2:
theta, xi = res
if ratios:
theta_th, xi_th = res_theory
a.plot(theta, xi / xi_th, label="TXPipe Data/CCL Theory")
a.axhline(y=1, color="k", ls=":")
else:
if xlogscale:
(l,) = a.loglog(theta, xi, "x", label=obs["name"])
a.loglog(theta, -xi, "s", color=l.get_color())
else:
(l,) = a.plot(theta, xi, "x", label=obs["name"])
a.plot(theta, -xi, "s", color=l.get_color())
a.set_yscale("log")
else:
theta, xi, cov = res
err = cov.diagonal() ** 0.5
if ratios:
theta_th, xi_th = res_theory
a.errorbar(
theta,
xi / xi_th,
err / abs(xi_th),
fmt=".",
label="TXPipe Data/CCL Theory",
)
a.axhline(y=1, color="k", ls=":")
else:
a.errorbar(theta, xi, err, fmt=".", label=obs["name"], capsize=5)
a.set_yscale("log")
if xlogscale:
a.set_xscale("log")
if not ratios:
# plot theory
for theory in obs_theory:
theta, xi = theory[(corr, i, j)]
if xlogscale:
a.loglog(theta, xi, "-", label=theory["name"])
else:
a.plot(theta, xi, "-", label=theory["name"])
if ratios:
axis_setup(a, i, j, ny, 0.6, 1.4, name)
else:
axis_setup(a, i, j, ny, ymin, ymax, name)
if corr in [EE, ED, DD]:
a.set_xlim(90, 1500)
f.suptitle(rf"TXPipe ${name}$")
# plt.tight_layout()
# f.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.subplots_adjust(wspace=0.05, hspace=0.05)
return plt.gcf()
[docs]
def smooth_nz(nz):
return np.convolve(nz, np.exp(-0.5 * np.arange(-4, 5) ** 2) / 2**2, mode="same")
"""
This whole function below can be removed now since it is not used anymore.
"""
[docs]
def make_theory_plot_data(data, cosmo, obs, label, smooth=True, xi=None, ratios=False):
import pyccl
theory = {"name": label}
xi = ("galaxy_density_xi" in data.get_data_types()) if xi is None else xi
nbin_source = obs["nbin_source"]
nbin_lens = obs["nbin_lens"]
ell = np.unique(np.logspace(np.log10(2), 5, 400).astype(int))
tracers = {}
print("Cosmological parameters used in plots:")
print(cosmo)
# Make the lensing tracers
for i in range(nbin_source):
name = f"source_{i}"
Ti = data.get_tracer(name)
nz = smooth_nz(Ti.nz) if smooth else Ti.nz
print("smooth:", smooth)
# Convert to CCL form
tracers[name] = pyccl.WeakLensingTracer(cosmo, (Ti.z, nz))
# And the clustering tracers
for i in range(nbin_lens):
name = f"lens_{i}"
Ti = data.get_tracer(name)
nz = smooth_nz(Ti.nz) if smooth else Ti.nz
# Convert to CCL form
tracers[name] = pyccl.NumberCountsTracer(cosmo, has_rsd=False, dndz=(Ti.z, nz), bias=(Ti.z, np.ones_like(Ti.z)))
for i in range(nbin_source):
for j in range(i + 1):
print(f"Computing theory lensing-lensing ({i},{j})")
# uncomment to take ratio in Fourier space to get the exact ells
if ratios and (not xi):
if i == 0 and j == 0:
ell, _, _ = obs[("galaxy_shearDensity_cl_e", i, j)] # to get the ratio between theory and data
# compute power spectra
cl = pyccl.angular_cl(cosmo, tracers[f"source_{i}"], tracers[f"source_{j}"], ell)
theory[(EE, i, j)] = ell, cl
# Optionally also compute correlation functions
if xi:
theta, *_ = obs[(XIP, i, j)]
theory[(XIP, i, j)] = theta, pyccl.correlation(cosmo, ell, cl, theta / 60, corr_type="L+")
theory[(XIM, i, j)] = theta, pyccl.correlation(cosmo, ell, cl, theta / 60, corr_type="L-")
for i in range(nbin_lens):
print(f"Computing theory density-density ({i},{i})")
# compute power spectra
cl = pyccl.angular_cl(cosmo, tracers[f"lens_{i}"], tracers[f"lens_{i}"], ell)
theory[(DD, i, i)] = ell, cl
# Optionally also compute correlation functions
if xi:
theta, *_ = obs[(W, i, i)]
theory[W, i, i] = theta, pyccl.correlation(cosmo, ell, cl, theta / 60, corr_type="GG")
for i in range(nbin_source):
for j in range(nbin_lens):
print(f"Computing theory lensing-density (S{i},L{j})")
# compute power spectra
cl = pyccl.angular_cl(cosmo, tracers[f"source_{i}"], tracers[f"lens_{j}"], ell)
theory[(ED, i, j)] = ell, cl
# Optionally also compute correlation functions
if xi:
theta, *_ = obs[(GAMMA, i, j)]
theory[GAMMA, i, j] = theta, pyccl.correlation(cosmo, ell, cl, theta / 60, corr_type="GL")
return theory