<jupyter_start><jupyter_text>LENSING ANGULAR POWER SPECTRAThis notebook tests ``corfu.py`` functions. It compares the convergence and shear angular power spectra against cosmoSIS, using DES-Y1 data.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
plt.rcParams.update({
    'text.usetex': False,
    'font.family': 'serif',
    'legend.frameon': False,
    'legend.handlelength': 1.5,
})<jupyter_output><empty_output><jupyter_text>1. Load CosmoSIS quantities * Metadata<jupyter_code>nbin = 4
bin_a, bin_b = np.tril_indices(nbin)
bin_a += 1
bin_b += 1<jupyter_output><empty_output><jupyter_text>* Cosmology<jupyter_code>cosmo = {}
with open('../data/des-y1-test/cosmological_parameters/values.txt') as cosmo_values:
    for line in cosmo_values:
        if line:
            key, val = line.partition('=')[::2]
            cosmo[key.strip()] = float(val)<jupyter_output><empty_output><jupyter_text>* Distance functions<jupyter_code>zdM = np.loadtxt('../data/des-y1-test/distances/z.txt')
dM = np.loadtxt('../data/des-y1-test/distances/d_m.txt')<jupyter_output><empty_output><jupyter_text>* Matter power spectrum<jupyter_code>zp = np.loadtxt('../data/des-y1-test/matter_power_nl/z.txt')
k_h = np.loadtxt('../data/des-y1-test/matter_power_nl/k_h.txt')
p_h = np.loadtxt('../data/des-y1-test/matter_power_nl/p_k.txt')
xp = np.interp(zp, zdM, dM)
k0, kf = k_h[0]*(cosmo['hubble']/100), k_h[-1]*(cosmo['hubble']/100)
k_h2 = np.logspace(np.log10(k0), np.log10(kf), 1024)
p_h2 = np.exp([np.interp(np.log(k_h2), np.log(k_h), np.log(p)) for p in p_h])
k = k_h2*cosmo['h0']
p = p_h2*cosmo['h0']**(-3)
p_geom = np.sqrt(p[:, None]*p[None, :])
plt.figure(figsize=(5,3))
plt.loglog(k, p[0], lw=1)
plt.loglog(k, p_geom[0][0], '--',lw=2)
plt.xlabel(r'$k\; (1/Mpc)$')
plt.ylabel(r'$P\; (Mpc^3)$')
plt.show()<jupyter_output><empty_output><jupyter_text>* Angula power spectrum<jupyter_code>l_csis_ = np.loadtxt('../data/des-y1-test/shear_cl/ell.txt')
cl_csis_ = [np.loadtxt('../data/des-y1-test/shear_cl/bin_%d_%d.txt' % bins) for bins in zip(bin_a, bin_b)]
l_csis = np.arange(1, 100000.1, 1)
cl_csis = [np.interp(np.log(l_csis), np.log(l_csis_), cl) for cl in cl_csis_]<jupyter_output><empty_output><jupyter_text>2. The correlation function<jupyter_code>import corfu
r_limb, xi_limb = corfu.ptoxi(k, p, q=0, limber=True)
r_geom, xi_geom = corfu.ptoxi(k, p_geom, q=0)
plt.figure(figsize=(6,4))
plt.loglog(r_limb, +xi_limb[0], '--b', label='limber linear')
plt.loglog(r_limb, -xi_limb[0], ':b')
plt.loglog(r_geom, +xi_geom[0,0], '--r', label='geometric linear')
plt.loglog(r_geom, -xi_geom[0,0], ':r')

plt.legend()
plt.xlabel('r')
plt.ylabel(r'$\xi(r)$')
plt.show()<jupyter_output><empty_output><jupyter_text>3. Lensing filters<jupyter_code>THREE_2C2 = 1.669e-11
INVERSE_C = 3.333e-6<jupyter_output><empty_output><jupyter_text>* Redshift distribution of galaxies<jupyter_code>zn = np.loadtxt('../data/des-y1-test/nz_source/z.txt')
nz = [np.loadtxt('../data/des-y1-test/nz_source/bin_%d.txt' % i) for i in range(1, nbin+1)]
xf = np.interp(zn, zdM, dM)<jupyter_output><empty_output><jupyter_text>* Lensing efficiency<jupyter_code>def lensing_efficiency(x, zx, nz):
    TINY = np.finfo(0.).tiny
    return np.trapz(np.clip(x - x[:,None], 0, None) / (x + TINY) * nz, zx)
q = [lensing_efficiency(xf, zn, n) for n in nz]<jupyter_output><empty_output><jupyter_text>* Convergence<jupyter_code>def filter_convergence(x, zx, lens_efficiency):
    factor = THREE_2C2 * np.square(cosmo['hubble']) * cosmo['omega_m']
    return factor * (1.0 + zx) * lens_efficiency * x
fc = [filter_convergence(xf, zn, qq) for qq in q]
for qq in fc:
    plt.plot(xf, qq)
plt.show()<jupyter_output><empty_output><jupyter_text>4. Angular correlation function<jupyter_code>theta = np.logspace(-3, np.log10(np.pi), 512)
theta_arcmin = np.degrees(theta)*60<jupyter_output><empty_output><jupyter_text>* Convergence<jupyter_code>w_limb = [corfu.eqt(theta, (xf, fc[a-1]*fc[b-1]), (xp, r_limb, xi_limb)) for a, b in zip(bin_a, bin_b)]
w_geom = [corfu.uneqt(theta, (xf, fc[a-1]), (xf, fc[b-1]), (xp, xp, r_geom, xi_geom), True) for a, b in zip(bin_a, bin_b)]
w_csis = [corfu.cltow(cl, theta) for cl in cl_csis]
fig, axes = plt.subplots(4, 4, figsize=(9.5, 9.5), sharex=True, sharey=True)

for ax in axes.ravel():
    ax.axis('off')

for i, (a, b) in enumerate(zip(bin_a, bin_b)):
    ax = axes[a-1, b-1]
    ax.axis('on')
    ax.loglog(theta_arcmin, +w_csis[i], 'b', label='CosmoSIS', lw=1)
    ax.loglog(theta_arcmin, -w_csis[i], '--b', lw=1)
    ax.loglog(theta_arcmin, +w_limb[i], 'k', label='Limber', lw=1)
    ax.loglog(theta_arcmin, -w_limb[i], '--k', lw=1)
    ax.loglog(theta_arcmin, +w_geom[i], 'r', label='Exact', lw=2)
    ax.loglog(theta_arcmin, -w_geom[i], '--r', lw=2)
    ax.set_xlim(5e0, 2e4)
#     ax.set_ylim(5e-11, 2e-5)
    ax.tick_params(axis='y', which='minor', labelcolor='none')
    string = '({0},{1})'.format(a,b)
    ax.text(3e3,5e-6,string)

axes[0, 0].legend(markerfirst=False, loc='lower left')

ax = fig.add_subplot(111, frameon=False)
ax.set_xlabel(r'Angular separation, $\theta$ [arcmin]', size=12)
ax.set_ylabel(r'Angular correlation, $w(\theta)$', size=12)
ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
ax.tick_params(axis='y', pad=15)

fig.tight_layout(pad=0.5)
# fig.savefig('plots/w_kk.pdf', bbox_inches='tight')

plt.show()<jupyter_output><empty_output><jupyter_text>5. Angular power spectra * Convergence<jupyter_code>l_limb, cl_limb = np.transpose([corfu.wtocl(theta, w, lmax=2000) for w in w_limb], (1, 0, 2))
l_geom, cl_geom = np.transpose([corfu.wtocl(theta, w, lmax=2000) for w in w_geom], (1, 0, 2))
fig, axes = plt.subplots(4, 4, figsize=(11, 11), sharex=True, sharey=True)

for ax in axes.ravel():
    ax.axis('off')

for i, (a, b) in enumerate(zip(bin_a, bin_b)):
    ax = axes[a-1, b-1]
    ax.axis('on')
    ax.loglog(l_csis, cl_csis[i], 'b', label='CosmoSIS', lw=1)
    ax.loglog(l_limb[i], cl_limb[i], 'k', label='Limber', lw=1)
    ax.loglog(l_geom[i], cl_geom[i], '--r', label='Exact', lw=2)
    ax.set_xlim(9e-1, 2e3)
    ax.set_ylim(2e-10, 5e-8)
    ax.set_xticks([1e0, 1e1, 1e2, 1e3])
    string = '({0},{1})'.format(a,b)
    ax.text(5e2,1e-8,string)


axes[0, 0].legend(markerfirst=False, loc='lower left')

ax = fig.add_subplot(111, frameon=False)
ax.set_xlabel(r'Angular mode, $\ell$', size=12)
ax.set_ylabel(r'Angular power, $C_{\ell}$', size=12)
ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
ax.tick_params(axis='y', pad=12)

fig.tight_layout(pad=0.5)
# fig.savefig('plots/cl_kk.pdf', bbox_inches='tight')

plt.show()<jupyter_output><empty_output><jupyter_text>* ShearNote: $$C_{\ell}^{\gamma \gamma} = \frac{(\ell - 1)(\ell + 2)}{\ell (\ell + 1)}C_{\ell}^{\kappa \kappa}$$<jupyter_code>def ell_factor(ell):
    return (ell - 1) * (ell + 2) / (ell * (ell + 1))
factor_limb = ell_factor(l_limb[0])
factor_geom = ell_factor(l_geom[0])
factor_csis = ell_factor(l_csis)
fig, axes = plt.subplots(4, 4, figsize=(11, 11), sharex=True, sharey=True)

for ax in axes.ravel():
    ax.axis('off')

for i, (a, b) in enumerate(zip(bin_a, bin_b)):
    ax = axes[a-1, b-1]
    ax.axis('on')
    ax.loglog(l_csis, factor_csis * cl_csis[i], 'b', label='CosmoSIS', lw=1)
    ax.loglog(l_limb[i], factor_limb * cl_limb[i], 'k', label='Limber', lw=1)
    ax.loglog(l_geom[i], factor_geom * cl_geom[i], '--r', label='Exact', lw=2)
    ax.set_xlim(2, 2e3)
    ax.set_ylim(2e-10, 5e-8)
    ax.set_xticks([1e1, 1e2, 1e3])
    string = '({0},{1})'.format(a,b)
    ax.text(5e2,1e-8,string)

axes[0, 0].legend(markerfirst=False, loc='lower left')

ax = fig.add_subplot(111, frameon=False)
ax.set_xlabel(r'Angular mode, $\ell$', size=12)
ax.set_ylabel(r'Angular power, $C_{\ell}$', size=12)
ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
ax.tick_params(axis='y', pad=12)

fig.tight_layout(pad=0.5)
# fig.savefig('plots/cl_shear.pdf', bbox_inches='tight')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20245779543197617, 'fraction_numerical': 0.03289473684210526, 'mean_word_length': 5.430167597765363, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '122360', 'license': 'MIT', 'n_tokens_mistral': 3879, 'n_tokens_neox': 3697, 'n_words': 598, 'path': 'tests/test_cosmic_lensing.ipynb', 'repo_name': 'Lucia-Fonseca/unequalpy', 'size': '8056'}