Weak Lensing Systematics

environment.yaml

name: weak-lensing-systematics
channels:
  - conda-forge
  - defaults
dependencies:
  - galsim
  - jupyter
  - matplotlib
  - numpy

Weak Lensing Systematics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "78f5a606",
   "metadata": {},
   "source": [
    "# Weak Lensing Systematics: the Art of Precisely Measuring Fuzzy Blobs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab59e6b7",
   "metadata": {},
   "source": [
    "This notebook will introduce the basic concepts necessary for understanding the systematic uncertainties affecting weak lensing analyses.\n",
    "We will make a number of oversimplifications, each of which becomes a critical aspect for doing weak lensing at the accuracy and precision required for a project like LSST.\n",
    "Additionally, this notebook will aim to introduce you to practices and software commonly used in the field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4ec53c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import galsim\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4b691e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will override some Matplotlib parameters for better plots\n",
    "import matplotlib as mpl\n",
    "\n",
    "mpl.rcParams[\"font.family\"] = \"serif\"\n",
    "mpl.rcParams[\"mathtext.fontset\"] = \"dejavuserif\"\n",
    "\n",
    "mpl.rcParams[\"figure.autolayout\"] = \"True\"\n",
    "# mpl.rcParams[\"figure.dpi\"] = 200  # to increase the size of the plots displayed in notebooks (default: 100)\n",
    "\n",
    "mpl.rcParams[\"image.origin\"] = \"lower\"\n",
    "mpl.rcParams[\"image.cmap\"] = \"gray_r\"\n",
    "\n",
    "mpl.rcParams[\"xtick.minor.visible\"] = True\n",
    "mpl.rcParams[\"ytick.minor.visible\"] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcfccffd",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5ebb869",
   "metadata": {},
   "source": [
    "## Stars are Points and Galaxies are Ellipses"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b56f927",
   "metadata": {},
   "source": [
    "The first two assumptions we make are as follows:\n",
    "- stars are small and far away; therefore, we can model them as points (Dirac delta functions in spatial coordinates)\n",
    "- the morphology of galaxies is sufficiently simple that, by the time that a given galaxy is imaged by our cameras, every galaxy will look like an ellipse with some size, axis ratio, and rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0308026",
   "metadata": {},
   "source": [
    "Mathematically, we can write the surface brightness profiles of each type of object:\n",
    "- a star centered at $\\vec{\\mu}$: $$I(\\vec{x}) = \\delta(\\vec{x} - \\vec{\\mu})$$\n",
    "- a galaxy with a Gaussian profile centered at $\\vec{\\mu}$ with covariance matrix $\\Sigma$: $$I(\\vec{x}) = \\frac{1}{2 \\pi |\\Sigma|} \\exp \\left( -\\frac{1}{2} (\\vec{x} - \\vec{\\mu})^\\top \\Sigma^{-1} (\\vec{x} - \\vec{\\mu}) \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7093e72c",
   "metadata": {},
   "source": [
    "Let's use GalSim to draw an example galaxy image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de45d953",
   "metadata": {},
   "outputs": [],
   "source": [
    "galaxy = galsim.Exponential(half_light_radius=0.7)\n",
    "\n",
    "image = galsim.Image(64, 64, scale=0.2)\n",
    "\n",
    "galaxy.drawImage(image=image)\n",
    "\n",
    "plt.imshow(image.array)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c744098",
   "metadata": {},
   "source": [
    "A few comments here:\n",
    "- the $x$ and $y$ axes are in units of pixels\n",
    "- each pixel corresponds to $0.2\\ \\rm{arcsec}$ in either dimension (this is the scale of the LSST camera)\n",
    "- the flux of the galaxy is integrated over each pixel (rather than directly sampling the profile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c24a2f0",
   "metadata": {},
   "source": [
    "We can make elliptical galaxies with GalSim by applying a \"shear\" to the galaxy object.\n",
    "This shear is a linear transformation and can be paramaterized by the distortion $(e_1, e_2)$.\n",
    "See https://galsim-developers.github.io/GalSim/_build/html/shear.html#galsim.Shear for more information about shear definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0483ffd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = galsim.Image(64, 64, scale=0.2)\n",
    "\n",
    "fig, axs = plt.subplots(3, 3)\n",
    "\n",
    "for j, e2 in enumerate([0.5, 0.0, -0.5]):\n",
    "    for i, e1 in enumerate([-0.5, 0.0, 0.5]):\n",
    "        galaxy = galsim.Exponential(half_light_radius=0.7)\n",
    "        ellipticity = galsim.Shear(e1=e1, e2=e2)\n",
    "        galaxy = galaxy.shear(ellipticity)\n",
    "        galaxy.drawImage(image=image)\n",
    "        axs[j, i].imshow(image.array)\n",
    "        axs[j, i].set_title(f\"$(e_1, e_2) = ({e1}, {e2})$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98429c4f",
   "metadata": {},
   "source": [
    "### Things to try\n",
    "\n",
    "The above uses a simple exponential galaxy as an example.\n",
    "Take a look at some of the [simple profiles](https://galsim-developers.github.io/GalSim/_build/html/simple.html) in the GalSim documentation; then, make some plots to explore the differences.\n",
    "The original [GalSim paper](https://arxiv.org/abs/1407.7676) is a great reference for these."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "960a1cc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# placeholder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77f8dd00",
   "metadata": {},
   "source": [
    "## The Atmosphere is Blurry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9262161e",
   "metadata": {},
   "source": [
    "The images shown above represent an ideal measurement of these galaxies as if from space.\n",
    "In practice, the light coming from space will be blurred by the atmosphere before reaching our cameras.\n",
    "We refer to this blurring as the \"point-spread function\" (PSF) because the light from a point is spread around in real space; note that other names for this process are Green's functions, impulse response functions, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18c82435",
   "metadata": {},
   "source": [
    "Mathematically, we model PSFs as real-space convolutions.\n",
    "For some source surface brightness profile $I^0$ and some PSF $P$, the observed surface brightness profile is given by $$I = I^0 \\otimes P$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f65f83",
   "metadata": {},
   "source": [
    "Atmospheric PSFs are often represented by something like a Gaussian profile (i.e., each point gets spread out into a Gaussian) and may have some ellipticity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c41411c",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf = galsim.Gaussian(half_light_radius=0.9)\n",
    "\n",
    "image = galsim.Image(64, 64, scale=0.2)\n",
    "\n",
    "psf.drawImage(image=image)\n",
    "\n",
    "plt.imshow(image.array)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38ed4b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "star = galsim.DeltaFunction()\n",
    "psf = galsim.Gaussian(half_light_radius=0.9)\n",
    "\n",
    "image = galsim.Image(64, 64, scale=0.2)\n",
    "\n",
    "observed_star = galsim.Convolve([star, psf])\n",
    "\n",
    "observed_star.drawImage(image=image)\n",
    "\n",
    "plt.imshow(image.array)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d61d918",
   "metadata": {},
   "source": [
    "Note that the image of the star is the same as the image of the PSF!\n",
    "This follows from our assumption that a star is a delta function in space: convolving a profile with a delta function will yield the profile at the position of the delta function.\n",
    "This property has a very important consequence for weak lensing science: we can measure the PSF of an image by measuring the profiles of the stars in that image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2462434f",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf = galsim.Gaussian(half_light_radius=0.9)\n",
    "\n",
    "image = galsim.Image(64, 64, scale=0.2)\n",
    "\n",
    "fig, axs = plt.subplots(3, 3)\n",
    "\n",
    "for j, e2 in enumerate([0.5, 0.0, -0.5]):\n",
    "    for i, e1 in enumerate([-0.5, 0.0, 0.5]):\n",
    "        galaxy = galsim.Exponential(half_light_radius=0.7)\n",
    "        ellipticity = galsim.Shear(e1=e1, e2=e2)\n",
    "        galaxy = galaxy.shear(ellipticity)\n",
    "        observed_galaxy = galsim.Convolve([galaxy, psf])\n",
    "        observed_galaxy.drawImage(image=image)\n",
    "        axs[j, i].imshow(image.array)\n",
    "        axs[j, i].set_title(f\"$(e_1, e_2) = ({e1}, {e2})$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2dc0d2c",
   "metadata": {},
   "source": [
    "### Things to try\n",
    "\n",
    "The above uses a simple Gaussian PSF as an example.\n",
    "Take a look at some of the [simple profiles](https://galsim-developers.github.io/GalSim/_build/html/simple.html) again and [PSF-specific profiles](https://galsim-developers.github.io/GalSim/_build/html/psf.html) in the GalSim documentation; then, make some plots to explore the differences.\n",
    "Note that some of the galaxy profiles are PSF profiles are the same!\n",
    "As before, use the [GalSim paper](https://arxiv.org/abs/1407.7676) as a reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e9f5f24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# placeholder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8180076c",
   "metadata": {},
   "source": [
    "## Space is Crowded"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6d261ed",
   "metadata": {},
   "source": [
    "With real data, we often don't get to work with images of isolated objects.\n",
    "Instead, we will find many stars and galaxies in a given image.\n",
    "Moreover, we are looking at the projection of the actual positions of the objects to the two-dimensional surface that is the sky---objects that are extremely far apart in space may look close together to us (i.e., if one is much closer to us, but their angular positions are similar).\n",
    "You will see that this introduces particular challenges to our science."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c229aa3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = galsim.Image(320, 320, scale=0.2)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "psf = galsim.Gaussian(half_light_radius=0.9)\n",
    "\n",
    "n_objects = 100  # try changing this!\n",
    "xs = rng.uniform(-320 // 2, 320 // 2, n_objects) * 0.2\n",
    "ys = rng.uniform(-320 // 2, 320 // 2, n_objects) * 0.2\n",
    "objs = rng.uniform(0, 1, n_objects)\n",
    "\n",
    "object_ratio = 0.5  # try changing this!\n",
    "\n",
    "for x, y in zip(xs, ys):\n",
    "    if rng.uniform() < object_ratio:\n",
    "        # draw a galaxy\n",
    "        e1 = rng.uniform(-0.5, 0.5)\n",
    "        e2 = rng.uniform(-0.5, 0.5)\n",
    "        ellipticity = galsim.Shear(e1=e1, e2=e2)\n",
    "        galaxy = galsim.Exponential(half_light_radius=0.7)\n",
    "        galaxy = galaxy.shear(ellipticity).shift(x, y)\n",
    "        observed = galsim.Convolve([galaxy, psf])\n",
    "    else:\n",
    "        # draw a star\n",
    "        star = galsim.DeltaFunction().shift(x, y)\n",
    "        observed = galsim.Convolve([star, psf])\n",
    "\n",
    "    observed.drawImage(image=image, add_to_image=\"True\")\n",
    "\n",
    "ax.imshow(image.array)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eadeeb3d",
   "metadata": {},
   "source": [
    "### Things to try\n",
    "\n",
    "- adjust the number of objects (`n_objects`) in the image\n",
    "- adjust the ratio of galaxies to stars (`object_ratio`) in the image\n",
    "- try different galaxy profiles\n",
    "    - can you make the galaxies each have different sizes, brightnesses, etc.?\n",
    "- try different PSFs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3ed8e6a",
   "metadata": {},
   "source": [
    "Some problems to think about here:\n",
    "- How can we tell which objects are stars and which objects are galaxies?\n",
    "- How can we measure the PSF of this image?\n",
    "- How can we disentangle \"blended\" objects (i.e., overlapping objects)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8405fa1c",
   "metadata": {},
   "source": [
    "## Weak Lensing Shears Galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "492e1283",
   "metadata": {},
   "source": [
    "Recall that our principal goal is to measure the weak lensing of galaxies to glean information about dark matter (the lens) and dark energy (the geometry).\n",
    "In the regime where lensing is sufficiently weak, it's primary effect is to induce a shear that distorts the shape of the galaxies at the percent level.\n",
    "Note that this distortion is degenerate with the intrinsic ellipticity of each galaxy (and the reason that we can use a GalSim `Shear` to make elliptical galaxies)!\n",
    "To measure a weak lensing signal, we must average the shape over a population of galaxies (we assume that the average intrinsic ellipcitity of galaxies is zero)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bc34a24",
   "metadata": {},
   "source": [
    "First, let's see sheared galaxies.\n",
    "We'll use a shear that is _ten times stronger_ than we see in real data to exaggerate the effects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "021fd786",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = galsim.Image(64, 64, scale=0.2)\n",
    "\n",
    "fig, axs = plt.subplots(3, 3)\n",
    "\n",
    "for j, g2 in enumerate([0.1, 0.0, -0.1]):\n",
    "    for i, g1 in enumerate([-0.1, 0.0, 0.1]):\n",
    "        shear = galsim.Shear(g1=g1, g2=g2)\n",
    "        galaxy = galsim.Exponential(half_light_radius=0.7)\n",
    "        sheared_galaxy = galaxy.shear(shear)\n",
    "        sheared_galaxy.drawImage(image=image)\n",
    "        axs[j, i].imshow(image.array)\n",
    "        axs[j, i].set_title(f\"$(g_1, g_2) = ({g1}, {g2})$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8d13592",
   "metadata": {},
   "source": [
    "Now let's draw several example images where a field of galaxies is being coherently sheared due to weak lensing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "851cb316",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n_objects = 100\n",
    "xs = rng.uniform(-320 // 2, 320 // 2, n_objects) * 0.2\n",
    "ys = rng.uniform(-320 // 2, 320 // 2, n_objects) * 0.2\n",
    "objs = rng.uniform(0, 1, n_objects)\n",
    "\n",
    "object_ratio = 0.5\n",
    "\n",
    "fig, axs = plt.subplots(3, 3)\n",
    "\n",
    "psf = galsim.Gaussian(half_light_radius=0.9)\n",
    "\n",
    "ellipticity_seed = 1\n",
    "\n",
    "for j, g2 in enumerate([0.1, 0.0, -0.1]):\n",
    "    for i, g1 in enumerate([-0.1, 0.0, 0.1]):\n",
    "        image = galsim.Image(320, 320, scale=0.2)\n",
    "        shear = galsim.Shear(g1=g1, g2=g2)\n",
    "        ellipticity_rng = np.random.default_rng(ellipticity_seed)\n",
    "        for obj, x, y in zip(objs, xs, ys):\n",
    "            if obj < object_ratio:\n",
    "                # draw a galaxy\n",
    "                e1 = ellipticity_rng.uniform(-0.5, 0.5)\n",
    "                e2 = ellipticity_rng.uniform(-0.5, 0.5)\n",
    "                ellipticity = galsim.Shear(e1=e1, e2=e2)\n",
    "                galaxy = galsim.Exponential(half_light_radius=0.7)\n",
    "                galaxy = galaxy.shear(ellipticity).shift(x, y)\n",
    "                sheared_galaxy = galaxy.shear(shear)\n",
    "                observed = galsim.Convolve([sheared_galaxy, psf])\n",
    "            else:\n",
    "                # draw a star\n",
    "                star = galsim.DeltaFunction().shift(x, y)\n",
    "                observed = galsim.Convolve([star, psf])\n",
    "            observed.drawImage(image=image, add_to_image=\"True\")\n",
    "        axs[j, i].imshow(image.array, origin=\"lower\")\n",
    "        axs[j, i].set_title(f\"$(g_1, g_2) = ({g1}, {g2})$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ab5ceea",
   "metadata": {},
   "source": [
    "This is the same field (stars and galaxies) being affected by different cosmic shears.\n",
    "What do you notice about the shapes and positions of each object between two different shears?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c88533b2",
   "metadata": {},
   "source": [
    "Shear measurements are affected by all of these systematics:\n",
    "- how do you identify blended objects?\n",
    "- how do you identify stars vs. galaxies?\n",
    "- how do you accurately measure PSFs?\n",
    "- how do you accurately measure galaxy shapes (especially in the presence of PSFs, blending)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70a275c5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

Download the ZIP of this gist.
Then create the conda environment with

conda env create -f environment.yaml

To run the notebook, activate the environment with

conda activate weak-lensing-systematics

then

jupyter notebook "Weak Lensing Systematics.ipynb"