workshops/lsst2018/intro-with-globular.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing a (fake) globular cluster with the DM Science Pipelines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook provides an introduction to using some of the most important pieces of the DM Science Pipelines codebase:\n",
    "\n",
    " - Geometry classes from `lsst.geom`, such as points and boxes.\n",
    " - Higher-level astronomical primitives from `lsst.afw`, such as the `Image`, `Exposure`, and `Psf` classes.\n",
    " - Our core algorithmic `Task` classes, including those for source detection, deblending, and measurement.\n",
    " \n",
    "We'll be working with coadded images made from Subaru Hyper Suprime-Cam (HSC) data in the COSMOS field, augmented with a simulated globular cluster.  We've taken a recent LSST reprocessing of the HSC-SSP UltraDeep COSMOS field (see [this page](https://confluence.lsstcorp.org/display/DM/S18+HSC+PDR1+reprocessing) for information on that reprocessing, and [this page](https://hsc-release.mtk.nao.ac.jp/doc/) for the data), and added simulated stars from a scaled [SDSS catalog](http://www.sdss.org/dr14/data_access/value-added-catalogs/?vac_id=photometry-of-crowded-fields-in-sdss-for-galactic-globular-and-open-clusters).  The result is a very deep image (deeper than the 10-year LSST Deep-Wide-Fast survey, though not as deep as LSST Deep Drilling fields will be) with both a large number of galaxies and region full of stars.  As we'll see, that'll present a challenge for the vanilla DM pipelines (at least today), and hence a good excuse to do some custom processing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and Custom Installs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start with some standard imports of both LSST and third-party packages.\n",
    "\n",
    "Most of these are already installed automatically (assuming you're using an \"LSST\" Jupyter Kernel, not a vanilla Python one - and you should be).  But the `tqdm` package, a nice collection of progress bar widgets for both notebooks and the command-line, is not.  That means we can use it to show how to install custom PyPI packages in the notebook environment.  Open a new Terminal tab (either via the Launcher tab, if you have one, or by the \"+\" button in the upper left), and type:\n",
    "\n",
    "    pip install tqdm --user\n",
    "\n",
    "You may need to then restart the kernel in this notebook afterwards (circle-with-arrow button, at the top of this tab).  If you can execute the next cell, it's properly installed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from tqdm import tqdm_notebook\n",
    "from lsst.daf.persistence import Butler\n",
    "from lsst.geom import Box2I, Box2D, Point2I, Point2D, Extent2I, Extent2D\n",
    "from lsst.afw.image import Exposure, Image, PARENT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be retrieving data using the `Butler` tool, which manages where various datasets are stored on the filesystem (and can in principle manage datasets that aren't even stored as files, though all of these are).\n",
    "\n",
    "We start by creating a `Butler` instance, pointing it at a *Data Repository* (which here is just a root directory).  If you're interesting in looked at the original HSC data without the simulated cluster, change the path below to just `/datasets/hsc/cosmos`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butler = Butler(\"/project/jbosch/tutorials/lsst2018/data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Datasets managed by a butler are identified by a dictionary *Data ID* (specifying things like the visit number or sky patch) and a string *DatasetType* (such as a particular image or catalog).  Different DatasetTypes have different keys, while different instances of the same Dataset Type have different values.  All of the datasets we use in this tutorial will correspond to the same patch of sky, so they'll have at least the keys in the dictionary in the next cell (they will also have `filter`, but with different values):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = {\"tract\": 9813, \"patch\": \"4,4\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use those to load a set of *griz* coadds, which we'll put directly in a dictionary.  The result of each `Butler.get` call is in this case an `lsst.afw.image.Exposure` object, an image that actually contains three \"planes\" (the main image, a bit mask, and a variance image) as well as many other objects that describe the image, such as its PSF and WCS.  Note that we (confusingly) use `Exposures` to hold coadd images as well as true single-exposure images.\n",
    "\n",
    "The DatasetType here is `deepCoadd_calexp` (a coadd on which we've already done some additional processing, such as subtracting the background and setting some mask values), and the extra `filter` argument gets appended to the Data ID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coadds = {b: butler.get(\"deepCoadd_calexp\", dataId, filter=\"HSC-{}\".format(b.upper())) for b in \"griz\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making and displaying color composite images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by just looking at the images, as 3-color composites.  We'll use astropy to build those as a nice way to demonstrate how to get NumPy arrays from the `Exposure` objects in the `coadds` dict.  (LSST also has code to make 3-color composites using the same algorithm, and in fact the Astropy implementation is based on ours, but now that it's in Astropy we'll probably retire ours.)\n",
    "\n",
    "We'll just use matplotlib to display the images themselves.  We'll use Firefly for other image display tasks later, but while Firefly itself supports color-composites, work on our preferred composition algorithm is still in progress, and we haven't quite finished connecting that functionality to the Python client we'll demonstrate here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.visualization import make_lupton_rgb\n",
    "from matplotlib import pyplot\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the following function a few times to display color images.  It's worth reading through the implementation carefully to see what's going on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def showRGB(exps, bgr=\"gri\"):\n",
    "    \"\"\"Display an RGB color composite image with matplotlib.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    exps : `dict`\n",
    "        Dictionary of `lsst.afw.image.Exposure` objects, keyed by filter name.\n",
    "    bgr : sequence\n",
    "        A 3-element sequence of filter names (i.e. keys of the exps dict) indicating what band\n",
    "        to use for each channel.\n",
    "    \"\"\"\n",
    "    # Extract the primary image component of each Exposure with the .image property, and use .array to get a NumPy array view.\n",
    "    rgb = make_lupton_rgb(image_r=exps[bgr[2]].image.array,  # numpy array for the r channel\n",
    "                          image_g=exps[bgr[1]].image.array,  # numpy array for the g channel\n",
    "                          image_b=exps[bgr[0]].image.array,  # numpy array for the b channel\n",
    "                          stretch=1, Q=10)  # parameters used to stretch and scale the pixel values\n",
    "    pyplot.figure(figsize=(20, 15))\n",
    "    # Exposure.getBBox() returns a Box2I, a box with integer pixel coordinates that correspond to the centers of pixels.\n",
    "    # Matplotlib's `extent` argument expects to receive the coordinates of the edges of pixels, which is what\n",
    "    # this Box2D (a box with floating-point coordinates) represents.\n",
    "    integerPixelBBox = exps[bgr[0]].getBBox()\n",
    "    bbox = Box2D(integerPixelBBox)\n",
    "    pyplot.imshow(rgb, interpolation='nearest', origin='lower', extent=(bbox.getMinX(), bbox.getMaxX(), bbox.getMinY(), bbox.getMaxY()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showRGB(coadds, bgr=\"gri\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showRGB(coadds, bgr=\"riz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those images are a full \"patch\", which is our usual unit of processing for coadds - it's about the same size as a single LSST sensor (exactly the same in pixels, smaller in terms of area because these use HSC's smaller pixel scale).  That's a bit unweildy (just because waiting for processing to happen isn't fun in a tutorial setting), so we'll reload our dict with sub-images centered on the cluster.  Note that we can load the sub-images directly with the `butler`, by appending `_sub` to the DatasetType and passing a `bbox` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coadds = {b: butler.get(\"deepCoadd_calexp_sub\", dataId, filter=\"HSC-{}\".format(b.upper()),\n",
    "                        bbox=Box2I(corner=Point2I(18325, 17725), dimensions=Extent2I(400, 350)))\n",
    "          for b in \"griz\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showRGB(coadds, \"gri\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll try the regular LSST processing tasks, with a simpler configuration than we usually use to process coadds, just to avoid being distracted by complexity.  This includes\n",
    "\n",
    " - Detection (`SourceDetectionTask`): given an `Exposure`, find above-threshold regions and peaks within them (`Footprints`), and create a *parent* source for each `Footprint`.\n",
    " - Deblending (`SourceDeblendTask`): given an `Exposure` and a catalog of parent sources, create a *child* source for each peak in every `Footprint` that contains more than one peak.  Each child source is given a `HeavyFootprint`, which contains both the pixel region that source covers and the fractional pixel values associated with that source.\n",
    " - Measurment (`SingleFrameMeasurementTask`): given an `Exposure` and a catalog of sources, run a set of \"measurement plugins\" on each source, using deblended pixel values if it is a child.\n",
    "\n",
    "We'll start by importing these, along with the `SourceCatalog` class we'll use to hold the outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SourceDetectionTask\n",
    "from lsst.meas.deblender import SourceDeblendTask\n",
    "from lsst.meas.base import SingleFrameMeasurementTask\n",
    "from lsst.afw.table import SourceCatalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now construct all of these `Tasks` before actually running any of them.  That's because `SourceDeblendTask` and `SingleFrameMeasurementTask` are constructed with a `Schema` object that records what fields they'll produce, and they modify that schema when they're constructed by adding columns to it.  When we run the tasks later, they'll need to be given a catalog that includes all of those columns, but we can't add columns to a catalog that already exists.\n",
    "\n",
    "To recap, the sequence looks like this:\n",
    "\n",
    " 1. Make a (mostly) empty schema.\n",
    " 2. Construct all of the `Task`s (in the order you plan to run them), which adds columns to the schema.\n",
    " 3. Make a `SourceCatalog` object from the *complete* schema.\n",
    " 4. Pass the same `SourceCatalog` object to each `Task` when you run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = SourceCatalog.Table.makeMinimalSchema()\n",
    "\n",
    "detectionTask = SourceDetectionTask(schema=schema)\n",
    "\n",
    "deblendTask = SourceDeblendTask(schema=schema)\n",
    "\n",
    "# We'll customize the configuration of measurement to just run a few plugins.\n",
    "# The default list of plugins is much longer (and hence slower).\n",
    "measureConfig = SingleFrameMeasurementTask.ConfigClass()\n",
    "measureConfig.plugins.names = [\"base_SdssCentroid\", \"base_PsfFlux\", \"base_SkyCoord\"]\n",
    "# \"Slots\" are aliases that provide easy access to certain plugins.\n",
    "# Because we're not running the plugin these slots refer to by default,\n",
    "# we need to disable them in the configuration.\n",
    "measureConfig.slots.apFlux = None\n",
    "measureConfig.slots.instFlux = None\n",
    "measureConfig.slots.shape = None\n",
    "measureConfig.slots.modelFlux = None\n",
    "measureConfig.slots.calibFlux = None\n",
    "measureTask = SingleFrameMeasurementTask(config=measureConfig, schema=schema)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step we'll run is detection, which actually returns a new `SourceCatalog` object rather than working on an existing one.\n",
    "\n",
    "Instead, it takes a `Table` object, which is sort of like a factory for records.  We won't use it directly after this, and it isn't actually necessary to make a new `Table` every time you run `SourceDetectionTask` (but you can only create one after you're done adding columns so the schema).\n",
    "\n",
    "`Task`s that return anything do so via a `lsst.pipe.base.Struct` object, which is just a simple collection of named attributes.  The only return values we're  interested is `sources`.  That's our new `SourceCatalog`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = SourceCatalog.Table.make(schema)\n",
    "detectionResult = detectionTask.run(table, coadds['r'])\n",
    "catalog = detectionResult.sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a quick look at what's in that catalog.  First off, we can look at its schema:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog.schema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this includes a lot of columns that were actually added by the deblend or measurement steps; those will all still be blank (`0` for integers or flags, `NaN` for floating-point columns).\n",
    "\n",
    "In fact, the only columns filled by `SourceDetectionTask` are the IDs.  But it also attaches `Footprint` objects, which don't appear in the schema.  You can retrieve the `Footprint` by calling `getFootprint()` on a row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "footprint = catalog[0].getFootprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Footprints` have two components:\n",
    " - a `SpanSet`, which represents an irregular region on an image via a list of (y, x0, x1) `Spans`;\n",
    " - a `PeakCatalog`, a slightly different kind of catalog whose rows represent peaks within that `Footprint`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(footprint.getSpans())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(footprint.getPeaks())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's worth noting that while the peaks *can* have both an integer-valued position and a floating-point position, they're the same right now; `SourceDetectionTask` currently just finds the pixels that are local minima and doesn't try to find their sub-pixel locations.  That's left to the centroider, which is part of the measurement stage.\n",
    "\n",
    "Before we can get to that point, we need to run the deblender:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deblendTask.run(coadds['r'], catalog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`SourceDeblendTask` doesn't actually return anything - all of its outputs are just modifications to the catalog that's passed in.  It both sets some columns (`parent` and all those whose names start with `deblend_`) and creates new rows (for the child sources).  It does *not* remove the parent rows it created those child rows from, and this is intentional, because we want to measure both \"interpretations\" of the blend family: one in which there is only one object (the parent version) and one in which there are several (the children). Before doing any science with the outputs of an LSST catalog, it's important to remove one of those interpretations (typically the parent one).  That can be done by looking at the `deblend_nChild` and `parent` fields:\n",
    "\n",
    " - `parent` is the ID of the source from which this was deblended, or `0` if the source is itself a parent.\n",
    " - `deblend_nChild` is the number of child sources this source has (so it's `0` for sources that are themselves children or were never blended).\n",
    " \n",
    "Together, these define two particularly useful filters:\n",
    "\n",
    " - `deblend_nChild == 0`: never-blended object or de-blended child\n",
    " - `deblend_nChild == 0 and parent == 0`: never-blended object\n",
    " \n",
    "The first is what you'll usually want to use; the second is what to use if you're willing to throw away some objects (possibly many) because you don't trust the deblender."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last processing step for our purposes is running measurement:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measureTask.run(catalog, coadds['r'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting processing outputs with Firefly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll look at the results using Firefly.  Because these are the low-level outputs of a particular configuration of some processing `Tasks`, Firefly doesn't know nearly as much about them as it will about the (more \"curated\") data products that will appear in actual LSST Data Releases.  The tooling to connect Firefly to Python is also pretty new, so you can expect a lot of improvement in the future in both what Firefly can do with the datasets we send it and how easy it is is to send them.\n",
    "\n",
    "First some imports and setup:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.display import setDefaultBackend, Display\n",
    "from IPython.display import IFrame\n",
    "import os\n",
    "setDefaultBackend(\"firefly\")\n",
    "channel = '{}_test_channel'.format(os.environ['USER'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In LSST Science Platform environments, we have prepopulated environment variables `FIREFLY_URL` and `FIREFLY_HTML` with the Firefly server URL and the landing page or html file on that server. These are picked up automatically when creating the display. To use an alternate server, you can supply a `url` argument to `Display` such as `url='https://lsst-demo.ncsa.illinois.edu/firefly'`.\n",
    "\n",
    "In the previous cell, we generated a user-specific channel name; if unspecified, a random channel name is generated, that will change every time the notebook is executed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display = Display(frame=1, name=channel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The link that appears above will open a different browser tab/window for the Firefly viewer. Execute the next cell if you want to use Firefly from within a notebook cell. There's nothing stopping you from using both, but it can be a little confusing since you might need to refresh one or wait a bit to see changes you make in the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IFrame(display._url, 1024, 768)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are using the IFrame in the cell above and are in Jupyterlab, you can move it to its own output tab. In the left edge of the cell, right-click and select `Create New View for Output`. This creates a new tab with the IFrame inside. Rearrange this inside Jupyterlab, to keep the IFrame visible while you are working through the rest of the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now display an `Exposure` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display.mtv(coadds[\"r\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After doing that, take some time to play around with the GUI. See [this help page](http://irsa.ipac.caltech.edu/onlinehelp/finderchart/visualization.html?bd=2018-07-12#imageoptions) for help with the Firefly toolbar. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a lot of transparent mask overlays that can make it hard to see the image, but the overlays button ![overlays button](http://irsa.ipac.caltech.edu/onlinehelp/finderchart/img/layers.png) on the toolbar gives you very detailed control over them. Mask transparency can also be controlled using `afw.display`, for individual planes or for all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display.setMaskTransparency(20)\n",
    "display.setMaskTransparency(10, 'DETECTED')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As of the `w_2018_32` build make these transparency settings will be \"sticky\" for each Display, and settable by default with `lsst.afw.display.setDefaultMaskTransparency`. For now, we can set the transparency after each image display."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scale can also be changed programmatically, or by `afw.display` commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display.scale(\"asinh\", \"zscale\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function below overplots the sources we detected on the image.  Firefly itself has much more sophisticated ways of interacting them with catalogs; an example of overlaying a catalog will be shown at the end of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def showCatalog(disp, cat, color='orange', sym='+'):\n",
    "    \"\"\"Display sources from a catalog.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    disp : `lsst.afw.display.Display`\n",
    "        Display interface object to send to.\n",
    "    cat : `lsst.afw.table.SourceCatalog`\n",
    "        Catalog containing sources to display.  Must have populated centroid columns.\n",
    "    color : `str`\n",
    "        Color to use for overlaid points (X11 color string).\n",
    "    sym : `str`\n",
    "        Symbol; one of \"+\", \"x\", or \"o\".\n",
    "    \"\"\"\n",
    "    with disp.Buffering():\n",
    "        for record in cat:\n",
    "            if record[\"deblend_nChild\"] != 0:  # don't show \"one object\" interpretation of blends\n",
    "                continue\n",
    "            disp.dot(sym, record.getX(), record.getY(), size=5, ctype=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showCatalog(display, catalog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Subtracting Stars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how well that measurement worked, we'll use the PSF model and the measured PSF fluxes to subtract all of the objects.  Some of those objects are actually well-resolved galaxies, so we don't expect them to subtract well, but we'll just ignore that for now and focus our attention on the stars when we look at the results.\n",
    "\n",
    "Once again, it's worth taking some time to read carefully through the function below before just running it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subtractStars(data, cat, fluxes=None):\n",
    "    \"\"\"Subtract point sources from an image.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : `lsst.afw.image.Exposure`\n",
    "        Image from which point sources will be subtracted (not modified).\n",
    "    cat : `lsst.afw.table.SourceCatalog`\n",
    "        Catalog providing centroids and possibly fluxes.\n",
    "    fluxes : `numpy.ndarray`, optional\n",
    "        If not None, an array of the same length as ``cat`` containing fluxes\n",
    "        to use.  If None, ``cat.getPsfFlux()` will be used instead.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    model : `lsst.afw.image.Exposure`\n",
    "        Image containing just the point source models.\n",
    "    residuals : `lsst.afw.image.Exposure`\n",
    "        ``data`` with ``model`` subtracted.\n",
    "    \"\"\"\n",
    "    # Make sure no one calls this still-blended sources in the catalog.\n",
    "    assert (cat[\"deblend_nChild\"] == 0).all()\n",
    "    # Get the PSF model from the given Exposure object.  The returned object\n",
    "    # can be evaluated anywhere in the image to obtain an image of the PSF at\n",
    "    # that point.\n",
    "    psf = data.getPsf()\n",
    "    # Make a new blank Exposure object with the same dimensions and pixel type.\n",
    "    model = Exposure(data.getBBox(), dtype=np.float32)\n",
    "    # Copy the WCS and PSF over from the original Exposure.\n",
    "    model.setWcs(data.getWcs())\n",
    "    model.setPsf(psf)\n",
    "    if fluxes is None:\n",
    "        fluxes = cat.getPsfFlux()\n",
    "    for flux, record in zip(fluxes, tqdm_notebook(cat)):\n",
    "        # Obtain a PSF model image at the position of this source\n",
    "        psfImage = psf.computeImage(record.getCentroid())\n",
    "        # Make sure the PSF model image fits within the larger image; if it doesn't, clip it so it does.\n",
    "        psfBBox = psfImage.getBBox()\n",
    "        if not data.getBBox().contains(psfBBox):\n",
    "            psfBBox.clip(data.getBBox())  # shrink the bounding box to the intersection\n",
    "            psfImage = psfImage[psfBBox, PARENT]   # obtain a subimage\n",
    "        # Make a subimage view of `model`, and subtract the PSF image, scaled by the flux.\n",
    "        # PARENT here sets the coordinate system to be the one shared by all patches in\n",
    "        # the tract rather than the one in which this patches' origin is (0, 0).\n",
    "        # PARENT is the coordinate system used by the PSF, and it will soon be the default\n",
    "        # here too (but isn't yet, so we need to make that explicit).\n",
    "        model.image[psfBBox, PARENT].scaledPlus(flux, psfImage.convertF())\n",
    "    # Now that we've made a model image, make a copy of the data and subtract the model from it.\n",
    "    residuals = data.clone()\n",
    "    residuals.image -= model.image\n",
    "    return model, residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can actually subtract the stars, we should remove the not-deblended sources; this function and the next one we write will assume they've been removed.\n",
    "\n",
    "We can do that with NumPy-style boolean indexing, with one catch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deblended = catalog[catalog[\"deblend_nChild\"] == 0].copy(deep=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Column arrays of `SourceCatalogs` can only be accessed when the catalog is stored in a single contiguous block of memory.  But unlike Numpy arrays, using boolean indexing on a catalog doesn't automatically make a copy to ensure memory is contiguous.  Instead it creates a view to the selected rows.  That can be useful or more efficient in some cases, but it also prevents us from accessing columns.  To fix that, we immediately make a deep copy of the catalog, which copies it into a new block of contiguous memory.\n",
    "\n",
    "We can now run our `subtractStars` function.  Note the nice `tqdm` progress bar in action!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model, residuals = subtractStars(coadds['r'], deblended)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To look at the results, we'll make two more display frames to show the residuals (frame 2) and model (frame 3).\n",
    "\n",
    "Use the \"WCS Match\" and single-frame view GUI options to blink between them and the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display2 = Display(frame=2, name=channel)\n",
    "display2.mtv(residuals, title=\"residuals\")\n",
    "display2.setMaskTransparency(20)\n",
    "display2.setMaskTransparency(10, 'DETECTED')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display3 = Display(frame=3, name=channel)\n",
    "display3.mtv(model, title=\"model\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The isolated objects seem to be subtracted well, and *some* of the deblended ones are too.  But others really didn't subtract well at all; the deblender is probably mangling those.  The current deblender depends a lot on having at least some sides of an object not blended with a neighbor, and that's manifestly untrue in the center of the cluster.  We're currently working on a new deblender ([Scarlet](https://github.com/fred3m/scarlet)) that we expect to handle this better.\n",
    "\n",
    "In any case, if we look at the overplotted positions, it looks like the centroids aren't too bad, even if the fluxes frequently are.  It seems our centroider is pretty robust to whatever failure modes the deblender is exhibiting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting Stars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For better photometry, let's try fitting the stars ourselves, simultaneously.  We can use the same PSF object we used to subtract stars to instead construct a model for the entire image, and then we can use basic linear least squares from NumPy to fit it.\n",
    "\n",
    "Once again, this will go poorly for galaxies, and we'll ignore that.  We could also imagine adding terms to fit small offsets in the positions, or using sparse matrices to make this scale better.  But this works well enough for this tutorial.\n",
    "\n",
    "This is another function worth reading through carefully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fitStars(data, cat):\n",
    "    \"\"\"\n",
    "    Fit all sources in the given catalog simultaneously, assuming they are all point sources.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : `lsst.afw.image.Exposure`\n",
    "        Image to fit and subtract model from.\n",
    "    cat : `lsst.afw.table.SourceCatalog`.\n",
    "        Catalog providing centroids.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    fluxes : `numpy.ndarray`\n",
    "        Array of best-fit flux values.\n",
    "    model : `lsst.afw.image.Exposure`\n",
    "        Image containing just the point source models.\n",
    "    residuals : `lsst.afw.image.Exposure`\n",
    "        ``data`` with ``model`` subtracted.\n",
    "    \"\"\"\n",
    "    assert (cat[\"deblend_nChild\"] == 0).all()\n",
    "    bbox = data.getBBox()\n",
    "    psf = data.getPsf()\n",
    "    # Dimensions of the problem:\n",
    "    N = len(cat)\n",
    "    W = bbox.getWidth()\n",
    "    H = bbox.getHeight()\n",
    "    M = W*H\n",
    "    # To solve a lineat-least squares problem, we need to construct a matrix that maps\n",
    "    # parameters (fluxes) to data (pixel values).  That'd have dimensions M x N.\n",
    "    # Instead, we'll construct an array with dimensions N x H x W, which we'll later\n",
    "    # transpose and reshape.  It's important that we order the dimensions this way\n",
    "    # because it makes each nested H x W array have the same memory layout as the\n",
    "    # LSST Image class, which lets us make image views into those subarrays.\n",
    "    matrix = np.zeros((N, H, W), dtype=float)\n",
    "    for n, record in enumerate(tqdm_notebook(cat)):\n",
    "        # Make an Image view to a nested sub-array.  Note that writing to this\n",
    "        # will modify the parent array.\n",
    "        matrixView = Image(matrix[n, :, :], xy0=bbox.getMin(), dtype=np.float64)\n",
    "        # Obtain a PSF image, and clip it to fit the larger image, just as we\n",
    "        # did in subtractStars().\n",
    "        psfImage = psf.computeImage(record.getCentroid())\n",
    "        psfBBox = psfImage.getBBox()\n",
    "        if not bbox.contains(psfBBox):\n",
    "            psfBBox.clip(bbox)\n",
    "            psfImage = psfImage[psfBBox, PARENT]\n",
    "        # Add the PSF image to the matrix sub-array.\n",
    "        matrixView[psfBBox, PARENT].scaledPlus(record.getPsfFlux(), psfImage)\n",
    "    # Reshape and transpose the matrix, as promised\n",
    "    A = matrix.reshape(N, M).transpose()\n",
    "    # Get an array view to the image we're fitting, and flatten it the same way.\n",
    "    b = data.image.array.reshape(M)\n",
    "    # Fit for the fluxes.  This does an SVD under the hood.\n",
    "    fluxes, _, _, _ = np.linalg.lstsq(A, b, rcond=None)\n",
    "    # Make a model image and subtract it from the data.  Because we've already\n",
    "    # evaluated the matrix, this is much more efficient than calling\n",
    "    # subtractStars with our flux array.\n",
    "    model = Exposure(bbox, dtype=np.float32)\n",
    "    model.setWcs(data.getWcs())\n",
    "    model.setPsf(psf)\n",
    "    model.image.array[:, :] = np.dot(A, fluxes).reshape(bbox.getHeight(), bbox.getWidth())\n",
    "    residuals = data.clone()\n",
    "    residuals.image -= model.image\n",
    "    return fluxes, model, residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes, model, residuals = fitStars(coadds['r'], deblended)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the new residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display2.mtv(residuals.image, title=\"residuals\")\n",
    "display2.setMaskTransparency(20)\n",
    "display2.setMaskTransparency(10, 'DETECTED')\n",
    "display3.mtv(model, title=\"model\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are better than the last subtraction, but still not great.\n",
    "\n",
    "Note that you may need to set the stretch manually (color stretch icon ![stretch button](http://irsa.ipac.caltech.edu/onlinehelp/finderchart/img/stretch.png)) to give the residuals image the same stretch as the others.  For this image, you can set the stretch using \"Data\" bounds of `(-0.05, 2)` with the original image selected, and then blinking to the other frames while the stretch window is open. The lock icon ![lock icon](http://irsa.ipac.caltech.edu/onlinehelp/finderchart/img/lockimages.png) on the Firefly toolbar can be used to lock the stretch for the displayed frames."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stretch is also settable programmatically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display.scale('asinh', -0.1, 0.5, Q=8)\n",
    "display2.scale('asinh', -0.1, 0.5, Q=8)\n",
    "display3.scale('asinh', -0.1, 0.5, Q=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can erase the overlays on the first display:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display.erase()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To improve things further, let's try running the detect/deblend/measure `Tasks` on the *residual* image, to pick up those peaks that just weren't included in the first round.\n",
    "\n",
    "Here's a function that just runs those (relying on us having constructed them all earlier; we don't have to re-construct them to run them again):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process(data):\n",
    "    \"\"\"Run detection, deblending, and measurement, returning a new SourceCatalog.\n",
    "    \"\"\"\n",
    "    result = detectionTask.run(table, data)\n",
    "    cat = result.sources\n",
    "    deblendTask.run(data, cat)\n",
    "    measureTask.run(cat, data)\n",
    "    return cat[cat[\"deblend_nChild\"] == 0].copy(deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "additional = process(residuals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make use of those additional sources, we'll have to merge them with the original ones.  We'll do that extremely naively - we won't even worry about whether we have any duplicates (from sources that were partially - but not completely - subtracted by the model).\n",
    "\n",
    "Here's a function to do that merge.  Note that we reserve space in the output catalog before actually concatenating the input catalogs into it.  That makes sure the result is contiguous in memory and hence we can get column arrays without doing a deep copy at the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def concatenate(*cats):\n",
    "    \"\"\"Concatenate multiple SourceCatalogs (assumed to have the same schema).\"\"\"\n",
    "    result = SourceCatalog(table)\n",
    "    result.reserve(sum(len(cat) for cat in cats))\n",
    "    for cat in cats:\n",
    "        result.extend(cat)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined = concatenate(deblended, additional)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll fit all the stars in the combine catalog together, and look at the residuals again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes, model, residuals = fitStars(coadds['r'], combined)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display2.mtv(residuals, title=\"residuals\")\n",
    "display2.setMaskTransparency(5)\n",
    "display2.setMaskTransparency(30, 'DETECTED')\n",
    "display3.mtv(model, title=\"model\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display2.scale('asinh', -0.1, 0.5, Q=8)\n",
    "display3.scale('asinh', -0.1, 0.5, Q=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even better!  One more iteration *might* make things better, but we're probably getting to the limits of this very simple algorithm (especially since there really are a lot of galaxies in this image)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uploading the catalog directly to Firefly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As noted earlier, Firefly has sophisticated capabilities for viewing tables and overlaying catalogs. The `firefly_client` Python package underlies the Firefly backend in `lsst.afw.display`. The `firefly_client.plot` module includes convenience functions that we'll demonstrate here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import firefly_client.plot as ffplt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `firefly_client.plot` module should default to the client that is used in the `afw.display` Displays we defined earlier, so long as a browser tab is connected to those displays. To be certain, the next line ensures we use the same setup as our Display instances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ffplt.use_client(display._client)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Upload our combined catalog to Firefly. By default, the catalog will be shown in an interactive table viewer. Since the catalog contains coordinate columns recognized by Firefly, the catalog entries will be overlaid on the images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_id = ffplt.upload_table(combined, title='combined catalog')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firefly includes plotting capabilities using the Plotly.js library. We'll end this section with a histogram of flux values. `ffplt.scatter` can be used to make a scatter plot of two columns for a table. These functions use the table ID of the last uploaded table by default; the `tbl_id` we saved in the last cell can be passed as an optional argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ffplt.hist('log10(base_PsfFlux_flux)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table, plot, and image cells can be resized and moved around in the Firefly viewer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra Credit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One More Time\n",
    "\n",
    "Try doing one more round of detection/deblend/measure, concatenate, and fit.  Does it help?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-Band Fitting\n",
    "\n",
    "While we started this tutorial by making color images, we've only processed the *r*-band data.  Run `fitStars` on each band separately, using the positions from detect/deblend/measure processing in *r*.  Make *gri* and *riz* color images of the model and residuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deblender Output Inspection\n",
    "\n",
    "Instead of just guessing that the deblender is mangling things, we could actually look.  If you call `getFootprint()` on a row corresponding to a child source, you'll get a `HeavyFootprint`.  Use its `getBBox()` and `insert` methods to make images of some `HeavyFootprints`, and display them with matplotlib or Firefly.\n",
    "\n",
    "Note that if you try this on a source that isn't a deblended child, you'll get a regular `Footprint`, which doesn't have an `insert` method (because it doesn't contain any pixel values).  What catalog filter can you apply to get only child objects?\n",
    "\n",
    "If you're ambitious, import the [ipywidgets](https://ipywidgets.readthedocs.io/en/latest/user_guide.html) package, and use the `interact` function and an `IntSlider` to interactively control which source's `HeavyFootprint` is displayed."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}