Skip to content

Commit

Permalink
add quantiles in tests/test_galaxycluster.py
Browse files Browse the repository at this point in the history
  • Loading branch information
@m-aguena
m-aguena committed Sep 27, 2024
1 parent 2d53891 commit 8fc36c8
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 37 deletions.
83 changes: 47 additions & 36 deletions examples/test_pz_integral.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
"from clmm import GalaxyCluster\n",
"from clmm.dataops import compute_galaxy_weights\n",
"from clmm.support import mock_data as mock\n",
"\n",
"clmm.__version__"
]
},
Expand All @@ -50,7 +51,7 @@
"metadata": {},
"outputs": [],
"source": [
"cosmo = Cosmology(H0 = 71.0, Omega_dm0 = 0.265 - 0.0448, Omega_b0 = 0.0448, Omega_k0 = 0.0)"
"cosmo = Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0)"
]
},
{
Expand Down Expand Up @@ -94,20 +95,26 @@
"np.random.seed(41363)\n",
"noisy_data_z = mock.generate_galaxy_catalog(*args, **kwargs, pzpdf_type=\"shared_bins\")\n",
"np.random.seed(41363)\n",
"noisy_data_z2 = mock.generate_galaxy_catalog(*args, **kwargs, pzpdf_type=\"individual_bins\")"
"noisy_data_z2 = mock.generate_galaxy_catalog(*args, **kwargs, pzpdf_type=\"individual_bins\")\n",
"np.random.seed(41363)\n",
"noisy_data_z3 = mock.generate_galaxy_catalog(*args, **kwargs, pzpdf_type=\"quantiles\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "596ee44c",
"id": "1c9abd54",
"metadata": {},
"outputs": [],
"source": [
"for i, data in enumerate(noisy_data_z[:10]):\n",
" plt.plot(noisy_data_z.pzpdf_info['zbins'], data['pzpdf'], lw=.5, color=f\"C{i}\")\n",
"for i, data in enumerate(noisy_data_z2[:10]):\n",
" plt.plot(data['pzbins'], data['pzpdf'], lw=.9, color=f\"C{i}\", ls=\"--\")"
"for i, data in enumerate(noisy_data_z[:5]):\n",
" plt.plot(noisy_data_z.pzpdf_info[\"zbins\"], data[\"pzpdf\"], lw=0.5, color=f\"C{i}\")\n",
"for i, data in enumerate(noisy_data_z2[:5]):\n",
" plt.plot(data[\"pzbins\"], data[\"pzpdf\"], lw=0.9, color=f\"C{i}\", ls=\"--\")\n",
"pzbins, pzpdfs = noisy_data_z3.get_pzpdfs()\n",
"for i, data in enumerate(pzpdfs[:5]):\n",
" plt.plot(pzbins, data, lw=3, color=f\"C{i}\", ls=\":\")\n",
"plt.xlim(0.4, 2.1)"
]
},
{
Expand Down Expand Up @@ -139,13 +146,14 @@
{
"cell_type": "code",
"execution_count": null,
"id": "21f3560f",
"id": "282d6fba",
"metadata": {},
"outputs": [],
"source": [
"integrals = {\n",
" \"clmm_shared\": _integ_pzfuncs(noisy_data_z['pzpdf'], noisy_data_z.pzpdf_info['zbins'], cluster_z),\n",
" \"clmm_individual\": _integ_pzfuncs(noisy_data_z2['pzpdf'], noisy_data_z2['pzbins'], cluster_z),\n",
" \"clmm_shared\": _integ_pzfuncs(*noisy_data_z.get_pzpdfs()[::-1], cluster_z),\n",
" \"clmm_individual\": _integ_pzfuncs(*noisy_data_z2.get_pzpdfs()[::-1], cluster_z),\n",
" \"clmm_quantiles\": _integ_pzfuncs(*noisy_data_z3.get_pzpdfs()[::-1], cluster_z),\n",
"}"
]
},
Expand Down Expand Up @@ -178,7 +186,7 @@
" qp.interp,\n",
" data={\"xvals\": noisy_data_z.pzpdf_info[\"zbins\"], \"yvals\": noisy_data_z[\"pzpdf\"]},\n",
")\n",
"integrals[\"qp_shared\"] = 1-qp_dat.cdf(cluster_z)[:,0]"
"integrals[\"qp_shared\"] = 1 - qp_dat.cdf(cluster_z)[:, 0]"
]
},
{
Expand All @@ -192,20 +200,34 @@
" qp.interp_irregular,\n",
" data={\"xvals\": noisy_data_z2[\"pzbins\"], \"yvals\": noisy_data_z2[\"pzpdf\"]},\n",
")\n",
"integrals[\"qp_individual\"] = 1-qp_dat2.cdf(cluster_z)[:,0]"
"integrals[\"qp_individual\"] = 1 - qp_dat2.cdf(cluster_z)[:, 0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6c6d950",
"metadata": {},
"outputs": [],
"source": [
"qp_dat3 = qp.Ensemble(\n",
" qp.quant,\n",
" data={\"locs\": noisy_data_z3[\"pzquantiles\"], \"quants\": noisy_data_z3.pzpdf_info[\"quantiles\"]},\n",
")\n",
"integrals[\"qp_quantiles\"] = 1 - qp_dat3.cdf(cluster_z)[:, 0]"
]
},
{
"cell_type": "markdown",
"id": "3d644338",
"id": "490d5786",
"metadata": {},
"source": [
"### True Values"
]
},
{
"cell_type": "markdown",
"id": "6b853399",
"id": "d3694d30",
"metadata": {},
"source": [
"For a gaussian distribution, the integral can be computed with the error function:\n",
Expand All @@ -220,7 +242,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "9e3783cc",
"id": "d7ab180f",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -230,13 +252,12 @@
{
"cell_type": "code",
"execution_count": null,
"id": "f02cec07",
"id": "51dfebbe",
"metadata": {},
"outputs": [],
"source": [
"true_integ = 0.5 * erfc(\n",
" (cluster_z - noisy_data_z2[\"z\"])\n",
" / (0.05 * (1 + noisy_data_z2[\"z\"][ind]) * np.sqrt(2))\n",
" (cluster_z - noisy_data_z2[\"z\"]) / (0.05 * (1 + noisy_data_z2[\"z\"]) * np.sqrt(2))\n",
")"
]
},
Expand All @@ -251,7 +272,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "a001a3dd",
"id": "d422f4ca",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -261,39 +282,29 @@
{
"cell_type": "code",
"execution_count": null,
"id": "1e614a9c",
"id": "93f46215",
"metadata": {},
"outputs": [],
"source": [
"fig, axes = plt.subplots(\n",
" 2, 2, sharex=True, sharey=True, figsize=(10, 5)\n",
")\n",
"fig, axes = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(10, 5))\n",
"bins = np.linspace(0, 1, 21)\n",
"for (label, integ), ax in zip(integrals.items(), axes.flatten()):\n",
" ax.errorbar(\n",
" 0.5 * (bins[1:] + bins[:-1]),\n",
" binned_statistic(integ, (integ / true_integ - 1) * 100, bins=bins)[0],\n",
" binned_statistic(integ, (integ / true_integ - 1) * 100, bins=bins, statistic='std')[0],\n",
" binned_statistic(integ, (integ / true_integ - 1) * 100, bins=bins, statistic=\"std\")[0],\n",
" )\n",
" ax.axhline(0, c=\"r\", ls=\"--\")\n",
" ax.minorticks_on()\n",
" ax.grid()\n",
" ax.grid(which='minor', lw=.3)\n",
" ax.set_ylim(-20, 15)\n",
" ax.grid(which=\"minor\", lw=0.3)\n",
" ax.set_ylim(-20, 30)\n",
" ax.set_title(label)\n",
"for ax in axes[1]:\n",
" ax.set_xlabel(\"True integral\")\n",
"for ax in axes[:,0]:\n",
" ax.set_ylabel('rel. diff [%]')"
"for ax in axes[:, 0]:\n",
" ax.set_ylabel(\"rel. diff [%]\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5c781ff5",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down
2 changes: 1 addition & 1 deletion tests/test_galaxycluster.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -337,7 +337,7 @@ def test_integrity_of_probfuncs():
1 + cluster.galcat["z"][:, None]
)
cluster.compute_background_probability(use_pdz=True, p_background_name="p_bkg_pz")
assert_allclose(cluster.galcat["p_bkg_pz"], expected, rtol=5e-4)
assert_allclose(cluster.galcat["p_bkg_pz"], expected, rtol=2e-3)


def test_integrity_of_weightfuncs():
Expand Down

0 comments on commit 8fc36c8

Please sign in to comment.