presentation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "from IPython.display import YouTubeVideo, Markdown, Code\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# NumPy: A look at the past, present, and future of array computation\n",
    "\n",
    "[Ross Barnowski](https://bids.berkeley.edu/people/ross-barnowski) `rossbar@berkeley.edu` | [@rossbar](https://github.com/rossbar) on GitHub\n",
    "\n",
    "University of Michigan EECS | 1/30/2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Part I - Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What is NumPy?\n",
    "\n",
    "> *NumPy is the fundamental package for scientific computing with Python*\n",
    "> \n",
    ">  [numpy.org](https://numpy.org/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Strong stuff."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<center>\n",
    "\n",
    "## The scientific Python ecosystem\n",
    "\n",
    "<img src=\"images/state_of_the_stack_2015.png\" alt=\"scientific_python_ecosystem\"/>\n",
    "\n",
    "</center>\n",
    "\n",
    "Image credit: [Jake VanderPlas](http://vanderplas.com/) **circa 2015**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A bit of history\n",
    "\n",
    " - **Mid 90's/Early 00's**: desire for high-performance numerical computation in Python eventually culminates in the [`Numeric`](https://numpy.org/_downloads/768fa66c250a0335ad3a6a30fae48e34/numeric-manual.pdf) library\n",
    " - Early adopters included the [Space Telescope Science Institute (STScI)](http://www.stsci.edu/) who developed another array computation package to better suit their needs: `NumArray`.\n",
    " - **2005** The best ideas from `Numeric` and `NumArray` were combined in the development of a new library, `NumPy`\n",
    "   * This work was largely done by [Travis Oliphant](https://github.com/teoliphant), then an assistant professor at BYU\n",
    " - **2006** NumPy v1.0 released in October\n",
    " \n",
    "[NumPy Development History](https://github.com/numpy/numpy/graphs/contributors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What does NumPy provide?\n",
    "\n",
    " - `ndarray`: A generic, n-dimensional array data structure\n",
    " - Sophisticated machinery for operating on array data\n",
    "   * Powerful indexing\n",
    "   * Built-in, array-aware operations\n",
    "   * Vectorization and broadcasting\n",
    " - All features exposed by a concise, expressive syntax\n",
    " - Language extension/integration (C-API, `f2py`) and interoperability\n",
    "   * [Array API](https://numpy.org/doc/1.17/reference/c-api.array.html) for accessing/extending array functionality\n",
    "   * Protocols for replicating the NumPy interface (stay tuned...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## What else?\n",
    "\n",
    "`numpy` also includes tools for common scientific/numerical tasks:\n",
    "   * Random number generation (`np.random`)\n",
    "   * Fourier analysis (`np.fft`)\n",
    "   * Linear algebra (`np.linalg`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The `scipy` package includes modules with the same name? What's the deal?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy, scipy.linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(scipy.random) # scipy.stats\n",
    "scipy.random is np.random "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(scipy.fft)\n",
    "scipy.fft is np.fft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(scipy.linalg)\n",
    "scipy.linalg is np.linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## A useful analogy...\n",
    "\n",
    "<center><img src=\"images/tool_analogy.png\" alt=\"socket set analogy\"/></center>\n",
    "\n",
    " - E.g. see [this quick comparison](https://numpy.org/devdocs/reference/routines.linalg.html) of the `numpy` and `scipy` `linalg` modules."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Where is NumPy used?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# At a glance\n",
    "\n",
    "<center><img src=\"images/NumPy_info3.jpg\" alt=\"NumPy Overview Infographic\" width=1152 height=658/></center>\n",
    "\n",
    "**See also**: [Anaconda 2019 Year in Review](https://www.anaconda.com/2019-year-in-review/),  [NumPy PyPI Stats page](https://pypistats.org/packages/numpy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Part II - NumPy in the wild"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Code example: github graphql query for top starred projects with numpy as a dependency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Neuroimaging Analysis\n",
    "\n",
    "Like much of the scientific python ecosystem, [nipy](https://nipy.org/) relies on `np.ndarray` as the fundamental structure for neuroimaging data.\n",
    "\n",
    "The following example is adapted from [Machine learning for neuroimaging with scikit learn](https://www.frontiersin.org/articles/10.3389/fninf.2014.00014/full). The dataset used comes from the [nilearn data](https://www.nitrc.org/frs/?group_id=728)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import nibabel   # package for loading/saving neuroimaging data\n",
    "bg_img = nibabel.load('data/bg.nii.gz')\n",
    "bg = bg_img.get_fdata()\n",
    "type(bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Create activation map by thresholding the data\n",
    "act_thresh = 6000\n",
    "act = bg.copy()\n",
    "# Set \"unactivated\" voxels to NaN for visualization\n",
    "act[act <= act_thresh] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# imshow kwargs\n",
    "imshow_opts = {\n",
    "    \"origin\" : \"lower\",\n",
    "    \"interpolation\" : \"nearest\"\n",
    "}\n",
    "\n",
    "# Axial slice of activation map overlay\n",
    "plt.imshow(bg[...,10].T, cmap=\"gray\");             # Background\n",
    "plt.imshow(act[...,10].T, cmap=\"plasma\");          # Activation map\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Interested in neuroimaging? Check out [openneuro.org](https://openneuro.org/) for curated data sets from published neuroimaging studies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Detecting gravitational wave signature of black hole and neutron star mergers\n",
    "\n",
    "<center><img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/d/db/LIGO_measurement_of_gravitational_waves.svg/710px-LIGO_measurement_of_gravitational_waves.svg.png\" alt=\"CBC Chirp\"></center>\n",
    "\n",
    "From [Wikipedia](https://en.wikipedia.org/wiki/First_observation_of_gravitational_waves)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "YouTubeVideo('I_88S8DWbcU', autoplay=1, loop=1, playlist='I_88S8DWbcU')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "[PyCBC](https://pycbc.org/) is the toolkit used to analyze data from gravitational wave observatories like [LIGO](https://www.ligo.caltech.edu/) and [Virgo](http://www.virgo-gw.eu/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The [PyCBC tutorials](https://github.com/gwastro/PyCBC-Tutorials) have some really cool examples - let's recreate the \"chirp\" from [first ever direct detection of gravitational waves](https://en.wikipedia.org/wiki/First_observation_of_gravitational_waves) that resulted from two black holes merging. For more info, see [the second PyCBC tutorial](https://colab.research.google.com/github/gwastro/pycbc-tutorials/blob/master/tutorial/2_VisualizationSignalProcessing.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import pycbc\n",
    "from pycbc import catalog\n",
    "\n",
    "merger_data = catalog.Merger('GW150914')\n",
    "# Though the catalog includes data from multiple observatories,\n",
    "# let's focus on just one\n",
    "ligo_data = merger_data.strain('L1')\n",
    "type(ligo_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "`pycbc` has its own (quite extensive) API that uses `numpy` and `scipy` under the hood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(type(ligo_data._data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To re-create the \"chirp\" we have to do some analysis on the raw data. \n",
    "\n",
    "Let's start by applying a simple band-pass filter. This is simpler than the analysis method [used in the official pycbc tutorial](https://colab.research.google.com/github/gwastro/pycbc-tutorials/blob/master/tutorial/2_VisualizationSignalProcessing.ipynb), but works suprisingly well!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Apply a bandpass filter to the data\n",
    "res = ligo_data.highpass_fir(20, 512).lowpass_fir(350, 512)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "`pycbc` relies on tools in `scipy.signal` to conduct the frequency analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "pycbc.filter.lowpass_fir??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "pycbc.filter.fir_zero_filter??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.signal import lfilter\n",
    "lfilter?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's take a look at the results of our filter analysis..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "time_of_merger = merger_data.time\n",
    "\n",
    "# Look 500 msec-worth of data around the merger time\n",
    "roi = res.time_slice(time_of_merger - 0.25, time_of_merger + 0.25)\n",
    "\n",
    "# Similar to a spectrogram with more sophisticated, irregular sampling\n",
    "times, freqs, power = roi.qtransform(\n",
    "    delta_t=0.001,\n",
    "    logfsteps=100,\n",
    "    qrange=(8, 8),\n",
    "    frange=(30, 512),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6,3))\n",
    "ax.pcolormesh(times, freqs, power**0.5)\n",
    "ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Generating the first ever direct image of a black hole\n",
    "\n",
    "On April 10th 2019, the [Event Horizon Telescope](https://eventhorizontelescope.org/) collaboration released the [first ever image of a black hole](https://eventhorizontelescope.org/press-release-april-10-2019-astronomers-capture-first-image-black-hole):\n",
    "\n",
    "<center><img src=\"https://static.projects.iq.harvard.edu/files/styles/os_files_xlarge/public/eht/files/20190410-78m-800x466.png?m=1554877319&itok=ryK319ed\" alt=\"EHT_M87_04-10-19\"/></center>\n",
    "\n",
    "Image source: The [official blog post](https://eventhorizontelescope.org/press-release-april-10-2019-astronomers-capture-first-image-black-hole) from the EHT collaboration announcing the result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The data and analysis pipeline are *way* too complicated to cover in a few slides. \n",
    "\n",
    "Instead, we'll just take advantage of the fact that the imaging pipeline is built on the tools of the scientific python ecosystem:\n",
    "\n",
    "<center><img src=\"images/ehtim_dependency_graphic.png\" alt=\"eht-imaging dependency graph\" height=648 width=1152/></center>\n",
    "\n",
    "Image credit: [Shaloo Shalini (@shaloo)](https://github.com/shaloo). For info on how this graphic was created, check out [shaloo's script](https://github.com/numpy/numpy.org/pull/23)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's run [the eht imaging pipeline](https://github.com/eventhorizontelescope/2019-D01-02) provided by the Event Horizons collaborators to help produce images from their [calibrated data](https://github.com/eventhorizontelescope/2019-D01-01). \n",
    "\n",
    "These repos with the helper-script for running the pipeline and the calibrated data have been included as submodules in `event_horizon_example/`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run event_horizons_example/2019-D01-02/eht-imaging/eht-imaging_pipeline.py  -i event_horizons_example/2019-D01-01/uvfits/SR1_M87_2017_101_lo_hops_netcal_StokesI.uvfits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "See? ...It's complicated. Here's the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = im_out.display(cfun=plt.cm.plasma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Parker Solar Probe\n",
    "\n",
    "The [Parker Solar Probe](https://www.nasa.gov/content/goddard/parker-solar-probe-humanity-s-first-visit-to-a-star) was launched in 2018 to study the solar atmosphere, coming nearer to the sun than any previous space probe. Data from the Parker probe is already yielding [unexpected results](https://news.engin.umich.edu/2019/12/were-missing-something-fundamental-about-the-sun/?utm_source=newsletter&utm_medium=email&utm_campaign=January_2020).\n",
    "\n",
    "One of the instruments on the probe is [SWEAP](http://sweap.cfa.harvard.edu/), a set of charged-particle detectors. Publicly available data from SWEAP can be found [here](http://sweap.cfa.harvard.edu/Data.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Full disclosure: the original analysis for the [results published in Nature](https://www.nature.com/articles/s41586-019-1813-z) were produced with IDL, not Python. However, NASA just released [the first batch of data from the probe](https://sppgway.jhuapl.edu/) to the public, so let's see if we can't replicate some results..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### A quick aside: data formats\n",
    "\n",
    "Data from the Parker probe is stored in NASA's [Common Data Format (CDF)](https://cdf.gsfc.nasa.gov/). Python libraries such as [spacepy](https://spacepy.github.io/) are used for I/O from the CDF format. As you might expect, `spacepy`'s `pycdf` module loads data from CDF files into NumPy arrays. Unfortunately, `spacepy`'s `pycdf` module depends on an external C-library, and there is not (yet) a `conda` recipe for installing it automatically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To get around this I've used `spacepy.pycdf` to save a small amount of SWEAP in the more Python-friendly `.npz` format. I took data from the [SPC instrument collected on 11-08-18](http://sweap.cfa.harvard.edu/pub/data/sci/sweap/spc/L2/2018/11/). If you'd like to work with the full dataset, you can [manually install CDF](https://spacepy.github.io/install_linux.html#cdf), download the data (or any other dataset), and use `devlogs/parker_probe_velocity_log.py` as an example of how to load and interact with the raw data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# See the full notebook to download this dataset\n",
    "fname = 'data/parker_probe_spcL2_data_11-08-18.npz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# The dataset is ~150MB, so it is not included in the repo.\n",
    "\n",
    "import os, requests, tqdm\n",
    "dsize = 152538956   # File size in bytes\n",
    "dlink = 'https://www.dropbox.com/s/z45tbkqwjpyu6tz/parker_probe_spcL2_data_11-08-18.npz?dl=0'\n",
    "if not os.path.exists(fname):\n",
    "    r = requests.get(\n",
    "        dlink,\n",
    "        headers={'user-agent':'Wget/1.20 (linux-gnu)'},\n",
    "        stream=True\n",
    "    )\n",
    "    with open(fname, 'wb') as fh:\n",
    "        for chunk in tqdm.tqdm(r.iter_content(chunk_size=1024), total=dsize/1024):\n",
    "            if chunk:\n",
    "                fh.write(chunk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now that we have the data, let's try to replicate the top pane of [this image](http://sweap.cfa.harvard.edu/Images/example_spc_ql.png) from the [SWEAP data page](http://sweap.cfa.harvard.edu/Data.html).\n",
    "\n",
    "We don't have time to discuss the data in detail, but the [Appendix 3 of the SWEAP Data User's Guide](http://sweap.cfa.harvard.edu/sweap_data_user_guide.pdf) outlines a procedure we can use to reproduce the desired figure. We start by loading the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Load data from the SPC instrument on the Parker probe\n",
    "data = np.load(fname)\n",
    "# Measurement time (x-axis of image)\n",
    "t = data['t']\n",
    "# Edges of Voltage bins (y-axis of image)\n",
    "mv_lo = data['mv_lo'].T\n",
    "mv_hi = data['mv_hi'].T\n",
    "# Differential charge flux density\n",
    "dcfd = data['diff_charge_flux_density'].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# The data include timestamps with microsecond resolution\n",
    "# and 128 channels per data point\n",
    "print(t.shape, mv_lo.shape, dcfd.shape)\n",
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# The CDF file uses a fill value (-1e31) to denote invalid data\n",
    "print(dcfd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Let's replace them so we can keep track of non-data in the arrays\n",
    "for arr in (mv_lo, mv_hi, dcfd):\n",
    "    arr[arr == -1e31] = np.nan\n",
    "print(dcfd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Upon closer inspection, only the first 30 of the 128 channels store valid data\n",
    "np.sum(np.isfinite(dcfd), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Limit computation to valid voltage bins\n",
    "mv_lo, mv_hi, dcfd = mv_lo[:31,:], mv_hi[:31,:], dcfd[:31,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "After removing the unused data channels, there are still individual measurements that resulted in invalid data. Let's remove these as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Mask out time samples that have bad data\n",
    "bad_data = np.any(~np.isfinite(dcfd), axis=0)\n",
    "t = t[~bad_data]\n",
    "mv_lo = mv_lo[:, ~bad_data]\n",
    "mv_hi = mv_hi[:, ~bad_data]\n",
    "dcfd = dcfd[:, ~bad_data]\n",
    "print(\"{} time samples out of {} discarded\".format(bad_data.sum(), bad_data.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "That takes care of the data munging, now on to the computation. The procedure in [Appendix 3 of the SWEAP data user's guide](http://sweap.cfa.harvard.edu/sweap_data_user_guide.pdf) boils down to a few straight-forward steps:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### 1. Compute the center of the voltage bins, $V$,  from `mv_lo` and `mv_hi`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### 2. Transform from particle *energy* to particle *velocity*\n",
    "\n",
    "[Appendix 3](http://sweap.cfa.harvard.edu/sweap_data_user_guide.pdf) provides some helpful formulae:\n",
    " \n",
    "$v_{p} = \\sqrt{\\frac{2qmv_{hi}}{m_{p}}} \\frac{2}{\\pi}E(\\frac{mv_{lo}}{mv_{hi}})$\n",
    "    \n",
    "$dv_{p} = \\sqrt{\\frac{4qV}{m_{p}}} - v_{p}^{2}$\n",
    "\n",
    "Where $v_{p}$ is the proton velocity, $dv_{p}$ is the bin width in velocity space, $q$ and $m_p$ are the charge and mass of the proton, respectively, and $E(x)$ represents an approximation to the elliptical integral."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### 3. Finally, compute the distribution of proton velocity, $F(v_{p})$, from the differential charge flux density measured by the instrument.\n",
    "\n",
    "Again, [Appendix 3](http://sweap.cfa.harvard.edu/sweap_data_user_guide.pdf) gives us everything we need in describing the relationship between the differential charge flux density and the distribution of proton velocity:\n",
    "\n",
    "$ dcfd = q v_{p} F(v_{p})dv_{p} \\cdot 10^{8} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Compute center and widths of voltage bins\n",
    "V = (mv_hi + mv_lo) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# We need the mass and charge of the proton for the next calculation\n",
    "from scipy.constants import m_p as mp   # Proton mass [kg]\n",
    "from scipy.constants import e as q      # Fundamental charge [C]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We'll also need an approximation to the elliptic integral equation, $E(x)$. We'll use one of the approximations provided in `scipy.special`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.special import ellipe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2: Convert from energy -> velocity\n",
    "v = np.sqrt(2 * q * mv_hi / mp) * (2 / np.pi) * ellipe(mv_lo/mv_hi)\n",
    "dv = np.sqrt((4 * (q/mp) * V) - v**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Finally, compute the proton distribution as a function of proton velocity\n",
    "Fv = (dcfd / (q* v* 10**8)) * (1 / dv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's see how we did..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "viz_kwargs = {\n",
    "    \"cmap\" : plt.cm.plasma,\n",
    "    \"norm\" : colors.LogNorm(vmin=1, vmax=200)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Visualize\n",
    "fig = plt.figure(figsize=(8,4)); ax = fig.add_subplot(111)\n",
    "ax.pcolormesh(t[np.newaxis,:], v/1000, Fv, **viz_kwargs)\n",
    "fig.colorbar(ax.collections[0])\n",
    "fig.autofmt_xdate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Part III - Developing NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# [Scope of NumPy](https://numpy.org/neps/scope.html)\n",
    "\n",
    "The NumPy execution engine currently targets:\n",
    "\n",
    " * in-memory, homogenously-typed array data\n",
    " * cpu-based operations\n",
    " \n",
    "Specialized hardware (e.g. GPUs), features for scalable computing (e.g. distributed arrays) are currently out of scope\n",
    " - Supporting libraries that provide these features *is in scope*: **extensibility** and **interoperability**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Important guiding principles:\n",
    " - **Stability**: Foundational component of the scientific python ecosystem for going-on 15 years\n",
    " - **Interoperability**\n",
    "   * NumPy is the standard array data structure within the scientific Python ecosystem\n",
    "   * What about all the new array libraries?\n",
    "     - [XArray](http://xarray.pydata.org/en/stable/)\n",
    "     - [Dask Arrays](https://docs.dask.org/en/latest/array.html)\n",
    "     - [Jax](https://jax.readthedocs.io/en/latest/)\n",
    "     - [pydata sparse](https://sparse.pydata.org/en/latest/)\n",
    "     - [PyTorch](https://pytorch.org/)\n",
    "     - [TensorFlow](https://www.tensorflow.org/api_docs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# How is NumPy Developed?\n",
    "\n",
    " - **Collaboratively** - https://github.com/numpy/numpy/\n",
    "\n",
    "Commitment to stability means proposed changes must go through extensive design and review:\n",
    " - [Numpy Enhancement Proposals (NEPs)](https://numpy.org/neps/) - analogous to PEPs, specific to NumPy\n",
    "   * Community-driven development and consensus among contributors/developers\n",
    "     - Mailing list\n",
    "     - PRs/Issues on GitHub\n",
    " - Steering council for high-level direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A case-study: `np.random`\n",
    " - Changes proposed in [NEP 19](https://numpy.org/neps/nep-0019-rng-policy.html), subsequently approved by the community via discussion on the mailing list and on GitHub.\n",
    " - Overhaul of `np.random` landed in version 1.17\n",
    " \n",
    "   * Improve *performance* and *flexibility* without sacrificing stability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Generate 1,000,000 random numbers the old way\n",
    "old_rands = np.random.random(int(1e6))\n",
    "print(\"Uniform random numbers from legacy np.random.random:\\n  {}\".format(old_rands))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# ... and the new way\n",
    "from numpy.random import default_rng\n",
    "rg = default_rng()\n",
    "new_rands = rg.random(int(1e6))\n",
    "print(\"Uniform random numbers with new tools:\\n  {}\".format(new_rands))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Compatibility\n",
    "\n",
    "There are many, many LOC (both in test suites and in production) that depend on the original `numpy.random`, so both the *interface* and the *results* must remain unchanged\n",
    " * <font color=\"green\">**Upside: Stability**</font> - output of `np.random` remains consistent with previous versions\n",
    " * <font color=\"orange\">**Downside: Discoverability**</font> - users need to know about new interface to access improvements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Choose a seed for generator\n",
    "seed = 1817\n",
    "\n",
    "# Random numbers generated by np.random in v1.15\n",
    "rands_from_v1_15 = np.load('data/npy_v1.15_random_seed1817_1000samples.npy')\n",
    "# Generate random numbers with legacy interface\n",
    "np.random.seed(seed)\n",
    "legacy_rands = np.random.random(1000)\n",
    "\n",
    "print(\"Arrays equivalent: \", np.allclose(rands_from_v1_15, legacy_rands))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "It is possible (though clunky) to replicate legacy behavior with new interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "seed = 1817\n",
    "\n",
    "from numpy.random import MT19937, RandomState, Generator\n",
    "# Set random state with legacy seeding\n",
    "rs = RandomState(seed)\n",
    "mt = MT19937()\n",
    "mt.state = rs.get_state()\n",
    "\n",
    "# New interface for generation\n",
    "rg = Generator(mt)\n",
    "mt_rands = rg.random(1000)\n",
    "print(\"Legacy: {}\\nGenerator: {}\".format(legacy_rands[:4], mt_rands[:4]))\n",
    "print(\"Arrays equivalent: \", np.allclose(legacy_rands, mt_rands))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Performance\n",
    "\n",
    "The [PCG64](https://docs.scipy.org/doc/numpy/reference/random/bit_generators/pcg64.html) BitGenerator is a \n",
    "[significant improvement](http://www.pcg-random.org/) over the legacy Marsenne Twister in many areas, including speed:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "#NOTE: PCG64 is the new default bit_generator, so default_rng() equivalent to Generator(PCG64())\n",
    "from numpy.random import default_rng\n",
    "rg = default_rng()\n",
    "num_samples = int(1e5) \n",
    "\n",
    "print(\"Uniform random numbers:\")\n",
    "%timeit np.random.random(num_samples)\n",
    "%timeit rg.random(num_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In addition, `Generator` includes improved methods for drawing samples from distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Standard Normal:\")\n",
    "%timeit -n 100 np.random.standard_normal(num_samples)\n",
    "%timeit -n 100 rg.standard_normal(num_samples) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Standard Exponential:\")\n",
    "%timeit -n 100 np.random.standard_exponential(num_samples)\n",
    "%timeit -n 100 rg.standard_exponential(num_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Standard Gamma:\")\n",
    "shape_param = 3.0\n",
    "%timeit -n 100 np.random.standard_gamma(shape_param, num_samples)\n",
    "%timeit -n 100 rg.standard_gamma(shape_param, num_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Parallel Generation\n",
    "\n",
    "`np.random` includes new functionality to produce high-quality initital states for multiple generators to produce reproducible random numbers accross multiple processes.\n",
    "\n",
    "For one example, let's take a look at `SeedSequence` and an [example from the documentation](https://numpy.org/devdocs/reference/random/multithreading.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Code(filename=\"mrng.py\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from mrng import MultithreadedRNG, default_rng\n",
    "num_workers = 4\n",
    "seed = 1817\n",
    "n = int(1e7)\n",
    "\n",
    "# Compare concurrent.futures multithreaded generation to single thread\n",
    "rg = default_rng()\n",
    "mg = MultithreadedRNG(n, seed=seed, threads=num_workers)\n",
    "\n",
    "%timeit rg.standard_normal(n)\n",
    "%timeit mg.fill()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Maintain reproducible random number generators\n",
    "ex1 = MultithreadedRNG(n, seed=seed, threads=num_workers)\n",
    "ex2 = MultithreadedRNG(n, seed=seed, threads=num_workers)\n",
    "\n",
    "# Generate numbers\n",
    "ex1.fill()\n",
    "ex2.fill()\n",
    "\n",
    "# Results are reproducible\n",
    "np.allclose(ex1.values, ex2.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Part IV - Looking Ahead"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The changing landscape\n",
    "\n",
    " - In the early days, many new NumPy users were converts from languages like Matlab and IDL\n",
    "   * See the [NumPy for Matlab users](https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html) article in the docs\n",
    "   \n",
    " - **Now**: The scientific Python ecosystem (including libraries for data science and ML) is incredibly feature-rich and powerful, and is attracting many new users.\n",
    "   * Users interested in specific applications (machine learning, image processing, geoscience, bioinformatics, etc.) end up interacting with NumPy indirectly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Google Trends"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Data downloaded from google trends\n",
    "!ls data/*.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "!head data/datascience.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "timeseries_dtype = np.dtype([\n",
    "    ('date', 'datetime64[M]'),\n",
    "    ('relpop', float)\n",
    "])\n",
    "\n",
    "parse_kwargs = {\n",
    "    \"skiprows\" : 3,\n",
    "    \"delimiter\" : \",\",\n",
    "    \"dtype\" : timeseries_dtype\n",
    "}\n",
    "\n",
    "fnames = (\"numpy\", \"datascience\", \"matlab\")\n",
    "\n",
    "data = {\n",
    "    fname : np.loadtxt(\"data/{}.csv\".format(fname), **parse_kwargs) for fname in fnames\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for name, vals in data.items():\n",
    "    plt.plot(vals['date'], vals['relpop'], label=name)\n",
    "ax.set_title('Google Trends (US): 2004 - Present')\n",
    "ax.set_ylabel('Relative Popularity of Search Term [arb]')\n",
    "fig.autofmt_xdate()\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def smooth(s, kernsize=21):\n",
    "    s_padded = np.hstack((s[kernsize-1:0:-1], s, s[-2:-kernsize-1:-1]))\n",
    "    kern = np.hamming(kernsize)\n",
    "    res_padded = np.convolve(kern/kern.sum(), s_padded, mode='valid')\n",
    "    # De-pad and renormalize\n",
    "    return 100 * res_padded[kernsize//2:-kernsize//2+1] / res_padded.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for name, vals in data.items():\n",
    "    plt.plot(vals['date'], smooth(vals['relpop']), label=name)\n",
    "ax.set_title('Google Trends (US): 2004 - Present')\n",
    "ax.set_ylabel('Relative Popularity of Search Term [arb]')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What's next for NumPy?\n",
    "\n",
    "<center><img src=\"images/numpy_roadmap_graphic.png\" alt=\"Numpy-near-future-graphic\" height=648 width=1152/></center>\n",
    "\n",
    "Image modified from [this PyData Amsterdam 2019 presentation](https://www.slideshare.net/RalfGommers/the-evolution-of-array-computing-in-python/14) by [Ralf Gommers](https://github.com/rgommers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Interoperability\n",
    "\n",
    "Separate NumPy API from NumPy *execution engine*\n",
    " - Allow other libraries ([Dask](https://dask.org/), [CuPy](https://cupy.chainer.org/), [PyTorch](https://pytorch.org/), etc.) to support NumPy API.\n",
    " - Mitigate ecosystem fragmentation\n",
    "   * E.g. don't want a re-implementation of `scipy` for each ML framework (`pytorch.scipy`, `tensorflow.scipy`, etc.)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Current n-dimensional array landscape\n",
    "\n",
    "<center><img src=\"images/array_landscape_now.png\" alt=\"Array-ecosystem-now\" height=648 width=1152/></center>\n",
    "\n",
    "Images from this [talk at PyData NY 2019](https://www.slideshare.net/RalfGommers/pydata-nyc-whatsnew-numpyscipy-2019?next_slideshow=1) by [Ralf Gommers](https://github.com/rgommers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Vision for the future\n",
    "\n",
    "<center><img src=\"images/array_landscape_vision.png\" alt=\"array-ecosystem-vision\" height=658 width=1152/></center>\n",
    "\n",
    "Images from this [talk at PyData NY 2019](https://www.slideshare.net/RalfGommers/pydata-nyc-whatsnew-numpyscipy-2019?next_slideshow=1) by [Ralf Gommers](https://github.com/rgommers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## One approach: `__array_function__` protocol\n",
    "\n",
    " - Proposed in [NEP 18](https://numpy.org/neps/nep-0018-array-function-protocol.html)\n",
    " - Array function protocol enabled by default as of version 1.17\n",
    " \n",
    "<center><img src=\"images/array_function_descr.png\" alt=\"array_function_protocol\"/></center>\n",
    " \n",
    "Image source: [this presentation](https://www.slideshare.net/RalfGommers/arrayfunction-conceptual-design-related-concepts?from_action=save) by [Ralf Gommers](https://github.com/rgommers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### `__array_function__` example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "rg = np.random.default_rng()\n",
    "x = rg.random((5000, 1000))\n",
    "\n",
    "# Factorize with np.linalg\n",
    "q, r = np.linalg.qr(x)\n",
    "type(q), type(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import dask.array as da\n",
    "\n",
    "d = da.from_array(x, chunks=(1000, 1000))\n",
    "\n",
    "# Same call signature!\n",
    "q, r = np.linalg.qr(d)\n",
    "type(q), type(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "da.core.Array??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lessons learned from `__array_function__`\n",
    "\n",
    " - The `__array_function__` protocol has been successful, but has fallen short of universal adoption.\n",
    " - Valuable feedback from the community has resulted in [NEP 37](https://numpy.org/neps/nep-0037-array-module.html)\n",
    "   * Defines `__array_module__` protocol\n",
    "   * Currently under development (interested?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The problem with data types...\n",
    "\n",
    "Current `dtype` system has some flexibility issues\n",
    " - Difficult to specify fully-featured types\n",
    " - Some mechanisms (e.g. casting rules) are difficult to extend to new types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Goal: Improve NumPy maintainability\n",
    " * Improve organization of dtype checking/comparison machinery\n",
    " * Use the same API for built-in and user-defined dtypes\n",
    " * Improve extensibility of API: facilitate future additions/modifications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### User impact\n",
    " - Easier-to-use mechanism for defining fully-feature dtypes (including from Python)\n",
    " - Host of new dtypes for the ecosystem:\n",
    "   * Physical units (cf. [astropy.unit](https://docs.astropy.org/en/stable/units/))\n",
    "   * `bfloat16`, `int24`, etc.\n",
    "   * Categorical types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "An approach to overhauling the dtype system is [currently being fleshed out in a new NEP](https://github.com/numpy/numpy/blob/a111b551ae940d7d5f8523fef1cf3589c6ba00a0/doc/neps/nep-0033-extensible-dtypes.rst)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Improved SIMD incorporation for `ufuncs`\n",
    "\n",
    "Strike a balance between **optimization** and **maintainability**\n",
    "\n",
    " - Define set of architecture-agnostic universal intrinsics\n",
    "   * At build time, build code paths based on features available for the host architecture\n",
    "   * At run time, detect which features are available and select which of available code paths to use\n",
    " - In the process of being formalized in a [draft NEP](https://github.com/mattip/numpy/blob/nep_simd/doc/neps/nep-XXXX-SIMD-optimizations.rst)\n",
    "   * Preliminary work in support of this proposed enhancement can be found [here](https://github.com/numpy/numpy/pull/13421/files) and [here](https://github.com/numpy/numpy/pull/13516)\n",
    "   \n",
    "**N.B.** this approach (i.e. using universal intrinsics) was adopted by OpenCV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supporting language features: type annotations\n",
    "\n",
    "Thinking about how best to support type annotations became especially important when they became an official core language feature in Python 3.7."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is currently the most cross-referenced issue in the NumPy GH repository:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Markdown(filename=\"data/top_issues_table.md\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Work on type annotations is located in the [numpy-stubs](https://github.com/numpy/numpy-stubs) repository. Basic type annotations are supported:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Code(filename=\"type_annotations.py\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "!mypy type_annotations.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# ... and beyond: NumPy 2.0?\n",
    "\n",
    " - Major revision -> opportunity for refactoring/enhancements that break API\n",
    "   * Weigh potential for improvements against the pain of breaking changes\n",
    "   * Commitment to stability still a central theme!\n",
    " - So much new functionality being developed in external libraries\n",
    "   * Changes that facilitate external development are priorities\n",
    " \n",
    "A bit of the history surrounding the idea of NumPy 2.0 can be found [here](https://github.com/numpy/numpy/issues/9066)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Getting involved\n",
    "\n",
    "NumPy presents an opportunity to work on a project that is depended on by tens of millions of users (and counting). Here's how you can get involved:\n",
    " 1. Where discussion happens:\n",
    "  - [Numpy discussion mailing list](https://www.scipy.org/scipylib/mailing-lists.html)\n",
    "  - Numpy community meetings - video conference every-other-week: [Community calendar link](https://calendar.google.com/calendar?cid=YmVya2VsZXkuZWR1X2lla2dwaWdtMjMyamJobGRzZmIyYzJqODFjQGdyb3VwLmNhbGVuZGFyLmdvb2dsZS5jb20)\n",
    "  - slack channel: numpy-team.slack.com\n",
    " 2. Contribute\n",
    "   - [GitHub Issues](https://github.com/numpy/numpy/issues) and [open PRs](https://github.com/numpy/numpy/pulls) are a great entry point\n",
    "     * If you want to get your hands dirty immediately, try starting with the [good first issue](https://github.com/numpy/numpy/issues?q=is%3Aopen+is%3Aissue+label%3A%22good+first+issue%22) label\n",
    "     * For challenges with a greater scope, try the [Enhancement](https://github.com/numpy/numpy/labels/01%20-%20Enhancement) or [Wish List](https://github.com/numpy/numpy/labels/23%20-%20Wish%20List) labels\n",
    "   - Check out the discussion revolving around accepted and proposed [NEPs](https://numpy.org/neps/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Thank you!\n",
    "\n",
    "If you have any questions, comments, or ideas please don't hesitate to contact me: rossbar@berkeley.edu\n",
    "\n",
    "Also feel free to ask about/use/modify/contribute to this presentation on GitHub!\n",
    "\n",
    "## Acknowledgements\n",
    "\n",
    "> This project is funded in part by the Gordon and Betty Moore Foundation through\n",
    "> [Grant GBMF5447](https://www.moore.org/grant-detail?grantId=GBMF5447f) and by\n",
    "> the Alfred P. Sloan Foundation through \n",
    "> [Grant G-2017-9960](https://sloan.org/grant-detail/8222)\n",
    "> to the University of California, Berkeley.\n",
    "\n",
    "And a special thanks to [@stefanv](https://github.com/stefanv), [@seberg](https://github.com/seberg), and [@mattip](https://github.com/mattip/) for their generous input in the creation of this presentation.\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}