notebooks/11_EfficientNumpy.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:1d479d0c5c65d1105d153a32be2ce657e26c643868f0ab6f916a3b28320fd76c"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<small><i>This notebook was put together by [Jake Vanderplas](http://www.vanderplas.com) for UW's [Astro 599](http://www.astro.washington.edu/users/vanderplas/Astr599_2014/) course. Source and license info is on [GitHub](https://github.com/jakevdp/2014_fall_ASTR599/).</i></small>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%run talktools.py"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<style>\n",
        "div.cell, div.text_cell_render{\n",
        "  max-width:750px;\n",
        "  margin-left:auto;\n",
        "  margin-right:auto;\n",
        "}\n",
        "\n",
        ".rendered_html\n",
        "{\n",
        "  font-size: 140%;\n",
        "  }\n",
        "\n",
        ".rendered_html li\n",
        "{\n",
        "  line-height: 1.8;\n",
        "  }\n",
        "\n",
        ".rendered_html h1, h2 {\n",
        "  text-align:center;\n",
        "  font-familly:\"Charis SIL\", serif;\n",
        "}\n",
        "</style>"
       ],
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.HTML at 0x103787810>"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Efficient Numerical Computing with Numpy\n",
      "\n",
      "In this session, we'll discuss how to make your programs as efficient as possible, mainly by taking advantage of *vectorization* in NumPy."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Aside: the \"Unladen Swallow\""
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import YouTubeVideo\n",
      "YouTubeVideo(\"y2R3FvS4xr4\")"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "html": [
        "\n",
        "        <iframe\n",
        "            width=\"400\"\n",
        "            height=300\"\n",
        "            src=\"https://www.youtube.com/embed/y2R3FvS4xr4\"\n",
        "            frameborder=\"0\"\n",
        "            allowfullscreen\n",
        "        ></iframe>\n",
        "        "
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "<IPython.lib.display.YouTubeVideo at 0x103787b90>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "That video was the original reference behind the now-defunct [unladen swallow](http://en.wikipedia.org/wiki/Unladen_Swallow) project, which had the goal of making Python's C implementation faster.\n",
      "\n",
      "Yes, Python is (unfortunately) a rather slow language.\n",
      "\n",
      "Here is an example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 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"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use IPython timeit magic to time the execution\n",
      "%timeit func_python(10000)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.63 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "To compare to a compiled language, let's write the same function in fortran and use the ``f2py`` tool (included in NumPy) to compile it"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%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"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting func_fortran.f\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# use f2py rather than f2py3 for Python 2\n",
      "!f2py3 -c func_fortran.f -m func_fortran > /dev/null"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from func_fortran import func_fort\n",
      "%timeit func_fort(10000)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 15.7 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "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",
      "- <font color=\"green\">C, Fortran, etc.</font>: **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",
      "- <font color=\"red\">Python, R, Matlab, IDL, etc.</font>: **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: Using UFuncs in Numpy\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",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "x = np.random.random(4)\n",
      "print(x)\n",
      "print(x + 1)  # add 1 to each element of x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.14407406  0.1104411   0.82417477  0.61817359]\n",
        "[ 1.14407406  1.1104411   1.82417477  1.61817359]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x * 2  # multiply each element of x by 2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "array([ 0.28814812,  0.22088221,  1.64834954,  1.23634718])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x * x  # multiply each element of x by itself"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "array([ 0.02075734,  0.01219724,  0.67926405,  0.38213859])"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[1:] - x[:-1]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "array([-0.03363296,  0.71373367, -0.20600118])"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "These are *binary ufuncs*: they take two arguments.\n",
      "\n",
      "There are also many *unary ufuncs:*"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "-x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([-0.14407406, -0.1104411 , -0.82417477, -0.61817359])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sin(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "array([ 0.14357615,  0.11021673,  0.73398757,  0.57954772])"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### The Speed of Ufuncs"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random(10000)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit\n",
      "# compute element-wise x + 1 via a ufunc \n",
      "y = np.zeros_like(x)\n",
      "y = x + 1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 21.1 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit\n",
      "# compute element-wise x + 1 via a loop\n",
      "y = np.zeros_like(x)\n",
      "for i in range(len(x)):\n",
      "    y[i] = x[i] + 1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 3.43 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Why is NumPy so much faster?\n",
      "\n",
      "Numpy UFuncs are faster than Python functions involving loops, because *the looping happens in compiled code*.  This is only possible when types are known beforehand, which is why numpy arrays must be typed."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Other Available Ufuncs\n",
      "\n",
      "- Trigonometric functions (``np.sin``, ``np.cos``, etc.)\n",
      "- Scipy special functions (``scipy.special.j0``, ``scipy.special.gammaln``, etc.)\n",
      "- Element-wise minimum/maximum (``np.minimum``, ``np.maximum``)\n",
      "- User-defined ufuncs (read more [here](http://scipy-lectures.github.io/advanced/advanced_numpy/#universal-functions))"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random(5)\n",
      "print(x)\n",
      "print(np.minimum(x, 0.5))\n",
      "print(np.maximum(x, 0.5))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.18846954  0.66824324  0.10376865  0.43608056  0.36835284]\n",
        "[ 0.18846954  0.5         0.10376865  0.43608056  0.36835284]\n",
        "[ 0.5         0.66824324  0.5         0.5         0.5       ]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# contrast this behavior with that of min() and max()\n",
      "print(np.min(x))\n",
      "print(np.max(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.103768651674\n",
        "0.668243237764\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "# On older IPython versions, use %pylab inline\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(0, 10, 1000)\n",
      "plt.plot(x, np.sin(x));"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VWW2P/DvIqFLEUE6NsBBRQFRKYNE2o1IgEAUsWBB\nxVGujk7F8f7keh11vDr4U3REsTDSpSi9G2EAUUA6SIdQpIog0kLe+8c6GSKmnbPP2e8u38/z8JCE\nfc5eJOesvHvt912vGGNAREThUMJ2AERE5B4mfSKiEGHSJyIKESZ9IqIQYdInIgoRJn0iohBxnPRF\n5AMR2Sciqws55g0R2SQiK0WkqdNzEhFRbOIx0v8QQGpB/yginQHUN8Y0APAIgH/E4ZxERBQDx0nf\nGLMAwPeFHNIVwLDIsUsAVBaR6k7PS0RE0XOjpl8bQFaez3cBqOPCeYmI6Dxu3ciV8z5n7wciIguS\nXTjHbgB183xeJ/K1nxER/iIgIoqBMeb8gXWB3BjpTwLQBwBEpAWAI8aYffkdaIzxzZ+9ew169jRo\n2NBgzBiD7OzCj9+3z+DPfza46CKDN94wyMkp+NjnnnvO+v/PK3/4vQj29+LsWYO//13fF88+a3Dg\nQOHHnzljMHKkQeXKz6FXL4P9++3/H2z/iZbjkb6IjALQFkBVEckC8ByAkpEkPsQYM01EOovIZgDH\nATzg9Jy2LV4MZGQA994LDB8OlClT9GMuvhh46SXg/vuBPn2AzEzgww+BihUTHS2RNx05ou+hw4f1\nPdWgQdGPSU4GevcG1qwBTp0CmjcHxo0Dbrgh8fEGheOkb4zpXYxj+js9j1dMmAA8+ijw/vtAWlr0\nj7/ySmD+fOCJJ4A2bYBZs4DqnMtEIbN3L9ChA3DLLcD48UCpUtE9vmRJ4K9/BVq3Bm67Lfb3Yxhx\nRW4Uxo4FHn9cE7WTF1jp0sA77wA9emji333eHY6UlBRHcQYJvxfnBOV7kZUFtG0L3H03MHhw9Akf\nOPe9SE8Hpk4FHn5Y359UNImlJpQIImK8Ekt+Jk0C+vXThN+4cfye96WXgFGjdPRfuXL8npfIiw4e\nBFq1Ah55BPj97+P3vKtWAR07arm1Y8f4Pa8fiAhMFDdymfSLYflyIDUVmDZNa4jxZAzw5JP6op0x\no3j3B4j86ORJLem0aaODnXhbsADo2VMHZk2axP/5vYpJP8727AFuugl44w29lEyEs2eB22/Xm73v\nvJOYcxDZZIxOYDh1Chg9GiiRoMLyJ5/oFcSyZUDVqok5h9dEm/RZ0y/E2bNad3z44cQlfABISgI+\n+giYNw/4+OPEnYfIlqFDgZUrgWHDEpfwAR089eqls4JychJ3Hj/jSL8Qzz8PfPGFXi4mJSX+fKtX\nA+3a6Tmvuirx5yNyQ+7resEC4Fe/Svz5zpwB2rcH/uM/gL/8JfHns43lnThZtEjrg8uWAbVquXfe\nd97R6WeLF+ucZCI/O3UKuP564A9/AO67z73z7t4NNG0KzJypfwcZyztxcPIk0Lcv8Oab7iZ8QGcI\nXXgh8Le/uXteokR48UWgfn2t57updm3gtdd0MeTp0+6e2+s40s/Hf/2XrvibMAGQYv/+jJ+sLKBZ\nM+Dzz4FrrnH//ETxkFvWWbFCk7DbjAG6ddOZPM8/7/753cLyjkOrVmk9cOVK90f5eb39NjBmjLZr\nsPGLh8iJs2eBli11EsTDD9uLY88e4NprtVzbsKG9OBKJ5R0HjAEeewx44QW7CR/QMs+xY7pwi8hv\nhg7VNScPPWQ3jlq1gAEDdC2MB8aUnsCRfh5jx+qikaVL3ZmtU5Tcxm7r17MxG/nHkSM6S2f6dG/c\nRD1zBrjuOuDll4GuXW1HE38s78ToxAmgUSOdL++lFicPPABUqwa88ortSIiK5/e/18Q/dKjtSM6Z\nM0dbP6xbF7xV70z6MXrxRR3hT5hgLYR87dmjvX5WrADq1i36eCKbNm3SWv7atd7rHtutmw7onnrK\ndiTxxaQfgwMH9HJ0yRKdXuY1zzwD7Nun8/eJvOyOO7SkM2CA7Uh+ac0anaSxaVOwyqVM+jH4wx+A\n48d1xowXHTmiMw8yM7lSl7xr5UpdBbtlC1C+vO1o8nf//UC9esGawsmkH6XvvtNEunq1nbnExfXq\nq8DChcDEibYjIcpf9+7aJ9/L5ZPt23WF8Lp13is/xYpJP0pPPqnz4F9/3fVTR+XECS09TZnijRkR\nRHktW6YzYzZvBsqWtR1N4Z54Qnfeeu0125HEB5N+FHbt0qlca9cCNWq4euqYDBqko/1x42xHQvRz\nt90GdO6sO8t53a5dumBr48ZgtF9m0o/CU09pm1e//MY/fhy4/HJtwXz11bajIVJff63NCTdt0q1A\n/eCRR7S88z//YzsS55j0i+nQIaBBA+/X8s/30ks6C2HECNuREKmMDN0N68knbUdSfFu26OZIW7YA\nlSrZjsYZJv1iev55YMcO/02DPHpUR/uLF+svLSKbNm3SPW+3b/fujJ2C3HOPXjF7cXppNJj0iyG3\nTPLFF+5s6hBvzz4LfP898NZbtiOhsHv0Ud3m049TINeu1S6g27YB5crZjiZ2TPrF8Oab2rbYa6tv\ni2vvXp1mumULUKWK7WgorHKnO3/7rbYK8aOuXYEuXbTG71fsslmE7Gy9cfunP9mOJHY1a+qS8iFD\nbEdCYfbmm0Dv3v5N+ADw29/qdG2PjH1dEbqkP2mS3ri96SbbkTjz1FPA4MHcFYjs+Okn4N13gaef\nth2JM7fconP2Z8+2HYl7Qpf0Bw8G+ve3HYVz112n9yPGjrUdCYXRqFFAixbAFVfYjsQZER3tDxpk\nOxL3hCrpr1kDbNigc4qD4Omngb//PVyXpmSfMcEZPAFaolq+XPetCINQJf3Bg3VHqlKlbEcSH7fe\nqrtrffml7UgoTBYt0hlwHTvajiQ+ypTRWUhvvGE7EneEZvbO99/rNM116/RGaFC8+qouMBs2zHYk\nFBZ33aX3xPy0GKsouTPiduzwX9tlTtkswKBBulx85MiEncKKgwe1EduWLcBFF9mOhoJu715d0LR1\nK1C5su1o4isjA+jQQUf9fsIpm/nIydGFTEGpQeZVtSqQlsaRPrnjvfeAXr2Cl/ABLf2+807w75GF\nIul//rkuEW/Z0nYkifHoo+F4sZJd2dm6NuSxx2xHkhjt2+s9sq++sh1JYoUi6Q8dCjz0kE7PCqJW\nrfRm1Lx5tiOhIJsxQ3edatzYdiSJUaKErswN+qLHwNf0Dx3SucRbtwa7ZcHbb+sVzSef2I6Egio9\nXVsW9O1rO5LE2b8fuPJK7cfjlxIWa/rnGTFCN3gIcsIHtGPgnDm6gTpRvH33ne7RfMcdtiNJrIsv\n1n1+hw+3HUniBDrpG6M3noI8MslVsaL242GffUqEYcN0UWOFCrYjSbx+/bTFhEeKIHEX6KT/9dfa\nIyQlxXYk7njgAeDDD4P7YiU7jNF9J8IweAJ0c/djx4AVK2xHkhiBTvq5L9QSgf5fntOmja6UXLbM\ndiQUJAsWaFOyFi1sR+KOEiWAPn2Ajz6yHUliBPZG7o8/AnXrar8dP22H6NTzz+vNqMGDbUdCQdGn\nD9Ckif87akZjyxad4r1rl/fbtvBGbsT48UDr1uFK+ABw333A6NHAyZO2I6EgOHZM25Hfe6/tSNx1\nxRXaxXbaNNuRxF9gk/7HH2sCDJtLLtFR2aRJtiOhIJg4UcuGft4oJVb33x/Mle6BTPq7d2ur1C5d\nbEdiR+4NXSKnPv44fKP8XBkZuvblwAHbkcRXIJP+yJFAjx5A2bK2I7EjPV3bLe/dazsS8rM9e4Cl\nS7W3UxhVrKgDx1GjbEcSX4FM+sOH62KlsCpXTufsjxljOxLys7APngAt8QRtFk/gkv6qVdo7/+ab\nbUdi1913c6EWORP2wROge+geOKCzAIMicEl/+HBNeGGZm1+Qdu10utnGjbYjIT9avVr7VrVtazsS\nu5KSgDvvDFaJx3FqFJFUEdkgIptE5E/5/HuKiPwgIt9E/jzr9JwFOXtWR7dhvfGUV1KS9j0P2qYx\n5A4Ons7p3VuTvkeWNDnm6EcqIkkABgNIBXAVgN4i0iifQ78wxjSN/HnByTkLk5kJVK+u257RuRJP\nUF6s5I6cHH3dhL20k6tpU12RvGSJ7Ujiw+nv8RsBbDbGbDfGnAEwGkC3fI5zpZM9a5A/17y5/r10\nqd04yF8yM3Ve/jXX2I7EG0R0X+CglHicJv3aALLyfL4r8rW8DIBWIrJSRKaJSELG4adOAZ9+qvU3\nUiK8oUvRGzFCXzd0Tu/ewNixWkL2u2SHjy9O4WA5gLrGmJ9E5FYAnwJomN+BAwcO/PfHKSkpSImi\nPebMmcC11wK1ahX7IaFw1106k+nVV4Fkpz9tCrzTp3XwlOetSAAaNtTckpmp2yralJmZiczMzJgf\n76jhmoi0ADDQGJMa+XwAgBxjzN8Kecw2ANcbYw6f93VHDdfuvlu3DXz88ZifIrBuuAF48UWgY0fb\nkZDXTZ8OvPACsHCh7Ui859VXgQ0bdPtVL3G74dpSAA1E5FIRKQWgF4CfdX0RkeoiujutiNwI/UVz\n+JdPFbsTJ4CpU3XZNP3SnXfqpSlRUcaMCf7uWLHq1Ut7EZ06ZTsSZxwlfWNMNoD+AGYCWAdgjDFm\nvYj0E5F+kcMyAKwWkRUAXgcQ96r79OnA9dfrzB36pYwMfbGeOWM7EvKyU6e0UR8HT/mrWxe4+mrd\nIN7PHFd5jTHTAUw/72tD8nz8FoC3nJ6nMGPHcnRSmEsuAerX1+ZRnTrZjoa8avZsnbETtnbk0cid\ns98tvzmKPuH7pRfHj+tv3p49bUfibbffDnzyie0oyMtY2ilaRoZWFn76yXYksfN90p86FbjpJqBq\nVduReBtLPFSYkyeBKVM4eCpKtWq6/mXmTNuRxM73SZ+lneLJW+IhOt/Mmbr5Ts2atiPxvowMYNw4\n21HEztdJ/9gxrUOmp9uOxB9Y4qGCsLRTfOnpWmHw65akvk76kycDv/41UKWK7Uj8ISNDF96wxEN5\nnTihe8GytFM8NWoA112nA04/8nXSZ2knOpdcohs+s8RDeU2frnXqiy+2HYl/+LnE49ukf/w4MG8e\n0LWr7Uj8hSUeOt8nn+jrgoqvRw+tNJw+bTuS6Pk26c+YAbRsCVx4oe1I/IUlHsrr5Ekd6XfvbjsS\nf6ldG2jUCJg713Yk0fNt0p8wgTdwY8ESD+U1d642KuRq9uj5tcTjy6R/6pTeePLzqjibevbUX5pE\nEydy8BSrnj2Bzz7z31WzL5P+vHnaA4NzimOTnq4v1pwc25GQTdnZ2muHST829erpVbODLsdW+DLp\nc3TiTP36uoL5yy9tR0I2LVwI1KkDXHqp7Uj8y48lHt8l/bNndZTKpO9Merr+8qTwmjBBZ6FQ7Hr2\n1PeRn3bU8l3SX7hQd7C5/HLbkfhbbtLnpunhZAyvmOPh8su1zLx4se1Iis93SZ+zduKjSROt6a5Z\nYzsSsmHZMqBsWeCqhOxYHS7p6ToN2i98lfRzRye8JHVOhCWeMMst7UixN9mjgnTvrknfL1fNvkr6\ny5cDpUvrzB1yjkk/vHjFHD/XXadXzWvX2o6keHyV9HNfqBydxEfr1sDu3cD27bYjITetX69tTJo3\ntx1JMIicG+37ge+SPks78ZOUBKSl+efFSvExYYImqRK+evd7G5N+Aqxfr/3zb7jBdiTBwhJP+HDw\nFH+//rVeMWdl2Y6kaL5J+rnTyzg6ia8OHYCVK4EDB2xHQm7YsQPYuRNo08Z2JMGSnAx06aJriLzO\nNyl04kR2AkyEMmWATp10OT4F32efAbfdpkmK4ssvJR5fJP09e4AtW4Cbb7YdSTCxxBMekyezUWGi\ndOoEfPUV8P33tiMpnC+S/pQpQGoqULKk7UiCqXNnYP58vWdCwfXDD8CSJUDHjrYjCaZy5YB27XT/\nXC/zRdKfNIk7ZCVSpUpAq1bArFm2I6FEmjFDa/kXXGA7kuDyQ4nH80n/+HEdhaam2o4k2NLSWNcP\nOg6eEq9LF90w/cQJ25EUzPNJf84cnaZZubLtSIItLU0vS/3ULZCK78wZ3RaxSxfbkQRb1apA06be\n3kbR80mfoxN31KsH1K3rr26BVHz/+pdu+FG7tu1Igs/rJR5PJ/2zZ/Umblqa7UjCoWtXlniCavJk\nvo/c0q2bfr+9etXs6aT/1VdAtWrsne8WJv1gMoZXzG667DLg4os1f3mRp5P+5Ml8obqpWTOdtvnt\nt7YjoXhavx44fVq7QZI7unbV/OVFnk76HJ24S0RLAF59sVJsct9H7E7rHi+/jzyb9LduBQ4eBG68\n0XYk4cIST/Bw8OS+G28E9u8Htm2zHckveTbpT56s08vYYM1d7dppA7ZDh2xHQvGwbx+wbh3Qtq3t\nSMKlRAntceTF0b5nU+qkSZxtYEOZMkD79sC0abYjoXiYOlV7wpQubTuS8PFqXd+TSf/774Gvv9a2\nv+Q+rs4NDk7VtKdjR+119MMPtiP5OU8m/Rkz9HK0fHnbkYTTbbfpUvJTp2xHQk6cOKErQzt3th1J\nOJUvr5urzJxpO5Kf82TS51RNuy6+WDefz8y0HQk5MW+etgS46CLbkYSXF2fxeC7pnzmjI332CLGL\ns3j8j7N27EtL055H2dm2IznHc0l/wQKgfn2gZk3bkYRb7k0oY2xHQrHIyeEVsxfUqaN9rRYtsh3J\nOZ5L+hydeMOvfgWUKqXTN8l/li3TfRIaNLAdCXmtxOOppM8eId4hwhKPn3GU7x1M+oVYt0470zVu\nbDsSApj0/YyDJ+/I7Wm1caPtSJSnkj57hHhL69a6jHz3btuRUDR27NCfWYsWtiMhQFfndunindG+\np5I+L0m9pWRJ3abSKy9WKp7Jk3WtRVKS7Ugol5dKPJ5K+uvXs0eI13h1KTkVjKUd72nfHvjmG+Dw\nYduReCzpd+yoM0bIO1JTdRrtjz/ajoSK44cfdMvLTp1sR0J5lS0LpKTonH3bHCd9EUkVkQ0isklE\n/lTAMW9E/n2liDQt6Lk4OvGeSpW0Tezs2bYjoeKYORNo0wa44ALbkdD5vFLicZT0RSQJwGAAqQCu\nAtBbRBqdd0xnAPWNMQ0APALgHwU9H3uEeBNLPP7BBmveddtt+kv5zBm7cTgd6d8IYLMxZrsx5gyA\n0QC6nXdMVwDDAMAYswRAZRGpnt+TVaniMBpKiLQ03aDeqxs9k8rO1pbYTPreVLMm0LChlkttcpr0\nawPIyvP5rsjXijqmjsPzkosuuwyoUcO7Gz2TWrQIuOQSXfpP3uSFtuXJDh9f3M4s58+8z/dxAwcO\n/PfHKSkpSElJiSkoir/cF2vLlrYjoYJw4yHvS0sDevQABg2KfT1SZmYmMh20wBXjoKOWiLQAMNAY\nkxr5fACAHGPM3/Ic8w6ATGPM6MjnGwC0NcbsO++5jJNYKLGWLAEefBBYu9Z2JFSQK68ERo4Err/e\ndiRUEGP0amzGDOCqq+LznCICY0yxf4U4Le8sBdBARC4VkVIAegE4/+JlEoA+keBaADhyfsIn77vh\nBt03d8sW25FQfjZu1Gm1zZrZjoQK44WeVo6SvjEmG0B/ADMBrAMwxhizXkT6iUi/yDHTAGwVkc0A\nhgB4zGHMZIHXlpLTz02erD8ftjDxPttTNx2Vd+KJ5R3vmzQJeP113ZGJvKVtW+CPf9RpgeRtp04B\n1asDmzYB1ao5fz63yzsUIh06AEuX6sb15B2HDwMrVgDt2tmOhIqjdGntPjB1qp3zM+lTsZUrpyPK\nGTNsR0J5TZ+uS/zLlrUdCRWXzambTPoUFds3oeiXOFXTfzp3BubOBU6edP/cTPoUlS5dvLGUnNTp\n08CsWfpzIf+oWhW47jo798eY9CkqNWvqxvW2l5KTWrBAl/bXqGE7EoqWratmJn2KGks83sHSjn/l\nNjLMyXH3vEz6FLXcpM8ZtnYZw93m/KxhQ6BCBWD5cnfPy6RPUWvcWDturltnO5JwW7dOfw6NG9uO\nhGJlo205kz5FzQtLyelc73yuwvUvG+8jJn2KCTdWsY/1fP9r2RLYtQvYudO9czLpU0zattXywj62\nzrNi/379/rP7uL8lJWnrDDcHUEz6FJNSpXTzbVtLycNu2jRti1G6tO1IyCm3SzxM+hQzlnjsYWkn\nODp1AhYvBo4eded87LJJMTt0SLdS3LePfV/cdPKkdmncvDk+XRrJvltv1U2Kbr89+seyyya55qKL\ngKZN2WrZbZmZOk2TCT843CzxMOmTIyzxuI+lneBJS9P7NNnZiT8Xkz45krsLkNtLycPKGGDKFCb9\noKlTB7j0UmDRosSfi0mfHLG1lDysVq7UmVONGtmOhOLNrRIPkz45xhKPez77jKtwg6prV/35Jno+\nC5M+OWZzF6Cw+fRTID3ddhSUCE2a6Mysb79N7HmY9Mmxli2BrCx3l5KH0fbtumS/VSvbkVAiuNXT\nikmfHEtO1qXkU6bYjiTYcks7ycm2I6FEYdIn32CJJ/E+/RTo3t12FJRIKSnA6tXAgQOJOwdX5FJc\nHDsG1K4N7N6ts3kovrj6OTwyMnTP4/vvL97xXJFLVlSooLXmWbNsRxJMU6ZogzUm/OBLdImHSZ/i\nhiWexGFpJzw6dwbmzgVOnEjM8zPpU9x07aqtlt1YSh4mP/2k/Y26dLEdCbmhalWgWTNg9uzEPD+T\nPsVN3brA5ZcD8+fbjiRYZs8GmjcHqlSxHQm5pUcPYMKExDw3kz7FVXp64l6sYcXSTvh07673cc6c\nif9zM+lTXPXoAUycyAZs8ZKdrW/+bt1sR0JuSuRVM5M+xdWVVwKVKwNffWU7kmBYuBCoV0//ULgk\nqsTDpE9xl8h6ZNiwtBNePXrozz/eV81M+hR3uUmfa+2cMYZJP8waNkzMVTOTPsVdkybA2bO6nJxi\nt2oVUKIEcM01tiMhWxJx1cykT3Enwlk88TBxoo7y2Ts/vHInRsTzqplJnxIi98VKsRs3Drj9dttR\nkE1Nmui0zTVr4vecTPqUEC1banOwzZttR+JP69cDP/wA3Hij7UjIJpH4l3iY9CkhkpK0NMHRfmzG\njwd69tSaPoUbkz75Bqduxu6TT7TFLlHuVfOWLfF5PiZ9SpiUFN3vc/du25H4y8aNuokGt0UkQK+a\nu3WL31Uzkz4lTKlSuo0iSzzRGT9er5JY2qFcPXvqjf144MuKEiojQ0sVVHzjxrG0Qz93yy1a3tm+\n3flzMelTQqWm6iKjPXtsR+IPW7cCu3YBbdrYjoS8pGRJXfsSjwEUkz4lVOnSurlKvC5Ng278eH1z\nJyXZjoS8plcvYOxY58/DpE8Jd8cdwJgxtqPwB5Z2qCBt2wI7dzqfxcOkTwnXsSOwYQOQlWU7Em/b\nsUPLO23b2o6EvCg5WW/oOi3xMOlTwpUqpQu1eEO3cOPH69S8kiVtR0Jedccdzks8MSd9EakiIrNF\nZKOIzBKRygUct11EVonINyLCrTVCqlcvlniKMmoUcOedtqMgL2vTBti7F9i0KfbncDLS/zOA2caY\nhgDmRj7PjwGQYoxpaoxhJ5GQatdOSxfxmHIWRJs2afnrlltsR0JelpSk93ycjPadJP2uAIZFPh4G\noLCtHtgcNuSSk3XBUTxmHwTRqFF66c5ZO1QUpxMjnCT96saYfZGP9wGoXsBxBsAcEVkqIg87OB/5\nHEs8+TNGk37v3rYjIT9o3Ro4dEg7scYiubB/FJHZAGrk809/yfuJMcaISEFt/lsbY/aKSDUAs0Vk\ngzFmQX4HDhw48N8fp6SkICUlpbDwyGfattU+PJs3A/Xr247GO1auBE6eBFq0sB0J+cH8+ZmoWzcT\nv/mN9reKlpgYt2QRkQ3QWv13IlITwOfGmF8V8ZjnAPxojHktn38zscZC/tG/P1CzJvCXvxR9bFj8\n8Y9a/nrxRduRkF98+SVw//062i9RQmCMKXYJ3Ul5ZxKA+yIf3wfg0/MPEJFyIlIh8nF5AJ0AcOfU\nELvrLmDECG6anisnBxg9mqUdis5NNwGnTwPffBP9Y50k/ZcBdBSRjQDaRT6HiNQSkamRY2oAWCAi\nKwAsATDFGDPLwTnJ51q2BE6dApYvtx2JNyxaBFSsCDRubDsS8hMRHUCNHBnDY71SUmF5Jzyeew44\nehQYNMh2JPY9/jhQqxbLXRS9det0tfuePdGVd5j0yXWbNukik127tJYdVtnZmvAXLwauuMJ2NORH\nH3wA9O3rXk2fKCYNGgCXXgrMmWM7ErtmztRkz4RPsXrwwegfw6RPVtx7L/Dxx7ajsOujj3QGBpGb\nWN4hKw4e1Ln6WVlAhQq2o3Hf4cPAZZdpZ83K+XatIioeEZZ3yAeqVgVuvjm8++eOHg3ceisTPrmP\nSZ+sCXOJh6UdsoXlHbLm5EmgTh1g2TLgkktsR+Oe3Kl2O3eywRo5x/IO+UaZMto//qOPbEfirmHD\ngHvuYcInOzjSJ6tWrNDdorZuDUcSzM4G6tXT6apXXWU7GgoCjvTJV5o00Zu6c+fajsQdc+ZoSYsJ\nn2xh0ifr+vYF3n/fdhTueO+92BbUEMULyztk3ZEjukJ382Yd9QfVnj3A1VfrDdwwrk2gxGB5h3yn\ncmUgLQ0YPtx2JIn1wQe61R0TPtnEkT55whdfaMfJ1au1bWzQnD0LXH65LkZr1sx2NBQkHOmTL918\nsybGBflupOl/M2YA1asz4ZN9TPrkCSI60h882HYkiTFkCNCvn+0oiFjeIQ85elRv6K5eDdSubTua\n+MnK0qmpO3cC5cvbjoaChuUd8q2KFXULuCFDbEcSX+++q3vgMuGTF3CkT56yfj3Qrp22HC5VynY0\nzp08qX2F5s8HrrzSdjQURBzpk681aqSrVcePtx1JfIwYATRvzoRP3sGkT57Tvz/w5pu2o3DOGOD1\n14Hf/tZ2JETnMOmT56SlAXv36obhfjZvnib+Dh1sR0J0DpM+eU5yMvC73wH/+7+2I3Emd5QfxMVm\n5F+8kUuedPy47iG7YIE/6+EbNuiCsx07gLJlbUdDQcYbuRQI5csDjz0GvPaa7Uhi8/LLwBNPMOGT\n93CkT57aXgtqAAAG8UlEQVR18CDQsKFuL1ijhu1oim/bNp2xs2ULNz6nxONInwKjalXg7ru1Nu4n\nr7yiLReY8MmLONInT9u5E2jaVGvk1arZjqZoe/YA11wDfPutP+Il/+NInwKlXj3dPP2VV2xHUjyv\nvgr06cOET97FkT553u7dwLXXAmvXeru2v3s30Lhx8BrGkbdFO9Jn0idfeOqpcytcvapfP6BSJf9c\nlVAwMOlTIH33ne4vu3IlUKeO7Wh+aeNGoFUr/btKFdvRUJiwpk+BVKOGjqSffdZ2JPl79lldRcyE\nT17HkT75xrFjujr3s8+AG26wHc05S5YAPXroKJ8988ltHOlTYFWoALzwgvaz8cr44OxZ3ebx5ZeZ\n8MkfmPTJV+67DzhxAhgzxnYk6oMPgDJlgHvusR0JUfGwvEO+869/Ab166RROm6teDx/WTV9mztQ9\ncIls4OwdCoXf/AbIzgbee89eDA89BJQuDbz1lr0YiJj0KRSOHtUpnP/8J3DLLe6ff8YM4NFHgVWr\ndEN3Ilt4I5dCoWJF4O23gYcf1lk9bjpyRM87dCgTPvkPR/rka337AqdP64jfjR2qjAHuugu48EL9\npUNkG0f6FCpvvAEsXw4MG+bO+f7xD+3v79fNXYg40iffW7sWSEnROvv11yfuPF99BXTpAixcCDRo\nkLjzEEWDI30KnauvBoYMAbp1A7KyEnOObduA9HTg3XeZ8Mnfkm0HQBQPPXoAW7fqSPzzz+PbA+fg\nQSA1FXjmGaB79/g9L5ENHOlTYPzud0CnTkD79pqo42H/fqBjRyAjQ9stEPldzElfRG4XkbUiclZE\nmhVyXKqIbBCRTSLyp1jPR1QUEe1lf+utOnd/xw5nz7dtG3DzzUDXrtrzhygInIz0VwNIBzC/oANE\nJAnAYACpAK4C0FtEGjk4ZyhkZmbaDsEzov1eiAB//SvwwANAixbA3LmxnXfGDH18//7Af/+3O9NB\ni8LXxTn8XsQu5qRvjNlgjNlYxGE3AthsjNlujDkDYDSAbrGeMyz4gj4nlu+FCPD008Dw4dqgrV8/\n4NCh4j32wAHgwQd18dW4cZr0vYKvi3P4vYhdomv6tQHknU+xK/I1ooRr316ncyYn64ybp58Gli79\nZVvmnBydjvmf/6n9+itW1Ln4bdrYiZsokQqdvSMiswHktxX1M8aYycV4fk68J6sqVdKGaAMG6N/3\n3KM3Zxs21OR+7JhuflKjhk7JXLMGqFXLdtREieN4cZaIfA7gd8aY5fn8WwsAA40xqZHPBwDIMcb8\nLZ9j+QuCiCgG0SzOitc8/YJOuBRAAxG5FMAeAL0A9M7vwGiCJiKi2DiZspkuIlkAWgCYKiLTI1+v\nJSJTAcAYkw2gP4CZANYBGGOMWe88bCIiioVneu8QEVHiWV+Ry8VbSkTqisjnkQVva0TkCdsx2SYi\nSSLyjYgUZ9JAYIlIZREZJyLrRWRd5F5ZKInIgMh7ZLWIjBSR0rZjcouIfCAi+0RkdZ6vVRGR2SKy\nUURmiUiRG4haTfpcvPUzZwA8ZYy5GloyezzE34tcT0LLgmG/HP3/AKYZYxoBuBZAKEukkXuDDwNo\nZoxpDCAJwJ02Y3LZh9BcmdefAcw2xjQEMDfyeaFsj/S5eCvCGPOdMWZF5OMfoW/s0E4eFJE6ADoD\nGIqCJwoEnohUAtDGGPMBoPfJjDE/WA7LlqPQwVE5EUkGUA7AbrshuccYswDA9+d9uSuA3N0khgEo\nsiWg7aTPxVv5iIxomgJYYjcSqwYB+AOAHNuBWHYZgAMi8qGILBeR90SknO2gbDDGHAbwGoCd0NmA\nR4wxc+xGZV11Y8y+yMf7AFQv6gG2k37YL9t/QUQuADAOwJOREX/oiEgXAPuNMd8gxKP8iGQAzQC8\nbYxpBuA4inEJH0QicgWA3wK4FHoVfIGI3G01KA+J7EJVZE61nfR3A6ib5/O60NF+KIlISQDjAQw3\nxnxqOx6LWgHoKiLbAIwC0E5E/mk5Jlt2AdhljPk68vk46C+BMGoOYJEx5lBkOvgE6GslzPaJSA0A\nEJGaAPYX9QDbSf/fi7dEpBR08dYkyzFZISIC4H0A64wxr9uOxyZjzDPGmLrGmMugN+rmGWP62I7L\nBmPMdwCyRKRh5EsdAKy1GJJNGwC0EJGykfdLB+iN/jCbBOC+yMf3AShysGh15yxjTLaI5C7eSgLw\nfogXb7UGcA+AVSLyTeRrA4wxMyzG5BVhLwP+J4ARkYHRFgAPWI7HCmPMysgV31LovZ7lAN61G5V7\nRGQUgLYAqkYWxv4/AC8DGCsifQFsB3BHkc/DxVlEROFhu7xDREQuYtInIgoRJn0iohBh0iciChEm\nfSKiEGHSJyIKESZ9IqIQYdInIgqR/wM+68FI5Ooc3wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106a02cd0>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.special import gammaln\n",
      "plt.plot(x, gammaln(x));"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHORJREFUeJzt3Xu8lXPe//HXhxxS5Dg5ZTqMQyWHbocSWmhIIzxoamqG\nyHkcMz8jjNrD/RPm5sYYZpAKiYimkIRWlCKiUko1DiXtjM1WOtn1uf/47rQ1HdZep2uta72fj8d6\nWIdrretjqbfv+l7fg7k7IiJS3LaKugAREcmcwlxEJAYU5iIiMaAwFxGJAYW5iEgMKMxFRGIgpTA3\ns0fNrNzMZmzktT+Y2Voz2zX75YmISCpSbZkPBDpu+KSZNQJ+CXyWzaJERKR2Ugpzd38T+GYjL90N\n/DGrFYmISK2l3WduZmcAC919ehbrERGRNNRJ501mtgNwI6GL5cens1KRiIjUWlphDjQDGgPTzAxg\nX+A9MzvK3ZfUPNDMtPiLiEga3D3lRnJa3SzuPsPdG7p7E3dvAiwEWm8Y5DWO182dfv36RV5Dodz0\nXei70Hex+VttpTo0cSjwFnCAmS0ws/M3zOtan1lERLImpW4Wd+++hdebZqccERFJh2aA5lEikYi6\nhIKh72I9fRfr6btIn6XTN1OrE5h5rs8hIhI3Zobn+gKoiIgUFoW5iEgMKMxFRGJAYS4iEgMKcxGR\nGFCYi4jEgMJcRCQGFOYiIjGgMBcRiQGFuYhIgZk9u/bvUZiLiBSQYcOgffvav09hLiJSIB5+GHr3\nhrFja//edHcaEhGRLLrzTnjwQRg/Hn7xi9q/X2EuIhIhd7jxRvjnP2HCBNhnn/Q+R2EuIhKRNWvg\n8sth6lR44w3Yfff0P0thLiISgdWr4dxzYckSeO012HHHzD5PYS4ikmfLl0OXLrDNNvDSS7D99pl/\nZqobOj9qZuVmNqPGc38xs4/MbJqZPWdmDTIvR0Qk3ior4ZRTQpfK8OHZCXJIfWjiQKDjBs+9ArR0\n90OBj4EbslOSiEg8lZfDCSfA4YfDoEFQJ4t9IymFubu/CXyzwXNj3X1t9cO3gX2zV5aISLx88gkc\neyyccQbcey9sleVZPtn6uF7AS1n6LBGRWJk+HY47LkwI6tcPLOVtmlOXcSPfzG4CVrv7k1moR0Qk\nViZMgLPPhr/+Fbp2zd15MgpzMzsP6ASctLnjysrKfryfSCRIJBKZnFZEpCiMGgUXXABDhsAvf7n5\nY5PJJMlkMu1zmbundqBZY2CUu7eqftwRuAto7+7/3sz7PNVziIjExeDBcP31MHIkHHVU7d9vZrh7\nyh0yKYW5mQ0F2gO7A+VAP8LolW2BiurDJrn77zfyXoW5iJSUu+6C++6DMWPgoIPS+4ychHkmFOYi\nUircoU+f0L0yZgw0apT+Z9U2zDUDVEQkC6qq4JJLYOZMePNN2G23/J5fYS4ikqEVK6B79/DPV1+F\n+vXzX4M2pxARyUBlJXTsCHXrhu6VKIIcFOYiImlbvDhs8XbIIWH44bbbRleLwlxEJA3z54fp+Wed\nFUauZHt6fm2pz1xEpJamTYNOneBPf4LLLou6mkBhLiJSC+PGQbducP/9uZ2eX1vqZhERSdHTT4cg\nf/rpwgpyUMtcRCQl994Lf/kLjB0Lhx4adTX/SWEuIrIZa9eGWZ0jR8LEifDzn0dd0cYpzEVENmH1\naujVC/71rxDk+Z7VWRsKcxGRjVi6NKxDXrdumNW5ww5RV7R5ugAqIrKBxYshkYAmTcKmy4Ue5KAw\nFxH5iblzoV27sFfn3/+e3U2Xc6lIyhQRyb133gkhfuutcOGFUVdTOwpzERFg9Gjo2RMGDIDOnaOu\npvbUzSIiJW/QIDj//DD8sBiDHNQyF5ES5g633QaPPALjx8OBB0ZdUfoU5iJSktasgauuCuPHJ06E\nvfeOuqLMpNTNYmaPmlm5mc2o8dyuZjbWzD42s1fMbOfclSkikj3Ll4cx5LNnhxZ5sQc5pN5nPhDo\nuMFzfYCx7n4A8Fr1YxGRgrZkCZxwAjRoEC56NmgQdUXZkVKYu/ubwDcbPH06MLj6/mDgzCzWJSKS\ndXPmQNu2YZu3QYOi3Rko2zLpM2/o7uXV98uBhlmoR0QkJyZMgC5dwgXPXr2irib7snIB1N3dzHxT\nr5eVlf14P5FIkEgksnFaEZGUPPMMXH45PP44nHJK1NVsXDKZJJlMpv1+c99kBv/0QLPGwCh3b1X9\neDaQcPfFZrYXMM7dD9rI+zzVc4iIZJM73HVXWIt81Cg47LCoK0qdmeHulurxmUwaGgn0rL7fExiR\nwWeJiGTVuqGHjz0Gb71VXEGejpRa5mY2FGgP7E7oH+8L/BMYBuwHfAp0dfdvN/JetcxFJK++/x56\n9Aj/HD68OEes1LZlnnI3S7oU5iKST+XlYUp+ixbw0EPFO2Iln90sIiIFZd3Qw06dYODA4g3ydGg6\nv4jEwrqhh/37h0WzSo3CXESK3rBhcMUV8MQTcPLJUVcTDYW5iBQtd7j9dnjwQRg7Fg49NOqKoqMw\nF5GitHo1XHopTJsGkyfHY7GsTCjMRaToVFSEVQ932gneeAPq1Yu6ouhpNIuIFJV58+CYY6B1a3ju\nOQX5OgpzESkaEybAscfCNdeEafpbbx11RYVD3SwiUhSGDIHevUt7xMrmKMxFpKC5w5//DIMHw+uv\nw8EHR11RYVKYi0jBWrkSLrgA5s8PI1YaateETVKfuYgUpK++gg4dwhDEceMU5FuiMBeRgjN7dlhj\n5bjj4OmnoW7dqCsqfApzESkor78O7dvDjTeGdVa2UkqlRF+TiBSMAQOge3d46ql47tOZS7oAKiKR\nq6qC666DF18MMzoPPDDqioqPwlxEIlVZCb/5TQj0t9+GXXaJuqLipG4WEYnMvHnQpg00awYvvaQg\nz4TCXEQiMW4ctGsXNl2+/37YZpuoKypuGYe5md1gZjPNbIaZPWlm22WjMBGJr3/8I3StDB0Kl10W\ndTXxkFGfuZk1Bi4Cmrv7KjN7GvgNMDjz0kQkbqqq4Npr4ZVXwqJZ++8fdUXxkekF0O+AH4AdzGwN\nsAPwRcZViUjsfPstdOsW7k+eDDvvHG09cZNRN4u7VwB3AZ8Di4Bv3f3VbBQmIvExd2640HnQQWH4\noYI8+zLtZmkGXAM0BiqBZ8zst+4+pOZxZWVlP95PJBIkEolMTisiReS116BHD7j1Vrj44qirKVzJ\nZJJkMpn2+83d03+zWTfgl+5+YfXjc4A27n55jWM8k3OISPF64AG45ZYwo1NtuNoxM9zdUj0+0z7z\n2cDNZlYXWAl0AN7J8DNFpMitXg1XXw3jx8PEiWEcueRWRmHu7tPM7DHgXWAtMBV4KBuFiUhxKi+H\nLl3CBKDJk8Omy5J7GXWzpHQCdbOIlIz33oOzzoKePaGsTCseZiLf3SwiIkDYo/Oaa8KEoLPOirqa\n0qMwF5GMrFkDffrAc8+FtchbtYq6otKkMBeRtFVUhPXH16yBd96B3XaLuqLSpR4tEUnLzJlw1FHQ\nsiW8/LKCPGoKcxGptREjwrjxvn3h7ruhjn7jR07/CUQkZWvXhklAjz4a1h8/8sioK5J1FOYikpKl\nS+Hcc2HJktA/vueeUVckNambRUS2aO5caNsW9tgjjFhRkBcehbmIbNbIkWFHoCuuCGPIt9P2MwVJ\n3SwislFr1oRZnIMGhUBv0ybqimRzFOYi8h8qKsKytStXwrvvQsOGUVckW6JuFhH5ifffhyOOCOPH\nX31VQV4sFOYi8qPHHoOTT4b+/eGuuzR+vJjoP5WIsHo19O4NY8fCuHFw8MFRVyS1pTAXKXGLFoX1\nx/fYA6ZMgQYNoq5I0qFuFpES9sYbYRbnr34Fzz+vIC9mapmLlCB3uO8+uO02GDwYOnaMuiLJlMJc\npMQsXQoXXggffwyTJkHTplFXJNmgbhaREjJjRhh22KABvPWWgjxOMg5zM9vZzJ41s4/MbJaZaZ6Y\nSAEaPBhOPBFuugkeegjq1o26IsmmbHSz3Au85O5dzKwOUC8LnykiWbJiBVx5JUyYoGGHcZZRy9zM\nGgDHufujAO5e5e6VWalMRDI2b15Y7XDZsjDsUEEeX5l2szQBvjKzgWY21cweNrMdslGYiGTmuefg\nmGPgootg6FDYcceoK5JcyrSbpQ7QGrjC3aeY2T1AH6BvzYPKysp+vJ9IJEgkEhmeVkQ25Ycf4Prr\nQ5i/8ELYp1MKXzKZJJlMpv1+c/f032y2JzDJ3ZtUPz4W6OPup9U4xjM5h4ikbuFC6NYNdtklrLOy\n665RVyTpMjPc3VI9PqNuFndfDCwwswOqn+oAzMzkM0UkPa+8EmZzdu4c1h9XkJeWjFrmAGZ2KPAI\nsC0wHzi/5kVQtcxFcquqCvr1C5tIDBkC6sWMh9q2zDMemuju0wDt0S0SgQULoHt3qFcvrEP+s59F\nXZFERTNARYrUyJFhNmfnzjB6tIK81GltFpEis2pVGK0yYkRY6fCYY6KuSAqBwlykiMybF0ar7Lcf\nTJ2qi5yynrpZRIrEU0+F2ZznnRfGkCvIpSa1zEUK3PLlcPXVkEzCmDHQunXUFUkhUstcpIDNmhVm\ncC5fHrpVFOSyKQpzkQLkDg8/DO3bw7XXwhNPaG0V2Tx1s4gUmK+/Dotj/etfMH48tGgRdUVSDNQy\nFykgr78Ohx0GjRvD228ryCV1apmLFIDVq+Hmm0N3ysCBcPLJUVckxUZhLhKxjz+GHj1gr73ggw9g\njz2irkiKkbpZRCLiDo88Au3aQa9eYXq+glzSpZa5SAQqKuDii0OrPJmEli2jrkiKnVrmInmWTIaL\nnPvuC++8oyCX7FDLXCRPVq2Cvn3h8cfh0UehY8eoK5I4UZiL5MG0aXDOOdCsWbivvnHJNnWziOTQ\nmjVw++3QoQP84Q9hgSwFueSCWuYiOTJ/PvTsCdtsA+++Cz//edQVSZypZS6SZe7w0ENw9NFw9tnw\n2msKcsm9rLTMzWxr4F1gobt3zsZnihSjxYvhwgth0aKwropGqki+ZKtlfjUwC/AsfZ5I0Rk+PAw5\nPOwwmDxZQS75lXHL3Mz2BToB/x+4NuOKRIrMN9+EzSMmTQr7crZpE3VFUoqy0TL/X+A6YG0WPkuk\nqLzwArRqBTvtFNZVUZBLVDJqmZvZacASd3/fzBKbOq6srOzH+4lEgkRik4eKFIV1rfGJE8NKh/oj\nLZlKJpMkk8m032/u6Xdzm9ltwDlAFbA9sBMw3N3PrXGMr1njbKVxMxITo0bBpZeGkSr9+0O9elFX\nJHFkZri7pXx8JmG+wYnbA/9vw9EsZuYVFc4uu2TlNCKRqahY3zc+YEDY0k0kV2ob5tluL2/0/wzf\nfJPls4jk2ciRoW98113DdHwFuRSarM0AdffxwPiNvVZRAU2bZutMIvlTUQFXXRWGGg4dCscfH3VF\nIhuXl55stcyl2LjDs8+G1vjuu4fWuIJcClle1mapqMjHWUSy44sv4PLLYc4cGDYs7AQkUujUMhep\ntnYtPPjg+lmcH3ygIJfikZeWucJcCt1HH8FFF4VA1zZuUozy0jJXN4sUqtWr4ZZbQn949+4wYYKC\nXIqTWuZSsiZNCq3xJk1g6lRo1CjqikTSl5cw/+qrfJxFJDWVlfCnP4VVDu+5B379a7CUp2aIFKa8\ndLMsWJCPs4hsnnsYK96iBaxcCR9+CF27KsglHvLSMv/883ycRWTT5swJww3//e8wfrxt26grEsmu\nvLTMly2D5cvzcSaRn1qxAvr2DUMMf/WrsBenglziKC9h3qiRulok/0aPhoMPDsMOP/gAeveGOtrC\nXGIqL3+0GzUKXS0HHpiPs0mpW7gQrrkG3n8f/vY36Ngx6opEci8vLfOmTWH+/HycSUrZqlVw551h\n9mbLluECp4JcSkVeWubNm4efuiK5Mnp0WGv8wAPD+PH994+6IpH8ykuYt2gBY8bk40xSaubNC33h\nc+aEMeOdOkVdkUg08tLN0rw5zJqVjzNJqVi2DG64IWygfPzxoUtFQS6lLC9hvt9+sHQpfP11Ps4m\nceYOQ4bAQQfBokUwYwZcdx1su23UlYlEKy/dLFttBUccAe+8A6eemo8zShxNmRK6VFauhGee0Xhx\nkZry0jKH8HN48uR8nU3i5LPP4Le/hTPPhPPPD40CBbnIT2UU5mbWyMzGmdlMM/vQzK7a1LEKc6mt\nysrQL966NfziF+Ei5wUXhF96IvJTmf61+AHo7e4tgTbA5WbWfGMHHn10aFGtXZvhGSX2qqrggQfC\nMMPFi2H6dPjzn6F+/agrEylcGfWZu/tiYHH1/WVm9hGwN/Afo8p/9jPYe2947z048shMzipx5Q4v\nvhguaO69N7z8cpgAJCJblrULoGbWGDgceHtTx3TsGCZ3KMxlQ1OmQJ8+8OWX8D//E4YZamlakdRl\nJczNrD7wLHC1uy/b8PWysjIgbFIxenSCvn0T2TitxMBHH4WNIt5+G26+OfSJazEsKUXJZJJkMpn2\n+83dMyrAzLYBXgBGu/s9G3nd151j1Spo2BBmz4Y998zotFLkPv8cyspg1Cj44x/hiiugbt2oqxIp\nHGaGu6f8+zTT0SwGDABmbSzIN7TddnDGGTBsWCZnlWK2ZElY0fDww0O/+Ny5oY9cQS6SmUxHs7QD\nfgecYGbvV982u05djx5hBp+UlspK6NcvLO2wZg3MnAn//d+w885RVyYSDxl3s2zxBDW6WSAMO9tv\nv7DwVqtWOT21FIDKSrjvvnA79dQwxLBJk6irEil8ee1mSUedOnDppfDXv+b7zJJPlZVw661hss/c\nuTBxIjz2mIJcJFfy3jIHKC8PE0Lmz4fddsvp6SXParbEO3WCm26CAw6IuiqR4lPwLXMII1rOPBPu\nvz+Ks0sufPvt+pb4vHnw1lsweLCCXCRfImmZQ/gL36YNfPwx7LprTkuQHPryy7ApxCOPwGmnhTHj\n2uVHJHNF0TKH0II7+2y4/faoKpBMzJsHl1wS9tpcsQKmTg0tcQW5SDQia5lDaNUdcgiMHx+2lpPC\nN3Uq3HEHvP46XHYZXHkl7LFH1FWJxE/RtMwB9torzAK85BKtpljI3MNQ0lNOgdNPD91jn3wCt9yi\nIBcpFJG2zCFMIDnmGOjZE37/+5yWIrW0fDk8/jjcey9ss02YudmjR5jJKyK5VduWeeRhDuEiaLt2\n8OqrcOihOS1HUrBwIfztb+GiZrt2IcTbt9cqhiL5VFTdLOsccEBo/XXtGoa4Sf65h4k93buH6xgr\nVoSdoUaMgERCQS5S6AqiZb5O797wwQdhUwL9lM+Pykp44gn4+9/DqpaXXQa9ekGDBlFXJlLairKb\nZZ01a6Bbt9A/+8QTsPXWOS2tpL37bgjw4cPh5JPDRegTTlALXKRQFHWYQ/h536kT7LMPDBqkjQqy\n6Ztv4Omn4eGHoaICLr447HavteVFCk/RhzmEURRnnQX16oUWuta6Tl9VVRhWOHjw+uGFvXqF1rh2\nuRcpXLEIcwj9t+edF2YaPv887Ltv9muLsxkzQoAPGQKNG4ehn926wS67RF2ZiKSiKEezbMx228GT\nT0KXLnD00ZDB1nglY86csNhVq1ahq2rbbcP3NmlSWHZYQS4SXwXbMq9pzJjQNdC1K9x2m7pdapo3\nL2zDN2xY2JLt178O31PbtupGESlmselm2dDXX4dNf6dMgbvvhs6dS3Pkxdq14TsYNSrcysvDr5eu\nXeHYYxXgInGR9zCv3vPzHmBr4BF3v2OD17MS5uu8/DJce23YDPiOO+C//itrH12wvv8+zI4dNQpe\nfDF0l3TuHG5t22oIp0gc5TXMzWxrYA7QAfgCmAJ0d/ePahyT1TCHMELjH/+A/v3DEqw33BCv6eY/\n/BBa36++Cq+9Bu+9F64brAvwZs2irlBEci3fYd4W6OfuHasf9wFw99trHJP1MF9n1aowdPHOO8N4\n9AsvhHPOgd13z8npcmblyrC07KRJ4YLlG29A06bQoQOcdBIcd1wYpikipSPfYd4FOMXdL6p+/Dvg\naHe/ssYxOQvzddxDAA4YACNHhiVazzgjLNe6zz45PXWtrV0Ln34aWtuTJoXb9OlhT9S2beH44+HE\nE7W0rEipq22YZzq/MrcpnSKz0M3Svj0sXRpGv4wYETYTbtgwBORxx4WuiqZN89PH7B4uTs6bBzNn\nwrRp4TZjRlj35PDDQ3j37w9HHqmWt4hkJtMw/wJoVONxI2DhhgeVlZX9eD+RSJBIJDI87abtuGMY\n3dGlS1jrZfr00Gp//vmwP+VXX0Hz5nDwwdCkCTRqFCYk7bMP7Lwz7LQT1K+/6f73qipYtizcKipC\nYC9eHG5ffhla3fPnh1vduqF/u0WLsLRv165hDPhuu+XsX19EilQymSSZwYSaTLtZ6hAugJ4ELALe\nIQ8XQDPx3Xehpfzhh/DZZ2Ht7gULYNGisILgd9+F9WHq1Vs/zG9dsK9YES5O1q8f/qfRoEHYLWnP\nPcOtYcMw27JZs3DTyoMikq4ohiaeyvqhiQPcvf8GrxdUmKeiqioMB3QPNwj/rFsXtt8+PqNmRKRw\nxXbSkIhIKYnN2iwiIpI6hbmISAwozEVEYkBhLiISAwpzEZEYUJiLiMSAwlxEJAYU5iIiMaAwFxGJ\nAYW5iEgMKMxFRGJAYS4iEgMKcxGRGFCYi4jEgMJcRCQGFOYiIjGgMBcRiQGFuYhIDCjMRURiIO0w\nN7O/mNlHZjbNzJ4zM+1FLyISkUxa5q8ALd39UOBj4IbslBRfyWQy6hIKhr6L9fRdrKfvIn1ph7m7\nj3X3tdUP3wb2zU5J8aU/qOvpu1hP38V6+i7Sl60+817AS1n6LBERqaU6m3vRzMYCe27kpRvdfVT1\nMTcBq939yRzUJyIiKTB3T//NZucBFwEnufvKTRyT/glEREqYu1uqx262Zb45ZtYRuA5ov6kgr20x\nIiKSnrRb5mY2F9gWqKh+apK7/z5bhYmISOoy6mYREZHCkLMZoGbW0cxmm9lcM7s+V+cpdGbWyMzG\nmdlMM/vQzK6KuqaomdnWZva+mY2KupYomdnOZvZs9eS7WWbWJuqaomJmN1T/HZlhZk+a2XZR15Qv\nZvaomZWb2Ywaz+1qZmPN7GMze8XMdt7S5+QkzM1sa+B+oCPQAuhuZs1zca4i8APQ291bAm2Ay0v4\nu1jnamAWUOo/C+8FXnL35sAhwEcR1xMJM2tMGEjR2t1bAVsDv4mypjwbSMjKmvoAY939AOC16seb\nlauW+VHAPHf/1N1/AJ4CzsjRuQqauy929w+q7y8j/IXdO9qqomNm+wKdgEeAkr04Xr38xXHu/iiA\nu1e5e2XEZUXlO0KjZwczqwPsAHwRbUn54+5vAt9s8PTpwODq+4OBM7f0ObkK832ABTUeL6x+rqRV\nt0AOJ8yYLVX/SxgFtXZLB8ZcE+ArMxtoZlPN7GEz2yHqoqLg7hXAXcDnwCLgW3d/NdqqItfQ3cur\n75cDDbf0hlyFean/fP4PZlYfeBa4urqFXnLM7DRgibu/Twm3yqvVAVoDD7h7a+B7UvgpHUdm1gy4\nBmhM+NVa38x+G2lRBcTDKJUtZmquwvwLoFGNx40IrfOSZGbbAMOBJ9x9RNT1ROgY4HQz+wQYCpxo\nZo9FXFNUFgIL3X1K9eNnCeFeio4A3nL3r929CniO8GellJWb2Z4AZrYXsGRLb8hVmL8L7G9mjc1s\nW6AbMDJH5ypoZmbAAGCWu98TdT1Rcvcb3b2RuzchXOB63d3PjbquKLj7YmCBmR1Q/VQHYGaEJUVp\nNtDGzOpW/33pQLhAXspGAj2r7/cEttgITHsG6Oa4e5WZXQGMIVyZHuDuJXmlHmgH/A6YbmbvVz93\ng7u/HGFNhaLUu+OuBIZUN3jmA+dHXE8k3H1a9S+0dwnXUqYCD0VbVf6Y2VCgPbC7mS0A+gK3A8PM\n7ALgU6DrFj9Hk4ZERIqfto0TEYkBhbmISAwozEVEYkBhLiISAwpzEZEYUJiLiMSAwlxEJAYU5iIi\nMfB/7DXjYjSZufkAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106af0290>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Some interesting properties of UFuncs\n",
      "\n",
      "UFuncs have some methods built-in, which allow for some very interesting, flexible, and fast operations:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(5)\n",
      "y = np.arange(1, 6)\n",
      "np.add(x, y)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "array([1, 3, 5, 7, 9])"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.accumulate(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "array([ 0,  1,  3,  6, 10])"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.accumulate(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "array([0, 0, 0, 0, 0])"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.accumulate(y)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "array([  1,   2,   6,  24, 120])"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.identity"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.identity"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "1"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.outer(x, y)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 28,
       "text": [
        "array([[1, 2, 3, 4, 5],\n",
        "       [2, 3, 4, 5, 6],\n",
        "       [3, 4, 5, 6, 7],\n",
        "       [4, 5, 6, 7, 8],\n",
        "       [5, 6, 7, 8, 9]])"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# make a times-table\n",
      "x = np.arange(1, 13)\n",
      "np.multiply.outer(x, x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 30,
       "text": [
        "array([[  1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12],\n",
        "       [  2,   4,   6,   8,  10,  12,  14,  16,  18,  20,  22,  24],\n",
        "       [  3,   6,   9,  12,  15,  18,  21,  24,  27,  30,  33,  36],\n",
        "       [  4,   8,  12,  16,  20,  24,  28,  32,  36,  40,  44,  48],\n",
        "       [  5,  10,  15,  20,  25,  30,  35,  40,  45,  50,  55,  60],\n",
        "       [  6,  12,  18,  24,  30,  36,  42,  48,  54,  60,  66,  72],\n",
        "       [  7,  14,  21,  28,  35,  42,  49,  56,  63,  70,  77,  84],\n",
        "       [  8,  16,  24,  32,  40,  48,  56,  64,  72,  80,  88,  96],\n",
        "       [  9,  18,  27,  36,  45,  54,  63,  72,  81,  90,  99, 108],\n",
        "       [ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100, 110, 120],\n",
        "       [ 11,  22,  33,  44,  55,  66,  77,  88,  99, 110, 121, 132],\n",
        "       [ 12,  24,  36,  48,  60,  72,  84,  96, 108, 120, 132, 144]])"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Ufunc mini-exercises\n",
      "\n",
      "Each of the following functions take an array as input, and return an array as output.  They are implemented using loops, which is not very efficient.\n",
      "\n",
      "1. For each function, implement a fast version which uses ufuncs to calculate the result more efficiently.  Double-check that you get the same result for several different arrays.\n",
      "\n",
      "2. use the %timeit magic to time the execution of the two implementations for a large array (say, 1000 elements).  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 1. computing the element-wise sine + cosine\n",
      "\n",
      "from math import sin, cos\n",
      "def slow_sincos(x):\n",
      "    \"\"\"x is a 1-dimensional array\"\"\"\n",
      "    y = np.zeros_like(x)\n",
      "    for i in range(len(x)):\n",
      "        y[i] = sin(x[i]) + cos(x[i])\n",
      "    return y\n",
      "\n",
      "x = np.random.random(5)\n",
      "print(slow_sincos(x))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1.38638594  1.41068192  1.41074488  1.39904614  1.39058355]\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# write a fast_sincos function\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 2. computing the difference between adjacent squares\n",
      "\n",
      "def slow_sqdiff(x):\n",
      "    \"\"\"x is a 1-dimensional array\"\"\"\n",
      "    y = np.zeros(len(x) - 1)\n",
      "    for i in range(len(y)):\n",
      "        y[i] = x[i + 1] ** 2 - x[i] ** 2\n",
      "    return y\n",
      "\n",
      "x = np.random.random(5)\n",
      "print(slow_sqdiff(x))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.40936747 -0.30871454  0.49088252 -0.1586789 ]\n"
       ]
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# write a fast_sqdiff function\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 3. computing the outer-product of each consecutive pair\n",
      "\n",
      "def slow_pairprod(x):\n",
      "    \"\"\"x is a 1-dimensional array\"\"\"\n",
      "    if len(x) % 2 != 0:\n",
      "        raise ValueError(\"length of x must be even\")\n",
      "    N = len(x) // 2\n",
      "    y = np.zeros((N, N))\n",
      "    for i in range(N):\n",
      "        for j in range(N):\n",
      "            y[i, j] = x[2 * i] * x[2 * j + 1]\n",
      "    return y\n",
      "\n",
      "x = np.arange(1, 9)\n",
      "print(slow_pairprod(x))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  2.   4.   6.   8.]\n",
        " [  6.  12.  18.  24.]\n",
        " [ 10.  20.  30.  40.]\n",
        " [ 14.  28.  42.  56.]]\n"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# write a fast_pairprod function\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Strategy 2. Using Numpy Aggregates\n",
      "\n",
      "Aggregates are functions over arrays which return smaller arrays.\n",
      "\n",
      "Numpy has several built-in"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 10 x 10 array drawn from a standard normal\n",
      "x = np.random.randn(10, 10)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.mean()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 40,
       "text": [
        "0.0086655841045156605"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.mean(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 41,
       "text": [
        "0.0086655841045156605"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.std()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 42,
       "text": [
        "0.99647690340540895"
       ]
      }
     ],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.var()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 43,
       "text": [
        "0.99296621902043269"
       ]
      }
     ],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.sum()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 44,
       "text": [
        "0.86655841045156601"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.prod()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 45,
       "text": [
        "4.788617672159386e-22"
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.median(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "0.12432664008018213"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.percentile(x, 50)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 47,
       "text": [
        "0.12432664008018213"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.percentile(x, (25, 75))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "array([-0.59558374,  0.6702099 ])"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Aggregates along certain dimensions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.rand(3, 5)\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "array([[ 0.90877938,  0.58853087,  0.3233339 ,  0.63680331,  0.64820305],\n",
        "       [ 0.37713506,  0.68006887,  0.68975908,  0.53422032,  0.3331116 ],\n",
        "       [ 0.93016689,  0.55381498,  0.35948829,  0.80807909,  0.34102259]])"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.sum(0)  # sum all rows"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 50,
       "text": [
        "array([ 2.21608134,  1.82241471,  1.37258128,  1.97910272,  1.32233723])"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.sum(1)  # sum all columns"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 51,
       "text": [
        "array([ 3.10565051,  2.61429493,  2.99257184])"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.median(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 52,
       "text": [
        "array([ 0.63680331,  0.53422032,  0.55381498])"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.mean(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 53,
       "text": [
        "array([ 0.6211301 ,  0.52285899,  0.59851437])"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Binary ufuncs as aggregates\n",
      "\n",
      "Any binary ufunc (a ufunc taking two arguments) can be turned into an aggregate using the ``reduce()`` method:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sum(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 54,
       "text": [
        "array([ 3.10565051,  2.61429493,  2.99257184])"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.reduce(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 55,
       "text": [
        "array([ 3.10565051,  2.61429493,  2.99257184])"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.prod(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 56,
       "text": [
        "array([ 0.07138309,  0.03148164,  0.05103256])"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.reduce(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 57,
       "text": [
        "array([ 0.07138309,  0.03148164,  0.05103256])"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.divide.reduce(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "array([ 11.56968674,   4.51789808,  16.95408724])"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "A caution: for ``reduce`` methods, the default axis is 0:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.reduce(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 59,
       "text": [
        "array([ 2.21608134,  1.82241471,  1.37258128,  1.97910272,  1.32233723])"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sum(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 60,
       "text": [
        "8.7125172829842636"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### A quick efficiency note:\n",
      "### Beware the built-in Python aggregates!\n",
      "\n",
      "Python has a ``min``, ``max``, and ``sum`` aggregate built-in.  These are much more general than the versions in NumPy:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random(10000)\n",
      "\n",
      "%timeit np.sum(x)\n",
      "%timeit sum(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 8.9 \u00b5s per loop\n",
        "1000 loops, best of 3: 920 \u00b5s per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "**Dynamic type-checking is slow**.\n",
      "\n",
      "Make sure to use Numpy's ``sum``, ``min``, and ``max``."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Aggregate Mini-exercises\n",
      "\n",
      "Take the following functions, and convert them into an efficient form using aggregates.  Each function expects a 1-dimensional array as input.  Double-check that your function returns the same result as the original"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def slow_cubesum(x):\n",
      "    \"\"\"x is a 1D array\"\"\"\n",
      "    result = 0\n",
      "    for i in range(len(x)):\n",
      "        result += x[i] ** 3\n",
      "    return result\n",
      "\n",
      "x = np.random.random(100)\n",
      "slow_cubesum(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 62,
       "text": [
        "22.195828759459477"
       ]
      }
     ],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# implement fast_cubesum\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def slow_rms(x):\n",
      "    \"\"\"x is a 1D array\"\"\"\n",
      "    m = np.mean(x)\n",
      "    rms = 0\n",
      "    for i in range(len(x)):\n",
      "        rms += (x[i] - m) ** 2\n",
      "    rms /= len(x)\n",
      "    return np.sqrt(rms)\n",
      "\n",
      "x = np.random.random(100)\n",
      "slow_rms(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 64,
       "text": [
        "0.30973927909541354"
       ]
      }
     ],
     "prompt_number": 64
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# implement fast_rms\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Now we return to our silly function from the beginning of this section.  Can you implement a fast version using ufuncs and aggregates?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def slow_sillyfunc(N):\n",
      "    \"\"\"N is an integer\"\"\"\n",
      "    d = 0.0\n",
      "    for i in range(N):\n",
      "        d += (i % 3 - 1) * i\n",
      "    return d\n",
      "\n",
      "slow_sillyfunc(100)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 66,
       "text": [
        "-33.0"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Implement fast_sillyfunc using ufuncs & aggragates\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 67
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Strategy 3: Using Numpy Broadcasting\n",
      "\n",
      "We've taken a look at broadcasting previously.  But it's important enough that we'll review it quickly here:\n",
      "\n",
      "\n",
      "<img src=\"http://www.astroml.org/_images/fig_broadcast_visual_1.png\">\n",
      "\n",
      "([image source](http://www.astroml.org/book_figures/appendix/fig_broadcast_visual.html))\n",
      "\n",
      "**Broadcasting rules:**\n",
      "\n",
      "1. If the two arrays differ in their number of dimensions, the shape of the array with fewer dimensions is *padded* with ones on its leading (left) side.\n",
      "\n",
      "2. If the shape of the two arrays does not match in any dimension, the array with shape equal to 1 in that dimension is *stretched* to match the other shape.\n",
      "\n",
      "3. If in any dimension the sizes disagree and neither is equal to 1, an error is raised."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Some Broadcasting examples..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(10)\n",
      "x ** 2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 68,
       "text": [
        "array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = x * x[:, np.newaxis]\n",
      "Y"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 69,
       "text": [
        "array([[ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0],\n",
        "       [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],\n",
        "       [ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18],\n",
        "       [ 0,  3,  6,  9, 12, 15, 18, 21, 24, 27],\n",
        "       [ 0,  4,  8, 12, 16, 20, 24, 28, 32, 36],\n",
        "       [ 0,  5, 10, 15, 20, 25, 30, 35, 40, 45],\n",
        "       [ 0,  6, 12, 18, 24, 30, 36, 42, 48, 54],\n",
        "       [ 0,  7, 14, 21, 28, 35, 42, 49, 56, 63],\n",
        "       [ 0,  8, 16, 24, 32, 40, 48, 56, 64, 72],\n",
        "       [ 0,  9, 18, 27, 36, 45, 54, 63, 72, 81]])"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y + 10 * x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 71,
       "text": [
        "array([[  0,  10,  20,  30,  40,  50,  60,  70,  80,  90],\n",
        "       [  0,  11,  22,  33,  44,  55,  66,  77,  88,  99],\n",
        "       [  0,  12,  24,  36,  48,  60,  72,  84,  96, 108],\n",
        "       [  0,  13,  26,  39,  52,  65,  78,  91, 104, 117],\n",
        "       [  0,  14,  28,  42,  56,  70,  84,  98, 112, 126],\n",
        "       [  0,  15,  30,  45,  60,  75,  90, 105, 120, 135],\n",
        "       [  0,  16,  32,  48,  64,  80,  96, 112, 128, 144],\n",
        "       [  0,  17,  34,  51,  68,  85, 102, 119, 136, 153],\n",
        "       [  0,  18,  36,  54,  72,  90, 108, 126, 144, 162],\n",
        "       [  0,  19,  38,  57,  76,  95, 114, 133, 152, 171]])"
       ]
      }
     ],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y + 10 * x[:, np.newaxis]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 72,
       "text": [
        "array([[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0],\n",
        "       [ 10,  11,  12,  13,  14,  15,  16,  17,  18,  19],\n",
        "       [ 20,  22,  24,  26,  28,  30,  32,  34,  36,  38],\n",
        "       [ 30,  33,  36,  39,  42,  45,  48,  51,  54,  57],\n",
        "       [ 40,  44,  48,  52,  56,  60,  64,  68,  72,  76],\n",
        "       [ 50,  55,  60,  65,  70,  75,  80,  85,  90,  95],\n",
        "       [ 60,  66,  72,  78,  84,  90,  96, 102, 108, 114],\n",
        "       [ 70,  77,  84,  91,  98, 105, 112, 119, 126, 133],\n",
        "       [ 80,  88,  96, 104, 112, 120, 128, 136, 144, 152],\n",
        "       [ 90,  99, 108, 117, 126, 135, 144, 153, 162, 171]])"
       ]
      }
     ],
     "prompt_number": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = np.random.random((2, 3, 4))\n",
      "x = 10 * np.arange(3)\n",
      "\n",
      "Y + x[:, np.newaxis]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 73,
       "text": [
        "array([[[  0.16013996,   0.61070178,   0.70907938,   0.52638606],\n",
        "        [ 10.74057435,  10.40082693,  10.48134161,  10.70561624],\n",
        "        [ 20.45622939,  20.68181045,  20.40559978,  20.0304505 ]],\n",
        "\n",
        "       [[  0.99607108,   0.93981662,   0.25316248,   0.41212159],\n",
        "        [ 10.06903578,  10.41168541,  10.48591661,  10.24673383],\n",
        "        [ 20.46827554,  20.32719887,  20.46975577,  20.39303078]]])"
       ]
      }
     ],
     "prompt_number": 73
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Quick Broadcasting Exercise\n",
      "\n",
      "Now, assume you have $N$ points in $D$ dimensions, represented by an array of shape ``[N, D]``.\n",
      "\n",
      "1. Compute the mean of the distribution of points efficiently using the built-in ``np.mean`` aggregate (that is, find the ``D``-dimensional point which is the mean of the rest of the points)\n",
      "2. Compute the mean of the distribution of points efficiently using the ``np.add`` ufunc.\n",
      "3. Compute the standard error of the mean $\\sigma_{mean} = \\sigma N^{-1/2}$, where $\\sigma$ is the standard-deviation, using the ``np.std`` aggregate.\n",
      "4. Compute this again using the ``np.add`` ufunc.\n",
      "5. Construct the matrix ``M``, the centered and normalized version of the ``X`` array: $$ M_{ij} = (X_{ij} - \\mu_j) / \\sigma_j $$  This is one version of *whitening* the array.\n",
      "    "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.random.random((1000, 5))  # 1000 points in 5 dimensions"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 74
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 1. Compute the mean of the 1000 points in X\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 2. Compute the standard deviation across the 1000 points\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 73
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 5. Compute the whitened version of the array\n",
      "#    (Each row centered, then divided by standard deviation)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 75
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Strategy 4: Fancy Indexing and Masking\n",
      "\n",
      "The last strategy we will cover is fancy indexing and masking.\n",
      "\n",
      "For example, imagine you have an array of data where negative values indicate some kind of error."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([1, 2, 3, -999, 2, 4, -999])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 75
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "How might you clean this array, setting all negative values to, say, zero?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(len(x)):\n",
      "    if x[i] < 0:\n",
      "        x[i] = 0\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 76,
       "text": [
        "array([1, 2, 3, 0, 2, 4, 0])"
       ]
      }
     ],
     "prompt_number": 76
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "A faster way is to construct a *boolean mask*:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([1, 2, 3, -999, 2, 4, -999])\n",
      "\n",
      "mask = (x < 0)\n",
      "mask"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 77,
       "text": [
        "array([False, False, False,  True, False, False,  True], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 77
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "And the mask can be used directly to set the value you desire:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[mask] = 0\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 78,
       "text": [
        "array([1, 2, 3, 0, 2, 4, 0])"
       ]
      }
     ],
     "prompt_number": 78
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Typically this is done directly:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([1, 2, 3, -999, 2, 4, -999])\n",
      "x[x < 0] = 0\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 79,
       "text": [
        "array([1, 2, 3, 0, 2, 4, 0])"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Useful masking functions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random(5)\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 80,
       "text": [
        "array([ 0.67649933,  0.96325803,  0.94898948,  0.44587652,  0.63444064])"
       ]
      }
     ],
     "prompt_number": 80
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[x > 0.5] = np.nan\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 81,
       "text": [
        "array([        nan,         nan,         nan,  0.44587652,         nan])"
       ]
      }
     ],
     "prompt_number": 81
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[np.isnan(x)] = np.inf\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 82,
       "text": [
        "array([        inf,         inf,         inf,  0.44587652,         inf])"
       ]
      }
     ],
     "prompt_number": 82
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.nan == np.nan"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 84,
       "text": [
        "False"
       ]
      }
     ],
     "prompt_number": 84
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[np.isinf(x)] = 0\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 83,
       "text": [
        "array([ 0.        ,  0.        ,  0.        ,  0.44587652,  0.        ])"
       ]
      }
     ],
     "prompt_number": 83
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([1, 0, -np.inf, np.inf, np.nan])\n",
      "print(\"input   \", x)\n",
      "print(\"x < 0   \", (x < 0))\n",
      "print(\"x > 0   \", (x > 0))\n",
      "print(\"isinf   \", np.isinf(x))\n",
      "print(\"isnan   \", np.isnan(x))\n",
      "print(\"isposinf\", np.isposinf(x))\n",
      "print(\"isneginf\", np.isneginf(x))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "input    [  1.   0. -inf  inf  nan]\n",
        "x < 0    [False False  True False False]\n",
        "x > 0    [ True False False  True False]\n",
        "isinf    [False False  True  True False]\n",
        "isnan    [False False False False  True]\n",
        "isposinf [False False False  True False]\n",
        "isneginf [False False  True False False]\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:3: RuntimeWarning: invalid value encountered in less\n",
        "-c:4: RuntimeWarning: invalid value encountered in greater\n"
       ]
      }
     ],
     "prompt_number": 84
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Boolean Operations on Masks\n",
      "\n",
      "Use bitwise operators (and make sure to use parentheses!)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(16).reshape((4, 4))\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 85,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [ 4,  5,  6,  7],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15]])"
       ]
      }
     ],
     "prompt_number": 85
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x < 5"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 86,
       "text": [
        "array([[ True,  True,  True,  True],\n",
        "       [ True, False, False, False],\n",
        "       [False, False, False, False],\n",
        "       [False, False, False, False]], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 86
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "~(x < 5)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 87,
       "text": [
        "array([[False, False, False, False],\n",
        "       [False,  True,  True,  True],\n",
        "       [ True,  True,  True,  True],\n",
        "       [ True,  True,  True,  True]], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 87
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(x < 10) & (x % 2 == 0)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 88,
       "text": [
        "array([[ True, False,  True, False],\n",
        "       [ True, False,  True, False],\n",
        "       [ True, False, False, False],\n",
        "       [False, False, False, False]], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 88
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(x > 3) & (x < 8)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 89,
       "text": [
        "array([[False, False, False, False],\n",
        "       [ True,  True,  True,  True],\n",
        "       [False, False, False, False],\n",
        "       [False, False, False, False]], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 89
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Counting elements with a mask\n",
      "\n",
      "Sum over a mask to find the number of ``True`` elements:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random(100)\n",
      "print(\"array is length\", len(x), \"and has\")\n",
      "print((x > 0.5).sum(), \"elements are greater than 0.5\")"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "array is length 100 and has\n",
        "61 elements are greater than 0.5\n"
       ]
      }
     ],
     "prompt_number": 90
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# clip is a useful function:\n",
      "x = np.clip(x, 0.3, 0.6)\n",
      "\n",
      "print(np.sum(x < 0.3))\n",
      "print(np.sum(x > 0.6))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0\n",
        "0\n"
       ]
      }
     ],
     "prompt_number": 91
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# works for 2D arrays as well\n",
      "X = np.random.random((10, 10))\n",
      "(X < 0.1).sum()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 92,
       "text": [
        "13"
       ]
      }
     ],
     "prompt_number": 92
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### ``where`` function: Turning a mask into indices"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random((3, 3))\n",
      "x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 93,
       "text": [
        "array([[ 0.05999021,  0.03303966,  0.89629243],\n",
        "       [ 0.87950144,  0.30331336,  0.17156985],\n",
        "       [ 0.56147418,  0.2192196 ,  0.90690265]])"
       ]
      }
     ],
     "prompt_number": 93
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.where(x < 0.3)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 94,
       "text": [
        "(array([0, 0, 1, 2]), array([0, 1, 2, 1]))"
       ]
      }
     ],
     "prompt_number": 94
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[x < 0.3]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 95,
       "text": [
        "array([ 0.05999021,  0.03303966,  0.17156985,  0.2192196 ])"
       ]
      }
     ],
     "prompt_number": 95
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[np.where(x < 0.3)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 96,
       "text": [
        "array([ 0.05999021,  0.03303966,  0.17156985,  0.2192196 ])"
       ]
      }
     ],
     "prompt_number": 96
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "When you index with the result of a ``where`` function, you are using what is called *fancy indexing*: indexing with tuples"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Fancy Indexing (indexing with sequences)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.arange(16).reshape((4, 4))\n",
      "X"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 97,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [ 4,  5,  6,  7],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15]])"
       ]
      }
     ],
     "prompt_number": 97
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[(0, 1), (1, 0)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 98,
       "text": [
        "array([1, 4])"
       ]
      }
     ],
     "prompt_number": 98
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[range(4), range(4)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 99,
       "text": [
        "array([ 0,  5, 10, 15])"
       ]
      }
     ],
     "prompt_number": 99
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X.diagonal()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 100,
       "text": [
        "array([ 0,  5, 10, 15])"
       ]
      }
     ],
     "prompt_number": 100
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X.diagonal() = 100"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "ename": "SyntaxError",
       "evalue": "can't assign to function call (<ipython-input-101-3064a6b8dbd8>, line 1)",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-101-3064a6b8dbd8>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    X.diagonal() = 100\u001b[0m\n\u001b[0m                      ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m can't assign to function call\n"
       ]
      }
     ],
     "prompt_number": 101
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[range(4), range(4)] = 100"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 102
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 103,
       "text": [
        "array([[100,   1,   2,   3],\n",
        "       [  4, 100,   6,   7],\n",
        "       [  8,   9, 100,  11],\n",
        "       [ 12,  13,  14, 100]])"
       ]
      }
     ],
     "prompt_number": 103
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "#### Randomizing the rows"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.arange(24).reshape((6, 4))\n",
      "X"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 104,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [ 4,  5,  6,  7],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15],\n",
        "       [16, 17, 18, 19],\n",
        "       [20, 21, 22, 23]])"
       ]
      }
     ],
     "prompt_number": 104
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "i = np.arange(6)\n",
      "np.random.shuffle(i)\n",
      "i"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 105,
       "text": [
        "array([0, 1, 5, 4, 2, 3])"
       ]
      }
     ],
     "prompt_number": 105
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[i]  # X[i, :] is identical"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 106,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [ 4,  5,  6,  7],\n",
        "       [20, 21, 22, 23],\n",
        "       [16, 17, 18, 19],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15]])"
       ]
      }
     ],
     "prompt_number": 106
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Fancy indexing also works for multi-dimensional index arrays"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "i2 = i.reshape(3, 2)\n",
      "X[i2]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 107,
       "text": [
        "array([[[ 0,  1,  2,  3],\n",
        "        [ 4,  5,  6,  7]],\n",
        "\n",
        "       [[20, 21, 22, 23],\n",
        "        [16, 17, 18, 19]],\n",
        "\n",
        "       [[ 8,  9, 10, 11],\n",
        "        [12, 13, 14, 15]]])"
       ]
      }
     ],
     "prompt_number": 107
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[i2].shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 108,
       "text": [
        "(3, 2, 4)"
       ]
      }
     ],
     "prompt_number": 108
    },
    {
     "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": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Breakout: asteroid data\n",
      "\n",
      "In the [github repository](https://github.com/jakevdp/2014_fall_ASTR599/), there is a file containing measurements of 5000 asteroid orbits, at ``notebooks/data/asteroids5000.csv``.\n",
      "\n",
      "These are compiled from a query at [http://ssd.jpl.nasa.gov/sbdb_query.cgi](http://ssd.jpl.nasa.gov/sbdb_query.cgi)\n",
      "\n",
      "### Part 1: loading and exploring the data\n",
      "\n",
      "1. Use ``np.genfromtxt`` to load the data from the file.  This is like ``loadtxt``, but can handle missing data.\n",
      "\n",
      "    - Remember to set the appropriate ``delimiter`` keyword.\n",
      "\n",
      "2. ``genfromtxt`` sets all missing values to ``np.nan``.  Use the operations we discussed here to answer these questions:\n",
      "\n",
      "     - How many values are missing in this data?\n",
      "     - How many complete *rows* are there?  i.e. how many objects have **no** missing values?\n",
      "\n",
      "3. Create a new array containing only the rows with no missing values.\n",
      "\n",
      "4. Compute the maximum, minimum, mean, and standard deviation of the values in each column."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Part 2: Plotting the data\n",
      "\n",
      "1. Use the bash ``head`` command to display the first line of the data file: this lists the names of the columns in the dataset. (remember that bash commands in the notebook are indicated by ``!``, and that ``head -n`` displays the first ``n`` lines of a file)\n",
      "\n",
      "2. Invoke the ``matplotlib inline`` magic to make figures appear inline in the notebook\n",
      "\n",
      "3. Use ``plt.scatter`` to plot the semi-major axis versus the sine of the inclination angle (note that the inclination angle is listed in degrees -- you'll have to convert it to radians to compute the sine).  What do you notice about the distribution?  What do you think could explain this?\n",
      "\n",
      "4. Use ``plt.scatter`` to plot a color-magnitude diagram of the asteroids (H vs B-V).  You should see two distinct \"families\" of asteroids indicated in this plot.  Over-plot a line that divides these.\n",
      "\n",
      "5. Repeat the orbital parameter plot from above, but plot the two \"families\" from #4 in different colors.  Note that this magnitude is undefined for many of the asteroids.  Do you see any correlation between color and orbit?\n",
      "\n",
      "6. Compare what you found to plots in [Parker et al. 2008](http://adsabs.harvard.edu/abs/2008Icar..198..138P).  Note that we're not using the same data here, but we're looking at a similar collection of objects."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}