Skip to content

Non-limber C_ells #42

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
damonge wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Non-limber C_ells #42

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
303 changes: 303 additions & 0 deletions NonLimberCells.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,303 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "ad3754f0-0e21-4641-8f3b-cd79e1ac9f16",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pyccl as ccl\n",
"import pyccl.nl_pt as pt\n",
"\n",
"import time\n",
"import pytest\n",
"\n",
"root = \"benchmarks/data/nonlimber/\"\n",
"\n",
"\n",
"def get_cosmological_parameters():\n",
" return {\n",
" \"Omega_m\": 0.3156,\n",
" \"Omega_b\": 0.0492,\n",
" \"w0\": -1.0,\n",
" \"h\": 0.6727,\n",
" \"A_s\": 2.12107e-9,\n",
" \"n_s\": 0.9645,\n",
" \"Neff\": 3.046,\n",
" \"T_CMB\": 2.725,\n",
" }\n",
"\n",
"\n",
"def get_tracer_parameters():\n",
" # Per-bin galaxy bias\n",
" b_g = np.array(\n",
" [\n",
" 1.376695,\n",
" 1.451179,\n",
" 1.528404,\n",
" 1.607983,\n",
" 1.689579,\n",
" 1.772899,\n",
" 1.857700,\n",
" 1.943754,\n",
" 2.030887,\n",
" 2.118943,\n",
" ]\n",
" )\n",
" return {\"b_g\": b_g}\n",
"\n",
"\n",
"def get_ells():\n",
" return np.unique(np.geomspace(2, 2000, 128).astype(int)).astype(float)\n",
"\n",
"\n",
"def get_nmodes(fsky=0.4):\n",
" \"\"\"Returns the number of modes in each ell bin\"\"\"\n",
" ls = get_ells()\n",
" nmodes = list(ls[1:] ** 2 - ls[:-1] ** 2)\n",
" lp = ls[-1] ** 2 / ls[-2]\n",
" nmodes.append(lp**2 - ls[-1] ** 2)\n",
" return np.array(nmodes) * 0.5 * fsky\n",
"\n",
"\n",
"def get_tracer_dNdzs():\n",
" filename = root + \"/dNdzs_fullwidth.npz\"\n",
" d = np.load(filename)\n",
" return {\n",
" \"z_sh\": d[\"z_sh\"],\n",
" \"dNdz_sh\": d[\"dNdz_sh\"],\n",
" \"z_cl\": d[\"z_cl\"],\n",
" \"dNdz_cl\": d[\"dNdz_cl\"],\n",
" }\n",
"\n",
"\n",
"def read_cls():\n",
" d = np.load(root + \"/benchmarks_nl_full_clgg.npz\")\n",
" ls = d[\"ls\"]\n",
" cls_gg = d[\"cls\"]\n",
" d = np.load(root + \"/benchmarks_nl_full_clgs.npz\")\n",
" cls_gs = d[\"cls\"]\n",
" d = np.load(root + \"/benchmarks_nl_full_clss.npz\")\n",
" cls_ss = d[\"cls\"]\n",
" return ls, cls_gg, cls_gs, cls_ss\n",
"\n",
"\n",
"@pytest.fixture(scope=\"module\")\n",
"def set_up():\n",
" par = get_cosmological_parameters()\n",
" cosmo = ccl.Cosmology(\n",
" Omega_c=par[\"Omega_m\"] - par[\"Omega_b\"],\n",
" Omega_b=par[\"Omega_b\"],\n",
" h=par[\"h\"],\n",
" n_s=par[\"n_s\"],\n",
" A_s=par[\"A_s\"],\n",
" w0=par[\"w0\"],\n",
" )\n",
" tpar = get_tracer_parameters()\n",
" Nzs = get_tracer_dNdzs()\n",
"\n",
" t_g = []\n",
" # we pass unity bias here and will multiply by the actual bias later\n",
" # this allows us to use the same setup\n",
" # for both direct and PTtracer approach\n",
" for Nz in Nzs[\"dNdz_cl\"].T:\n",
" t = ccl.NumberCountsTracer(\n",
" cosmo,\n",
" has_rsd=False,\n",
" dndz=(Nzs[\"z_cl\"], Nz),\n",
" bias=(Nzs[\"z_cl\"], np.ones(len(Nzs[\"z_cl\"]))))\n",
" t_g.append(t)\n",
" t_s = []\n",
" for Nz in Nzs[\"dNdz_sh\"].T:\n",
" t = ccl.WeakLensingTracer(cosmo, dndz=(Nzs[\"z_sh\"], Nz))\n",
" t_s.append(t)\n",
" ells = get_ells()\n",
" raw_truth = read_cls()\n",
" indices_gg = []\n",
" indices_gs = []\n",
" indices_ss = []\n",
" rind_gg = {}\n",
" rind_gs = {}\n",
" rind_ss = {}\n",
" Ng, Ns = len(t_g), len(t_s)\n",
" for i1 in range(Ng):\n",
" for i2 in range(i1, Ng):\n",
" rind_gg[(i1, i2)] = len(indices_gg)\n",
" rind_gg[(i2, i1)] = len(indices_gg)\n",
" indices_gg.append((i1, i2))\n",
"\n",
" for i2 in range(Ns):\n",
" rind_gs[(i1, i2)] = len(indices_gs)\n",
" rind_gs[(i2, i1)] = len(indices_gs)\n",
" indices_gs.append((i1, i2))\n",
"\n",
" for i1 in range(Ns):\n",
" for i2 in range(i1, Ns):\n",
" rind_ss[(i1, i2)] = len(indices_ss)\n",
" rind_ss[(i2, i1)] = len(indices_ss)\n",
" indices_ss.append((i1, i2))\n",
"\n",
" # Sanity checks\n",
" assert np.allclose(raw_truth[0], ells)\n",
" Nell = len(ells)\n",
" tgg, tgs, tss = raw_truth[1:]\n",
" assert tgg.shape == (len(indices_gg), Nell)\n",
" assert tgs.shape == (len(indices_gs), Nell)\n",
" assert tss.shape == (len(indices_ss), Nell)\n",
"\n",
" # now generate errors\n",
" err_gg = []\n",
" err_gs = []\n",
" err_ss = []\n",
" nmodes = get_nmodes()\n",
" for i1, i2 in indices_gg:\n",
" err_gg.append(\n",
" np.sqrt(\n",
" (\n",
" tgg[rind_gg[(i1, i1)]] * tgg[rind_gg[(i2, i2)]]\n",
" + tgg[rind_gg[(i1, i2)]] ** 2\n",
" )\n",
" / nmodes\n",
" )\n",
" )\n",
" for i1, i2 in indices_gs:\n",
" err_gs.append(\n",
" np.sqrt(\n",
" (\n",
" tgg[rind_gg[(i1, i1)]] * tss[rind_ss[(i2, i2)]]\n",
" + tgs[rind_gs[(i1, i2)]] ** 2\n",
" )\n",
" / nmodes\n",
" )\n",
" )\n",
" for i1, i2 in indices_ss:\n",
" err_ss.append(\n",
" np.sqrt(\n",
" (\n",
" tss[rind_ss[(i1, i1)]] * tss[rind_ss[(i2, i2)]]\n",
" + tss[rind_ss[(i1, i2)]] ** 2\n",
" )\n",
" / nmodes\n",
" )\n",
" )\n",
" # Finally get the PT thing going\n",
" ptc = pt.EulerianPTCalculator(with_NC=True, with_IA=False,\n",
" log10k_min=-4, log10k_max=4,\n",
" nk_per_decade=200)\n",
" ptc.update_ingredients(cosmo)\n",
" ptc_lin = pt.EulerianPTCalculator(with_NC=True, with_IA=False,\n",
" log10k_min=-4, log10k_max=4,\n",
" nk_per_decade=200,\n",
" b1_pk_kind='linear',\n",
" bk2_pk_kind='linear')\n",
" ptc_lin.update_ingredients(cosmo)\n",
" ptt_g = [pt.PTNumberCountsTracer(b1=bias) for bias in tpar[\"b_g\"]]\n",
" ptt_m = pt.PTMatterTracer()\n",
"\n",
" tracers1 = {\"gg\": t_g, \"gs\": t_g, \"ss\": t_s}\n",
" tracers2 = {\"gg\": t_g, \"gs\": t_s, \"ss\": t_s}\n",
" truth = {\"gg\": tgg, \"gs\": tgs, \"ss\": tss}\n",
" errors = {\"gg\": err_gg, \"gs\": err_gs, \"ss\": err_ss}\n",
" indices = {\"gg\": indices_gg, \"gs\": indices_gs, \"ss\": indices_ss}\n",
" ptobj = {\"ptc\": ptc, \"ptc_lin\": ptc_lin, \"ptt_g\": ptt_g, \"ptt_m\": ptt_m}\n",
" return cosmo, ells, tracers1, tracers2, truth, errors, indices, ptobj\n",
"\n",
"\n",
"@pytest.mark.parametrize(\"method\", [None, \"FKEM\"])\n",
"@pytest.mark.parametrize(\"cross_type\", [\"gg\", \"gs\", \"ss\"])\n",
"@pytest.mark.parametrize(\"pt_path\", [False, True])\n",
"def test_cells(set_up, method, cross_type, pt_path):\n",
" cosmo, ells, tracers1, tracers2, truth, errors, indices, ptobj = set_up\n",
" t0 = time.time()\n",
" chi2max = 0\n",
" biases = get_tracer_parameters()[\"b_g\"]\n",
" for pair_index, (i1, i2) in enumerate(indices[cross_type]):\n",
" bias_fact = 1.0\n",
" p_of_k_a = p_of_k_a_lin = ccl.DEFAULT_POWER_SPECTRUM\n",
" if cross_type == \"gg\":\n",
" if pt_path:\n",
" p_of_k_a = ptobj['ptc'].get_biased_pk2d(\n",
" ptobj['ptt_g'][i1],\n",
" tracer2=ptobj['ptt_g'][i2])\n",
" p_of_k_a_lin = ptobj['ptc_lin'].get_biased_pk2d(\n",
" ptobj['ptt_g'][i1],\n",
" tracer2=ptobj['ptt_g'][i2])\n",
" else:\n",
" bias_fact = biases[i1] * biases[i2]\n",
" elif cross_type == \"gs\":\n",
" if pt_path:\n",
" p_of_k_a = ptobj['ptc'].get_biased_pk2d(\n",
" ptobj['ptt_g'][i1],\n",
" tracer2=ptobj['ptt_m'])\n",
" p_of_k_a_lin = ptobj['ptc_lin'].get_biased_pk2d(\n",
" ptobj['ptt_g'][i1],\n",
" tracer2=ptobj['ptt_m'])\n",
" else:\n",
" bias_fact = biases[i1]\n",
"\n",
" if method is None:\n",
" cls = ccl.angular_cl(\n",
" cosmo,\n",
" tracers1[cross_type][i1],\n",
" tracers2[cross_type][i2],\n",
" ells,\n",
" l_limber=-1,\n",
" p_of_k_a=p_of_k_a,\n",
" p_of_k_a_lin=p_of_k_a_lin\n",
" )\n",
" l_limber = 0\n",
" else:\n",
" cls, meta = ccl.angular_cl(\n",
" cosmo,\n",
" tracers1[cross_type][i1],\n",
" tracers2[cross_type][i2],\n",
" ells,\n",
" p_of_k_a=p_of_k_a,\n",
" p_of_k_a_lin=p_of_k_a_lin,\n",
" l_limber='auto',\n",
" non_limber_integration_method=method,\n",
" return_meta=True\n",
" )\n",
" l_limber = meta['l_limber']\n",
" cls *= bias_fact\n",
" chi2 = (cls - truth[cross_type][pair_index, :]) ** 2 / errors[\n",
" cross_type\n",
" ][pair_index] ** 2\n",
" chi2max = max(chi2.max(), chi2max)\n",
" if method is not None: # Limber is going to fail by default\n",
" assert np.all(chi2 < 0.3)\n",
" t1 = time.time()\n",
" method_name = \"Limber\"\n",
" if method is not None:\n",
" method_name = method\n",
" print(\n",
" f'Time taken for {method_name} on {cross_type} = {(t1-t0):3.2f};\\\n",
" worst chi2 = {chi2max:5.3f} l_limber = {l_limber}'\n",
" )"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}