diff --git a/README.md b/README.md index 677d7f5..f432622 100644 --- a/README.md +++ b/README.md @@ -30,7 +30,7 @@ $ conda activate shear ## Install the `shrbk` library `shrbk` is the library made for this book. -Install it ither in development mode +Install it either in development mode ```bash $ python -m pip install -e . diff --git a/environment.yml b/environment.yml index 3d0e635..ca98dfc 100644 --- a/environment.yml +++ b/environment.yml @@ -10,3 +10,4 @@ dependencies: - galsim - pip: - jupyter-book + - shear_bias diff --git a/shearbook/_bibliography/z_bias.bib b/shearbook/_bibliography/z_bias.bib index 531eb20..faf4631 100644 --- a/shearbook/_bibliography/z_bias.bib +++ b/shearbook/_bibliography/z_bias.bib @@ -36,6 +36,50 @@ @ARTICLE{2012MNRAS.424.2757M adsnote = {Provided by the SAO/NASA Astrophysics Data System} } +@ARTICLE{PKSB17, + author = {{Pujol}, A. and {Kilbinger}, M. and {Sureau}, F. and {Bobin}, J. + }, + title = "{A highly precise shape-noise-free shear bias estimator}", + journal = {\aap}, +archivePrefix = "arXiv", + volume = {621}, + pages = {A21}, + keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics}, + year = 2019, + month = jan, + adsurl = {http://cdsads.u-strasbg.fr/abs/2018arXiv180610537P}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@ARTICLE{2017arXiv170202600H, + author = {{Huff}, E. and {Mandelbaum}, R.}, + title = "{Metacalibration: Direct Self-Calibration of Biases in Shear Measurement}", + journal = {arXiv}, +archivePrefix = "arXiv", + volume = {1702.02600}, + keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics}, + year = 2017, + month = feb, + adsurl = {http://adsabs.harvard.edu/abs/2017arXiv170202600H}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@ARTICLE{SH17, + author = {{Sheldon}, E.~S. and {Huff}, E.~M.}, + title = "{Practical Weak-lensing Shear Measurement with Metacalibration}", + journal = {\apj}, +archivePrefix = "arXiv", + eprint = {1702.02601}, + keywords = {cosmology: observations, gravitational lensing: weak, methods: observational}, + year = 2017, + month = may, + volume = 841, + eid = {24}, + pages = {24}, + doi = {10.3847/1538-4357/aa704b}, + adsurl = {http://adsabs.harvard.edu/abs/2017ApJ...841...24S}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} @ARTICLE{Bosch2018, author = {{Bosch}, James and {Armstrong}, Robert and {Bickerton}, Steven and {Furusawa}, Hisanori and {Ikeda}, Hiroyuki and {Koike}, Michitaro and {Lupton}, Robert and {Mineo}, Sogo and {Price}, Paul and {Takata}, Tadafumi and {Tanaka}, Masayuki and {Yasuda}, Naoki and {AlSayyad}, Yusra and {Becker}, Andrew C. and {Coulton}, William and {Coupon}, Jean and {Garmilla}, Jose and {Huang}, Song and {Krughoff}, K. Simon and {Lang}, Dustin and {Leauthaud}, Alexie and {Lim}, Kian-Tat and {Lust}, Nate B. and {MacArthur}, Lauren A. and {Mandelbaum}, Rachel and {Miyatake}, Hironao and {Miyazaki}, Satoshi and {Murata}, Ryoma and {More}, Surhud and {Okura}, Yuki and {Owen}, Russell and {Swinbank}, John D. and {Strauss}, Michael A. and {Yamada}, Yoshihiko and {Yamanoi}, Hitomi}, title = "{The Hyper Suprime-Cam software pipeline}", @@ -151,3 +195,21 @@ @article{Barbary2016 title = {SEP: Source Extractor as a library}, journal = {Journal of Open Source Software} } + +@ARTICLE{pujol_shear_bias_20, + author = {{Pujol}, Arnau and {Bobin}, Jerome and {Sureau}, Florent and {Guinot}, Axel and {Kilbinger}, Martin}, + title = "{Shear measurement bias. II. A fast machine-learning calibration method}", + journal = {\aap}, + keywords = {gravitational lensing: weak, methods: numerical, methods: data analysis, methods: observational, methods: statistical, cosmology: observations, Astrophysics - Cosmology and Nongalactic Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics}, + year = 2020, + month = nov, + volume = {643}, + eid = {A158}, + pages = {A158}, + doi = {10.1051/0004-6361/202038658}, +archivePrefix = {arXiv}, + eprint = {2006.07011}, + primaryClass = {astro-ph.CO}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2020A&A...643A.158P}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} diff --git a/shearbook/_toc.yml b/shearbook/_toc.yml index 9f9829c..0c89ec4 100644 --- a/shearbook/_toc.yml +++ b/shearbook/_toc.yml @@ -43,8 +43,10 @@ chapters: - file: bias/bias-intro sections: + - file: bias/bias-definition - file: bias/model-bias - file: bias/noise-bias + - file: bias/measuring-bias - file: bias/blending-bias - file: bias/bias-ref - file: calibration/calibration-intro diff --git a/shearbook/bias/bias-definition.ipynb b/shearbook/bias/bias-definition.ipynb new file mode 100644 index 0000000..968e36f --- /dev/null +++ b/shearbook/bias/bias-definition.ipynb @@ -0,0 +1,76 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " (sec:bias-definition)=\n", + "# Shear bias definition\n", + "\n", + "In most cases, shear bias can be decomposed in a multiplicative and an additive part.\n", + "\n", + "## Population shear bias\n", + "\n", + "First, we define the shear bias for a population of galaxies. Let us recall the equation that,\n", + "in the weak-lensing regime, relates observed and intrincic ellipticity to reduced shear,\n", + "\n", + "$$\\varepsilon^{\\textrm{obs}}_\\alpha = \\varepsilon^{\\textrm{I}}_\\alpha + g_\\alpha,$$\n", + "\n", + "for components $\\alpha = 1, 2$. If the mean intrinsic ellipticity is zero, $\\langle \\varepsilon^{\\textrm{I}}_\\alpha \\rangle = 0$, the observed ellipticity is an unbiased estimator of the shear, $\\langle \\varepsilon^{\\textrm{obs}}_\\alpha \\rangle = g_\\alpha$ in this basic case.\n", + "\n", + "Next, we introduce the ensemble-average additive bias $c_\\alpha$, and multiplicative bias $m_\\alpha$, and write\n", + "\n", + "$$\n", + "\\langle \\varepsilon^{\\textrm{obs}}_\\alpha \\rangle \\equiv g^{\\textrm{obs}} = c_\\alpha + (1 + m_\\alpha) g_\\alpha.\n", + "$$ (shear-estim-biased)\n", + "\n", + "A common calibration method is to estimate the population biases $c_\\alpha$, $m_\\alpha$, by simulating a large number of galaxy images with different properties and shear. These measured values are then applied to the observed galaxy ellipticities. The quantity $(\\varepsilon^{\\textrm{obs}}_\\alpha - c_\\alpha) / (1 + m_\\alpha)$ is then an unbiased estimator of the reduced shear.\n", + "\n", + "## Individual galaxy shear bias\n", + "\n", + "We now define the shear bias for an individual galaxy, and show that the ensemble average of this individual bias\n", + "can be identified with the population shear bias defined above.\n", + "\n", + "The multiplicative bias in {eq}`shear-estim-biased` can be interpreted as derivative, or response, of the average observed ellipticity with respect to the shear. Instead of using the average ellipticity however, we instead compute this derivative of the ellipticity of an individual galaxy. In addition, we generalise this derivative and \n", + "account for the two-component nature of ellipticity. This defines the response matrix $R$ ({cite}`2017arXiv170202600H`, {cite}`SH17`) by\n", + "\n", + "$$\n", + "R_{\\alpha\\beta} = \\frac{\\partial \\varepsilon^{\\textrm{obs}}_\\alpha}{\\partial g_\\beta} .\n", + "$$ (R)\n", + "\n", + "The additive shear bias for an individual galaxy can be defined as well,\n", + "\n", + "$$a_\\alpha = \\varepsilon^{\\textrm{obs}}_\\alpha - R_{\\alpha\\alpha} g_\\alpha - \\varepsilon^{\\textrm{I}}_\\alpha .$$\n", + "\n", + "This quantity can be computed in the case of image simulations where the intrinsic ellipticity $\varepsilon^{\\textrm{I}}$ is known.\n", + "\n", + "By applying the ensemble average to {eq}`R` and comparing to {eq}`shear-estim-biased`, we find the equalities\n", + "$\\langle R_{\\alpha\\alpha} \\rangle = m_\\alpha$, and $\\langle a_\\alpha \\rangle = c_\\alpha$.\n", + "\n", + "With calibration schemes where individual shear bias are measured, for example using metacalibration ({cite}`2017arXiv170202600H`), or deep learning ({cite}`pujol_shear_bias_20`), one applies the ensemble-averaged response matrix and additive shear vector o the observed galaxy ellipticities as\n", + "$\\langle R \\rangle^{-1} \\left( \\vec \\varepsilon^{\\textrm{obs}} - \\langle \\vec a \\rangle \\right)$." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/shearbook/bias/measuring-bias.ipynb b/shearbook/bias/measuring-bias.ipynb new file mode 100644 index 0000000..9085705 --- /dev/null +++ b/shearbook/bias/measuring-bias.ipynb @@ -0,0 +1,479 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{code-block}\n", + ":class: thebe, thebe-init\n", + "# Automatic import for live code\n", + "import numpy as np\n", + "from shrbk.plot import *\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " (sec:measuring-bias)=\n", + "# Measuring Bias\n", + "\n", + "In this section we will discuss several methods to measure shear bias.\n", + "\n", + "\n", + "## Shear bias from individual galaxies\n", + "\n", + "{cite}`PKSB17` present a shape-noise-free shear bias estimator from individual simulated galaxies. The principle is similar to the metacalibration technique: A small shear is applied to a galaxy image, and the response matrix of the measured ellipticity to the shear is determined by finite differences.\n", + "\n", + "In general, the shear estimate for a single galaxy is dominated by shape noise. To beat down this noise, traditional methods use a very large number of galaxies. Some additional improvement can be gained by using rotated versions of the same galaxy, such that the sum of the intrinsic ellipticities is zero. The resulting shear estimate has a reduced variance.\n", + "\n", + "However, this does not account for the variance from the measured ellipticity, which does not only depend on the intrinsic one, but on the PSF, pixel noise, etc. In addition, all versions of the rotated galaxy need to be detected and have an measured shape, and therefore, selection biases are difficult to quantify with this method.\n", + "\n", + "To reduce shape noise (to virtually zero), {cite}`PKSB17` uses the same noise realisation is used for each sheared, and the unsheared image. The resulting estimate of the response per galaxy is highly precise. It is still not very accurate, since the measured value depends on the noise realisation. To increase the accuracy, a large number of simulated galaxies with different noise realisations is required. These galaxies are however anyway simulated to cover the large parameter space of galaxy properties.\n", + "\n", + "## Example\n", + "Here we compute the individual shear bias from image simulations. This method is implemented in the library `shear_bias`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the libraries\n", + "import shear_bias as sb\n", + "import galsim\n", + "import matplotlib.pylab as plt\n", + "import numpy as np\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Set up" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Shear values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, set the value of the small shear $\\Delta g$ to add to the galaxy images." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dg = 0.15" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the five steps, one in each positive and negative direction along the two coordinate axis, plus the original unchanged image (0, 0)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "g_steps = [(-1, 0), (0, -1), (1, 0), (0, 1), (0, 0)]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create dictionary of shear values with step tuples as keys." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "g_values = {}\n", + "for step in g_steps:\n", + " g_values[step] = (step[0] * dg, step[1] * dg)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Number of galaxy images" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_gal = 100" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Downloading HST COSMOS galaxy template images" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_dir = 'data'\n", + "cat_file_name = 'real_galaxy_catalog_23.5_example.fits'\n", + "\n", + "sky_level = 1e3 # ADU / arcsec^2\n", + "gal_flux = 1e5 # arbitrary, choose large value for not too noisy images\n", + "\n", + "pixel_scale = 0.16 # arcsec\n", + "random_seed = 5693562491" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sb.download_HST_images(dest_dir=data_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "real_galaxy_catalog = galsim.RealGalaxyCatalog(cat_file_name, dir=data_dir)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create the PSF image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fwhm_psf = [0.6, 2.3]\n", + "flux_frac_psf = [0.8, 0.2]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create PSF\n", + "psf = galsim.Gaussian(fwhm=fwhm_psf[0], flux=flux_frac_psf[0])\n", + "for i in range(1, len(fwhm_psf)):\n", + " psf = psf + galsim.Gaussian(fwhm=fwhm_psf[i], flux=flux_frac_psf[i])\n", + "\n", + "# Draw PSF\n", + "psf_image = psf.drawImage(scale=pixel_scale)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Galaxies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Lensing magnification. Set to 1 for no magnification.\n", + "magnification = 1\n", + "\n", + "# Constant background\n", + "background = sky_level * pixel_scale**2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Create (sheared) galaxy images" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "all_images = {}\n", + "im_test = {}\n", + "k_test = 3\n", + "\n", + "for step in g_steps:\n", + " all_images[step] = []\n", + "\n", + "for k in range(n_gal):\n", + "\n", + " print(k, end=' ')\n", + " rng = galsim.UniformDeviate(random_seed+k+1)\n", + " gal = galsim.RealGalaxy(real_galaxy_catalog, index=k, flux=gal_flux)\n", + " \n", + " for step in g_steps:\n", + " gal_sh = gal.shear(g1=g_values[step][0], g2=g_values[step][1])\n", + "\n", + " if magnification != 1:\n", + " gal_sh = gal_sh.magnify(magnification)\n", + " \n", + " final = galsim.Convolve([psf, gal_sh])\n", + " \n", + " dx = rng() - 0.5\n", + " dy = rng() - 0.5\n", + "\n", + " if k == 0:\n", + " im = final.drawImage(scale=pixel_scale, offset=(dx, dy))\n", + " xsize, ysize = im.array.shape\n", + " else:\n", + " im = galsim.ImageF(xsize, ysize)\n", + " final.drawImage(im, scale=pixel_scale, offset=(dx, dy))\n", + " \n", + " if k == k_test and step == (0, 0):\n", + " im_test['no_noise'] = np.copy(im.array)\n", + "\n", + " im.addNoise(galsim.PoissonNoise(rng=rng, sky_level=sky_level))\n", + "\n", + " #ccd_noise = galsim.CCDNoise(rng, sky_level=0., gain=1., read_noise=0.)\n", + " #im.addNoise(ccd_noise)\n", + " if k == k_test:\n", + " if step == (0, 0):\n", + " im_test['no_shear'] = np.copy(im.array)\n", + " if step == (1, 0):\n", + " im_test['+g1'] = np.copy(im.array)\n", + " if step == (0, 1):\n", + " im_test['+g2'] = np.copy(im.array)\n", + "\n", + " all_images[step].append(im)\n", + "print()\n", + "\n", + "im_test['+g1-(-g1)'] = all_images[(1,0)][k_test].array - all_images[(-1, 0)][k_test].array\n", + "im_test['+g2-0'] = all_images[(0,1)][k_test].array - all_images[(0, 0)][k_test].array" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plot some galaxies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nim = len(im_test)\n", + "fig, (axes) = plt.subplots(nrows=1, ncols=nim, figsize=(20, 10), tight_layout=True)\n", + "\n", + "for ax, key in zip(axes, im_test):\n", + " ax.imshow(im_test[key])\n", + " ax.set_title(key)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Shape measurement" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "e1 = {}\n", + "e2 = {}\n", + "scale = {}\n", + "sn = {}\n", + "beta = {}\n", + "q = {}\n", + "ep = {}\n", + "ex = {}\n", + "\n", + "# Initialise dictionaries as arrays\n", + "for x in [e1, e2, scale, sn, beta, q, ep, ex]:\n", + " for step in g_steps:\n", + " x[step] = []\n", + "\n", + "for k in range(n_gal):\n", + "\n", + " print(k, end=' ')\n", + " \n", + " for step in g_steps:\n", + " result = galsim.hsm.EstimateShear(all_images[step][k], psf_image, shear_est=\"KSB\", strict=False)\n", + " e1[step].append(result.corrected_g1)\n", + " e2[step].append(result.corrected_g2)\n", + " scale[step].append(result.moments_sigma)\n", + " \n", + " # TODO: compute the following parameters\n", + " sn[step].append(0)\n", + " beta[step].append(0)\n", + " q[step].append(0)\n", + " ep[step].append(0)\n", + " ex[step].append(0)\n", + "\n", + "results = {}\n", + "for step in g_steps:\n", + " results[step] = sb.gal_par.from_values(e1[step], e2[step], scale[step], sn[step],\n", + " beta[step], q[step], ep[step], ex[step])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = np.zeros(2)\n", + "s = np.zeros(2)\n", + "print('{:6s}: {:5s} +- {:3s}'.format('dg', 'mean', 'std'))\n", + "for step in g_steps:\n", + " m[0], s[0] = np.mean(e1[step]), np.std(e1[step])\n", + " m[1], s[1] = np.mean(e2[step]), np.std(e2[step])\n", + " for i in (0, 1):\n", + " print('{: 6.3g}: {: .3f} +- {: .3f}'.format(dg * step[i], m[i], s[i]))\n", + " print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Shear response" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = sb.shear_response(results, dg, output_dir='.')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = sb.read_R('.')\n", + "\n", + "print('mean(R)')\n", + "print(R.mean(axis=2))\n", + "\n", + "print('std(R)')\n", + "print(R.std(axis=2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plot shear response" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Variable on x-axis\n", + "xvar = results[(0, 0)].scale\n", + "\n", + "# Variable name = xlabel\n", + "xname = 'scale'\n", + "\n", + "# Number of x-bins\n", + "nbins = 3\n", + "\n", + "# Variables on y-axis, can be array for multiple curves in plot\n", + "yvar = [sb.shear_bias_m(R, i) for i in [0, 1]]\n", + "yvar.append(R[0,1])\n", + "\n", + "# Variable names for legend\n", + "yname = ['$m_1$', '$m_2$', '$R_{12}$']\n", + "\n", + "# Color of points and error bars\n", + "color = ['g', 'r', 'b']\n", + "\n", + "# Point types\n", + "marker = ['o', 's', 'd']\n", + "\n", + "# Create plot and save to file\n", + "x_mean, y_mean, y_std = sb.plot_mean_per_bin_several(xvar, xname, yvar, yname, nbins, error_mode='std',\n", + " color=color, lw=2, marker=marker,\n", + " out_name='{}_m1.pdf'.format(xname))\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}