day1/EfficientNumpy.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:2e3455a4efe8212a1ffd9adc03e22d2566cec40f55fe29733324593a5e7c7c68"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "\n",
      "from __future__ import print_function"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Efficient Numerical Computing with Numpy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*Jake VanderPlas*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "In this session, we'll discuss how to make your programs as efficient as possible, mainly by taking advantage of *vectorization* in NumPy.\n",
      "\n",
      "If you use Python for any amount of time, you'll quickly find that there are some things it is not so good at.\n",
      "In particular, performing repeated operations via loops is one of its weaknesses.\n",
      "\n",
      "For 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": 2
    },
    {
     "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": [
        "100 loops, best of 3: 1.97 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "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": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!f2py -c func_fortran.f -m func_fortran > /dev/null"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 5
    },
    {
     "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: 16.1 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Fortran is about 100 times faster for this task!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Why is Python so slow?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Languages tend to compromise between convenience and performance.\n",
      "\n",
      "- **High level languages** (Python, R, IDL, Matlab, etc.) tend to use **High-level objects**, **dynamic typing**, and **interpreted execution** for convenience in writing code\n",
      "\n",
      "  + But: have generally large execution overhead for type-checking, etc.\n",
      "  + JIT (just in time) compilation technologies, like LLVM (e.g. Julia, Numba) may be changing this, though\n",
      "  \n",
      "- **Low-level languages** (C, Fortran, etc.) tend to use **static types** and **compiled execution** for performance\n",
      "\n",
      "  + But: no interactive prompt, lots of development time spent on details like typing, etc.\n",
      "    \n",
      "We like high-level languages because our **development time** is generally more valuable than **execution time**.  But sometimes speed can be an issue."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Strategies for making Python fast"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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, masking, and fancy indexing** 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 strategies."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Strategy 1: Using UFuncs in Numpy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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": [
      "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.1231047   0.34071966  0.97596715  0.38643665]\n",
        "[ 1.1231047   1.34071966  1.97596715  1.38643665]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(x * 2)  # multiply each element of x by 2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.24620939  0.68143932  1.95193429  0.77287329]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(x * x)  # multiply each element of x by itself"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.01515477  0.11608989  0.95251187  0.14933328]\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(x[1:] - x[:-1])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.21761497  0.63524749 -0.5895305 ]\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "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": 11,
       "text": [
        "array([-0.1231047 , -0.34071966, -0.97596715, -0.38643665])"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sin(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([ 0.122794  ,  0.33416547,  0.82824423,  0.37689023])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "The Speed of Ufuncs"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.random.random(100000)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit\n",
      "# compute element-wise x + 1 via a ufunc \n",
      "y = x + 1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 57.1 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit\n",
      "# compute element-wise x + 1 via a loop\n",
      "y = np.empty_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": [
        "10 loops, best of 3: 118 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Why is NumPy so much faster?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "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": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Other Available Ufuncs"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "- 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": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "``np.min`` vs ``np.minimum``"
     ]
    },
    {
     "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.64960338  0.44854455  0.93193523  0.60049152  0.0451027 ]\n",
        "[ 0.5         0.44854455  0.5         0.5         0.0451027 ]\n",
        "[ 0.64960338  0.5         0.93193523  0.60049152  0.5       ]\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "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.0451027007923\n",
        "0.931935234484\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Sines, Gamma functions, etc."
     ]
    },
    {
     "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+/AAAIABJREFUeJzt3Xl4ldW1x/FvSBhVQESZFWWwDggoKmCRCMJFJEAgijiP\nxSpXqx0sbe8j9Vq1vVp8lFqnWqkoggzKICCDUQqIAoKMAmEeRAZRRKaQc/9YSYkQIDnDu9/h93me\nPDkJ55x3kZyzst/17r02iIiIiIiIiIiIiIiIiIiIiIiIiIhIALwGbAUWHec+zwErgYVASy+CEhGR\n1GiHJfJjJf2uwPuFty8HPvEiKBERSZ2GHDvpvwj0Kfb1cqBWqgMSEZGjlfPgGPWADcW+3gjU9+C4\nIiJyBC+SPkDaEV/HPDquiIgUk+HBMTYBDYp9Xb/wez/SqFGjWF5engfhiIiESh7QuLR39mKkPxa4\ntfB2a2AXNtvnR/Ly8ojFYoH52LIlRu/eMZo2jTF8eIz8/OPff+vWGL/9bYzTTovx3HMxCgqOfd9H\nH33U+f/PLx/6WYT7Z3HoUIy//tXeF3/4Q4xt245//4MHY7z1Vozq1R+lT58YX3/t/v/g+gNoVJaE\nnIyR/jCgPVATq90/CpQv/LeXsJk7XYFVwB7gjiQc06nZsyEnB265BYYOhUqVTvyYM86AJ5+E22+H\nW2+F3Fz45z+hatVURyviT7t22Xto5057TzVpcuLHZGRA376weDHs3w+tWsHIkXDppamPNyySkfT7\nluI+/ZNwHF8YPRruvRf+8Q/Iyir74889Fz7+GB54ANq1gw8+gFqayyQRs2ULXH01XHUVjBoFFSqU\n7fHly8Of/gRXXAHXXhv/+zGKvLqQGwojRsD991uiTuQFVrEivPgi9OpliX/TEVc4MjMzE4ozTPSz\nOCwsP4sNG6B9e7jpJhg8uOwJHw7/LLKzYcIEuOcee3/KiR05q8alWGF9ypfGjoV+/SzhN2uWvOd9\n8kkYNsxG/9WrJ+95Rfxo+3Zo2xZ+9jP41a+S97xffAGdOlm5tVOn5D1vEKSlpUEZcrmSfinMnw9d\nusD771sNMZliMXjwQXvRTppUuusDIkG0b5+VdNq1s8FOss2YAb1728CsRYvkP79fKekn2ebNcPnl\n8NxzdiqZCocOwXXX2cXeF19MzTFEXIrFbALD/v3w9ttQLkWF5XfesTOIefOgZs3UHMNvypr0VdM/\njkOHrO54zz2pS/gA6enw+uswfTq88UbqjiPiyquvwsKFMGRI6hI+2OCpTx+bFVRQkLrjBJlG+sfx\n2GPw0Ud2upienvrjLVoEHTrYMc8/P/XHE/FC0et6xgz4yU9Sf7yDB6FjR/iv/4Lf/z71x3NN5Z0k\nmTXL6oPz5kHdut4d98UXbfrZ7Nk2J1kkyPbvh0sugV//Gm67zbvjbtoELVvC5Mn2OcxU3kmCffvg\nrrvg+ee9TfhgM4ROPRX+/GdvjyuSCk88AY0bWz3fS/XqwTPP2GLIAwe8PbbfaaRfgv/5H1vxN3o0\npDn4CW3YABdfDB9+CBde6P3xRZKhqKyzYIElYa/FYtCjh83keewx74/vFZV3EvTFF1YPXLjQ+1F+\ncS+8AMOHW7sGF394RBJx6BC0aWOTIO65x10cmzfDRRdZubZpU3dxpJLKOwmIxeC+++Dxx90mfLAy\nz+7dtnBLJGhefdXWnNx9t9s46taFAQNsLYwPxpS+4KcxpPOR/ogRtmhk7lxvZuucSFFjt2XL1JhN\ngmPXLpulM3GiPy6iHjwIzZvDU09B9+6uo0k+lXfitHcvnHeezZf3U4uTO+6A00+Hv/zFdSQipfOr\nX1nif/VV15EcNnWqtX5YujR8q96V9OP0xBM2wh892lkIJdq82Xr9LFgADRqc+P4iLq1cabX8JUv8\n1z22Rw8b0D30kOtIkktJPw7bttnp6Jw5Nr3Mb373O9i61ebvi/jZ9ddbSWfAANeRHG3xYpuksXJl\nuMqlSvpx+PWvYc8emzHjR7t22cyD3Fyt1BX/WrjQVsHm5cFJJ7mOpmS33w5nnhmuKZxK+mX01VeW\nSBctcjOXuLSefhpmzoQxY1xHIlKynj2tT76fyydr19oK4aVL/Vd+ipeSfhk9+KDNg3/2Wc8PXSZ7\n91rpafx4f8yIEClu3jybGbNqFVSu7Dqa43vgAdt565lnXEeSHEr6ZbBxo03lWrIEatf29NBxGTTI\nRvsjR7qOROTHrr0Wuna1neX8buNGW7C1YkU42i8r6ZfBQw9Zm9eg/MXfswfOOcdaMF9wgetoRMxn\nn1lzwpUrbSvQIPjZz6y887//6zqSxCnpl9KOHdCkif9r+Ud68kmbhfDmm64jETE5ObYb1oMPuo6k\n9PLybHOkvDyoVs11NIlR0i+lxx6DdeuCNw3yu+9stD97tv3REnFp5Urb83btWv/O2DmWm2+2M2Y/\nTi8tCyX9Uigqk3z0kTebOiTbH/4A33wDf/ub60gk6u6917b5DOIUyCVLrAvomjVQpYrraOKnpF8K\nzz9vbYv9tvq2tLZssWmmeXlQo4braCSqiqY7f/mltQoJou7doVs3q/EHlbpsnkB+vl24feQR15HE\nr04dW1L+0kuuI5Eoe/556Ns3uAkf4Be/sOnaUerAGbmkP3asXbi9/HLXkSTmoYdg8GDtCiRu/PAD\nvPwyPPyw60gSc9VVNmd/yhTXkXgnckl/8GDo3991FIlr3tyuR4wY4ToSiaJhw6B1a2jUyHUkiUlL\ns9H+oEGuI/FOpJL+4sWwfLnNKQ6Dhx+Gv/41Wqem4l4sFp7BE1iJav5827ciCiKV9AcPth2pKlRw\nHUlyXHON7a71ySeuI5EomTXLZsB16uQ6kuSoVMlmIT33nOtIvBGZ2TvffGPTNJcutQuhYfH007bA\nbMgQ15FIVNx4o10TC9JirBMpmhG3bl3w2i5ryuYxDBpky8Xfeitlh3Bi+3ZrxJaXB6ed5joaCbst\nW2xB0+rVUL2662iSKycHrr7aRv1BoimbJSgosIVMYalBFlezJmRlaaQv3njlFejTJ3wJH6z0++KL\n4b9GFomk/+GHtkS8TRvXkaTGvfdG48UqbuXn29qQ++5zHUlqdOxo18g+/dR1JKkViaT/6qtw9902\nPSuM2ra1i1HTp7uORMJs0iTbdapZM9eRpEa5crYyN+yLHv2UBlNS09+xw+YSr14d7pYFL7xgZzTv\nvOM6Egmr7GxrWXDXXa4jSZ2vv4Zzz7V+PEEpYammf4Q337QNHsKc8ME6Bk6dahuoiyTbV1/ZHs3X\nX+86ktQ64wzb53foUNeRpE6ok34sZheewjwyKVK1qvXjUZ99SYUhQ2xR4ymnuI4k9fr1sxYTYb1G\nFuqk/9ln1iMkM9N1JN644w745z/D+2IVN2Ix23ciCoMnsM3dd++GBQtcR5IaoU76RS/UcqH+Xx7W\nrp2tlJw3z3UkEiYzZlhTstatXUfijXLl4NZb4fXXXUeSGqG9kPv999CggfXbCdJ2iIl67DG7GDV4\nsOtIJCxuvRVatAh+R82yyMuzKd4bN/q/bYsu5BYaNQquuCJaCR/gttvg7bdh3z7XkUgY7N5t7chv\nucV1JN5q1Mi62L7/vutIki+0Sf+NNywBRs1ZZ9mobOxY15FIGIwZY2XDIG+UEq/bbw/nSvdQJv1N\nm6xVarduriNxo+iCrkii3ngjeqP8Ijk5tvZl2zbXkSRXKJP+W29Br15QubLrSNzIzrZ2y1u2uI5E\ngmzzZpg713o7RVHVqjZwHDbMdSTJFcqkP3SoLVaKqipVbM7+8OGuI5Egi/rgCazEE7ZZPKFL+l98\nYb3zr7zSdSRu3XSTFmpJYqI+eALbQ3fbNpsFGBahS/pDh1rCi8rc/GPp0MGmm61Y4ToSCaJFi6xv\nVfv2riNxKz0dbrghXCWeZKTGLsByYCXwSAn/ngl8C3xe+PGHJByzRIcO2eg2qheeiktPt77nYds0\nRryhwdNhffta0g/LSvdEf6XpwGAs8Z8P9AXOK+F+HwEtCz8eT/CYx5SbC7Vq2bZncrjEE5YXq3ij\noMBeN1Ev7RRp2dJWJM+Z4zqS5Eg06V8GrALWAgeBt4EeJdzPk5W/qkH+WKtW9nnuXLdxSLDk5tq8\n/AsvdB2JP6Sl2b7AYSnxJJr06wEbin29sfB7xcWAtsBC4H3sjCDp9u+Hd9+1+puYtDRd0JWye/NN\ne93IYX37wogRVkIOuowEH1+awsF8oAHwA3AN8C7QtKQ7Dhw48D+3MzMzySxDe8zJk+Gii6Bu3VI/\nJBJuvNFmMj39NGQk+tuW0DtwwAZPxd6KAjRtarklN9e2VXQpNzeX3NzcuB+faNmlNTAQq+kDDAAK\ngD8f5zFrgEuAnUd8P6GGazfdZNsG3n9/3E8RWpdeCk88AZ06uY5E/G7iRHj8cZg503Uk/vP007B8\nuW2/6ideN1ybCzQBGgIVgD7AkV1fahUL6LLC20cm/ITs3QsTJtiyaTnaDTfYqanIiQwfHv7dseLV\np4/1Itq/33UkiUk06ecD/YHJwFJgOLAM6Ff4AZADLAIWAM8CSa+6T5wIl1xiM3fkaDk59mI9eNB1\nJOJn+/dboz4NnkrWoAFccIFtEB9kyajyTiz8KK74fvJ/K/xImREjNDo5nrPOgsaNrXlU586uoxG/\nmjLFZuxErR15WRTN2e9R0hzFgAj80os9e+wvb+/eriPxt+uug3fecR2F+JlKOyeWk2OVhR9+cB1J\n/AKf9CdMgMsvh5o1XUfibyrxyPHs2wfjx2vwdCKnn27rXyZPdh1J/AKf9FXaKZ3iJR6RI02ebJvv\n1KnjOhL/y8mBkSNdRxG/QCf93butDpmd7TqSYFCJR45FpZ3Sy862CkNQtyQNdNIfNw5++lOoUcN1\nJMGQk2MLb1TikeL27rW9YFXaKZ3ataF5cxtwBlGgk75KO2Vz1lm24bNKPFLcxIlWpz7jDNeRBEeQ\nSzyBTfp79sD06dC9u+tIgkUlHjnSO+/Y60JKr1cvqzQcOOA6krILbNKfNAnatIFTT3UdSbCoxCPF\n7dtnI/2ePV1HEiz16sF558G0aa4jKbvAJv3Ro3UBNx4q8Uhx06ZZo0KtZi+7oJZ4Apn09++3C09B\nXhXnUu/e9kdTZMwYDZ7i1bs3vPde8M6aA5n0p0+3HhiaUxyf7Gx7sRYUuI5EXMrPt147SvrxOfNM\nO2tOoMuxE4FM+hqdJKZxY1vB/MknriMRl2bOhPr1oWFD15EEVxBLPIFL+ocO2ShVST8x2dn2x1Oi\na/Rom4Ui8evd295HQdpRK3BJf+ZM28HmnHNcRxJsRUlfm6ZHUyymM+ZkOOccKzPPnu06ktILXNLX\nrJ3kaNHCarqLF7uORFyYNw8qV4bzU7JjdbRkZ9s06KAIVNIvGp3olDRxaWkq8URZUWknLdENU4We\nPS3pB+WsOVBJf/58qFjRZu5I4pT0o0tnzMnTvLmdNS9Z4jqS0glU0i96oWp0khxXXAGbNsHata4j\nES8tW2ZtTFq1ch1JOKSlHR7tB0Hgkr5KO8mTng5ZWcF5sUpyjB5tSapcoN79/qaknwLLlln//Esv\ndR1JuKjEEz0aPCXfT39qZ8wbNriO5MQCk/SLppdpdJJcV18NCxfCtm2uIxEvrFsH69dDu3auIwmX\njAzo1s3WEPldYFLomDHqBJgKlSpB5862HF/C77334NprLUlJcgWlxBOIpL95M+TlwZVXuo4knFTi\niY5x49SoMFU6d4ZPP4VvvnEdyfEFIumPHw9dukD58q4jCaeuXeHjj+2aiYTXt9/CnDnQqZPrSMKp\nShXo0MH2z/WzQCT9sWO1Q1YqVasGbdvCBx+4jkRSadIkq+WffLLrSMIrCCUe3yf9PXtsFNqli+tI\nwi0rS3X9sNPgKfW6dbMN0/fudR3Jsfk+6U+datM0q1d3HUm4ZWXZaWmQugVK6R08aNsiduvmOpJw\nq1kTWrb09zaKvk/6Gp1448wzoUGDYHULlNL7979tw4969VxHEn5+L/H4OukfOmQXcbOyXEcSDd27\nq8QTVuPG6X3klR497Oft17NmXyf9Tz+F009X73yvKOmHUyymM2YvnX02nHGG5S8/8nXSHzdOL1Qv\nXXyxTdv88kvXkUgyLVsGBw5YN0jxRvfulr/8yNdJX6MTb6WlWQnAry9WiU/R+0jdab3j5/eRb5P+\n6tWwfTtcdpnrSKJFJZ7w0eDJe5ddBl9/DWvWuI7kaL5N+uPG2fQyNVjzVocO1oBtxw7XkUgybN0K\nS5dC+/auI4mWcuWsx5EfR/u+Taljx2q2gQuVKkHHjvD++64jkWSYMMF6wlSs6DqS6PFrXd+XSf+b\nb+Czz6ztr3hPq3PDQ1M13enUyXodffut60h+zJdJf9IkOx096STXkUTTtdfaUvL9+11HIonYu9dW\nhnbt6jqSaDrpJNtcZfJk15H8mC+TvqZqunXGGbb5fG6u60gkEdOnW0uA005zHUl0+XEWj++S/sGD\nNtJXjxC3NIsn+DRrx72sLOt5lJ/vOpLDfJf0Z8yAxo2hTh3XkURb0UWoWMx1JBKPggKdMftB/frW\n12rWLNeRHOa7pK/RiT/85CdQoYJN35TgmTfP9klo0sR1JOK3Eo+vkr56hPhHWppKPEGmUb5/KOkf\nx9Kl1pmuWTPXkQgo6QeZBk/+UdTTasUK15EYXyV99QjxlyuusGXkmza5jkTKYt06+521bu06EgFb\nndutm39G+75K+jol9Zfy5W2bSr+8WKV0xo2ztRbp6a4jkSJ+KvH4KukvW6YeIX7j16Xkcmwq7fhP\nx47w+eewc6frSHyW9Dt1shkj4h9dutg02u+/dx2JlMa339qWl507u45EiqtcGTIzbc6+a8lI+l2A\n5cBK4JFj3Oe5wn9fCLQ81hNpdOI/1apZm9gpU1xHIqUxeTK0awcnn+w6EjmSX0o8iSb9dGAwlvjP\nB/oC5x1xn65AY6AJ8DPg78d6MvUI8SeVeIJDDdb869pr7Y/ywYNu40g06V8GrALWAgeBt4EeR9yn\nOzCk8PYcoDpQq6Qnq1EjwWgkJbKybIN6v270LCY/31piK+n7U5060LSplUtdSjTp1wM2FPt6Y+H3\nTnSf+gkeVzx09tlQu7Z/N3oWM2sWnHWWLf0Xf/JD2/KMBB9f2s4sR868L/FxAwcO/M/tzMxMMjMz\n4wpKkq/oxdqmjetI5Fi08ZD/ZWVBr14waFD865Fyc3PJTaAFbqLLoFoDA7GaPsAAoAD4c7H7vAjk\nYqUfsIu+7YGtRzxXLKbuXr41Zw7ceScsWeI6EjmWc8+Ft96CSy5xHYkcSyxmZ2OTJsH55yfnOdPs\nr0epc3mi5Z252AXahkAFoA9w5MnLWODWwtutgV0cnfDF5y691PbNzctzHYmUZMUKm1Z78cWuI5Hj\n8UNPq0STfj7QH5gMLAWGA8uAfoUfAO8Dq7ELvi8B9yV4THHAb0vJ5cfGjbPfj1qY+J/rqZt+eomo\nvONzY8fCs8/ajkziL+3bw29+Y9MCxd/274datWDlSjj99MSfz+vyjkTI1VfD3Lm2cb34x86dsGAB\ndOjgOhIpjYoVrfvAhAlujq+kL6VWpYqNKCdNch2JFDdxoi3xr1zZdSRSWi6nbirpS5m4vgglR9NU\nzeDp2hWmTYN9+7w/tpK+lEm3bv5YSi7mwAH44AP7vUhw1KwJzZu7uT6mpC9lUqeObVzveim5mBkz\nbGl/7dquI5GycnXWrKQvZaYSj3+otBNcRY0MCwq8Pa6SvpRZUdLXDFu3YjHtNhdkTZvCKafA/Pne\nHldJX8qsWTPruLl0qetIom3pUvs9NGvmOhKJl4u25Ur6UmZ+WEouh3vnaxVucLl4HynpS1y0sYp7\nqucHX5s2sHEjrF/v3TGV9CUu7dtbeWGrWuc58fXX9vNX9/FgS0+31hleDqCU9CUuFSrY5tuulpJH\n3fvvW1uMihVdRyKJ8rrEo6QvcVOJxx2VdsKjc2eYPRu++86b4/npEpC6bAbMjh22leLWrer74qV9\n+6xL46pVyenSKO5dc41tUnTddWV/rLpsimdOOw1atlSrZa/l5to0TSX88PCyxKOkLwlRicd7Ku2E\nT1aWXafJz0/9sZT0JSFFuwB5vZQ8qmIxGD9eST9s6teHhg1h1qzUH0tJXxLiail5VC1caDOnzjvP\ndSSSbF6VeJT0JWEq8Xjnvfe0Cjesune332+q57Mo6UvCXO4CFDXvvgvZ2a6jkFRo0cJmZn35ZWqP\no6QvCWvTBjZs8HYpeRStXWtL9tu2dR2JpIJXPa2U9CVhGRm2lHz8eNeRhFtRaScjw3UkkipK+hIY\nKvGk3rvvQs+erqOQVMrMhEWLYNu21B3DT5eDtCI3wHbvhnr1YNMmm80jyaXVz9GRk2N7Ht9+e+nu\nrxW54sQpp1it+YMPXEcSTuPHW4M1JfzwS3WJR0lfkkYlntRRaSc6unaFadNg797UPL+SviRN9+7W\natmLpeRR8sMP1t+oWzfXkYgXataEiy+GKVNS8/xK+pI0DRrAOefAxx+7jiRcpkyBVq2gRg3XkYhX\nevWC0aNT89xK+pJU2dmpe7FGlUo70dOzp13HOXgw+c+tpC9J1asXjBmjBmzJkp9vb/4ePVxHIl5K\n5Vmzkr4k1bnnQvXq8OmnriMJh5kz4cwz7UOiJVUlHiV9SbpU1iOjRqWd6OrVy37/yT5rVtKXpCtK\n+lprl5hYTEk/ypo2Tc1Zs5K+JF2LFnDokC0nl/h98QWUKwcXXug6EnElFWfNSvqSdGlpmsWTDGPG\n2ChfvfOjq2hiRDLPmpX0JSWKXqwSv5Ej4brrXEchLrVoYdM2Fy9O3nMq6UtKtGljzcFWrXIdSTAt\nWwbffguXXeY6EnEpLS35JR4lfUmJ9HQrTWi0H59Ro6B3b6vpS7Qp6UtgaOpm/N55x1rsihSdNefl\nJef5lPQlZTIzbb/PTZtcRxIsK1bYJhraFlHAzpp79EjeWbOSvqRMhQq2jaJKPGUzapSdJam0I0V6\n97YL+8mgl5WkVE6OlSqk9EaOVGlHfuyqq6y8s3Zt4s+lpC8p1aWLLTLavNl1JMGwejVs3Ajt2rmO\nRPykfHlb+5KMAZSSvqRUxYq2uUqyTk3DbtQoe3Onp7uORPymTx8YMSLx51HSl5S7/noYPtx1FMGg\n0o4cS/v2sH594rN4lPQl5Tp1guXLYcMG15H427p1Vt5p3951JOJHGRl2QTfREo+SvqRchQq2UEsX\ndI9v1Cibmle+vOtIxK+uvz7xEk8iSb8GMAVYAXwAVD/G/dYCXwCfA9paI6L69FGJ50SGDYMbbnAd\nhfhZu3awZQusXBn/cySS9H+LJf2mwLTCr0sSAzKBloA6iURUhw5WukjGlLMwWrnSyl9XXeU6EvGz\n9HS75pPIaD+RpN8dGFJ4ewhwvK0e1Bw24jIybMFRMmYfhNGwYXbqrlk7ciKJToxIJOnXArYW3t5a\n+HVJYsBUYC5wTwLHk4BTiadksZgl/b59XUciQXDFFbBjh3VijUfGCf59ClC7hO///oivY4UfJbkC\n2AKcXvh8y4EZJd1x4MCB/7mdmZlJZmbmCcKTIGnf3vrwrFoFjRu7jsY/Fi6EffugdWvXkUgQfPxx\nLg0a5PLzn1t/q7JKpOyyHKvVfwXUAT4EfnKCxzwKfA88U8K/xWLaVDX0+veHOnXg90cOGyLsN7+x\n8tcTT7iORILik0/g9ttttF+uXBqUIZcnUt4ZC9xWePs24N0S7lMFOKXw9klAZ0A7p0bYjTfCm29q\n0/QiBQXw9tsq7UjZXH45HDgAn39e9scmkvSfAjphUzY7FH4NUBeYUHi7NlbKWQDMAcZj0zslotq0\ngf37Yf5815H4w6xZULUqNGvmOhIJkrQ0G0C99VYcj01+OHFTeSciHn0UvvsOBg1yHYl7998Pdeuq\n3CVlt3SprXbfvLls5R0lffHcypW2yGTjRqtlR1V+viX82bOhUSPX0UgQvfYa3HWXdzV9kbg0aQIN\nG8LUqa4jcWvyZEv2SvgSrzvvLPtjlPTFiVtugTfecB2FW6+/bjMwRLyk8o44sX27zdXfsAFOOeXE\n9w+bnTvh7LOts2b1Y3WtEimFtDSVdyQAataEK6+M7v65b78N11yjhC/eU9IXZ6Jc4lFpR1xReUec\n2bcP6teHefPgrLNcR+Odoql269erwZokTuUdCYxKlax//Ouvu47EW0OGwM03K+GLGxrpi1MLFthu\nUatXRyMJ5ufDmWfadNXzz3cdjYSBRvoSKC1a2EXdadNcR+KNqVOtpKWEL64o6Ytzd90F//iH6yi8\n8cor8S2oEUkWlXfEuV27bIXuqlU26g+rzZvhggvsAm4U1yZIaqi8I4FTvTpkZcHQoa4jSa3XXrOt\n7pTwxSWN9MUXPvrIOk4uWmRtY8Pm0CE45xxbjHbxxa6jkTDRSF8C6corLTHOKHEjzeCbNAlq1VLC\nF/eU9MUX0tJspD94sOtIUuOll6BfP9dRiKi8Iz7y3Xd2QXfRIqhXz3U0ybNhg01NXb8eTjrJdTQS\nNirvSGBVrWpbwL30kutIkuvll20PXCV88QON9MVXli2DDh2s5XCFCq6jSdy+fdZX6OOP4dxzXUcj\nYaSRvgTaeefZatVRo1xHkhxvvgmtWinhi38o6Yvv9O8Pzz/vOorExWLw7LPwi1+4jkTkMCV98Z2s\nLNiyxTYMD7Lp0y3xX32160hEDlPSF9/JyIBf/hL+7/9cR5KYolF+GBebSXD56eWoC7nyH3v22B6y\nM2YEsx6+fLktOFu3DipXdh2NhJku5EoonHQS3HcfPPOM60ji89RT8MADSvjiPxrpi29t3w5Nm9r2\ngrVru46m9NassRk7eXna+FxSTyN9CY2aNeGmm6w2HiR/+Yu1XFDCFz/SSF98bf16aNnSauSnn+46\nmhPbvBkuvBC+/DIY8UrwlXWkr6Qvvnf//VClSjBm8zz8MBQUBO/sRIJLSV9CZ9MmuOgiWLLE37X9\nTZugWbPwNYwTf1PSl1B66KHDK1z9ql8/qFbNavoiXlHSl1D66ivbX3bhQqhf33U0R1uxAtq2tc81\nariORqJ1CMEQAAAFQklEQVRESV9C63e/swulr7/uOpKjXX+9XXAeMMB1JBI1SvoSWrt32+rc996D\nSy91Hc1hc+ZAr142ylfPfPGa5ulLaJ1yCjz+uPWz8cv44NAhm1301FNK+BIMSvoSKLfdBnv3wvDh\nriMxr70GlSrBzTe7jkSkdFTekcD597+hTx+bwuly1evOnbbpy+TJtgeuiAuq6Usk/PznkJ8Pr7zi\nLoa774aKFeFvf3MXg4iSvkTCd9/ZFM5//Quuusr740+aBPfeC198YRu6i7iiC7kSCVWrwgsvwD33\n2KweL+3aZcd99VUlfAkejfQl0O66Cw4csBG/FztUxWJw441w6qn2R0fENY30JVKeew7mz4chQ7w5\n3t//bv39g7q5i4hG+hJ4S5ZAZqbV2S+5JHXH+fRT6NYNZs6EJk1SdxyRstBIXyLnggvgpZegRw/Y\nsCE1x1izBrKz4eWXlfAl2DJcByCSDL16werVNhL/8MPkNj3bvh26dLHePz17Ju95RVzQSF9C45e/\nhM6doWNHS9TJ8PXX0KkT5ORYuwWRoEsk6V8HLAEOARcf535dgOXASuCRBI4nclxpadbL/pprbO7+\nunWJPd+aNXDlldC9u/X8EQmDRJL+IiAb+Pg490kHBmOJ/3ygL3BeAseMhNzcXNch+EZZfxZpafCn\nP8Edd0Dr1jBtWnzHnTTJHt+/P/zxj95MBz0RvS4O088ifokk/eXAihPc5zJgFbAWOAi8DfRI4JiR\noBf0YfH8LNLSbK/aoUOtQVu/frBjR+keu20b3HmnLb4aOdKSvl/odXGYfhbxS3VNvx5QfD7FxsLv\niaRcx442nTMjw2bcPPwwzJ17dFvmggKbjvnf/239+qtWtbn47dq5iVsklU40e2cKUNJW1L8DxpXi\n+TXxXpyqVs0aog0YYJ9vvtkuzjZtasl9927b/KR2bZuSuXgx1K3rOmqR1ElGpfJD4JfA/BL+rTUw\nEKvpAwwACoA/l3DfVUCjJMQjIhIleUBjLw/4IXCsdZAZWEANgQrAAnQhV0QkkLKxev1e4CtgYuH3\n6wITit3vGuBLbCSvbaNFRERERKJCi7dMA6xUtgRYDDzgNhxfSAc+p3STBsKsOjASWAYsxa6VRdUA\n7D2yCHgLqOg2HE+9BmzF/u9FamATblYAH2CvFV9Lx8o+DYHyRLvmXxso2mn1ZKwkFtWfRZGHgTeB\nsa4DcWwIcGfh7QygmsNYXGoIrOZwoh8O3OYsGu+1A1ry46T/F+A3hbcfAZ7yOqiyagNMKvb1bws/\nBN4FOroOwqH6wFTgKqI90q+GJTqxUe2XwKnYH79xwNVOI/JeQ36c9JcDtQpv1y78+rhcN1zT4q2S\nNcT+os9xHIdLg4BfY1N8o+xsYBvwT2xa9CtAFacRubMTeAZYD2wGdmEDgyirhZV8KPxc6zj3Bdwn\nfS3eOtrJWP32QeB7x7G40g34Gqvn+6DrjVMZWEPDFwo/7yG6Z8ONgF9gg6K62HvlJpcB+UyMUuRU\n10l/E3YBs0gDbLQfVeWBUcBQrLwTVW2B7sAaYBjQAfiX04jc2Vj48Vnh1yM5flfbMGsFzAJ2APnA\naOy1EmVbOdw1oQ42WPI1Ld46LA1LbINcB+Iz7Yl2TR+sk23TwtsDKXlFexQ0x2a2VcbeL0OAqO1y\n0JCjL+QWzXr8LQG4kAtavFXkp1j9egFW1vicw+0roqw9mr3THBvpL8RGt1GdvQM2U6VoyuYQ7Ow4\nKoZh1zIOYNdC78Aubk8lQFM2RUREREREREREREREREREREREREREREREREQC6f8BSu5xhUIF2d0A\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10631ce50>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "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+/AAAGxBJREFUeJzt3Xl8VOW9x/FPFFAEZRELClgWqwIiynVhERmVKlKpvpTC\nhV5FEbfihr1W1Aqp3itqL221VluVJQpSUJQSFZFtoCCbIovIlmgVRMAajSAgEub+8ZuQgAFn5pw5\nz8yZ7/v1mldmkjPn/Bzhy5PnPAuIiIiIiIiIiIiIiIiIiIiIiIiIiARuFLAFWFnFz34N7AXqB1qR\niIjsc1iCx40Gulfx/abAT4GPfatIRETSqhnfb5m/BJwOfIRa5iIiziTaMq/K5cBGYIVPtYiISIqq\npfi+o4D7sC6WcnneyxERkVSkGuYtsW6X5fHXTYB3gXOArfsd2LJlrLi4ONX6RERyVTFwUqIHp9rN\nshJoCDSPPzYC7TkgyAGKi4uJxWJ6xGIMGzbMeQ2Z8tBnoc9Cn8WhH1ijOWGJhvl44G3gZGADcN0B\nP48lc1EREfFXot0sfX/g5y28FiIiIqnzMppFkhSJRFyXkDH0WVTQZ1FBn0XqghiBEov3/4iISILy\n8vIgiYxWy1xEJAQU5iIiIaAwFxEJAYW5iEgIKMxFREJAYS4iEgIKcxGREFCYi4iEgMJcRCQEFOYi\nIhlmzZrk36MwFxHJIBMnQteuyb9PYS4ikiGefRYGD4bp05N/b6o7DYmIiI8eewyefhrmzIGTEt5f\nqILCXETEoVgM7rsP/vEPmDcPGjdO7TwKcxERR8rKYNAgWLoU5s6FBg1SP5fCXETEgd274ZprYOtW\nmDkTjj7a2/kU5iIiAduxA3r1gurV4Y034MgjvZ8z0dEso4AtwMpK3/s9sBpYDrwC1PFejohIuJWW\nwiWXWJfKpEn+BDkkHuajge4HfO8toA3QDlgH3OtPSSIi4bRlC1xwAZx5JowZA9V87BtJNMz/CXx5\nwPemA3vjzxcBTfwqSkQkbD76CM47Dy6/HB5/HA7zeZaPX6cbALzh07lEREJlxQro0sUmBA0bBnkJ\nb9OcOD8a+fcDu4EXfTiXiEiozJsHV10Ff/4z9O6dvut4DfNrgR7ARYc6KD8/f9/zSCRCJBLxeFkR\nkcxXWAjXXw/jxsFPf3roY6PRKNFoNOVrJdPYbwYUAm3jr7sDI4CuwL8P8b5YLBZLqTgRkWxVUAD3\n3ANTpsA55yT//jzri0k4oxM9cDwW2g2wIYrDsNErNYCS+DELgF9V8V6FuYjklBEj4IknYNo0OPXU\n1M6RrjD3QmEuIjkhFoMhQ6x7Zdo0aNo09XMlG+aaASoi4oM9e+Cmm2DVKvjnP+HYY4O9vsJcRMSj\nnTuhb1/7OmMG1K4dfA3anEJExIPSUujeHWrWtO4VF0EOCnMRkZRt3mxbvJ1+ug0/rFHDXS0KcxGR\nFBQX2/T8K6+0kSt+T89PlvrMRUSStHw59OgBv/0t3HKL62qMwlxEJAmzZ0OfPvDkk+mdnp8sdbOI\niCRowgQL8gkTMivIQS1zEZGEPP44/P73MH06tGvnuprvU5iLiBzC3r02q3PKFJg/H378Y9cVVU1h\nLiJyELt3w4AB8OGHFuRBz+pMhsJcRKQK27bZOuQ1a9qszqOOcl3RoekGqIjIATZvhkgEmje3TZcz\nPchBYS4isp/166FzZ9ur869/9XfT5XTKkjJFRNJv8WIL8YcegoEDXVeTHIW5iAgwdSr07w8jR0LP\nnq6rSZ66WUQk540ZA9ddZ8MPszHIQS1zEclhsRg8/DA89xzMmQOnnOK6otQpzEUkJ5WVwe232/jx\n+fPhhBNcV+RNot0so7CNnFdW+l59YDqwDngLqOtvaSIi6bFjh40hX7PGWuTZHuSQeJiPBrof8L0h\nWJifDMyMvxYRyWhbt8IFF0CdOnbTs04d1xX5I9Ew/yfw5QHf+zlQEH9eAFzhV1EiIumwdi107Gjb\nvI0Z43ZnIL956TNviHW9EP/a0Hs5IiLpMW8e9OplNzwHDHBdjf/8ugEaiz+qlJ+fv+95JBIhEon4\ndFkRkR/20kswaBC88AJcconraqoWjUaJRqMpvz8viWObAYVA2/jrNUAE2AwcD8wGTq3ifbFY7KA5\nLyKSNrEYjBhha5EXFsIZZ7iuKHF5eXmQREZ7mTQ0Begff94fmOzhXCIiviofevj88/D229kV5KlI\nNPXHA12BBlj/+FDgH8BE4ETgX0Bv4Ksq3quWuYgE6ptvoF8/+zppUnaOWEm2ZZ5MN0uqFOYiEpgt\nW2xKfuvW8Mwz2TtiJchuFhGRjFI+9LBHDxg9OnuDPBWazi8ioVA+9HD4cFs0K9cozEUk602cCLfe\nCmPHwsUXu67GDYW5iGStWAweeQSefhqmT4d27VxX5I7CXESy0u7dcPPNsHw5LFwYjsWyvFCYi0jW\nKSmxVQ+POQbmzoVatVxX5J5Gs4hIVikqgk6doH17eOUVBXk5hbmIZI158+C88+DOO22a/uGHu64o\nc6ibRUSywrhxMHhwbo9YORSFuYhktFgMfvc7KCiAWbPgtNNcV5SZFOYikrF27YLrr4fiYhux0lC7\nJhyU+sxFJCN9/jl062ZDEGfPVpD/EIW5iGScNWtsjZUuXWDCBKhZ03VFmU9hLiIZZdYs6NoV7rvP\n1lk5TCmVEH1MIpIxRo6Evn3h738P5z6d6aQboCLi3J49cPfd8PrrNqPzlFNcV5R9FOYi4lRpKfzn\nf1qgL1oE9eq5rig7qZtFRJwpKoIOHaBlS3jjDQW5FwpzEXFi9mzo3Nk2XX7ySahe3XVF2c2PML8X\nWAWsBF4EjvDhnCISYn/7m3WtjB8Pt9ziuppw8LqhczNgFtAK+BaYALwBFFQ6Rhs6iwhg/eJ33QVv\nvQWFhfCTn7iuKHMlu6Gz1xugXwPfAUcBZfGvn3o8p4iE0FdfQZ8+9nzhQqhb1209YeO1m6UEGAF8\nAmwCvgJmeC1KRMJl/Xq70XnqqTb8UEHuP68t85bAnVh3SynwEvBLYFzlg/Lz8/c9j0QiRCIRj5cV\nkWwxcyb06wcPPQQ33ui6mswVjUaJRqMpv99rn3kf4KfAwPjrq4EOwKBKx6jPXCRHPfUUPPigzehU\nGy45QfeZrwEeAGoCu4BuwGKP5xSRLLd7N9xxB8yZA/Pn2zhySS+vYb4ceB54B9gLLAWe8VqUiGSv\nLVugVy+bALRwoW26LOnntZslEepmEckR774LV14J/ftDfr5WPPQi6G4WERHA9ui8806bEHTlla6r\nyT0KcxHxpKwMhgyBV16xtcjbtnVdUW5SmItIykpKbP3xsjJYvBiOPdZ1RblLPVoikpJVq+Ccc6BN\nG3jzTQW5awpzEUna5Mk2bnzoUPjDH6Cafsd3Tv8LRCRhe/faJKBRo2z98bPPdl2RlFOYi0hCtm2D\na66BrVutf7xRI9cVSWXqZhGRH7R+PXTsCMcdZyNWFOSZR2EuIoc0ZYrtCHTrrTaG/AhtP5OR1M0i\nIlUqK7NZnGPGWKB36OC6IjkUhbmIfE9JiS1bu2sXvPMONGzouiL5IepmEZH9vPcenHWWjR+fMUNB\nni0U5iKyz/PPw8UXw/DhMGKExo9nE/2vEhF274bBg2H6dJg9G047zXVFkiyFuUiO27TJ1h8/7jhY\nsgTq1HFdkaRC3SwiOWzuXJvF+bOfwauvKsizmVrmIjkoFoMnnoCHH4aCAuje3XVF4pXCXCTHbNsG\nAwfCunWwYAG0aOG6IvGDullEcsjKlTbssE4dePttBXmY+BHmdYGXgdXAB4DmiYlkoIICuPBCuP9+\neOYZqFnTdUXiJz+6WR4H3gB6xc9Xy4dziohPdu6E226DefM07DDMvLbM6wBdgFHx13uAUo/nFBGf\nFBXZaofbt9uwQwV5eHkN8+bA58BoYCnwLHCU16JExLtXXoFOneCGG2D8eDj6aNcVSTp57WapBrQH\nbgWWAH8ChgBDKx+Un5+/73kkEiESiXi8rIgczHffwT33WJi/9prt0ymZLxqNEo1GU35/nsfrNwIW\nYC10gPOwML+s0jGxWCzm8TIikoiNG6FPH6hXz9ZZqV/fdUWSqry8PEgio712s2wGNgAnx193A1Z5\nPKeIpOCtt2w2Z8+etv64gjy3eG2ZA7QDngNqAMXAdex/E1Qtc5E02rMHhg2zTSTGjQP1YoZDsi1z\nP8L8hyjMRdJkwwbo2xdq1YIXXoAf/ch1ReKXoLtZRMSRKVNsNmfPnjB1qoI812ltFpEs8+23Nlpl\n8mRb6bBTJ9cVSSZQmItkkaIiG61y4omwdKluckoFdbOIZIm//91mc157rY0hV5BLZWqZi2S4HTvg\njjsgGoVp06B9e9cVSSZSy1wkg33wgc3g3LHDulUU5HIwCnORDBSLwbPPQteucNddMHas1laRQ1M3\ni0iG+eILWxzrww9hzhxo3dp1RZIN1DIXySCzZsEZZ0CzZrBokYJcEqeWuUgG2L0bHnjAulNGj4aL\nL3ZdkWQbhbmIY+vWQb9+cPzxsGwZHHec64okG6mbRcSRWAyeew46d4YBA2x6voJcUqWWuYgDJSVw\n443WKo9GoU0b1xVJtlPLXCRg0ajd5GzSBBYvVpCLP9QyFwnIt9/C0KG2VO2oUdC9u+uKJEwU5iIB\nWL4crr4aWra05+obF7+pm0UkjcrK4JFHoFs3+PWvbYEsBbmkg1rmImlSXAz9+0P16vDOO/DjH7uu\nSMJMLXMRn8Vi8MwzcO65cNVVMHOmglzSz6+W+eHAO8BGoKdP5xTJOps3w8CBsGmTrauikSoSFL9a\n5ncAHwDauVly1qRJNuTwjDNg4UIFuQTLj5Z5E6AH8L/AXT6cTySrfPmlbR6xYIHty9mhg+uKJBf5\n0TL/I3A3sNeHc4lklddeg7Zt4ZhjbF0VBbm44rVlfhmwFXgPiBzsoPz8/H3PI5EIkchBDxXJCuWt\n8fnzbaVD/ZEWr6LRKNFoNOX353m8/sPA1cAe4EjgGGAScE2lY2JlZTEO07gZCYnCQrj5ZhupMnw4\n1KrluiIJo7y8PEgio72GeWVdgf/m+6NZYiUlMerV8/FKIg6UlFT0jY8caVu6iaRLsmHud3u5ytEs\nX37p81VEAjZlivWN169v0/EV5JJp/JwBOif++J6SEmjRwscriQSkpARuv92GGo4fD+ef77oikaoF\n0pOtlrlkm1gMXn7ZWuMNGlhrXEEumSyQtVlKSoK4iog/Pv0UBg2CtWth4kTbCUgk06llLhK3dy88\n/XTFLM5lyxTkkj0CaZkrzCXTrV4NN9xgga5t3CQbBdIyVzeLZKrdu+HBB60/vG9fmDdPQS7ZSS1z\nyVkLFlhrvHlzWLoUmjZ1XZFI6gIJ888/D+IqIokpLYXf/tZWOfzTn+AXv4A8P6fPiTgQSDfLhg1B\nXEXk0GIxGyveujXs2gXvvw+9eyvIJRwCaZl/8kkQVxE5uLVrbbjhv/9t48c7dnRdkYi/AmmZb98O\nO3YEcSWR/e3cCUOH2hDDn/3M9uJUkEsYBRLmTZuqq0WCN3UqnHaaDTtctgwGD4Zq2sJcQiqQP9pN\nm1pXyymnBHE1yXUbN8Kdd8J778Ff/gLdu7uuSCT9AmmZt2gBxcVBXEly2bffwmOP2ezNNm3sBqeC\nXHJFIC3zVq3sV12RdJk61dYaP+UUGz/+k5+4rkgkWIGEeevWMG1aEFeSXFNUZH3ha9famPEePVxX\nJOJGIN0srVrBBx8EcSXJFdu3w7332gbK559vXSoKcsllgYT5iSfCtm3wxRdBXE3CLBaDcePg1FNh\n0yZYuRLuvhtq1HBdmYhbgXSzHHYYnHUWLF4Ml14axBUljJYssS6VXbvgpZc0XlykskBa5mC/Di9c\nGNTVJEw+/hh++Uu44gq47jprFCjIRfbnNcybArOBVcD7wO0HO1BhLskqLbV+8fbt4aST7Cbn9dfb\nb3oisj+vfy2+AwYDbYAOwCCgVVUHnnuutaj27vV4RQm9PXvgqadsmOHmzbBiBfzud1C7tuvKRDKX\n1z7zzfEHwHZgNXBC/Ot+fvQjOOEEePddOPtsj1eVUIrF4PXX7YbmCSfAm2/aBCAR+WF+3gBtBpwJ\nLDrYAd272+QOhbkcaMkSGDIEPvsM/u//bJihlqYVSZxfYV4beBm4A2uh7yc/Px+wTSqmTo0wdGjE\np8tKtlu92jaKWLQIHnjA+sS1GJbkomg0SjQaTfn9frR9qgOvAVOBP1Xx81gsFgNs7YyGDWHNGmjU\nyIcrS9b65BPIz4fCQvjNb+DWW6FmTddViWSOPPvVNOGM9noDNA8YCXxA1UG+nyOOgMsvh4kTPV5V\nstbWrbai4ZlnWr/4+vXWR64gF/HGa5h3Bv4LuAB4L/445Dp1/frZDD7JLaWlMGyYLe1QVgarVsH/\n/A/Ureu6MpFwCOIW075uFrBhZyeeaAtvtW0bwNXFqdJSeOIJe1x6qQ0xbN7cdVUimS/obpakVasG\nN98Mf/5z0FeWIJWWwkMP2WSf9eth/nx4/nkFuUi6BN4yB9iyxSaEFBfDsccGUIEEpnJLvEcPuP9+\nOPlk11WJZJ+Mb5mDjWi54gp48kkXV5d0+OqripZ4URG8/TYUFCjIRYLipGUO9he+QwdYtw7q1w+g\nCkmLzz6zTSGeew4uu8zGjGuXHxHvsqJlDtaCu+oqeOQRVxWIF0VFcNNNttfmzp2wdKm1xBXkIm44\na5mDtepOPx3mzLGt5STzLV0Kjz4Ks2bBLbfAbbfBcce5rkokfLKmZQ5w/PE2C/Cmm7SaYiaLxWwo\n6SWXwM9/bt1jH30EDz6oIBfJFE5b5mATSDp1gv794Ve/CqAaSdiOHfDCC/D441C9us3c7NfPZvKK\nSHol2zJ3HuZgN0E7d4YZM6BduwAqkkPauBH+8he7qdm5s4V4165axVAkSFnVzVLu5JOt9de7tw1x\nk+DFYjaxp29fu4+xc6ftDDV5MkQiCnKRTJcRLfNygwfDsmW2KYF+lQ9GaSmMHQt//autannLLTBg\nANSp47oykdyWld0s5crKoE8f658dOxYOPzzNleWwd96xAJ80CS6+2G5CX3CBWuAimSKrwxzs1/se\nPaBxYxgzRhsV+OnLL2HCBHj2WSgpgRtvtN3utba8SObJ+jAHG0Vx5ZVQq5a10LXWder27LFhhQUF\nFcMLBwyw1rh2uRfJXKEIc7D+22uvtZmGr74KTZr4X1iYrVxpAT5uHDRrZkM/+/SBevVcVyYiicjK\n0SxVOeIIePFF6NULzj0XPGyNlzPWrrXFrtq2ta6qGjXsc1uwwJYdVpCLhFfGtswrmzbNugZ694aH\nH1a3S2VFRbYN38SJtiXbL35hn1PHjupGEclmoelmOdAXX9imv0uWwB/+AD175ubIi7177TMoLLTH\nli3220vv3nDeeQpwkbBwEebdsc2cDweeAx494Oe+hHm5N9+Eu+6yzYAffRT+4z98O3XG+uYbmx1b\nWAivv27dJT172qNjRw3hFAmjoMP8cGAt0A34FFgC9AVWVzrG1zAHG6Hxt7/B8OG2BOu994Zruvl3\n31nre8YMmDkT3n3X7huUB3jLlq4rFJF0CzrMOwLDsNY5wJD418qrlPse5uW+/daGLj72mI1HHzgQ\nrr4aGjRIy+XSZtcuW1p2wQK7YTl3LrRoAd26wUUXQZcuNkxTRHJH0GHeC7gEuCH++r+Ac4HbKh2T\ntjCvuIAF4MiRMGWKLdF6+eW2XGvjxmm9dNL27oV//cta2wsW2GPFCtsTtWNHOP98uPBCLS0rkuuS\nDXOv8yvTm9IJysuzbpauXWHbNhv9MnmybSbcsKEFZJcu1lXRokUwfcyxmN2cLCqCVatg+XJ7rFxp\n656ceaaF9/DhcPbZanmLiDdew/xToGml102BjQcelJ+fv+95JBIhEol4vOzBHX20je7o1cvWelmx\nwlrtr75q+1N+/jm0agWnnQbNm0PTpjYhqXFjqFsXjjkGatc+eP/7nj2wfbs9SkossDdvtsdnn1mr\nu7jYHjVrWv9269a2tG/v3jYG/Nhj0/afLyJZKhqNEvUwocZrN0s17AboRcAmYDEB3AD14uuvraX8\n/vvw8ce2dveGDbBpk60g+PXXtj5MrVoVw/zKg33nTrs5Wbu2/aNRp47tltSokT0aNrTZli1b2kMr\nD4pIqlwMTbyUiqGJI4HhB/w8o8I8EXv22HDAWMweYF9r1oQjjwzPqBkRyVyhnTQkIpJLQrM2i4iI\nJE5hLiISAgpzEZEQUJiLiISAwlxEJAQU5iIiIaAwFxEJAYW5iEgIKMxFREJAYS4iEgIKcxGREFCY\ni4iEgMJcRCQEFOYiIiGgMBcRCQGFuYhICCjMRURCQGEuIhICCnMRkRDwEua/B1YDy4FXAO1FLyLi\niJcwfwtoA7QD1gH3+lJRiEWjUdclZAx9FhX0WVTQZ5E6L2E+Hdgbf74IaOK9nHDTH9QK+iwq6LOo\noM8idX71mQ8A3vDpXCIikqRqP/Dz6UCjKr5/H1AYf34/sBt40ce6REQkCXke338tcANwEbDrIMcU\nAS09XkdEJNcUAycFcaHuwCqgQRAXExGRg/PSMl8P1ABK4q8XAL/yXJGIiIiIiPirO7AGa8Hf47gW\nl5oCs7EuqfeB292WkxEOB96j4iZ6rqoLvIxNvvsA6OC2HKfuxf6OrMQGUxzhtpxAjQK2YP/t5epj\nA1DWYXN66jqoC7C/rEVAM6A6sAxo5aoYxxoBZ8Sf1wbWkrufRbm7gHHAFNeFOFaADesFG1mWq7Oo\nmwEfUhHgE4D+zqoJXhfgTPYP88eA38Sf3wM8EnRR5ToCb1Z6PST+EJiMjf7JVU2AGcAF5HbLvA4W\nYGKt0LVAPewftUKgm9OKgteM/cN8DdAw/rxR/PUhpWuhrcbAhkqvN8a/l+uaYf8CL3Jch0t/BO6m\nYvZwrmoOfA6MBpYCzwJHOa3InRJgBPAJsAn4CvsHP5c1xLpeiH9teIhjgfSFeSxN581mtbH+0TuA\n7Y5rceUyYCvWX+51jkO2qwa0B56Kf/2G3P3ttSVwJ9bYOQH7u/JLlwVlmBgJZGq6wvxT7MZfuaZY\n6zxXVQcmAWOxbpZc1Qn4OfARMB64EHjeaUXubIw/lsRfv4yFei46C3gb+ALYg63C2slpRe5toWL2\n/fFYI8iJatjspWbYWPRcvgGahwXWH10XkmG6ktt95gBzgZPjz/OBR92V4lQ7bKRXTezvSwEwyGlF\nwWvG92+Alo8CHILDG6AAl2I3NYrI7eVxz8P6h5dh3QvvYcM2c11XNJqlHdYy154ANnKjfGhiAfbb\nbK4Yj90r2I3da7wOuyk8gwwYmigiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIijvw/V+nHwhAzfXgA\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106994650>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "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",
      "x + y"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "array([1, 3, 5, 7, 9])"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All operators (like ``+``) actuall call an underlying numpy function: in this case ``np.add``:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add(x, y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "array([1, 3, 5, 7, 9])"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These ufuncs have some interesting and useful properties:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.accumulate(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "array([ 0,  1,  3,  6, 10])"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.accumulate(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "array([0, 0, 0, 0, 0])"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.accumulate(y)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "array([  1,   2,   6,  24, 120])"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.identity"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.multiply.identity"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "1"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "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": 27,
       "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": 27
    },
    {
     "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": 28,
       "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": 28
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Ufunc mini-exercises"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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",
      "slow_sincos(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "array([ 1.10954934,  1.41280477,  1.11025905,  1.1235107 ,  1.19782107])"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# write a fast_sincos function\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 30
    },
    {
     "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",
      "slow_sqdiff(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 31,
       "text": [
        "array([ 0.02843657,  0.22826022, -0.33808932,  0.45751827])"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# write a fast_sqdiff function\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 32
    },
    {
     "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",
      "slow_pairprod(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "array([[  2.,   4.,   6.,   8.],\n",
        "       [  6.,  12.,  18.,  24.],\n",
        "       [ 10.,  20.,  30.,  40.],\n",
        "       [ 14.,  28.,  42.,  56.]])"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# write a fast_pairprod function\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Strategy 2. Using Numpy Aggregates"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.mean()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 36,
       "text": [
        "0.005926372912100737"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.mean(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 37,
       "text": [
        "0.005926372912100737"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.std()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 38,
       "text": [
        "1.0681813993312947"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.var()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 39,
       "text": [
        "1.1410115018773628"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.sum()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 40,
       "text": [
        "0.59263729121007369"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.prod()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 41,
       "text": [
        "2.3263845450932399e-26"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.median(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 42,
       "text": [
        "-0.080455675624693163"
       ]
      }
     ],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.percentile(x, 50)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 43,
       "text": [
        "-0.080455675624693163"
       ]
      }
     ],
     "prompt_number": 43
    },
    {
     "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": 44,
       "text": [
        "[-0.63705411377136534, 0.66386400880543395]"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "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": 45,
       "text": [
        "array([[ 0.37521212,  0.16805226,  0.61283898,  0.21395026,  0.36733656],\n",
        "       [ 0.91795723,  0.02566346,  0.16991705,  0.62027897,  0.3057437 ],\n",
        "       [ 0.751662  ,  0.1830548 ,  0.31503027,  0.85763625,  0.90561931]])"
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.sum(0)  # sum along rows"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "array([ 2.04483135,  0.37677051,  1.0977863 ,  1.69186548,  1.57869957])"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x.sum(1)  # sum along columns"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 47,
       "text": [
        "array([ 1.73739018,  2.0395604 ,  3.01300262])"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.median(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "array([ 0.36733656,  0.3057437 ,  0.751662  ])"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.mean(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "array([ 0.34747804,  0.40791208,  0.60260052])"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Binary ufuncs as aggregates"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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": 50,
       "text": [
        "array([ 1.73739018,  2.0395604 ,  3.01300262])"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "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": 51,
       "text": [
        "array([ 1.73739018,  2.0395604 ,  3.01300262])"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.prod(x, 1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 52,
       "text": [
        "array([ 0.003037  ,  0.00075914,  0.03366703])"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "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": 53,
       "text": [
        "array([ 0.003037  ,  0.00075914,  0.03366703])"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Some of these reduce methods are a bit silly:"
     ]
    },
    {
     "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": 54,
       "text": [
        "array([   46.35634202,  1110.00709617,    16.78187137])"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "A caution: for ``reduce`` methods, the default axis is 0, while for ``sum`` the default axis is ``None``:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.add.reduce(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 55,
       "text": [
        "array([ 2.04483135,  0.37677051,  1.0977863 ,  1.69186548,  1.57869957])"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sum(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 56,
       "text": [
        "6.7899532114887027"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "A quick efficiency note: Beware the built-in Python aggregates!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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: 15.8 \u00b5s per loop\n",
        "100 loops, best of 3: 3.66 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "This is slow for the same reason that loops are slow: built-in ``min``, ``sum``, and ``max`` functions basically have to type-check every element.\n",
      "\n",
      "Make sure to use Numpy's ``np.sum``, ``np.min``, and ``np.max``."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Aggregate Mini-exercises"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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": 58,
       "text": [
        "22.986847893700467"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# implement fast_cubesum\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 59
    },
    {
     "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": 60,
       "text": [
        "0.28895814643282985"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# implement fast_rms\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 61
    },
    {
     "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": 62,
       "text": [
        "-33.0"
       ]
      }
     ],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Implement fast_sillyfunc using ufuncs & aggragates\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Strategy 3: Using Numpy Broadcasting"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*Broadcasting* is a way to extend the usefulness of ufuncs. For example, it allows you to subtract a vector from a matrix without writing a loop:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.random.rand(3, 4)\n",
      "mean = X.mean(0)\n",
      "print(X)\n",
      "print(mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.13033367  0.39084497  0.48740275  0.46706801]\n",
        " [ 0.54222919  0.71130202  0.53520395  0.27169094]\n",
        " [ 0.91013207  0.05874528  0.35727916  0.30515685]]\n",
        "[ 0.52756498  0.38696409  0.45996195  0.34797193]\n"
       ]
      }
     ],
     "prompt_number": 64
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(X - mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[-0.3972313   0.00388088  0.02744079  0.11909608]\n",
        " [ 0.01466421  0.32433793  0.075242   -0.076281  ]\n",
        " [ 0.38256709 -0.32821881 -0.10268279 -0.04281509]]\n"
       ]
      }
     ],
     "prompt_number": 65
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We've subtracted the mean vector from each row of the matrix with a single operation. This is an example of **broadcasting**."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Rules of Broadcasting"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "<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": "heading",
     "level": 3,
     "metadata": {},
     "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": 66,
       "text": [
        "array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "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": 67,
       "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": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y + 10 * x"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 68,
       "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": 68
    },
    {
     "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": 69,
       "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": 69
    },
    {
     "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": 70,
       "text": [
        "array([[[  0.83250853,   0.68313399,   0.37255257,   0.90963716],\n",
        "        [ 10.86953565,  10.49053538,  10.40165644,  10.06197872],\n",
        "        [ 20.95098578,  20.31752505,  20.18094159,  20.78663043]],\n",
        "\n",
        "       [[  0.13724524,   0.31820215,   0.8721424 ,   0.37733029],\n",
        "        [ 10.22608683,  10.18324573,  10.89380149,  10.85604921],\n",
        "        [ 20.36238462,  20.93689935,  20.53204387,  20.63015677]]])"
       ]
      }
     ],
     "prompt_number": 70
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Quick Broadcasting/Ufunc Exercises"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Now, assume you have $N$ points in $D$ dimensions, represented by an array of shape ``[N, D]``, try the following:\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": 71
    },
    {
     "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": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 2. Compute the mean using np.add\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 73
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 3. Compute the standard deviation across the 1000 points\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 74
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 4. Compute the standard deviation using np.add only\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 5. Compute the whitened version of the array\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 76
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Strategy 4: Fancy Indexing and Masking"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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": 77
    },
    {
     "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": 78,
       "text": [
        "array([1, 2, 3, 0, 2, 4, 0])"
       ]
      }
     ],
     "prompt_number": 78
    },
    {
     "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": 79,
       "text": [
        "array([False, False, False,  True, False, False,  True], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "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": 80,
       "text": [
        "array([1, 2, 3, 0, 2, 4, 0])"
       ]
      }
     ],
     "prompt_number": 80
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Often you'll see this done in a single step:"
     ]
    },
    {
     "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": 81,
       "text": [
        "array([1, 2, 3, 0, 2, 4, 0])"
       ]
      }
     ],
     "prompt_number": 81
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "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": 82,
       "text": [
        "array([ 0.38048517,  0.21218316,  0.72218801,  0.12978412,  0.73458475])"
       ]
      }
     ],
     "prompt_number": 82
    },
    {
     "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": 83,
       "text": [
        "array([ 0.38048517,  0.21218316,         nan,  0.12978412,         nan])"
       ]
      }
     ],
     "prompt_number": 83
    },
    {
     "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": 84,
       "text": [
        "array([ 0.38048517,  0.21218316,         inf,  0.12978412,         inf])"
       ]
      }
     ],
     "prompt_number": 84
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# warning: don't replace np.isnan(x) with (x == np.nan):\n",
      "np.nan == np.nan"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 85,
       "text": [
        "False"
       ]
      }
     ],
     "prompt_number": 85
    },
    {
     "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": 86,
       "text": [
        "array([ 0.38048517,  0.21218316,  0.        ,  0.12978412,  0.        ])"
       ]
      }
     ],
     "prompt_number": 86
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([1, 0, -1, -np.inf, np.inf, np.nan])\n",
      "print(\"input   \", x)\n",
      "print(\"x < 0   \", x[x < 0])\n",
      "print(\"x > 0   \", x[x > 0])\n",
      "print(\"isinf   \", x[np.isinf(x)])\n",
      "print(\"isnan   \", x[np.isnan(x)])\n",
      "print(\"isposinf\", x[np.isposinf(x)])\n",
      "print(\"isneginf\", x[np.isneginf(x)])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "input    [  1.   0.  -1. -inf  inf  nan]\n",
        "x < 0    [ -1. -inf]\n",
        "x > 0    [  1.  inf]\n",
        "isinf    [-inf  inf]\n",
        "isnan    [ nan]\n",
        "isposinf [ inf]\n",
        "isneginf [-inf]\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": 87
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Boolean Operations on Masks"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Boolean operators can be used on masks: order of operations can be a gotcha \u2013 be 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": 88,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [ 4,  5,  6,  7],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15]])"
       ]
      }
     ],
     "prompt_number": 88
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[x < 5]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 89,
       "text": [
        "array([0, 1, 2, 3, 4])"
       ]
      }
     ],
     "prompt_number": 89
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[~(x < 5)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 90,
       "text": [
        "array([ 5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15])"
       ]
      }
     ],
     "prompt_number": 90
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[(x < 10) & (x % 2 == 0)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 91,
       "text": [
        "array([0, 2, 4, 6, 8])"
       ]
      }
     ],
     "prompt_number": 91
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[(x > 3) & (x < 8)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 92,
       "text": [
        "array([4, 5, 6, 7])"
       ]
      }
     ],
     "prompt_number": 92
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Counting elements with a mask"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "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",
        "60 elements are greater than 0.5\n"
       ]
      }
     ],
     "prompt_number": 93
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# clip is a useful function:\n",
      "x = np.clip(x, 0.3, 0.6)\n",
      "\n",
      "np.sum(x < 0.3), np.sum(x > 0.6)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 94,
       "text": [
        "(0, 0)"
       ]
      }
     ],
     "prompt_number": 94
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# works for 2D arrays as well\n",
      "X = np.random.random((10, 10))\n",
      "print((X < 0.1).sum())"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9\n"
       ]
      }
     ],
     "prompt_number": 95
    },
    {
     "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": 96,
       "text": [
        "array([[ 0.21305713,  0.58388143,  0.43704044],\n",
        "       [ 0.1925056 ,  0.82110282,  0.77627151],\n",
        "       [ 0.79434441,  0.44246991,  0.48757473]])"
       ]
      }
     ],
     "prompt_number": 96
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.where(x < 0.5)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 97,
       "text": [
        "(array([0, 0, 1, 2, 2]), array([0, 2, 0, 1, 2]))"
       ]
      }
     ],
     "prompt_number": 97
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[x < 0.5]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 98,
       "text": [
        "array([ 0.21305713,  0.43704044,  0.1925056 ,  0.44246991,  0.48757473])"
       ]
      }
     ],
     "prompt_number": 98
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[np.where(x < 0.5)]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 99,
       "text": [
        "array([ 0.21305713,  0.43704044,  0.1925056 ,  0.44246991,  0.48757473])"
       ]
      }
     ],
     "prompt_number": 99
    },
    {
     "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": 100,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [ 4,  5,  6,  7],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15]])"
       ]
      }
     ],
     "prompt_number": 100
    },
    {
     "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": 101,
       "text": [
        "array([1, 4])"
       ]
      }
     ],
     "prompt_number": 101
    },
    {
     "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": 102,
       "text": [
        "array([ 0,  5, 10, 15])"
       ]
      }
     ],
     "prompt_number": 102
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X.diagonal()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 103,
       "text": [
        "array([ 0,  5, 10, 15])"
       ]
      }
     ],
     "prompt_number": 103
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Note: diagonal() is read-only (though this will change in numpy 1.9)\n",
      "X.diagonal() = 100"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "ename": "SyntaxError",
       "evalue": "can't assign to function call (<ipython-input-104-d4b0e796ca47>, line 2)",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-104-d4b0e796ca47>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m    X.diagonal() = 100\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m can't assign to function call\n"
       ]
      }
     ],
     "prompt_number": 104
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Instead, try this:\n",
      "X[range(4), range(4)] = 100\n",
      "X"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 105,
       "text": [
        "array([[100,   1,   2,   3],\n",
        "       [  4, 100,   6,   7],\n",
        "       [  8,   9, 100,  11],\n",
        "       [ 12,  13,  14, 100]])"
       ]
      }
     ],
     "prompt_number": 105
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "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": 106,
       "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": 106
    },
    {
     "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": 107,
       "text": [
        "array([0, 4, 2, 3, 5, 1])"
       ]
      }
     ],
     "prompt_number": 107
    },
    {
     "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": 108,
       "text": [
        "array([[ 0,  1,  2,  3],\n",
        "       [16, 17, 18, 19],\n",
        "       [ 8,  9, 10, 11],\n",
        "       [12, 13, 14, 15],\n",
        "       [20, 21, 22, 23],\n",
        "       [ 4,  5,  6,  7]])"
       ]
      }
     ],
     "prompt_number": 108
    },
    {
     "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": 109,
       "text": [
        "array([[[ 0,  1,  2,  3],\n",
        "        [16, 17, 18, 19]],\n",
        "\n",
        "       [[ 8,  9, 10, 11],\n",
        "        [12, 13, 14, 15]],\n",
        "\n",
        "       [[20, 21, 22, 23],\n",
        "        [ 4,  5,  6,  7]]])"
       ]
      }
     ],
     "prompt_number": 109
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[i2].shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 110,
       "text": [
        "(3, 2, 4)"
       ]
      }
     ],
     "prompt_number": 110
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary: Speeding up NumPy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "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."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}