{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Losing your Loops\n", "## Fast Numerical Computing with NumPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python is Fast\n", "For Writing, Testing, and Developing of Code" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``` python\n", "# Hello World in Python\n", "print(\"hello world\")\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``` java\n", "/* Hello World in Java */\n", "public class HelloWorld {\n", " public static void main(String[] args) {\n", " System.out.println(\"Hello, World\");\n", " }\n", "}\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import seaborn as sns\n", "data = sns.load_dataset(\"iris\")\n", "sns.pairplot(data, hue=\"species\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python is Slow\n", "Compared to compiled languages" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# A silly function implemented in Python\n", "\n", "def func_python(N):\n", " d = 0.0\n", " for i in range(N):\n", " d += (i % 3 - 1) * i\n", " return d" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Use IPython timeit magic to time the execution\n", "%timeit func_python(10000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fortran version:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%load_ext fortranmagic" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%fortran\n", "subroutine func_fort(n, d)\n", " integer, intent(in) :: n\n", " double precision, intent(out) :: d\n", " integer :: i\n", " d = 0\n", " do i = 0, n - 1\n", " d = d + (mod(i, 3) - 1) * i\n", " end do\n", "end subroutine func_fort" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "%%file func_fortran.f\n", "\n", " subroutine func_fort(n, d)\n", " integer, intent(in) :: n\n", " double precision, intent(out) :: d\n", " integer :: i\n", " d = 0\n", " do i = 0, n - 1\n", " d = d + (mod(i, 3) - 1) * i\n", " end do\n", " end subroutine func_fort\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# use f2py rather than f2py3 for Python 2\n", "!f2py3 -c func_fortran.f -m func_fortran > /dev/null" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from func_fortran import func_fort" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "%timeit func_fort(10000)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Outline\n", "\n", " 1. Use Numpy **ufuncs** to your advantage\n", "\n", " 2. Use Numpy **aggregates** to your advantage\n", "\n", " 3. Use Numpy **broadcasting** to your advantage\n", "\n", " 4. Use Numpy **slicing and masking** to your advantage\n", "\n", " 5. Use a tool like *SWIG*, *Weave*, *cython*, *f2py*, *Numba*, etc.\n", " to compile Python or to interface to compiled code." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Fortran is about 100 times faster for this task!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Why is Python so slow?\n", "\n", "We alluded to this yesterday, but languages tend to have a compromise between convenience and performance.\n", "\n", "- C, Fortran, etc.: **static typing** and **compiled code** leads to fast execution\n", "\n", " + But: lots of development overhead in declaring variables, no interactive prompt, etc.\n", "\n", "- Python, R, Matlab, IDL, etc.: **dynamic typing** and **interpreted excecution** leads to fast development\n", "\n", " + But: lots of execution overhead in dynamic type-checking, etc.\n", " \n", "We like Python because our **development time** is generally more valuable than **execution time**. But sometimes speed can be an issue." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategies for making Python fast\n", "\n", "1. Use Numpy **ufuncs** to your advantage\n", "\n", "2. Use Numpy **aggregates** to your advantage\n", "\n", "3. Use Numpy **broadcasting** to your advantage\n", "\n", "4. Use Numpy **slicing and masking** to your advantage\n", "\n", "5. Use a tool like *SWIG*, *cython* or *f2py* to interface to compiled code.\n", "\n", "Here we'll cover the first four, and leave the fifth strategy for a later session." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategy 1: Use ufuncs to your advantage\n", "\n", "A **ufunc** in numpy is a *Universal Function*. This is a function which operates element-wise on an array. We've already seen examples of these in the various arithmetic operations:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = [1, 3, 2, 4, 3, 1, 4, 2]\n", "b = [val + 5 for val in a]\n", "print(b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "a = np.array(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "b = a + 5 # element-wise\n", "print(b)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The speed of ufuncs" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = list(range(100000))\n", "%timeit [val + 5 for val in a]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = np.array(a)\n", "%timeit a + 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other ufuncs\n", "\n", "There are many, many ufuncs available:\n", "\n", "- Arithmetic:\t\t\t\t``+ - * / // % **``\n", "- Bitwise Operations:\t\t``& | ~ ^ >> <<``\n", "- Comparisons:\t\t\t``< > <= >= == !=``\n", "- Trig Functions:\t\t\t``np.sin np.cos np.tan`` ...etc.\n", "- Exponential Family:\t\t``np.exp np.log np.log10`` ...etc.\n", "- Special Functions:\t\t``scipy.special.*``\n", "\n", "and many, many more.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategy 2. use aggregations to your advantage\n", "\n", "Aggregations are functions over arrays which return smaller arrays.\n", "\n", "Suppose you want to compute the minimum of an array" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from random import random\n", "c = [random() for i in range(100000)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%timeit min(c)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "c = np.array(c)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%timeit c.min()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aggregates along axes" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M = np.random.randint(0, 10, (3, 5))\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M.sum()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M.sum(axis=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M.sum(axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other Aggregation Functions\n", "\n", "Numpy has many useful aggregation functions:\n", "\n", "``np.min np.max np.sum np.prod np.mean np.std np.var np.any np.all np.median np.percentile np.argmin np.argmax``\n", "\n", "Most also have a NaN-aware equivalent:\n", "\n", "``np.nanmin np.nanmax np.nansum`` ...etc." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Strategy 3: Use Broadcasting to your advantage\n", "\n", "Broadcasting in NumPy is the set of rules for applying ufuncs on arrays of different sizes and/or dimensions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "np.arange(3) + 5" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "np.ones((3, 3)) + np.arange(3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "np.arange(3).reshape((3, 1)) + np.arange(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualizing Broadcasting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "([image source](http://www.astroml.org/book_figures/appendix/fig_broadcast_visual.html))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Rules of Broadcasting\n", "\n", "1. If array dimensions differ, left-pad the smaller shape with 1s\n", "\n", "2. If any dimension does not match, stretch the dimension with size=1\n", "\n", "3. If neither non-matching dimension is 1, raise an error." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Examples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M = np.ones((2, 3))\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = np.arange(3)\n", "a" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M + a" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = np.arange(3).reshape((3, 1))\n", "a" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "b = np.arange(3)\n", "b" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a + b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example 3" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M = np.ones((3, 2))\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = np.arange(3)\n", "a" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M + a" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Strategy 4: Use slicing and masking to your advantage\n", "\n", "The last strategy we will cover is slicing and masking." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python lists can be indexed with integers or slices:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "L = [2, 3, 5, 7, 11]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L[0] # integer index" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L[1:3] # slice for multiple elements" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### NumPy arrays are like lists" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L = np.array(L)\n", "L" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L[0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L[1:3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Masking" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mask = np.array([False, True, True,\n", " False, True])\n", "L[mask]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mask = (L < 4) | (L > 8) # \"|\" = \"bitwise OR\"\n", "L[mask]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fancy Indexing" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "L" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ind = [0, 4, 2]\n", "L[ind]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Multiple Dimensions" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M = np.arange(6).reshape(2, 3)\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# multiple indices separated by comma\n", "M[0, 1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# mixing slices and indices\n", "M[:, 1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# masking the full array\n", "M[abs(M - 3) < 2]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# mixing fancy indexing and slicing\n", "M[[1, 0], :2]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# mixing masking and slicing \n", "M[M.sum(axis=1) > 4, 1:]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Putting it All Together\n", "Nearest Neighbors of some data" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 1000 points in 3 dimensions\n", "X = np.random.random((1000, 3))\n", "X.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Broadcasting to find pairwise differences\n", "diff = X.reshape(1000, 1, 3) - X\n", "diff.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Aggregate to find pairwise distances\n", "D = (diff ** 2).sum(2)\n", "D.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set diagonal to infinity to skip self-neighbors\n", "i = np.arange(1000)\n", "D[i, i] = np.inf" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# print the indices of the nearest neighbor\n", "i = np.argmin(D, 1)\n", "print(i[:10])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# double-check with scikit-learn\n", "from sklearn.neighbors import NearestNeighbors\n", "d, i = NearestNeighbors().fit(X).kneighbors(X, 2)\n", "print(i[:10, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary: Speeding up NumPy\n", "\n", "It's all about **moving loops into compiled code:**\n", "\n", "1. Use Numpy **ufuncs** to your advantage (eliminate loops!)\n", "\n", "2. Use Numpy **aggregates** to your advantage (eliminate loops!)\n", "\n", "3. Use Numpy **broadcasting** to your advantage (eliminate loops!)\n", "\n", "4. Use Numpy **slicing and masking** to your advantage (eliminate loops!)\n", "\n", "5. Use a tool like *SWIG*, *cython* or *f2py* to interface to compiled code." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "", "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }