finding_qsos.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Finding QSOs\n",
    "\n",
    "In this notebook we will attempt to mine for [Quasi-Stellar Objects](https://en.wikipedia.org/wiki/Quasar)\n",
    "(aka Quasars) in photometric catalogues. The objects mimic stars, since they are point-like, but are actually\n",
    "extra-galactic objects. Gaia provides equisite proper motion constraints and QSOs should appear as stationary objects.\n",
    "A cut on having ~zero proper motion is thus expected to be a powerful method to initially find QSO candidates.\n",
    "\n",
    "This notebook follows in part the methodology of [Heintz et al. 2018](https://arxiv.org/pdf/1805.03394.pdf), with\n",
    "additional scope by further investigating optical colours of QSOs and building a Random Forest classifier to find new\n",
    "QSO candidates.\n",
    "\n",
    "The notebook builds on concepts discussed in the previous session\n",
    "\n",
    "![concepts](concepts.png)\n",
    "\n",
    "The basic workflow is:\n",
    "\n",
    "1. Perform a high-galactic latitude 1-degree cone search in Gaia\n",
    "2. Calculate quantites and perform cuts on the data directly in Python\n",
    "3. Cross-match Gaia sources with SDSS photometry\n",
    "4. Find spectroscopically confirmed QSOs within the cross-matched sources\n",
    "5. Investigate features of QSOs compared to other field sources to find QSO candidates\n",
    "6. Train a Random Forest classifier to classify sources as QSOs based on proper-motion and colour\n",
    "\n",
    "Throughout this we will employ data-mining skills such as retrieving, exploring and analysing data efficiently, as well\n",
    "as employing simple machine learning techniques to our findings.\n",
    "\n",
    "There are example exercises contained within this notebook that you can choose from (or use as inspiration) to fulfill\n",
    "the assignment if you are required to do so for MPAGS credit."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Preamble\n",
    "\n",
    "Imports etc."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy import table\n",
    "from astroquery.gaia import Gaia\n",
    "from astroquery.sdss import SDSS\n",
    "\n",
    "# supress some numpy warnings about converting masked elements\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# do not limit the number of results returned by Gaia queries (default is 50 otherwise)\n",
    "Gaia.ROW_LIMIT = -1\n",
    "\n",
    "# plot styling\n",
    "plt.style.use(\"seaborn-darkgrid\")\n",
    "plt.rcParams['figure.figsize'] = [12, 8]"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 1. Perform a high-galactic latitude 1-degree cone search in Gaia\n",
    "\n",
    "We will use [Astroquery's Gaia interface](https://astroquery.readthedocs.io/en/latest/gaia/gaia.html)\n",
    "to programmatically query for sources in a 1-degree cone at a high galactic latitude.\n",
    "\n",
    "First define our cone centre. You can change the ra and dec values to make your own unique results, but\n",
    "keep within the [SDSS footprint](https://classic.sdss.org/dr7/coverage/) and at high galactic latitude - shifting\n",
    "either anywhere within +/-20 degrees is safe."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Try changing ra and dec to get your own results\n",
    "ra_cone = 180 * u.degree\n",
    "dec_cone = 30 * u.degree\n",
    "\n",
    "# These are used to define our cone centre, and we set the radius\n",
    "coord = SkyCoord(ra=ra_cone, dec=dec_cone)\n",
    "radius = 1 * u.degree"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Perform a Gaia cone search and get results\n",
    "job = Gaia.cone_search_async(coord, radius)  # use async if retrieving large results (>2000 rows)\n",
    "gaia_results = job.get_results()\n",
    "gaia_results  # the notebook will display the table in a nice format if we end the cell with it"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Our results are stored as an [Astropy Table](https://docs.astropy.org/en/stable/table/).\n",
    "Let's check what data we have in our table:\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "gaia_results.info"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "In order to further understand what the columns mean, we should always refer to the **schema** documentation\n",
    "for the database, which describes the data and its structure. For example, for the main `gaia_source`\n",
    "table we queried above, we can refer to\n",
    "[this page](https://gea.esac.esa.int/archive/documentation/GEDR3/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2. Calculate quantites and perform cuts on the data directly in Python\n",
    "\n",
    "Our data are held as an [Astropy Table](https://docs.astropy.org/en/stable/table/) - this allows a lot of quick\n",
    "numpy-like (array-based) analysis, with the benefit of a human-readable structure (column names etc.).\n",
    "\n",
    "We will investigate aspects of the data and to show examples of **array-based** operations and filtering data.\n",
    "\n",
    "### Calculate total proper motion of sources\n",
    "\n",
    "In our hyopthesis, a selection for sources consistent with zero proper motion should be a good first cut to find QSOs.\n",
    "Let's calculate the total proper motion (as the proper motion in ra and dec added in quadrature),\n",
    "and along the way demonstrate why we don't use loops whenever possible:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def get_pmtot_loop(tbl):\n",
    "    \"\"\"Don't do this!\"\"\"\n",
    "    pmtot = []\n",
    "    for row in tbl:\n",
    "        pmtot.append((row[\"pmra\"]**2 + row[\"pmdec\"]**2)**0.5)\n",
    "    return pmtot\n",
    "\n",
    "print(\"Timing for calculating pmtot with for loop:\")\n",
    "%time pmtot_loop = get_pmtot_loop(gaia_results)\n",
    "\n",
    "def get_pmtot(tbl):\n",
    "    \"\"\"Do this!\"\"\"\n",
    "    return (tbl[\"pmra\"]**2 + tbl[\"pmdec\"]**2)**0.5\n",
    "\n",
    "print(\"\\nTiming for calculating pmtot with array operation:\")\n",
    "%time pmtot = get_pmtot(gaia_results)\n",
    "\n",
    "assert np.all(pmtot_loop == pmtot), \"results don't match\"\n",
    "# this would raise an error if methods gave different results\n",
    "# (there is a slight detail in that the for loop returns a list, whereas the array operation returns a column\n",
    "# - practically they are the same results, but the column has additional astropy.table.Table methods.)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "A factor of ~hundreds speed up for this example!\n",
    "\n",
    "We can now assign these calculated total proper motion values to our main results.\n",
    "We also populate the corresponding proper motion uncertainty (taking care of the correlation coefficient), as well\n",
    "as the significance of the proper motion (`=pm/sigma_pm`)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "gaia_results[\"pmtot\"] = pmtot\n",
    "gaia_results[\"pmtot_error\"] = (\n",
    "                        (gaia_results[\"pmra_error\"]**2 + gaia_results[\"pmdec_error\"]**2\n",
    "                        + 2 * gaia_results[\"pmra_pmdec_corr\"] * gaia_results[\"pmra_error\"] * gaia_results[\"pmdec_error\"]\n",
    "                            )** 0.5)\n",
    "gaia_results[\"pm_significance\"] = gaia_results[\"pmtot\"]/gaia_results[\"pmtot_error\"]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Filter the sources based on magnitude\n",
    "\n",
    "Our aim is to find QSO candidates. We have astrophysical reasons to expect these are going to be relatively\n",
    "faint (cf. foreground Galactic stars, which are at a full range of brightness). Because of this, we\n",
    "will apply a simple cut on apparant magnitude to limit things at the faint end of Gaia's detection distribution.\n",
    "\n",
    "To decide, let's plot the apparant magnitude distribution in Gaia's `G` band."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "plt.hist(gaia_results[\"phot_g_mean_mag\"], bins=np.arange(10, 23, 0.5))\n",
    "plt.xlabel(\"G [mag]\")\n",
    "plt.ylabel(\"N\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see that Gaia detection efficiency falls off rapidly around 21 mag. We set our lower cut to 20 mag because we\n",
    "want reasonably robust photometry (and require that there is a full proper motion solution). Bright sources are very\n",
    "likely to be foreground stars, so we also cut anything brighter than 18 mag."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Produce a boolean mask used to select rows meeting our criteria\n",
    "mag_mask = (gaia_results[\"phot_g_mean_mag\"] > 18) & (gaia_results[\"phot_g_mean_mag\"] < 20)\n",
    "print(f\"mag_mask = {mag_mask}\")\n",
    "print(f\"There are {np.sum(mag_mask)} sources that pass magnitude cuts\")\n",
    "gaia_results_faint = gaia_results[mag_mask]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "(In reality you would perform such cuts prior to calculating any quantities on the data, for efficiency.)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 3. Cross-match Gaia sources with SDSS photometry\n",
    "\n",
    "Gaia contains pre-baked cross-matches to external catalogues we could use, but instead we will perform manual\n",
    "cross-matching in Python to demonstrate it.\n",
    "\n",
    "\n",
    "### Query SDSS for sources\n",
    "\n",
    "Firstly, we need to grab SDSS sources in our cone. The Astroquery SDSS interface is a bit different from the Gaia one,\n",
    "unfortunately, and we will also need to specify the columns we want explicitly. Again, we'll refer to the\n",
    "[schema](https://skyserver.sdss.org/dr14/en/help/browser/browser.aspx#&&history=description+PhotoObjAll+U) to find the\n",
    "columns we want - we're interested in the position, type, and psf photometry of each source"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "photoobj_fields=[\"ra\", \"dec\", \"type\", \"psfMag_u\", \"psfMag_g\", \"psfMag_r\", \"psfMag_i\", \"psfMag_z\"]\n",
    "region_radius = radius / np.cos(np.deg2rad(coord.dec.degree))  # this is a correction since query_region doesn't account for cos(dec) factor in ra.\n",
    "sdss_results = SDSS.query_region(coord, region_radius, photoobj_fields=photoobj_fields)\n",
    "# immediately cut out any that are not of type=6 (i.e. star)\n",
    "sdss_results = sdss_results[sdss_results[\"type\"] == 6]\n",
    "sdss_results\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Cross match Gaia and SDSS\n",
    "\n",
    "To perform a cross-match, we will use\n",
    "[Astropy's catalogue matching](https://docs.astropy.org/en/stable/coordinates/matchsep.html#matching-catalogs) tools.\n",
    "For this we need [`SkyCoord`](https://docs.astropy.org/en/stable/coordinates/) objects of our positions in each\n",
    "catalogue before performing sky-matching"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "gaia_coord = SkyCoord(gaia_results_faint[\"ra\"], gaia_results_faint[\"dec\"], unit=u.degree)\n",
    "sdss_coord = SkyCoord(sdss_results[\"ra\"], sdss_results[\"dec\"], unit=u.degree)\n",
    "idx, d2d, _ = gaia_coord.match_to_catalog_sky(sdss_coord)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "`idx` is now an array of indices in `sdss_coord` that are closest to each `gaia_coord`. We can use `d2d`, the 2D (on-sky)\n",
    "distance to make a cut on the matches to be within 1 arcsecond.\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "close_match = d2d < 1 * u.arcsec  # find where the separations are within our threshold\n",
    "print(f\"Found {np.sum(close_match)} close matches\")\n",
    "\n",
    "# limit our gaia results to those with a close match\n",
    "gaia_matched = gaia_results_faint[close_match]\n",
    "# select the close matched SDSS results\n",
    "sdss_matched = sdss_results[idx[close_match]]\n",
    "\n",
    "# a boolean mask denoting no proper motion matches\n",
    "no_pm_mask = gaia_matched[\"pm_significance\"] < 2"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "As a sanity check, plot the positions of the matched sources, as well as the no proper motion subset."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "plt.scatter(sdss_matched[\"ra\"], sdss_matched[\"dec\"], s=20, facecolor=None, label=\"Matched SDSS sources\")\n",
    "plt.scatter(gaia_matched[\"ra\"], gaia_matched[\"dec\"], s=5, label=\"Matched Gaia sources (18 < G < 20)\")\n",
    "plt.scatter(gaia_matched[no_pm_mask][\"ra\"], gaia_matched[no_pm_mask][\"dec\"], s=100, c=\"None\", marker=\"s\", edgecolors=\"C3\", linewidths=2, label=\"No PM sources\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.xlabel(\"RA [deg]\")\n",
    "plt.ylabel(\"Dec [deg]\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 4. Find spectroscopically confirmed QSOs within the cross-matched sources\n",
    "\n",
    "We can exploit spectroscopic datasets to find known QSOs in our cone. For this we will appeal to SDSS spectra. (Another\n",
    "general extragalactic tool to use could be the [NASA Extragalactic Database](https://ned.ipac.caltech.edu/).)\n",
    "\n",
    "Here we want sources with a spectrum associated (`spectro=True`), and we define key columns to return, again referring to\n",
    "the appropriate [schema](https://skyserver.sdss.org/dr14/en/help/browser/browser.aspx#&&history=description+SpecObjAll+U)\n",
    "information.\n",
    "\n",
    "*(The `sciencePrimary` is used since we query all spectra, and we only want the classification from the best\n",
    "science spectrum at this location. Otherwise some objects can have conflicting spectral classes.)*"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "specobj_fields=[\"ra\", \"dec\", \"class\", \"sciencePrimary\"]\n",
    "sdss_spec_results = SDSS.query_region(coord, region_radius, spectro=True, specobj_fields=specobj_fields)\n",
    "# immediately cut out any that are not sciencePrimary spectra\n",
    "sdss_spec_results = sdss_spec_results[sdss_spec_results[\"sciencePrimary\"] == 1]\n",
    "qso_mask = sdss_spec_results[\"class\"] == \"QSO\"  # a boolean mask denoting QSO spectra\n",
    "print(f\"Found {len(sdss_spec_results)} spectroscopic matches, with {np.sum(qso_mask)} QSOs\")\n",
    "sdss_spec_results"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's cross-match these QSO labels with our Gaia-SDSS sources so we have a record of what a known QSO looks like:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "sdss_qso_coord = SkyCoord(sdss_spec_results[qso_mask][\"ra\"], sdss_spec_results[qso_mask][\"dec\"], unit=u.degree)\n",
    "sdss_matched_coord = SkyCoord(sdss_matched[\"ra\"], sdss_matched[\"dec\"], unit=u.degree)\n",
    "idx, d2d, _ = sdss_qso_coord.match_to_catalog_sky(sdss_matched_coord)\n",
    "close_match = d2d < 1 * u.arcsec  # boolean mask for finding close matches\n",
    "print(f\"Found {np.sum(close_match)} close matches between QSO spectra and Gaia-SDSS sources\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now create a boolean column for the Gaia-SDSS sources to denote if they are a QSO"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "sdss_matched[\"known_qso\"] = False  # initially set as False for all\n",
    "sdss_matched[\"known_qso\"][idx[close_match]] = True  # and set to True for those indices with a close match\n",
    "sdss_matched"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 5. Investigate features of QSOs compared to other field sources to find QSO candidates\n",
    "\n",
    "We will now look at some properies of our sources, making use of the `\"known_qso\"` column to investigate what a QSO\n",
    "\"should\" look like in these properties.\n",
    "\n",
    "Firstly, we stack our two tables, `gaia_matched` and `sdss_matched`, in a column-wise manner\n",
    "([`hstack`](https://docs.astropy.org/en/stable/api/astropy.table.hstack.html#astropy.table.hstack)).\n",
    "This produces our final `sources` table"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "sources = table.hstack([gaia_matched, sdss_matched], table_names=[\"gaia\", \"sdss\"])\n",
    "qso_mask = sources[\"known_qso\"]  # boolean mask giving known QSO sources (we can use `~qso_mask` to mean `not qso_mask`)\n",
    "print(f\"There are {len(sources)} sources, with {len(sources[qso_mask])} known QSOs.\")\n",
    "sources.info\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Proper motions\n",
    "\n",
    "We hypothesised that QSOs could be found as zero proper motion sources within Gaia. We can determine some simple\n",
    "statistics on proper motions and plot them for known QSOs, and other sources."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "print(f\"Proper motions range from {sources[qso_mask]['pmtot'].min():.3f} to {sources[qso_mask]['pmtot'].max():.3f} mas/yr for known QSOs\")\n",
    "print(f\"Proper motions range from {sources[~qso_mask]['pmtot'].min():.3f} to {sources[~qso_mask]['pmtot'].max():.3f} mas/yr for other sources\")\n",
    "print(f\"The mean proper motion is {sources[qso_mask]['pmtot'].mean():.3f} mas/yr for known QSOs\")\n",
    "print(f\"The mean proper motion is {sources[~qso_mask]['pmtot'].mean():.3f} mas/yr for other sources\")\n",
    "print(f\"There are {np.sum(sources[qso_mask]['pm_significance'] < 2)} known QSOs with a < 2-sigma proper motion\")\n",
    "print(f\"There are {np.sum(sources[~qso_mask]['pm_significance'] < 2)} other sources with a < 2-sigma proper motion\")\n",
    "\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(1,3)\n",
    "\n",
    "ax0.errorbar(sources[~qso_mask][\"pmra\"], sources[~qso_mask][\"pmdec\"], xerr=sources[~qso_mask][\"pmra_error\"], yerr=sources[~qso_mask][\"pmdec_error\"], marker=\"o\", markersize=1, ls=\"\", alpha=0.1, label=\"other sources\")\n",
    "ax0.errorbar(sources[qso_mask][\"pmra\"], sources[qso_mask][\"pmdec\"], xerr=sources[qso_mask][\"pmra_error\"], yerr=sources[qso_mask][\"pmdec_error\"], marker=\"d\", markersize=3, ls=\"\", label=\"known QSOs\")\n",
    "ax0.set_xlabel(\"Proper Motion in RA [mas/yr]\")\n",
    "ax0.set_ylabel(\"Proper Motion in Dec [mas/yr]\")\n",
    "ax0.legend()\n",
    "\n",
    "ax1.errorbar(sources[qso_mask][\"pmra\"], sources[qso_mask][\"pmdec\"], xerr=sources[qso_mask][\"pmra_error\"], yerr=sources[qso_mask][\"pmdec_error\"], marker=\"d\", markersize=3, ls=\"\", label=\"known QSOs\", c=\"C1\")\n",
    "ax1.set_xlabel(\"Proper Motion in RA [mas/yr]\")\n",
    "ax1.set_ylabel(\"Proper Motion in Dec [mas/yr]\")\n",
    "ax1.legend()\n",
    "\n",
    "ax2.hist(sources[~qso_mask]['pm_significance'], bins=np.arange(41), label=\"other sources\")\n",
    "ax2.hist(sources[qso_mask]['pm_significance'], bins=np.arange(41), label=\"known QSOs\")\n",
    "ax2.set_xlabel(\"Proper motion significance (sigma)\")\n",
    "ax2.set_ylabel(\"N\")\n",
    "ax2.legend()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "It looks like QSOs, as expected, cluster around zero proper motion in Gaia! We also have a number of other\n",
    "sources which do the same, but are not known QSOs. We want to assess how this subset of sources look in other\n",
    "properties with respect to QSOs, in order to select from it targets for a hypothetical new QSO-confirmation spectroscopy\n",
    "proposal."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Colours\n",
    "\n",
    "Colours of sources can be revealing of their nature, albeit with significant overlap, particularly in optical colours.\n",
    "\n",
    "\n",
    "First we need to calculate colour columns for our sources. However, we need to be aware that SDSS can fill in dummy\n",
    "`-9999.0` values for its photometry, when the data is missing/bad. To avoid this we explicitly set those values to\n",
    "`np.nan` (i.e. Not A Number) so the are not used in any statistics, or plotting."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "for f in [\"u\", \"g\", \"r\", \"i\", \"z\"]:\n",
    "    sources[f\"psfMag_{f}\"][sources[f\"psfMag_{f}\"] == -9999.0] = np.nan"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "source": [
    "sources[\"u_g\"] = sources[\"psfMag_u\"] - sources[\"psfMag_g\"]\n",
    "sources[\"g_r\"] = sources[\"psfMag_g\"] - sources[\"psfMag_r\"]\n",
    "sources[\"r_i\"] = sources[\"psfMag_r\"] - sources[\"psfMag_i\"]\n",
    "sources[\"i_z\"] = sources[\"psfMag_i\"] - sources[\"psfMag_z\"]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "We also need a make a row mask to give us sources which are consistent with zero proper motion, but are **not** known QSOs:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "zero_pm = sources[\"pm_significance\"] < 2  # a boolean mask giving rows that are consistent with zero proper motion\n",
    "zero_pm_not_qso_mask = zero_pm & ~qso_mask  # combine with our \"not qso\" mask\n",
    "print(f\"There are {np.sum(zero_pm_not_qso_mask)} zero proper motion sources not classified as QSOs\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "And now we can plot colours for our three source types:\n",
    "* Known QSOs (`sources[qso_mask]`)\n",
    "* Sources consistent with zero proper motion that are not known to be QSOs (`sources[zero_pm_not_qso_mask]`)\n",
    "* Sources consistent with significant proper motion not known to be QSOs (`sources[~zero_pm_not_qso_mask]`)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 2)\n",
    "axes = axes.ravel()\n",
    "colour_columns = [[\"u_g\", \"g_r\"], [\"u_g\", \"r_i\"], [\"u_g\", \"i_z\"], [\"g_r\", \"r_i\"], [\"g_r\", \"i_z\"], [\"r_i\", \"i_z\"]]\n",
    "\n",
    "for i, (colour_column1, colour_column2) in enumerate(colour_columns):\n",
    "    ax = axes[i]\n",
    "    ax.scatter(sources[~zero_pm_not_qso_mask][colour_column1], sources[~zero_pm_not_qso_mask][colour_column2], s=5, alpha=0.1, c=\"C7\", label=\"Non-zero PM, not QSO\")\n",
    "    ax.scatter(sources[zero_pm_not_qso_mask][colour_column1], sources[zero_pm_not_qso_mask][colour_column2], s=10, c=\"C2\", label=\"Zero PM, not QSO\")\n",
    "    ax.scatter(sources[qso_mask][colour_column1], sources[qso_mask][colour_column2], s=10, c=\"C1\", marker=\"d\", label=\"Known QSOs\")\n",
    "    ax.set_xlabel(colour_column1)\n",
    "    ax.set_ylabel(colour_column2)\n",
    "    ax.set_xlim(-1, 4)\n",
    "    ax.set_ylim(-1, 4)\n",
    "    ax.legend()\n",
    "plt.tight_layout()\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Even by eye we can see that there is good agreement between a significant fraction of our ~zero proper motion sources\n",
    "and known QSOs in colour-colour space. One can imagine placing some cuts on colour to obtain a good list of QSO\n",
    "candidates.\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Write the sources data file here, to be used for cross-matching exercises\n",
    "# sources.write(\"data/sources_xm.csv\")\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "-----------\n",
    "\n",
    "## Example Exercises (cross-matching)\n",
    "\n",
    "You can pick an example exercise from here for the assessment if you wish. (Or perform your own project-relevant\n",
    "data-mining, as per the assessment instructions on the MPAGS page.)\n",
    "\n",
    "You can use the link to the csv file in the repo `sources_xm.csv` to import the source data easily into your program of\n",
    "choice. The sources are the result of this notebook's workings - sections 1 to 4 - using a 1 degree radius cone centred\n",
    "on RA=180 deg, Dec=30deg.\n",
    "\n",
    "### a. Cross match QSO candidates with WISE\n",
    "\n",
    "Cross-match any sources consistent with zero proper motion that are not known QSOs with the\n",
    "[WISE](https://irsa.ipac.caltech.edu/Missions/wise.html) mid-infrared survey. You can use any method or programme you\n",
    "like, or you can recreate somethign similar to the matching method above for Gaia and SDSS cross-matching.\n",
    "\n",
    "How many sources get a cross-match with WISE? Try varying your cross-matching radius (`close_match` above) since WISE\n",
    "has much poorer spatial resolution. What do you see in the WISE photometric colours of the cross-matched sources?\n",
    "See Section 3 and Figures 3 and 4 of [Heintz et al. 2018](https://arxiv.org/pdf/1805.03394.pdf) for inspiration.\n",
    "\n",
    "The catalogue within WISE to use is the AllWISE Source Catalogue. You could query it via the online gui\n",
    "[here](https://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-scan?mission=irsa&submit=Select&projshort=WISE) or use\n",
    "[Astroquery's Irsa](https://astroquery.readthedocs.io/en/latest/ipac/irsa/irsa.html) interface, similarly to how the SDSS\n",
    "interface is used in this notebook. Note to additionally pass `catalog=\"allwise_p3as_psd\"` argument to query the\n",
    "correct table.\n",
    "\n",
    " > For information, can obtain a full list of catalogs available via the Irsa interface in `astroquery` with:\n",
    " >\n",
    " > `from astroquery.ipac import irsa`\n",
    " > `irsa.Irsa.print_catalogs()`"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### b. Check other databases for variability and classification of QSO candidates\n",
    "\n",
    "QSOs are a subset of [Active Galactic Nuclei](https://en.wikipedia.org/wiki/Active_galactic_nucleus) and they should\n",
    "typically exhibit variability.\n",
    "\n",
    "Select zero proper-motion sources that are not known QSOs and apply some colour cuts to get a small subset of QSO\n",
    "candidates. Check other variability/transient databases (e.g. [Gaia Alerts](http://gsaweb.ast.cam.ac.uk/alerts/home) or\n",
    "the [Lasair](https://lasair.roe.ac.uk/) or [ALeRCE](https://alerce.online/) ZTF brokers) to try obtain light curves for\n",
    "any of them. Comment on the variability (or lack of) you find for any of the sources. Search other general data\n",
    "resources, such as [NED](https://ned.ipac.caltech.edu/) and [SIMBAD](http://simbad.u-strasbg.fr/simbad/), to find out\n",
    "more information on the objects.\n",
    "\n",
    "-----------"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 6. Train a Random Forest classifier to classify sources as QSOs based on proper-motion and colour\n",
    "\n",
    "We will train a (really simple) [Random Forest](https://en.wikipedia.org/wiki/Random_forest) classifier using\n",
    "[Scikit-learn](https://scikit-learn.org/stable/) that, for a given source, will predict a label (the probability of\n",
    "being a QSO) based on a set of features (its SDSS colours and proper motion).\n",
    "\n",
    "For this we will need larger numbers of confirmed QSO objects, since we need to show the RF the full extent of what a\n",
    "QSO can look like in our selected features. If we have only a few, as in the above working, our classifier will\n",
    "perform poorly. We can use this larger data set to train, and test a [`RandomForectClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html)\n",
    "from the Scikit-learn package.\n",
    "\n",
    "We will concisely perform a slimmed-down version of the above analysis, here selecting over a wider cone radius and\n",
    "cross-matching Gaia only to SDSS spectroscopic sources (such that we can robustly label all objects as QSO or not)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Define our new, larger cone\n",
    "ra_cone = 180 * u.degree\n",
    "dec_cone = 40 * u.degree\n",
    "coord = SkyCoord(ra=ra_cone, dec=dec_cone, frame='icrs')\n",
    "radius = 3 * u.degree\n",
    "region_radius = radius / np.cos(np.deg2rad(coord.dec.degree))  # correction for cos(dec) factor in ra"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Cone search in Gaia\n",
    "job = Gaia.cone_search_async(coord, radius)\n",
    "gaia_results = job.get_results()\n",
    "gaia_results = gaia_results[(gaia_results[\"phot_g_mean_mag\"] > 16) & (gaia_results[\"phot_g_mean_mag\"] < 20.5) & (gaia_results[\"astrometric_params_solved\"] == 31)]\n",
    "print(f\"Found {len(gaia_results)} Gaia sources\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Cone search in SDSS for spectroscopic objects\n",
    "specobj_fields=[\"class\", \"sciencePrimary\"]\n",
    "photoobj_fields=[\"ra\", \"dec\", \"type\", \"psfMag_u\", \"psfMag_g\", \"psfMag_r\", \"psfMag_i\", \"psfMag_z\"]\n",
    "sdss_spec_results = SDSS.query_region(coord, region_radius, spectro=True, specobj_fields=specobj_fields, photoobj_fields=photoobj_fields)\n",
    "sdss_spec_results = sdss_spec_results[(sdss_spec_results[\"sciencePrimary\"] == 1) & (sdss_spec_results[\"type\"] == 6)]\n",
    "print(f\"Found {len(sdss_spec_results)} SDSS spectroscopic sources, with {np.sum(sdss_spec_results['class'] == 'QSO')} QSOs\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Cross match Gaia and SDSS and stack into single sources table\n",
    "gaia_coord = SkyCoord(gaia_results[\"ra\"], gaia_results[\"dec\"], unit=u.degree)\n",
    "sdss_spec_coord = SkyCoord(sdss_spec_results[\"ra\"], sdss_spec_results[\"dec\"], unit=u.degree)\n",
    "idx, d2d, _ = gaia_coord.match_to_catalog_sky(sdss_spec_coord)\n",
    "close_match = d2d < 1 * u.arcsec\n",
    "gaia_matched = gaia_results[close_match]\n",
    "sdss_spec_matched = sdss_spec_results[idx[close_match]]\n",
    "sources = table.hstack([gaia_matched, sdss_spec_matched], table_names=[\"gaia\", \"sdss\"])\n",
    "print(f\"Found {len(sources)} cross-matches between Gaia and SDSS spectroscopic sources, with {np.sum(sources['class'] == 'QSO')} QSOs\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Add feature columns to sources\n",
    "pmtot = (sources[\"pmra\"]**2 + sources[\"pmdec\"]**2)**0.5\n",
    "pmtot_error = (sources[\"pmra_error\"]**2 + sources[\"pmdec_error\"]**2 + 2 * sources[\"pmra_pmdec_corr\"] * sources[\"pmra_error\"] * sources[\"pmdec_error\"])**0.5\n",
    "sources[\"pm_significance\"] = pmtot / pmtot_error\n",
    "sources[\"u_g\"] = sources[\"psfMag_u\"] - sources[\"psfMag_g\"]\n",
    "sources[\"g_r\"] = sources[\"psfMag_g\"] - sources[\"psfMag_r\"]\n",
    "sources[\"r_i\"] = sources[\"psfMag_r\"] - sources[\"psfMag_i\"]\n",
    "sources[\"i_z\"] = sources[\"psfMag_i\"] - sources[\"psfMag_z\"]\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Typically in Random Forests \"labels\" (`1` = QSO, `0`= not QSO) are denoted as `y`, with \"features\" (`pm_significance, u_g, g_r...`)\n",
    "being denoted as `X`."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "y = (sources[\"class\"] == \"QSO\").astype(np.int)  # convert the boolean condition into 1/0 labels\n",
    "# Now we have to do a slightly awkward way of making our features array,\n",
    "# since we cannot pass an astropy table - the input to sklearn must have a single dtype\n",
    "features = [\"pm_significance\", \"u_g\", \"g_r\", \"r_i\", \"i_z\", \"phot_g_mean_mag\"]\n",
    "X = np.empty((len(y), len(features)))  # just makes an array of shape (y, n_features), populated with random (empty) values\n",
    "for i, feat in enumerate(features):\n",
    "    X[:, i] = sources[feat] # fills the empty array with real values"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to split these into two subsets. We **train** the model on one subset, and **test** that it performs well on data it\n",
    "has never seen before. This can help preventing overfitting in the model (i.e. when model memorises its training data\n",
    "and is not generalisable to the problem at hand).\n",
    "\n",
    "Scikit learn has a function to do this for us, by default 25% will be reserved for testing."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0, stratify=y)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we train our model on the **train** subset:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "clf = RandomForestClassifier(random_state=0, n_jobs=1, class_weight=\"balanced_subsample\")\n",
    "model = clf.fit(X_train, y_train)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "And that's it! We have a trained random forest. But importantly we must now test the model with the **test** subset -\n",
    "data it has never seen before. We calculate some high-level assessment of how well it\n",
    "performs: defined as the accuracy (balanced for the mis-matched sample sizes), essentially \"what\n",
    "percentage of the time does the model correctly label unseen test data\"?"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "scores = cross_val_score(model, X_test, y_test, scoring=\"balanced_accuracy\", cv=10)\n",
    "print(f\"Random Forest gives a {scores.mean()*100:.1f} +/- {scores.std()*100:.1f}% accuracy\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Not bad! This *extremely* important step of testing is typically called [cross-validation](https://scikit-learn.org/stable/modules/cross_validation.html).\n",
    "\n",
    "The below code cells give hints on how to get the model to make a prediction, and also get the probabilities for each class,\n",
    "given a set of features. See the [`RandomForestClassifier`docs](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html) for information\n",
    "on how to use the model further."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "predicted_class = model.predict(X_test)  # remember 0 = not QSO, 1 = QSO\n",
    "# predicted_class is now an array of (2, len(X_test)), with the first and second columns being\n",
    "# the probabilities of class = 0 and class = 1, respectively.\n",
    "# We will grab the 2nd column ([:, 1]) of the probabilities to get the probability of\n",
    "# class = 1 (i.e. QSO) for our objects\n",
    "predicted_prob_qso = model.predict_proba(X_test)[:, 1]\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Now we can plot these probabilities, split by their y_test class value (0/1)\n",
    "# The y_test array holds the TRUE classes of the objects, so we expect the ones\n",
    "# where y_test == 0 to have low probabilities of being a QSO, and where y_test == 1\n",
    "# to have high probabilities.\n",
    "plt.hist(predicted_prob_qso[y_test == 0], histtype=\"step\", label=\"Not QSO\", bins=np.linspace(0,1,40), lw=2)\n",
    "plt.hist(predicted_prob_qso[y_test == 1], histtype=\"step\", label=\"QSO\", bins=np.linspace(0,1,40), lw=2)\n",
    "plt.xlabel(\"Model probability of being QSO\")\n",
    "plt.ylabel(\"N\")\n",
    "plt.legend()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Write the sources data file here, to be used for machine learning exercises\n",
    "# sources.write(\"data/sources_ml.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "-------\n",
    "\n",
    "## Example Exercises (machine learning)\n",
    "\n",
    "You can pick an example exercise from here for the assessment if you wish. (Or perform your own project-relevant\n",
    "data-mining, as per the assessment instructions on the MPAGS page.)\n",
    "\n",
    "You can use the link to the csv file in the repo `data/sources_ml.csv` to import the source data easily into your\n",
    "program of choice. The sources are the result of this notebook's workings - section 5 - using a 3 degree radius cone\n",
    "centred on RA=180 deg, Dec=40deg.\n",
    "\n",
    "### c. Use the model to find new QSOs\n",
    "\n",
    "Using the trained Random Forest model above, search another (small, ~0.5-1 deg radius) region of sky for Gaia-SDSS cross\n",
    "matches (getting only photometric sources) and run the model to\n",
    "[predict](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html#sklearn.ensemble.RandomForestClassifier.predict)\n",
    "the labels for your new sources. Manually inspect a few sources with high-probability of being a QSO and find out about\n",
    "them. Was the model right?\n",
    "\n",
    "What happens if you try a lower Galactic latitude SDSS-field? *(Be careful, the number of sources returned by searches\n",
    "over low Galactic latitude regions can shoot up very quickly!)*\n",
    "\n",
    "### d. Re-train the model using different features\n",
    "\n",
    "Explore removing or adding features for training the model. What is the impact of the model's accuracy? Are there any\n",
    "features that do no add any improvement in performance? Are there others in the SDSS or Gaia tables that do improve\n",
    "the model?\n",
    "\n",
    "For more advanced investigation, you could look at\n",
    "[feature selection](https://scikit-learn.org/stable/modules/feature_selection.html) in general for such\n",
    "feature-based models.\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}