docs/tutorials/ctapipe_overview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "# Analyzing Events Using ctapipe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div align=\"center\" style=\"font-size: 2rem\">\n",
    "\n",
    "<img heoght=\"300px\" src=\"https://cta-observatory.github.io/ctapipe/_images/ctapipe_logo.png\" alt=\"ctapipe\"/>\n",
    "\n",
    "\n",
    "<p style=\"text-align: center;\">Initially presented @ LST Analysis Bootcamp</p>\n",
    "\n",
    "<p style=\"text-align: center\">Padova, 26.11.2018</p>\n",
    "\n",
    "<p style=\"text-align: center\">Maximilian Nöthe (@maxnoe) & Kai A. Brügge (@mackaiver)</p>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (12, 8)\n",
    "plt.rcParams[\"font.size\"] = 14\n",
    "plt.rcParams[\"figure.figsize\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<h1 id=\"tocheading\">Table of Contents</h1>\n",
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## General Information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Design\n",
    "\n",
    "* DL0 → DL3 analysis\n",
    "\n",
    "* Currently some R0 → DL2 code to be able to analyze simtel files\n",
    "\n",
    "* ctapipe is built upon the Scientific Python Stack, core dependencies are\n",
    "  * numpy\n",
    "  * scipy\n",
    "  * astropy\n",
    "  * numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Developement\n",
    "\n",
    "* ctapipe is developed as Open Source Software (BSD 3-Clause License) at <https://github.com/cta-observatory/ctapipe>\n",
    "\n",
    "* We use the \"Github-Workflow\": \n",
    "  * Few people (e.g. @kosack, @maxnoe) have write access to the main repository\n",
    "  * Contributors fork the main repository and work on branches\n",
    "  * Pull Requests are merged after Code Review and automatic execution of the test suite\n",
    "\n",
    "* Early developement stage ⇒ backwards-incompatible API changes might and will happen \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What's there?\n",
    "\n",
    "* Reading simtel simulation files\n",
    "* Simple calibration, cleaning and feature extraction functions\n",
    "* Camera and Array plotting\n",
    "* Coordinate frames and transformations \n",
    "* Stereo-reconstruction using line intersections\n",
    "  \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What's still missing?\n",
    "\n",
    "* Good integration with machine learning techniques\n",
    "* IRF calculation \n",
    "* Documentation, e.g. formal definitions of coordinate frames  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What can you do?\n",
    "\n",
    "* Report issues\n",
    "  * Hard to get started? Tell us where you are stuck\n",
    "  * Tell user stories\n",
    "  * Missing features\n",
    "\n",
    "* Start contributing\n",
    "  * ctapipe needs more workpower\n",
    "  * Implement new reconstruction features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {},
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A simple hillas analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Reading in simtel files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.io import EventSource\n",
    "from ctapipe.utils.datasets import get_dataset_path\n",
    "\n",
    "input_url = get_dataset_path(\"gamma_prod5.simtel.zst\")\n",
    "\n",
    "# EventSource() automatically detects what kind of file we are giving it,\n",
    "# if already supported by ctapipe\n",
    "source = EventSource(input_url, max_events=5)\n",
    "\n",
    "print(type(source))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "for event in source:\n",
    "    print(\n",
    "        \"Id: {}, E = {:1.3f}, Telescopes: {}\".format(\n",
    "            event.count, event.simulation.shower.energy, len(event.r0.tel)\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "Each event is a `DataContainer` holding several `Field`s of data, which can be containers or just numbers.\n",
    "Let's look a one event:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "source.subarray.camera_types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "len(event.r0.tel), len(event.r1.tel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Data calibration\n",
    "\n",
    "The `CameraCalibrator` calibrates the event (obtaining the `dl1` images)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.calib import CameraCalibrator\n",
    "\n",
    "calibrator = CameraCalibrator(subarray=source.subarray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "calibrator(event)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Event displays\n",
    "\n",
    "Let's use ctapipe's plotting facilities to plot the telescope images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "event.dl1.tel.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "tel_id = 130"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "geometry = source.subarray.tel[tel_id].camera.geometry\n",
    "dl1 = event.dl1.tel[tel_id]\n",
    "\n",
    "geometry, dl1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "dl1.image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.visualization import CameraDisplay\n",
    "\n",
    "display = CameraDisplay(geometry)\n",
    "\n",
    "# right now, there might be one image per gain channel.\n",
    "# This will change as soon as\n",
    "display.image = dl1.image\n",
    "display.add_colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Image Cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.image.cleaning import tailcuts_clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "# unoptimized cleaning levels\n",
    "cleaning_level = {\n",
    "    \"CHEC\": (2, 4, 2),\n",
    "    \"LSTCam\": (3.5, 7, 2),\n",
    "    \"FlashCam\": (3.5, 7, 2),\n",
    "    \"NectarCam\": (4, 8, 2),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "boundary, picture, min_neighbors = cleaning_level[geometry.name]\n",
    "\n",
    "clean = tailcuts_clean(\n",
    "    geometry,\n",
    "    dl1.image,\n",
    "    boundary_thresh=boundary,\n",
    "    picture_thresh=picture,\n",
    "    min_number_picture_neighbors=min_neighbors,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "d1 = CameraDisplay(geometry, ax=ax1)\n",
    "d2 = CameraDisplay(geometry, ax=ax2)\n",
    "\n",
    "ax1.set_title(\"Image\")\n",
    "d1.image = dl1.image\n",
    "d1.add_colorbar(ax=ax1)\n",
    "\n",
    "ax2.set_title(\"Pulse Time\")\n",
    "d2.image = dl1.peak_time - np.average(dl1.peak_time, weights=dl1.image)\n",
    "d2.cmap = \"RdBu_r\"\n",
    "d2.add_colorbar(ax=ax2)\n",
    "d2.set_limits_minmax(-20, 20)\n",
    "\n",
    "d1.highlight_pixels(clean, color=\"red\", linewidth=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Image Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.image import (\n",
    "    hillas_parameters,\n",
    "    leakage_parameters,\n",
    "    concentration_parameters,\n",
    ")\n",
    "from ctapipe.image import timing_parameters\n",
    "from ctapipe.image import number_of_islands\n",
    "from ctapipe.image import camera_to_shower_coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "hillas = hillas_parameters(geometry[clean], dl1.image[clean])\n",
    "\n",
    "print(hillas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "display = CameraDisplay(geometry)\n",
    "\n",
    "# set \"unclean\" pixels to 0\n",
    "cleaned = dl1.image.copy()\n",
    "cleaned[~clean] = 0.0\n",
    "\n",
    "display.image = cleaned\n",
    "display.add_colorbar()\n",
    "\n",
    "display.overlay_moments(hillas, color=\"xkcd:red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "timing = timing_parameters(geometry, dl1.image, dl1.peak_time, hillas, clean)\n",
    "\n",
    "print(timing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "long, trans = camera_to_shower_coordinates(\n",
    "    geometry.pix_x, geometry.pix_y, hillas.x, hillas.y, hillas.psi\n",
    ")\n",
    "\n",
    "plt.plot(long[clean], dl1.peak_time[clean], \"o\")\n",
    "plt.plot(long[clean], timing.slope * long[clean] + timing.intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "l = leakage_parameters(geometry, dl1.image, clean)\n",
    "print(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "disp = CameraDisplay(geometry)\n",
    "disp.image = dl1.image\n",
    "disp.highlight_pixels(geometry.get_border_pixel_mask(1), linewidth=2, color=\"xkcd:red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "n_islands, island_id = number_of_islands(geometry, clean)\n",
    "\n",
    "print(n_islands)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "conc = concentration_parameters(geometry, dl1.image, hillas)\n",
    "print(conc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Putting it all together / Stereo reconstruction\n",
    "\n",
    "\n",
    "All these steps are now unified in several components configurable through the config system, mainly:\n",
    "\n",
    "* CameraCalibrator for DL0 → DL1 (Images)\n",
    "* ImageProcessor for DL1 (Images) → DL1 (Parameters)\n",
    "* ShowerProcessor for stereo reconstruction of the shower geometry\n",
    "* DataWriter for writing data into HDF5\n",
    "\n",
    "A command line tool doing these steps and writing out data in HDF5 format is available as `ctapipe-process`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {},
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord, AltAz\n",
    "\n",
    "from ctapipe.containers import ImageParametersContainer\n",
    "from ctapipe.io import EventSource\n",
    "from ctapipe.utils.datasets import get_dataset_path\n",
    "\n",
    "from ctapipe.calib import CameraCalibrator\n",
    "\n",
    "from ctapipe.image import ImageProcessor\n",
    "from ctapipe.reco import ShowerProcessor\n",
    "\n",
    "from ctapipe.io import DataWriter\n",
    "\n",
    "from copy import deepcopy\n",
    "import tempfile\n",
    "\n",
    "from traitlets.config import Config\n",
    "\n",
    "\n",
    "image_processor_config = Config(\n",
    "    {\n",
    "        \"ImageProcessor\": {\n",
    "            \"image_cleaner_type\": \"TailcutsImageCleaner\",\n",
    "            \"TailcutsImageCleaner\": {\n",
    "                \"picture_threshold_pe\": [\n",
    "                    (\"type\", \"LST_LST_LSTCam\", 7.5),\n",
    "                    (\"type\", \"MST_MST_FlashCam\", 8),\n",
    "                    (\"type\", \"MST_MST_NectarCam\", 8),\n",
    "                    (\"type\", \"SST_ASTRI_CHEC\", 7),\n",
    "                ],\n",
    "                \"boundary_threshold_pe\": [\n",
    "                    (\"type\", \"LST_LST_LSTCam\", 5),\n",
    "                    (\"type\", \"MST_MST_FlashCam\", 4),\n",
    "                    (\"type\", \"MST_MST_NectarCam\", 4),\n",
    "                    (\"type\", \"SST_ASTRI_CHEC\", 4),\n",
    "                ],\n",
    "            },\n",
    "        }\n",
    "    }\n",
    ")\n",
    "\n",
    "input_url = get_dataset_path(\"gamma_prod5.simtel.zst\")\n",
    "source = EventSource(input_url)\n",
    "\n",
    "calibrator = CameraCalibrator(subarray=source.subarray)\n",
    "image_processor = ImageProcessor(\n",
    "    subarray=source.subarray, config=image_processor_config\n",
    ")\n",
    "shower_processor = ShowerProcessor(subarray=source.subarray)\n",
    "horizon_frame = AltAz()\n",
    "\n",
    "f = tempfile.NamedTemporaryFile(suffix=\".hdf5\")\n",
    "\n",
    "with DataWriter(\n",
    "    source, output_path=f.name, overwrite=True, write_showers=True\n",
    ") as writer:\n",
    "\n",
    "    for event in source:\n",
    "        energy = event.simulation.shower.energy\n",
    "        n_telescopes_r1 = len(event.r1.tel)\n",
    "        event_id = event.index.event_id\n",
    "        print(f\"Id: {event_id}, E = {energy:1.3f}, Telescopes (R1): {n_telescopes_r1}\")\n",
    "\n",
    "        calibrator(event)\n",
    "        image_processor(event)\n",
    "        shower_processor(event)\n",
    "\n",
    "        stereo = event.dl2.stereo.geometry[\"HillasReconstructor\"]\n",
    "        if stereo.is_valid:\n",
    "            print(\"  Alt: {:.2f}°\".format(stereo.alt.deg))\n",
    "            print(\"  Az: {:.2f}°\".format(stereo.az.deg))\n",
    "            print(\"  Hmax: {:.0f}\".format(stereo.h_max))\n",
    "            print(\"  CoreX: {:.1f}\".format(stereo.core_x))\n",
    "            print(\"  CoreY: {:.1f}\".format(stereo.core_y))\n",
    "            print(\"  Multiplicity: {:d}\".format(len(stereo.telescopes)))\n",
    "\n",
    "        # save a nice event for plotting later\n",
    "        if event.count == 3:\n",
    "            plotting_event = deepcopy(event)\n",
    "\n",
    "        writer(event)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from astropy.coordinates.angle_utilities import angular_separation\n",
    "import pandas as pd\n",
    "\n",
    "from ctapipe.io import TableLoader\n",
    "\n",
    "loader = TableLoader(f.name, load_dl2=True, load_simulated=True)\n",
    "\n",
    "events = loader.read_subarray_events()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "theta = angular_separation(\n",
    "    events[\"HillasReconstructor_az\"].quantity,\n",
    "    events[\"HillasReconstructor_alt\"].quantity,\n",
    "    events[\"true_az\"].quantity,\n",
    "    events[\"true_alt\"].quantity,\n",
    ")\n",
    "\n",
    "plt.hist(theta.to_value(u.deg) ** 2, bins=25, range=[0, 0.3])\n",
    "plt.xlabel(r\"$\\theta² / deg²$\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "## ArrayDisplay\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.visualization import ArrayDisplay\n",
    "\n",
    "\n",
    "angle_offset = plotting_event.pointing.array_azimuth\n",
    "\n",
    "plotting_hillas = {\n",
    "    tel_id: dl1.parameters.hillas for tel_id, dl1 in plotting_event.dl1.tel.items()\n",
    "}\n",
    "\n",
    "plotting_core = {\n",
    "    tel_id: dl1.parameters.core.psi for tel_id, dl1 in plotting_event.dl1.tel.items()\n",
    "}\n",
    "\n",
    "\n",
    "disp = ArrayDisplay(source.subarray)\n",
    "\n",
    "disp.set_line_hillas(plotting_hillas, plotting_core, 500)\n",
    "\n",
    "plt.scatter(\n",
    "    plotting_event.simulation.shower.core_x,\n",
    "    plotting_event.simulation.shower.core_y,\n",
    "    s=200,\n",
    "    c=\"k\",\n",
    "    marker=\"x\",\n",
    "    label=\"True Impact\",\n",
    ")\n",
    "plt.scatter(\n",
    "    plotting_event.dl2.stereo.geometry[\"HillasReconstructor\"].core_x,\n",
    "    plotting_event.dl2.stereo.geometry[\"HillasReconstructor\"].core_y,\n",
    "    s=200,\n",
    "    c=\"r\",\n",
    "    marker=\"x\",\n",
    "    label=\"Estimated Impact\",\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "# plt.xlim(-400, 400)\n",
    "# plt.ylim(-400, 400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "### Reading the LST dl1 data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loader = TableLoader(f.name, load_simulated=True, load_dl1_parameters=True)\n",
    "\n",
    "dl1_table = loader.read_telescope_events([\"LST_LST_LSTCam\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "plt.scatter(\n",
    "    np.log10(dl1_table[\"true_energy\"].quantity / u.TeV),\n",
    "    np.log10(dl1_table[\"hillas_intensity\"]),\n",
    ")\n",
    "plt.xlabel(\"log10(E / TeV)\")\n",
    "plt.ylabel(\"log10(intensity)\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "## Isn't python slow?\n",
    "\n",
    "* Many of you might have heard: \"Python is slow\".\n",
    "* That's trueish.\n",
    "* All python objects are classes living on the heap, even integers.\n",
    "* Looping over lots of \"primitives\" is quite slow compared to other languages.\n",
    "\n",
    "⇒ Vectorize as much as possible using numpy  \n",
    "⇒ Use existing interfaces to fast C / C++ / Fortran code  \n",
    "⇒ Optimize using numba  \n",
    "\n",
    "**But: \"Premature Optimization is the root of all evil\" — Donald Knuth**\n",
    "\n",
    "So profile to find exactly what is slow.\n",
    "\n",
    "### Why use python then?\n",
    "\n",
    "* Python works very well as *glue* for libraries of all kinds of languages\n",
    "* Python has a rich ecosystem for data science, physics, algorithms, astronomy\n",
    "\n",
    "### Example: Number of Islands\n",
    "\n",
    "Find all groups of pixels, that survived the cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from ctapipe.image import toymodel\n",
    "from ctapipe.instrument import SubarrayDescription\n",
    "\n",
    "\n",
    "geometry = loader.subarray.tel[1].camera.geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "Let's create a toy images with several islands;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "image = np.zeros(geometry.n_pixels)\n",
    "\n",
    "\n",
    "for i in range(9):\n",
    "\n",
    "    model = toymodel.Gaussian(\n",
    "        x=np.random.uniform(-0.8, 0.8) * u.m,\n",
    "        y=np.random.uniform(-0.8, 0.8) * u.m,\n",
    "        width=np.random.uniform(0.05, 0.075) * u.m,\n",
    "        length=np.random.uniform(0.1, 0.15) * u.m,\n",
    "        psi=np.random.uniform(0, 2 * np.pi) * u.rad,\n",
    "    )\n",
    "\n",
    "    new_image, sig, bg = model.generate_image(\n",
    "        geometry, intensity=np.random.uniform(1000, 3000), nsb_level_pe=5\n",
    "    )\n",
    "    image += new_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "clean = tailcuts_clean(\n",
    "    geometry,\n",
    "    image,\n",
    "    picture_thresh=10,\n",
    "    boundary_thresh=5,\n",
    "    min_number_picture_neighbors=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "disp = CameraDisplay(geometry)\n",
    "disp.image = image\n",
    "disp.highlight_pixels(clean, color=\"xkcd:red\", linewidth=1.5)\n",
    "disp.add_colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "def num_islands_python(camera, clean):\n",
    "    \"\"\"A breadth first search to find connected islands of neighboring pixels in the cleaning set\"\"\"\n",
    "\n",
    "    # the camera geometry has a [n_pixel, n_pixel] boolean array\n",
    "    # that is True where two pixels are neighbors\n",
    "    neighbors = camera.neighbor_matrix\n",
    "\n",
    "    island_ids = np.zeros(camera.n_pixels)\n",
    "    current_island = 0\n",
    "\n",
    "    # a set to remember which pixels we already visited\n",
    "    visited = set()\n",
    "\n",
    "    # go only through the pixels, that survived cleaning\n",
    "    for pix_id in np.where(clean)[0]:\n",
    "        if pix_id not in visited:\n",
    "            # remember that we already checked this pixel\n",
    "            visited.add(pix_id)\n",
    "\n",
    "            # if we land in the outer loop again, we found a new island\n",
    "            current_island += 1\n",
    "            island_ids[pix_id] = current_island\n",
    "\n",
    "            # now check all neighbors of the current pixel recursively\n",
    "            to_check = set(np.where(neighbors[pix_id] & clean)[0])\n",
    "            while to_check:\n",
    "                pix_id = to_check.pop()\n",
    "\n",
    "                if pix_id not in visited:\n",
    "                    visited.add(pix_id)\n",
    "                    island_ids[pix_id] = current_island\n",
    "\n",
    "                    to_check.update(np.where(neighbors[pix_id] & clean)[0])\n",
    "\n",
    "    n_islands = current_island\n",
    "    return n_islands, island_ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "n_islands, island_ids = num_islands_python(geometry, clean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "cmap = plt.get_cmap(\"Paired\")\n",
    "cmap = ListedColormap(cmap.colors[:n_islands])\n",
    "cmap.set_under(\"k\")\n",
    "\n",
    "disp = CameraDisplay(geometry)\n",
    "disp.image = island_ids\n",
    "disp.cmap = cmap\n",
    "disp.set_limits_minmax(0.5, n_islands + 0.5)\n",
    "disp.add_colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "%timeit num_islands_python(geometry, clean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "from scipy.sparse.csgraph import connected_components\n",
    "\n",
    "\n",
    "def num_islands_scipy(geometry, clean):\n",
    "    neighbors = geometry.neighbor_matrix_sparse\n",
    "\n",
    "    clean_neighbors = neighbors[clean][:, clean]\n",
    "    num_islands, labels = connected_components(clean_neighbors, directed=False)\n",
    "\n",
    "    island_ids = np.zeros(geometry.n_pixels)\n",
    "    island_ids[clean] = labels + 1\n",
    "\n",
    "    return num_islands, island_ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "n_islands_s, island_ids_s = num_islands_scipy(geometry, clean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "disp = CameraDisplay(geometry)\n",
    "disp.image = island_ids_s\n",
    "disp.cmap = cmap\n",
    "disp.set_limits_minmax(0.5, n_islands_s + 0.5)\n",
    "disp.add_colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "%timeit num_islands_scipy(geometry, clean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "**A lot less code, and a factor 3 speed improvement**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, current ctapipe implementation is using numba:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit number_of_islands(geometry, clean)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}