shape_measurement_102.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "030f1314",
   "metadata": {},
   "source": [
    "# Shape Measurement 102\n",
    "\n",
    "This notebook has the first part of a small-ish tutorial on galaxy shape measurement. The goals of this tutorial are as follows\n",
    "\n",
    "1. Become familiar with how weak lensing shear changes the shapes of objects.\n",
    "2. Understand the basic principles behind `metacalibration`-based shape measurements.\n",
    "3. Learn how to use `metadetection` and `ngmix` to measure the shapes of objects.\n",
    "\n",
    "**Please do the 101 tutorial before doing this one!**\n",
    "\n",
    "In order to run the code in this tutorial, you will need the following packages installed locally\n",
    "\n",
    " - `galsim`\n",
    " - `numpy`\n",
    " - `matplotlib`\n",
    " - `metadetect`\n",
    " - `ngmix`\n",
    " - `meds`\n",
    " - `joblib`\n",
    " - `schwimmbad`\n",
    " - `tqdm`\n",
    " - `joblib`\n",
    " - `fitsio`\n",
    " - `pyyaml`\n",
    " - `des-sxdes`\n",
    " \n",
    "I suggest using `conda`. You can run the command\n",
    "\n",
    "```\n",
    "conda install galsim numpy metaplotlib metadetect ngmix meds joblib schwimmbad tqdm fitsio pyyaml des-sxdes\n",
    "```\n",
    "\n",
    "in your environment to get things going.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4178111e",
   "metadata": {},
   "source": [
    "## A. How does a weak lensing shear change the shape of an object?\n",
    "\n",
    "Remember last time we learned about the moments of a surface brightness profile $I(x,y$\n",
    "\n",
    "$$\n",
    "\\langle M_x\\rangle = \\frac{\\int I(x,y) x}{\\int I(x,y)}\\\\\n",
    "\\langle M_y\\rangle = \\frac{\\int I(x,y) y}{\\int I(x,y)}\\\\\n",
    "\\langle M_{xx}\\rangle = \\frac{\\int I(x,y) (x - M_x)^2}{\\int I(x,y)}\\\\\n",
    "\\langle M_{xy}\\rangle = \\frac{\\int I(x,y) (x - M_x)(y-M_y)}{\\int I(x,y)}\\\\\n",
    "\\langle M_{yy}\\rangle = \\frac{\\int I(x,y) (y - M_y)^2}{\\int I(x,y)}\n",
    "$$\n",
    "\n",
    "It turns out, one can compute specific combinations of these moments that are particularly sensitive to the shearing operation we discussed last time. These combinations are\n",
    "\n",
    "$$\n",
    "\\langle M_{r}\\rangle = M_{xx} + M_{yy}\\\\\n",
    "\\langle M_{+}\\rangle = M_{xx} - M_{yy}\\\\\n",
    "\\langle M_{\\times}\\rangle = 2 M_{xy}\\\\\n",
    "$$\n",
    "\n",
    "which can then be combined into \n",
    "\n",
    "$$\n",
    "e_1 = \\frac{M_{+}}{M_{r}}\\\\\n",
    "e_2 = \\frac{M_{\\times}}{M_{r}}\\\\\n",
    "$$\n",
    "\n",
    "I have not made the mathematical origin of these particular combinations clear at all, so this should be confusing. The key thing is to remember that one can compute moments and then use the equations above to compute the shape of the object $(e_1, e_2)$. Below we'll explore how shears change these moments.\n",
    "\n",
    "### Exercise 1: How does it change?\n",
    "\n",
    "In this exercise, we are going to draw objects with several different shears, make plots of them, and then explore how the object's moments change with shear. \n",
    "\n",
    "To do this, we'll want to use `galsim` to shear an object. This can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01bba66a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import galsim\n",
    "\n",
    "obj = galsim.Gaussian(sigma=1)\n",
    "sheared_obj = obj.shear(g1=0.01, g2=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e70f0f44",
   "metadata": {},
   "source": [
    "We have called the `.shear` command which returns a new object that has been sheared by the amount given. If we print the underlying representation of this object out, we can see that `galsim` is applying a transformation to the original object, which hints at what a shear actually means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c633144c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(repr(sheared_obj))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c99dec7",
   "metadata": {},
   "source": [
    "With these preliminaries out of the way, do the following.\n",
    "\n",
    "1. For each of the shears $(g_1, g_2) = \\{(0, 0), (0.1, 0), (-0.1, 0), (0, 0.1), (0, -0.1)\\}$, use `galsim` to make the sheared object and draw an image of it. Describe to yourself how the object changes in relationship to the shear you applied. Why might we have used $+$ and $\\times$ when describing the moments above?\n",
    "2. For a grid of shears on the 1-axis with $g_1$ going from -0.2 to 0.2 in steps of 0.025, measure the shapes `(e_1, e_2)` as defined above. Set $g_2=0$. What do you notice about the relationship between the shear and the shape when the shear is small?\n",
    "3. Repeat 2 but with the roles of $g_1$ and $g_2$ swapped so that $g_1=0$ and the grid of shears is in $g_2$.\n",
    "4. What is your guess for the generic relationship between shapes and shears?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69add774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# do your work here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c92bc4b4",
   "metadata": {},
   "source": [
    "## B. Counterfactual Images and Metacalibration\n",
    "\n",
    "In the last exercise, you should have discovered how the measured shape of an object changes as it is sheared by different amounts. In particular, for small or *weak* gravitational shears, we expect that\n",
    "\n",
    "$$\n",
    "e_i \\approx e^{0} + R_i g_i\n",
    "$$\n",
    "\n",
    "where $R_i$ is called the *response*. The response is formally defined as a derivative $\\left.de/dg\\right.|_{g=0}$, but this is a detail that is not needed to complete this tutorial. They key bit of information here is that the response is the slope of the linear relationship between the observed shape and the underlying shear, for small shears.\n",
    "\n",
    "### Exercise 2: Compute the response numerically with `galsim`!\n",
    "\n",
    "In this exercise, we will use galsim to simulate the same object with different shears. Then by measuring the moments of those simulated profiles, we will estimate the response.\n",
    "\n",
    "1. Using `galsim`, create a Gaussian object with `sigma=1` and a shear of `g1=0.01` and `g2=0`. Draw this object onto a grid of pixels with `scale=0.2` and measure the shape $(e_1, e_2)$.\n",
    "2. Repeat 2 above but with `g1=-0.01` and `g2=0.0`. \n",
    "3. Combine the results from 1 and 2 to compute $R_1$ via a finite-difference derivative\n",
    "$$\n",
    "R_1 = \\frac{e_1^{g1=+h} - e_1^{g1=-h}}{2h}\n",
    "$$\n",
    "This is called a central difference derivative is typically an accurate but efficient method for this kind of numerical work.\n",
    "\n",
    "4. Do the same steps 1-3 but with the roles of the shears/shapes on the 1-axis swapped with those on the 2-axis. In other words, use shears with `g1=0`, `g2=+/-0.01` and compute $R_2$ by using a central difference derivative with $e_2$.\n",
    "5. Now make the same object but apply the shear `g1=0.02` and `g2=0.0`. Measure its shape with the same pixel scale. Estimate the underlying shear applied to the object by inverting the linear relationship above\n",
    "$$\n",
    "g_{i}^{\\rm obs} \\approx e_i/R_i\n",
    "$$\n",
    "How close did you get?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56f1e184",
   "metadata": {},
   "outputs": [],
   "source": [
    "# do your work here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83083d00",
   "metadata": {},
   "source": [
    "## C. Metadetection/Metacalibration in real data with `ngmix`\n",
    "\n",
    "Congrats! By completing exercise 2 of this tutorial you have coded up an extremely simple, but non-trivial form of the `metacalibration` weak lensing shape measurement technique!\n",
    "\n",
    "Moving from here to applying this to real observed galaxies involves a lot of complications. Namely\n",
    "\n",
    "1. We don't know the true underlying surface brightness profile of real galaxies. Thus we cannot directly compute the response from the unsheared surface brightness profile like we did above.\n",
    "2. Real galaxies are observed through telescopes and we need to account for this effect when computing the response.\n",
    "3. Images of real galaxies have noise in them. If we take an image with noise in it and shear the image, we will get correlated noise which can cause bias in our shape measurements.\n",
    "\n",
    "The solutions to these problems are the subject of years of research, so I'll give you the answers.\n",
    "\n",
    "1. Instead of using the true surface brightness profile of an object, we use the observed image of the object itself. We can do this because the shear applied by nature is small, so that the difference between the response computed using the observed object and that using the true object is even smaller. This process generates a *counter-factual* image of the object at another shear, hence the title above.\n",
    "2. We measure the PSF of the telescope and atmosphere using stars around the galaxy. We can then use this observed PSF to remove these effects via a process called deconvolution. This process cannot fully recover the galaxy as it would look without the telescope and atmosphere, but it can remove enough of these effects for our purposes.\n",
    "3. We can apply special noise corrections to our images which account for the correlated noise. The details of this process are unimportant except to say that if we don't we see 5% to 10% biases our measurements.\n",
    "\n",
    "All of these algorithms have been coded up into two packages, `metadetect` and `ngmix`. These are the main tools you will be using. The next exercise introduces them to you and generates the first kinds of simulations you'll run to test shape measurement.\n",
    "\n",
    "### Exercise 3: Making Sims w/ `ngmix` and `galsim`\n",
    "\n",
    "In this exercise, we will write a function to create a simple simulation of a single object using `ngmix` and `galsim`. Your goal tasks are to\n",
    "\n",
    "1. Fill in the parts of the function below to make the sim.\n",
    "2. Run the simulation with different shears and different psf sizes, measuring the output shape using the moments code above. Try psf sizes from 0.8 to 1.2 in steps of 0.1. What is the relationship between the slope of output shape vs input shape and the size of the PSF?\n",
    "\n",
    "The simulation code will output an `ngmix.Observation` object. You can read more about them in the docs: https://github.com/esheldon/ngmix/blob/master/ngmix/observation.py#L66. Here is a short set of things you can do:\n",
    "\n",
    "```\n",
    "# get the image in this observation as a numpy array\n",
    "obs.image \n",
    "\n",
    "# get the image of the PSF as a numpy array\n",
    "obs.psf.image\n",
    "```\n",
    "\n",
    "Note that these objects are read only, so you cannot modify the images etc once they are made without special syntax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33df84cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngmix\n",
    "import numpy as np\n",
    "\n",
    "def sim_func(*, g1, g2, seed, psf_fwhm):\n",
    "    # this is an RNG you can use to draw random numbers as needed\n",
    "    # always use this RNG and not np.random directly\n",
    "    # doing this makes sure the code is reproducible\n",
    "    rng = np.random.RandomState(seed=seed)\n",
    "\n",
    "    # make an Exponential object in galsim with a half light radius of 0.5\n",
    "    gal =  # fill this in\n",
    "    \n",
    "    # make a Gaussian object in galsim with a fwhm of `psf_fwhm`\n",
    "    psf =  # fill this in\n",
    "    \n",
    "    # apply the input shear `g1`, `g2` to the galaxy `gal`\n",
    "    sheared_gal =  # fill this in\n",
    "    \n",
    "    # here we are going to apply a random shift to the object's center\n",
    "    dx, dy = 2.0 * (rng.uniform(size=2)-0.5) * 0.2\n",
    "    sheared_gal = sheared_gal.shift(dx, dy)\n",
    "    \n",
    "    # convolve the sheared galaxy with the psf\n",
    "    obj =  # fill this in\n",
    "    \n",
    "    # render the object and the PSF on odd sized images of 53 pixels on a side with \n",
    "    # a pixel scale of 0.2\n",
    "    obj_im =  # fill this in\n",
    "    psf_im =  # fill this in\n",
    "    \n",
    "    # now we are going to add noise to the object image and setup the ngmix data\n",
    "    cen = (53-1)/2\n",
    "    nse_sd = 1e-5\n",
    "    nse = rng.normal(size=obj_im.array.shape, scale=nse_sd)\n",
    "    nse_im = rng.normal(size=obj_im.array.shape, scale=nse_sd)\n",
    "    jac = ngmix.jacobian.DiagonalJacobian(scale=0.2, row=cen+dy/0.2, col=cen+dx/0.2)\n",
    "    psf_jac = ngmix.jacobian.DiagonalJacobian(scale=0.2, row=cen, col=cen)\n",
    "    \n",
    "    # we have to add a little noise to the PSf to make things stable\n",
    "    target_psf_s2n = 500.0\n",
    "    target_psf_noise = np.sqrt(np.sum(psf_im.array ** 2)) / target_psf_s2n    \n",
    "    psf_obs = ngmix.Observation(\n",
    "        image=psf_im.array,\n",
    "        weight=np.ones_like(psf_im.array)/target_psf_noise**2,\n",
    "        jacobian=psf_jac,\n",
    "    )\n",
    "    \n",
    "    # here we build the final observation\n",
    "    obj_obs = ngmix.Observation(\n",
    "        image=obj_im.array + nse,\n",
    "        noise=nse_im,\n",
    "        weight=np.ones_like(nse_im) / nse_sd**2,\n",
    "        jacobian=psf_jac,\n",
    "        bmask=np.zeros_like(nse_im, dtype=np.int32),\n",
    "        ormask=np.zeros_like(nse_im, dtype=np.int32),\n",
    "        psf=psf_obs,\n",
    "    )\n",
    "\n",
    "    return obj_obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab3377f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# do your work here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa58a7a3",
   "metadata": {},
   "source": [
    "### Exercise 4: Run metadetect!\n",
    "\n",
    "In this exercise you will run metadetection for real! Yay! Before we do that, we need to understand a few things about the results. We are aiming to measure the shear and are using the estimate $g^{obs}$ as defined above. There is a standard way to characterize the bias in this estimate, which in our notation is\n",
    "\n",
    "$$\n",
    "g^{obs} = (1+m)g^{true} + c\n",
    "$$\n",
    "\n",
    "Here $m$ is called the *multiplicative bias* and is the fractional error in our shear estimate. The quantity $c$ is called the additive bias and denotes any mean bias in our shear estimate. For metadetect, we expect \n",
    "\n",
    "$$\n",
    "m\\approx 0.4\\times10^{-3}\\\\\n",
    "c\\approx 0\n",
    "$$\n",
    "\n",
    "Due to measurement noise in the simulations, we usually check that these conditions are true to within the $3\\sigma$ error bar.\n",
    "\n",
    "In the python source file `mdet_meas_tools.py`, I have put a bunch of code that you will use to run the simulation you made above through metadetection. The main functions to worry about are\n",
    "\n",
    " - `run_mdet_sims` - used to make the measurements\n",
    " - `write_sim_data` - used to write the simulation data to disk\n",
    " - `read_sim_data` - used to read the simulation data from disk\n",
    " - `estimate_m_and_c` - used to estimate $m$ and $c$ from a sim\n",
    " \n",
    "Read the doc strings on these functions for instructions on how to use them. Note that you need to pass a python dictionary for the `sim_kwargs` input and it should have the form `{\"psf_fwhm\": 0.9}`.\n",
    "\n",
    "Do the following:\n",
    "\n",
    "1. Use `run_mdet_sims` to run metadetection on your sim and verify that $m$ and $c$ meet the requirements above.\n",
    "2. Use the utilities to write your data to and then read your data from disk.\n",
    "3. Run your data through `estimate_m_and_c` and confirm you get the same answer as in 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "294148d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# do your work here!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd43baa4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:anl] *",
   "language": "python",
   "name": "conda-env-anl-py"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}