Data_Validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "35799286-1fa3-40ef-81d9-49afe47c79fe",
   "metadata": {},
   "source": [
    "# Data Validation Tutorial\n",
    "\n",
    "Have you ever seen a cosmology group unblind their data and then be left confused by the results? Underestimate a systematic effect? Needing to re-calibrate the data or fix uncovered issues? Well these things happen! The data has more statistical power than ever before and it's hard to foresee all issues that might come up in modern surveys. Running quality checks and validations as early as possible by as many groups as possible helps to minimize the risk to your science, but truly understanding the data is a never-ending goal."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4498668f-6c83-4c33-8aa9-064403f96767",
   "metadata": {},
   "source": [
    "<div markdown=\"1\" style=\"background-color:rgba(200, 230, 250, 0.4); text-align:left; vertical-align: middle; padding:10px; width:900px;\">\n",
    "<p>\n",
    "\n",
    "### Background: DESCQA / analysis_tools\n",
    "\n",
    "[DESCQA](https://github.com/LSSTDESC/descqa) is a DESC project for running validation tests on simulated catalogs. These tests were developed primarily for the cosmoDC2 catalog and go all the way from looking at whether there are outliers in the data to reproducing cosmology plots. When doing early-stage validation for the real data we want to be a little more careful to avoid peeking at the cosmology ahead of time or replicating analysis pipelines, and instead focus on understanding the quality of the data. For this array of early data tests we have an SRV branch.\n",
    "\n",
    "[analysis_tools](https://pipelines.lsst.io/v/daily/modules/lsst.analysis.tools/index.html) is the project's QA software. There's a lot of great functionality and pretty plots in there, it is installed anywhere with the lsst pipelines (There is a great description of the rho statistics test on the RSP [here](https://github.com/rubin-dp0/tutorial-notebooks/blob/main/DP02_12b_PSF_Science_Demo.ipynb) and I also have some notes on how to run analysis_tools at NERSC [here](https://lsstdesc.org/srv-dp02/tutorials/Rubin/tutorial3_Rubin_analysistools.html)) and there's a bunch of tests that you can see run by the operations team nightly here https://sitcomtn-124.lsst.io/. If you want to use or contribute to analysis_tools we are more than happy for tests to be developed in this format and run either by DESC or the project depending on the scope. \n",
    "\n",
    "\n",
    "</p></span>\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7160c701-2436-4af3-bfeb-e46dc53f1ecf",
   "metadata": {},
   "source": [
    "## Set-up\n",
    "This is much the same as the last tutorial except this time let's also add a descqa import\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34348e2e-fb7c-4050-be94-adc97e33ec55",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import time\n",
    "sys.path.insert(0,'/global/cfs/projectdirs/lsst/groups/SRV/gcr-catalogs-test')\n",
    "import GCRCatalogs\n",
    "sys.path.insert(0,'/global/cfs/projectdirs/lsst/groups/SRV/descqa')\n",
    "import descqa\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8880131b-356e-44b9-97c1-e69f9d17b012",
   "metadata": {},
   "source": [
    "## DESCQA - DESC-Quality Assurance Framework\n",
    "\n",
    "https://arxiv.org/abs/1709.09665\n",
    "https://arxiv.org/abs/2110.03769\n",
    "\n",
    "Designed to allow us to validate and inspect the results of a range of catalogs, read in by the generic catalog reader. We're expanding this in the SRV branch to scale up in data size and run validation on the DP0.2 data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65606acf-2868-4a24-b6f4-768e508c2dd0",
   "metadata": {},
   "source": [
    "<img src=\"./pictures/descqa.png\"  width=\"1100\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "922b28f0-6809-4e39-a89b-9766e637e8c8",
   "metadata": {
    "tags": []
   },
   "source": [
    "This has traditionally been run on the terminal at NERSC, this can be as easy as cloning the repository and within the main folder typing \n",
    "```bash\n",
    "./run_master.sh -c catalogs_to_run -t test_to_run\n",
    "```\n",
    "The test will then run, and output a url linking to the output of the test in the web browser. \n",
    "\n",
    "I'm going to show you how to run this in a jupyter notebook, as many people seem more comfortable developing interactively. To do this similarly the only two things we really need to know are which catalog to load and which test to run. Let's run one I made ahead of time for this called srv_tutorial on the catalog we currently have in hand, lsst_object. \n",
    "\n",
    "We're missing a bunch of exception catching here because I know these exist, if you want to see what those checks are we would typically run have a look [here](https://lsstdesc.org/srv-dp02/tutorials/Rubin/tutorial4_using_DESCQA.html)   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79c6a201-26cd-4f69-b43a-2e9dbdc3e4f0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "catalog = 'lsst_object'\n",
    "validation = 'srv_tutorial'\n",
    "\n",
    "catalog_instance =  GCRCatalogs.load_catalog(catalog)   \n",
    "validation_instance = descqa.load_validation(validation)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "081efffd-9a64-4f8c-8646-e33a3ea74928",
   "metadata": {},
   "source": [
    "We need to know where to output the result - for a proper terminal-based DESCQA run it would create an output directory in the DESC project space for you and show you where this is, but for developing I've created a sub-directory of this repo, and we can write the output to there. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3ceac3d-6464-4036-a371-e9f2d2f102b9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "output_dir_this = './test_direc'  # point to a directory where you want to store any plots or results\n",
    "test_result = validation_instance.run_on_single_catalog(catalog_instance, catalog, output_dir_this)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd754dbc-d150-4a20-992e-f3a700c4dbfd",
   "metadata": {},
   "source": [
    "And let's use the interpretation script I've made earlier to inspect the returned test result (look in utils.py if you're curious)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a0fa1df-2aab-4379-a848-6714a8b1f07a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from utils import interpret_result\n",
    "interpret_result(test_result)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3781e801-06c8-45af-b4b6-29f8acacaeb7",
   "metadata": {},
   "source": [
    "We can then look at any plots output into the directory (You can directly look through the file system or use a display within the notebook like I have here)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "825960c2-0e24-4699-be75-398a74ebad18",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import glob\n",
    "from PIL import Image\n",
    "from IPython.display import display\n",
    "\n",
    "for image in glob.glob( output_dir_this + '/*.png'):\n",
    "    print(image)\n",
    "    display(Image.open(image))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b27da7-3f1f-492c-bcf8-2c4b3011a60a",
   "metadata": {},
   "source": [
    "<div markdown=\"1\" style=\"background-color:rgba(200, 230, 250, 0.4); text-align:left; vertical-align: middle; padding:10px; width:900px;\">\n",
    "<p>\n",
    "\n",
    "### Background: SRV tests\n",
    "\n",
    "If you check out the SRV branch of DESCQA you'll see it's actually a little more complicated than that. In the run_master script you have options to do parallel reads, change the kernel to one with the lsst pipelines and set the environment. To run on the log-in node on Perlmutter you cannot import MPI or work with large datasets so for now we're developing things in a more simple environment. If you'd like to contribute a test and make it work in parallel, with GPUs or with lsst software, just let me know!\n",
    "\n",
    "</p></span>\n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adb37c9d-a09b-493b-8ab4-21c5e8f4a7e9",
   "metadata": {},
   "source": [
    "## How do we develop these tests?\n",
    "\n",
    "You need roughly three things to develop a test \n",
    "* Input arguments \n",
    "* Some operation to perform on the data \n",
    "* An output to display the result\n",
    "\n",
    "This is the exact code for what we just ran, except this time I thought we'd make it prettier by changing the color map and adding tight_layout to fix that pesky y-label issue. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc61286f-13f3-4a1c-8453-fcdee3625458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from descqa.base import BaseValidationTest, TestResult\n",
    "import numpy as np\n",
    "import healpy as hp\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "\n",
    "class Tutorial2(BaseValidationTest):\n",
    "    \"\"\"\n",
    "    validation test to plot RA/DEC\n",
    "    \"\"\"\n",
    "    def __init__(self, **kwargs): # initialize and get the inputs \n",
    "        self.ra = kwargs.get('ra')\n",
    "        self.dec = kwargs.get('dec')\n",
    "        self.catalog_filters = kwargs.get('catalog_filters', [])\n",
    "\n",
    "\n",
    "    def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir): # run the operation\n",
    "\n",
    "        if not catalog_instance.has_quantities([self.ra, self.dec]):\n",
    "            return TestResult(skipped=True)\n",
    "\n",
    "        # read in the data \n",
    "        data = catalog_instance.get_quantities([self.ra,self.dec],filters=self.catalog_filters, return_iterator=False)\n",
    "\n",
    "        xbins = np.linspace(np.min(data[self.ra]),np.max(data[self.ra]),50)\n",
    "        ybins = np.linspace(np.min(data[self.dec]),np.max(data[self.dec]),50)\n",
    "        plt.figure()\n",
    "        plt.hist2d(data[self.ra],data[self.dec], bins=(xbins,ybins),cmap='Blues')\n",
    "        plt.colorbar()\n",
    "        plt.xlabel('RA (degrees)')\n",
    "        plt.ylabel('dec (degrees)')\n",
    "        plt.tight_layout()\n",
    "        plt.savefig(os.path.join(output_dir, 'test_histogram.png')) # create output #1 \n",
    "        plt.close()\n",
    "\n",
    "        return TestResult(inspect_only=True, summary='made a histogram to inspect!') # add output information\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89c94ff8-609c-4a00-af0b-84490ea01c02",
   "metadata": {},
   "source": [
    "We need to create the inputs to pass in, we're just going to create a dictionary with the right key words, then all we need to do is initialize a validation with this class and the inputs and we're ready to run the test!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e946911c-bdd0-4ccc-b6a8-77e0e123a7fc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dict_config = {}\n",
    "dict_config['ra'] = 'coord_ra'\n",
    "dict_config['dec'] = 'coord_dec'\n",
    "dict_config['catalog_filters'] =  ['refExtendedness==1','detect_isPrimary']\n",
    "\n",
    "validation2 = Tutorial2(**dict_config)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bb18ae8-321c-45d2-9e93-537d95928ca0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "test_result = validation2.run_on_single_catalog(catalog_instance, catalog, output_dir_this)\n",
    "\n",
    "for image in glob.glob( output_dir_this + '/*.png'):\n",
    "    print(image)\n",
    "    display(Image.open(image))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b20097e-a0d1-4949-9f90-5f826e90d2fa",
   "metadata": {},
   "source": [
    "## Task! Make a validation plot\n",
    "\n",
    "Compute a plot of the position-dependent PSF size (we'll define this as the trace we computed in the last notebook)\n",
    "* Make the band an input argument, and run plots for the g-band and the i-band\n",
    "\n",
    "Hints \n",
    "* scipy.binned_statistic_2d is the function I like to use for this\n",
    "* You can use the code for the data read you used last challenge\n",
    "* If you were reading large amounts of data and don't want to repeat the read, just pass the data in through your input dictionary\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2f25810-152c-4947-97e1-dc4cd4938040",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-python",
   "language": "python",
   "name": "desc-python"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}