Skip to content

Commit

Permalink
Update examples/Paper_v1.0/validation_tests.ipynb
Browse files Browse the repository at this point in the history
  • Loading branch information
@hsinfan1996
hsinfan1996 committed Sep 7, 2023
1 parent 907b1bd commit 2b2946f
Showing 1 changed file with 49 additions and 25 deletions.
74 changes: 49 additions & 25 deletions examples/Paper_v1.0/validation_tests.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,9 @@
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"# Validation tests for the CLMM backends\n",
"\n",
Expand Down Expand Up @@ -161,28 +163,34 @@
"\n",
"# CCL\n",
"MDEF = \"matter\"\n",
"conc = ccl.halos.ConcentrationConstant(cvir)\n",
"mdef = ccl.halos.MassDef(Delta, \"matter\", c_m_relation=conc)\n",
"mdef = ccl.halos.MassDef(Delta, MDEF)\n",
"conc = ccl.halos.ConcentrationConstant(cvir, mass_def=mdef)\n",
"# mdef.concentration = conc\n",
"\n",
"ccl_nfw_num_opt = ccl.halos.HaloProfileNFW(\n",
" conc, truncated=False, projected_analytic=False, cumul2d_analytic=False, fourier_analytic=False\n",
" mass_def=mdef, concentration=conc,\n",
" truncated=False, projected_analytic=False, cumul2d_analytic=False, fourier_analytic=False\n",
")\n",
"ccl_nfw_num = ccl.halos.HaloProfileNFW(\n",
" conc, truncated=False, projected_analytic=False, cumul2d_analytic=False\n",
" mass_def=mdef, concentration=conc,\n",
" truncated=False, projected_analytic=False, cumul2d_analytic=False\n",
")\n",
"ccl_nfw_ana = ccl.halos.HaloProfileNFW(\n",
" conc, truncated=False, projected_analytic=True, cumul2d_analytic=True\n",
" mass_def=mdef, concentration=conc,\n",
" truncated=False, projected_analytic=True, cumul2d_analytic=True\n",
")\n",
"\n",
"# ccl_nfw_num.update_precision_fftlog (n_per_decade = 10000)\n",
"# ccl_nfw_num.update_precision_fftlog (plaw_fourier = -2)\n",
"\n",
"ccl_ein = ccl.halos.HaloProfileEinasto(conc, truncated=False)\n",
"ccl_her = ccl.halos.HaloProfileHernquist(conc, truncated=False)\n",
"ccl_ein = ccl.halos.HaloProfileEinasto(mass_def=mdef, concentration=conc,\n",
" truncated=False)\n",
"ccl_ein_quad = ccl.halos.HaloProfileEinasto(mass_def=mdef, concentration=conc,\n",
" truncated=False, projected_quad=True)\n",
"ccl_her = ccl.halos.HaloProfileHernquist(mass_def=mdef, concentration=conc, truncated=False)\n",
"\n",
"\n",
"alpha = ccl_ein._get_alpha(cosmo_ccl, Mvir, a, mdef)\n",
"alpha = ccl_ein._get_alpha(cosmo_ccl, Mvir, a)\n",
"\n",
"# Colossus\n",
"col_nfw = profile_nfw.NFWProfile(M=(Mvir * cosmo_col.h), c=cvir, z=z, mdef=\"200m\")\n",
Expand Down Expand Up @@ -277,26 +285,27 @@
"nc_Sigma_her = smd.sigma_array(nc_her, cosmo, r, 1.0, 1.0, z)\n",
"\n",
"# CCL\n",
"ccl_Sigma_nfw_ana = ccl_nfw_ana.projected(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_Sigma_nfw_num = ccl_nfw_num.projected(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_Sigma_ein = ccl_ein.projected(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_Sigma_her = ccl_her.projected(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_Sigma_nfw_ana = ccl_nfw_ana.projected(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_Sigma_nfw_num = ccl_nfw_num.projected(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_Sigma_ein = ccl_ein.projected(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_Sigma_ein_quad = ccl_ein.projected(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_Sigma_her = ccl_her.projected(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"\n",
"\n",
"# CCL numerical NFW, using optimised setup\n",
"# When using fourier_analytic=False in CCL profile definition, CCL performs\n",
"# better by first evaluating the profile on a wider range and then\n",
"# interpolating to the desired radii\n",
"rtmp = np.geomspace(1.0e-4, 100, 1000)\n",
"tmp = ccl_nfw_num_opt.projected(cosmo_ccl, rtmp / a, Mvir, a, mdef) / a**2\n",
"tmp = ccl_nfw_num_opt.projected(cosmo_ccl, rtmp / a, Mvir, a) / a**2\n",
"ptf = interp1d(np.log(rtmp), np.log(tmp), bounds_error=False, fill_value=-100)\n",
"ccl_Sigma_nfw_num_opt = np.exp(ptf(np.log(r)))\n",
"\n",
"tmp = ccl_ein.projected(cosmo_ccl, rtmp / a, Mvir, a, mdef) / a**2\n",
"tmp = ccl_ein.projected(cosmo_ccl, rtmp / a, Mvir, a) / a**2\n",
"ptf = interp1d(np.log(rtmp), np.log(tmp), bounds_error=False, fill_value=-100)\n",
"ccl_Sigma_ein_opt = np.exp(ptf(np.log(r)))\n",
"\n",
"tmp = ccl_her.projected(cosmo_ccl, rtmp / a, Mvir, a, mdef) / a**2\n",
"tmp = ccl_her.projected(cosmo_ccl, rtmp / a, Mvir, a) / a**2\n",
"ptf = interp1d(np.log(rtmp), np.log(tmp), bounds_error=False, fill_value=-100)\n",
"ccl_Sigma_her_opt = np.exp(ptf(np.log(r)))\n",
"\n",
Expand Down Expand Up @@ -335,15 +344,15 @@
"nc_DeltaSigma_her = np.array(smd.sigma_excess_array(nc_her, cosmo, r, 1.0, 1.0, z))\n",
"\n",
"# CCL\n",
"ccl_BarSigma_nfw_ana = ccl_nfw_ana.cumul2d(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_BarSigma_nfw_ana = ccl_nfw_ana.cumul2d(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_DeltaSigma_nfw_ana = ccl_BarSigma_nfw_ana - ccl_Sigma_nfw_ana\n",
"\n",
"# CCL numerical NFW, using default setup\n",
"ccl_BarSigma_nfw_num = ccl_nfw_num.cumul2d(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_BarSigma_nfw_num = ccl_nfw_num.cumul2d(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_DeltaSigma_nfw_num = ccl_BarSigma_nfw_num - ccl_Sigma_nfw_num\n",
"ccl_BarSigma_ein = ccl_ein.cumul2d(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_BarSigma_ein = ccl_ein.cumul2d(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_DeltaSigma_ein = ccl_BarSigma_ein - ccl_Sigma_ein\n",
"ccl_BarSigma_her = ccl_her.cumul2d(cosmo_ccl, r / a, Mvir, a, mdef) / a**2\n",
"ccl_BarSigma_her = ccl_her.cumul2d(cosmo_ccl, r / a, Mvir, a) / a**2\n",
"ccl_DeltaSigma_her = ccl_BarSigma_her - ccl_Sigma_her\n",
"\n",
"\n",
Expand All @@ -353,17 +362,17 @@
"# interpolating to the desired radii\n",
"\n",
"rtmp = np.geomspace(1.0e-4, 100, 1000) # extended radial range\n",
"tmp = ccl_nfw_num_opt.cumul2d(cosmo_ccl, rtmp / a, Mvir, a, mdef) / a**2 # CCL estimation\n",
"tmp = ccl_nfw_num_opt.cumul2d(cosmo_ccl, rtmp / a, Mvir, a) / a**2 # CCL estimation\n",
"ptf = interp1d(np.log(rtmp), np.log(tmp), bounds_error=False, fill_value=-100) # interpolation\n",
"ccl_BarSigma_nfw_num_opt = np.exp(ptf(np.log(r))) # evaluation on the desired radius array\n",
"ccl_DeltaSigma_nfw_num_opt = ccl_BarSigma_nfw_num_opt - ccl_Sigma_nfw_num_opt\n",
"\n",
"tmp = ccl_ein.cumul2d(cosmo_ccl, rtmp / a, Mvir, a, mdef) / a**2\n",
"tmp = ccl_ein.cumul2d(cosmo_ccl, rtmp / a, Mvir, a) / a**2\n",
"ptf = interp1d(np.log(rtmp), np.log(tmp), bounds_error=False, fill_value=-100)\n",
"ccl_BarSigma_ein_opt = np.exp(ptf(np.log(r)))\n",
"ccl_DeltaSigma_ein_opt = ccl_BarSigma_ein_opt - ccl_Sigma_ein_opt\n",
"\n",
"tmp = ccl_her.cumul2d(cosmo_ccl, rtmp / a, Mvir, a, mdef) / a**2\n",
"tmp = ccl_her.cumul2d(cosmo_ccl, rtmp / a, Mvir, a) / a**2\n",
"ptf = interp1d(np.log(rtmp), np.log(tmp), bounds_error=False, fill_value=-100)\n",
"ccl_BarSigma_her_opt = np.exp(ptf(np.log(r)))\n",
"ccl_DeltaSigma_her_opt = ccl_BarSigma_her_opt - ccl_Sigma_her_opt\n",
Expand Down Expand Up @@ -662,6 +671,14 @@
")\n",
"axs[0].plot(\n",
" r,\n",
" np.abs(ccl_Sigma_ein_quad / nc_Sigma_ein - 1.0),\n",
" label=\"CCL - EIN (quad)\",\n",
" linestyle=\":\",\n",
" color=\"darkorange\",\n",
" lw=0.5,\n",
")\n",
"axs[0].plot(\n",
" r,\n",
" np.abs(ccl_Sigma_ein_opt / nc_Sigma_ein - 1.0),\n",
" label=\"CCL - EIN (opt)\",\n",
" linestyle=\"--\",\n",
Expand Down Expand Up @@ -857,6 +874,13 @@
")\n",
"axs[1].plot(\n",
" r,\n",
" np.abs(ccl_Sigma_ein_quad / col_Sigma_ein - 1.0),\n",
" label=\"EIN - CCL (quad)\",\n",
" color=\"cadetblue\",\n",
" linestyle=\"--\",\n",
")\n",
"axs[1].plot(\n",
" r,\n",
" np.abs(ct_Sigma_ein / col_Sigma_ein - 1.0),\n",
" label=\"EIN - CT\",\n",
" color=\"cadetblue\",\n",
Expand Down Expand Up @@ -991,7 +1015,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -1005,7 +1029,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.8"
"version": "3.11.5"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 2b2946f

Please sign in to comment.