{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook was put together by [Jake Vanderplas](http://www.vanderplas.com) for UW's [Astro 599](http://www.astro.washington.edu/users/vanderplas/Astr599/) course. Source and license info is on [GitHub](https://github.com/jakevdp/2013_fall_ASTR599/)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Efficient Numerical Computing with Numpy\n", "\n", "In this session, we'll discuss how to make your programs as efficient as possible, mainly by taking advantage of *vectorization* in NumPy." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Aside: the \"Unladen Swallow\"" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.display import YouTubeVideo\n", "YouTubeVideo(\"y2R3FvS4xr4\")" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "html": [ "\n", " \n", " " ], "metadata": {}, "output_type": "pyout", "prompt_number": 2, "text": [ "" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "That video was the original reference behind the now-defunct [unladen swallow](http://en.wikipedia.org/wiki/Unladen_Swallow) project, which had the goal of making Python's C implementation faster.\n", "\n", "Yes, Python is (unfortunately) a rather slow language.\n", "\n", "Here is an example:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# A silly function implemented in Python\n", "\n", "def func_python(N):\n", " d = 0.0\n", " for i in range(N):\n", " d += (i % 3 - 1) * i\n", " return d" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "# Use IPython timeit magic to time the execution\n", "%timeit func_python(10000)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "1000 loops, best of 3: 1.81 ms per loop\n" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "To compare to a compiled language, let's write the same function in fortran and use the ``f2py`` tool (included in NumPy) to compile it" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%%file func_fortran.f\n", "\n", " subroutine func_fort(n, d)\n", " integer, intent(in) :: n\n", " double precision, intent(out) :: d\n", " integer :: i\n", " d = 0\n", " do i = 0, n - 1\n", " d = d + (mod(i, 3) - 1) * i\n", " end do\n", " end subroutine func_fort\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Overwriting func_fortran.f\n" ] } ], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "!f2py -c func_fortran.f -m func_fortran > /dev/null" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "from func_fortran import func_fort\n", "%timeit func_fort(10000)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "100000 loops, best of 3: 20.3 \u00b5s per loop\n" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Fortran is about 100 times faster for this task!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Why is Python so slow?\n", "\n", "We alluded to this yesterday, but languages tend to have a compromise between convenience and performance.\n", "\n", "- C, Fortran, etc.: **static typing** and **compiled code** leads to fast execution\n", "\n", " + But: lots of development overhead in declaring variables, no interactive prompt, etc.\n", "\n", "- Python, R, Matlab, IDL, etc.: **dynamic typing** and **interpreted excecution** leads to fast development\n", "\n", " + But: lots of execution overhead in dynamic type-checking, etc.\n", " \n", "We like Python because our **development time** is generally more valuable than **execution time**. But sometimes speed can be an issue." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategies for making Python fast\n", "\n", "1. Use Numpy **ufuncs** to your advantage\n", "\n", "2. Use Numpy **aggregates** to your advantage\n", "\n", "3. Use Numpy **broadcasting** to your advantage\n", "\n", "4. Use Numpy **slicing and masking** to your advantage\n", "\n", "5. Use a tool like *SWIG*, *cython* or *f2py* to interface to compiled code.\n", "\n", "Here we'll cover the first four, and leave the fifth strategy for a later session." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategy 1: Using UFuncs in Numpy\n", "\n", "A **ufunc** in numpy is a *Universal Function*. This is a function which operates element-wise on an array. We've already seen examples of these in the various arithmetic operations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "x = np.random.random(4)\n", "print x\n", "print x + 1 # add 1 to each element of x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.90839173 0.6614708 0.74607367 0.66114385]\n", "[ 1.90839173 1.6614708 1.74607367 1.66114385]\n" ] } ], "prompt_number": 8 }, { "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": [ "[ 1.81678346 1.32294159 1.49214734 1.32228771]\n" ] } ], "prompt_number": 9 }, { "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.82517553 0.43754362 0.55662592 0.43711119]\n" ] } ], "prompt_number": 10 }, { "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.24692093 0.08460287 -0.08492982]\n" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "These are *binary ufuncs*: they take two arguments.\n", "\n", "There are also many *unary ufuncs:*" ] }, { "cell_type": "code", "collapsed": false, "input": [ "-x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": [ "array([-0.90839173, -0.6614708 , -0.74607367, -0.66114385])" ] } ], "prompt_number": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "np.sin(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": [ "array([ 0.78851565, 0.61427811, 0.67876066, 0.61402009])" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Speed of Ufuncs" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.random(10000)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "%%timeit\n", "# compute element-wise x + 1 via a ufunc \n", "y = np.zeros_like(x)\n", "y = x + 1" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "10000 loops, best of 3: 23.7 \u00b5s per loop\n" ] } ], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "%%timeit\n", "# compute element-wise x + 1 via a loop\n", "y = np.zeros_like(x)\n", "for i in range(len(x)):\n", " y[i] = x[i] + 1" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "100 loops, best of 3: 13.7 ms per loop\n" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "#### Why is NumPy so much faster?\n", "\n", "Numpy UFuncs are faster than Python functions involving loops, because *the looping happens in compiled code*. This is only possible when types are known beforehand, which is why numpy arrays must be typed." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Other Available Ufuncs\n", "\n", "- Trigonometric functions (``np.sin``, ``np.cos``, etc.)\n", "- Scipy special functions (``scipy.special.j0``, ``scipy.special.gammaln``, etc.)\n", "- Element-wise minimum/maximum (``np.minimum``, ``np.maximum``)\n", "- User-defined ufuncs (read more [here](http://scipy-lectures.github.io/advanced/advanced_numpy/#universal-functions))" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.random(5)\n", "print x\n", "print np.minimum(x, 0.5)\n", "print np.maximum(x, 0.5)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.03731972 0.94058919 0.79808336 0.45500933 0.03967752]\n", "[ 0.03731972 0.5 0.5 0.45500933 0.03967752]\n", "[ 0.5 0.94058919 0.79808336 0.5 0.5 ]\n" ] } ], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "# contrast this behavior with that of min() and max()\n", "print np.min(x)\n", "print np.max(x)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.0373197161493\n", "0.940589193882\n" ] } ], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "# On older IPython versions, use %pylab inline\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 19 }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.linspace(0, 10, 1000)\n", "plt.plot(x, np.sin(x))" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 20, "text": [ "[]" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD9CAYAAABQvqc9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYlmXaBvDzFVDT3BdUwMyABDc0zbRBScVcAEXJpUVc\nx2ocy5r5vmrmO0abNG1sbLEmrWmyNCV3NMQ1ysmFFHOlRBRlk1zLVGTx+f64wh3kXe9nOX/H0TGj\nvbzPGb5c3s/13ItN0zQNRERkCVVUByAiIs9h0ScishAWfSIiC2HRJyKyEBZ9IiILYdEnIrIQp4v+\nmDFj4Ovri7Zt25b7mkmTJiEoKAjt27fH7t27nb0kERE5yOmiP3r0aCQnJ5f775OSknD48GFkZGRg\n3rx5eOaZZ5y9JBEROcjpoh8eHo569eqV++8TExMRHx8PAOjSpQvOnTuHgoICZy9LREQOcHtPPzc3\nFwEBAVd/7e/vj5ycHHdfloiIbsPbExe5eacHm812y2tu93tERHRn9uym4/aRvp+fH7Kzs6/+Oicn\nB35+frd9raZphvknP1/DkCEagoM1JCRoKCmp+PUFBRpeeklDgwYa3nlHw5Ur5b/2b3/7m/L/Pr38\nw++Fub8XpaUa/vlP+bn46181nDxZ8euLizV8/rmGevX+hmHDNPz0k/r/BtX/2MvtRT8mJgaffvop\nAGD79u2oW7cufH193X1Zt9q2DXjgASAwENizBxg6FPDyqvhrGjcGXn8d+PZbYMECIC4O+OUXz+Ql\n0qNz54CBA4GlS+Vn6u9/Bxo2rPhrvL2BESOAZ54B/P2BTp2A777zTF6zcLq9M2LECHz99dc4deoU\nAgICMHXqVBQXFwMAJkyYgP79+yMpKQmBgYGoWbMm/vOf/zgdWqXly4Gnnwb+/W8gOtr+r7//fuCb\nb4BJk4DwcGD9esDgfwcS2S0/H+jdG3jkEWDZMqBqVfu+3scHmDYNePhhYMAAx38eLUnTCR1FKVdC\ngqY1aaJpu3c7/15XrmjalCmaFhSkaTk5N/67r776yvkLmAS/F9eY5Xtx/Lh87qdNc/w9rv9epKZq\nmq+v/Hxakb210/bbFylns9kc6k95SmIiMGGCjMwrWIdmt9dfBxYtktF/3bque18iPTp1CujWDfj9\n74E//cl177t3LxAZKa3TyEjXva8R2Fs7WfQrIS0N6NsXSEqSHqIraRrw3HPyoU1OBqpXd+37E+lF\nYaG0dMLDZbDjalu2AEOGyMAsLMz1769XLPoulpcHdOkCvPMOEBvrnmuUlgKPPSYPez/4wD3XIFJJ\n04CRI4HLl4HFi4EqbppCsmSJ3EHs2nXnh8JmYW/t5IZrFSgtBZ54Ahg/3n0FH5CZP598AmzeDHz2\nmfuuQ6TKRx/JTLf5891X8AEZPA0bBjz1FHDlivuuY2Qc6Vfg1VeBr7+W28U7Tcl0hX37gJ495Zqh\noe6/HpEnlH2ut2wBWrVy//WKi4FevYBHHwX+8hf3X081tndcZOtW6Q/u2gU0a+a5637wgUw/27ZN\n5iQTGdnly7Km5c9/Bn7bgssjcnOBDh2Adevkf82M7R0XKCwExo4F3n3XswUfkBlC9eoBM2d69rpE\n7jB9uixiHDnSs9f18wPefBMYNQooKvLstfWOI/3b+L//A/bvl4VYKrYEys4GOnYEvvoKaNPG89cn\ncoWyts7330sR9jRNkxW/YWHSqjUrtnectHev9AP37PH8KP96778PJCQAKSlq/uIhckZpKdC1q0yC\nGD9eXY68PKBdO2nXBgery+FObO84QdOAZ58FXntNbcEHpM1z/rws3CIymo8+kjUn48apzdGsGfDy\ny7IWRgdjSl3gSP86X3whi0Z27vTMbJ072bZNNmZLTwdq11adhqhyzp2TWTpr1+rjIWpxMdC+PTBj\nBhATozqN67G946BLl4CQEJkvHxGhLMYtRo8GGjUC3nhDdRKiyvnTn6Twf/SR6iTXbNwoWz8cPGi+\nVe8s+g6aPl1G+MuXK4twW3l5stfP998D1x1ARqRLGRnSyz9wQH+7xw4cKAO6yZNVJ3EtFn0HnDwp\nt6M7dsj0Mr155RWgoEDm7xPp2dCh0tJ5+WXVSW61f79M0sjIMFe7lEXfAX/+M3DhgsyY0aNz52Tm\nQUoKV+qSfu3ZI6tgMzOBmjVVp7m9UaOA5s3NNYWTRd9OJ05IId23T81c4sqaNUtO3VqxQnUSotsb\nNAjo0UPf7ZOsLFkhfPCg/tpPjmLRt9Nzz8k8+Lfe8vil7XLpkrSe1qzRx4wIouvt2iUzYw4fBu66\nS3Waik2aJCdvvfmm6iSuwaJvh5wcmcp14ADQpIlHL+2Q2bNltL90qeokRDcaMADo3x/4wx9UJ7mz\nnBxZsHXokDm2X2bRt8PkybLNq1H+xr9wAWjZUrZgbt1adRoi8d13sjlhRgZQrZrqNJXz+99Le+fv\nf1edxHks+pV0+jQQFKT/Xv7NXn9dZiEsXKg6CZGIi5PTsJ57TnWSysvMlMORMjOBOnVUp3EOi34l\nvfoqcOyY8aZB/vKLjPa3bZO/tIhUysiQM2+zsvQ7Y6c8Tz4pd8x6nF5qDxb9Sihrk3z9tWcOdXC1\nv/4VOHsWeO891UnI6p5+Wo75NOIUyAMHZBfQo0eBGjVUp3Eci34lvPuubFust9W3lZWfL9NMMzOB\n+vVVpyGrKpvu/OOPslWIEcXEAFFR0uM3Ku6yeQclJfLg9n//V3USxzVtKkvK585VnYSs7N13gREj\njFvwAeD552W6tj6Gvp5huaKfmCgPbrt0UZ3EOZMnA3Pm8FQgUuPiRWDePOCFF1Qncc4jj8ic/Q0b\nVCfxHMsV/TlzgIkTVadwXvv28jziiy9UJyErWrQIeOgh4L77VCdxjs0mo/3Zs1Un8RxLFf39+4Ef\nfpA5xWbwwgvAP/9prVtTUk/TzDN4AqRFlZYm51ZYgaWK/pw5ciJV1aqqk7hGv35yutb27aqTkJVs\n3Soz4CIjVSdxjerVZRbSO++oTuIZlpm9c/asTNM8eFAehJrFrFmywGz+fNVJyCoef1yeiRlpMdad\nlM2IO3bMeNsuc8pmOWbPluXin3/utksoceqUbMSWmQk0aKA6DZldfr4saDpyBKhbV3Ua14qLA3r3\nllG/kXDK5m1cuSILmczSg7xew4ZAdDRH+uQZH34IDBtmvoIPSOv3gw/M/4zMEkX/q69kiXjXrqqT\nuMfTT1vjw0pqlZTI2pBnn1WdxD169ZJnZKmpqpO4lyWK/kcfAePGyfQsM+rWTR5Gbd6sOgmZWXKy\nnDrVtq3qJO5RpYqszDX7okfT9/RPn5a5xEeOmHvLgvfflzuaJUtUJyGzio2VLQvGjlWdxH1++gm4\n/37Zj8coLSz29G+ycKEc8GDmgg/IjoEbN8oB6kSuduKEnNE8dKjqJO7VuLGc87tggeok7mPqoq9p\n8uDJzCOTMrVry3483Gef3GH+fFnUWKuW6iTuN2GCbDGhjx6I65m66H/3newREhGhOolnjB4N/Oc/\n5v2wkhqaJudOWGHwBMjh7ufPA99/rzqJe5i66Jd9UKuY+r/ymvBwWSm5a5fqJGQmW7bIpmQPPaQ6\niWdUqQKMHAl88onqJO5h2ge5v/4KBATIfjtGOg7RWa++Kg+j5sxRnYTMYuRIICzM+Dtq2iMzU6Z4\n5+Tof9sWPsj9zbJlwMMPW6vgA0B8PLB4MVBYqDoJmcH587Id+VNPqU7iWffdJ7vYJiWpTuJ6pi36\nn30mBdBq7rlHRmWJiaqTkBmsWCFtQyMflOKoUaPMudLdlEU/N1e2So2KUp1EjbIHukTO+uwz643y\ny8TFydqXkydVJ3EtUxb9zz8HBg8G7rpLdRI1YmNlu+X8fNVJyMjy8oCdO2VvJyuqXVsGjosWqU7i\nWqYs+gsWyGIlq6pRQ+bsJySoTkJGZvXBEyAtHrPN4jFd0d+7V/bO795ddRK1nniCC7XIOVYfPAFy\nhu7JkzIL0CxMV/QXLJCCZ5W5+eXp2VOmmx06pDoJGdG+fbJvVY8eqpOo5eUFDB9urhaPqUpjaamM\nbq364Ol6Xl6y77nZDo0hz+Dg6ZoRI6To62NFk/Oc/iNNTk5Gq1atEBQUhJkzZ97y71NSUlCnTh10\n6NABHTp0wGuvvebsJcuVkgL4+sqxZ3StxWOWDyt5xpUr8rmxemunTIcOsiJ5xw7VSVzD25kvLi0t\nxcSJE7Fx40b4+fmhc+fOiImJQUhIyA2v69GjBxI9MHGcPcgbdeok/7tzJ9C5s9osZBwpKTIvv00b\n1Un0wWaTc4EXLTLHVhROjfRTU1MRGBiIFi1awMfHB8OHD8eqVatueZ0ndnq4fBlYuVL6byRsNj7Q\nJfstXCifG7pmxAjgiy+khWx0To30c3NzERAQcPXX/v7+2HHTPZDNZsPWrVvRvn17+Pn5YdasWQgt\np/8yZcqUq/8/IiICEXZsj7luHdCuHdCsmV3/Cab3+OMyk2nWLMDbqT9tsoKiIhk8XfejSACCg6W2\npKTIsYoqpaSkICUlxeGvd6oM2Cpx/mDHjh2RnZ2NGjVqYO3atRg0aBAOlTOlZIoTn7SEBPMf8OCI\n4GDZeO6rr4DISNVpSO82bZI9Z64by9Fvyh7oqi76Nw+Ip06datfXO9Xe8fPzQ3Z29tVfZ2dnw9/f\n/4bX1KpVCzVq1AAA9OvXD8XFxThz5owzl73FpUvAl1/Ksmm61fDhcmtKdCccPJVv2DDZi+jyZdVJ\nnONU0e/UqRMyMjKQlZWFoqIiJCQkICYm5obXFBQUXO3pp6amQtM01Hfx2YVr1wIPPCAzd+hWcXHy\nYS0uVp2E9OzyZdmoj4On2wsIAFq3lgPijcyp9o63tzfmzJmDRx99FKWlpRg7dixCQkIw97fj5CdM\nmIClS5fiX//6F7y9vVGjRg0sXrzYJcGv98UXHJ1U5J57gMBAafH06aM6DenVhg0yY8dq25Hbo6zF\nM3Cg6iSOM/whKhcuyIf08GGgYUM3BDOJN98EfvhBzgwmup2nngK6dAEmTlSdRL9OnpQBVH6+7HGl\nB5Y7ROXLL+WDyoJfMbZ4qCKFhcCaNXL4OZWvUSNZ/7JuneokjjN80Wdrp3Kub/EQ3WzdOjl8p2lT\n1Un0Ly4OWLpUdQrHGbronz8vfcjYWNVJjOGxx4AlS1SnID3irJ3Ki42VDoNRjyQ1dNFfvRr43e8A\nF08GMq24OFl4wxYPXe/SJTkLlq2dymnSBGjfXgacRmToos/Wjn3uuUcOfGaLh663dq30qRs3Vp3E\nOIzc4jFs0b9wAdi8GbhpWQDdAVs8dLMlS+RzQZU3eLB0GoqKVCexn2GLfnIy0LUrUK+e6iTGwhYP\nXa+wUEb6gwapTmIsfn5ASIhsW2E0hi36y5fzAa4j2OKh623aJBsVcjW7/Yza4jFk0b98WR48GXlV\nnEpDhshfmkQrVnDw5KghQ4BVq4x312zIor95s+yBwTnFjomNlQ/rlSuqk5BKJSWy1w6LvmOaN5e7\nZid2OVbCkEWfoxPnBAbKCubt21UnIZW+/Rbw9wdatFCdxLiM2OIxXNEvLZVRKou+c2Jj5S9Psq7l\ny2UWCjluyBD5OTLSiVqGK/rffisn2LRsqTqJsZUVfX1st0eepmm8Y3aFli2lzbxtm+oklWe4os9Z\nO64RFiY93f37VSchFXbtAu66Cyjn5FKyQ2ysTIM2CkMV/bLRCW9JnWezscVjZWWtnUqceEp3MGiQ\nFH2j3DUbquinpQHVqsnMHXIei7518Y7Zddq3l7vmAwdUJ6kcQxX9sg8qRyeu8fDDQG4ukJWlOgl5\nUnq6bGPSqZPqJOZgs10b7RuB4Yo+Wzuu4+UFREcb58NKrrF8uRSpKob66dc3Fn03SE+X/fM7d1ad\nxFzY4rEeDp5c73e/kzvm7GzVSe7MMEW/bHoZRyeu1bs3sGePnP1J5nfsGHD8OBAerjqJuXh7A1FR\nsoZI7wxTQles4E6A7lC9OtCnjyzHJ/NbtQoYMECKFLmWUVo8hij6eXlAZibQvbvqJObEFo91rF7N\njQrdpU8fIDUVOHtWdZKKGaLor1kD9O0L+PioTmJO/fsD33wjz0zIvH7+GdixA4iMVJ3EnGrUAHr2\nlPNz9cwQRT8xkSdkuVOdOkC3bsD69aqTkDslJ0sv/+67VScxLyO0eHRf9C9ckFFo376qk5hbdDT7\n+mbHwZP7RUXJgemXLqlOUj7dF/2NG2WaZt26qpOYW3S03JYaabdAqrziYjkWMSpKdRJza9gQ6NBB\n38co6r7oc3TiGc2bAwEBxtotkCrvv/+VAz/8/FQnMT+9t3h0XfRLS+UhbnS06iTWEBPDFo9ZrV7N\nnyNPGThQvt96vWvWddFPTQUaNeLe+Z7Com9OmsY7Zk+6916gcWOpX3qk66K/ejU/qJ7UsaNM2/zx\nR9VJyJXS04GiItkNkjwjJkbqlx7puuhzdOJZNpu0APT6YSXHlP0ccXdaz9Hzz5Fui/6RI8CpU8CD\nD6pOYi1s8ZgPB0+e9+CDwE8/AUePqk5yK90W/dWrZXoZN1jzrJ49ZQO206dVJyFXKCgADh4EevRQ\nncRaqlSRPY70ONrXbUlNTORsAxWqVwd69QKSklQnIVf48kvZE6ZaNdVJrEevfX1dFv2zZ4HvvpNt\nf8nzuDrXPDhVU53ISNnr6OefVSe5kS6LfnKy3I7WrKk6iTUNGCBLyS9fVp2EnHHpkqwM7d9fdRJr\nqllTDldZt051khvpsuhzqqZajRvL4fMpKaqTkDM2b5YtARo0UJ3EuvQ4i0d3Rb+4WEb63CNELc7i\nMT7O2lEvOlr2PCopUZ3kGt0V/S1bgMBAoGlT1UmsrewhlKapTkKOuHKFd8x64O8v+1pt3ao6yTW6\nK/ocnehDq1ZA1aoyfZOMZ9cuOSchKEh1EtJbi0dXRZ97hOiHzcYWj5FxlK8fLPoVOHhQdqZr21Z1\nEgJY9I2Mgyf9KNvT6tAh1UmEroo+9wjRl4cflmXkubmqk5A9jh2TP7OHHlKdhABZnRsVpZ/Rvq6K\nPm9J9cXHR46p1MuHlSpn9WpZa+HlpToJldFTi0dXRT89nXuE6I1el5JT+dja0Z9evYDdu4EzZ1Qn\n0VnRj4yUGSOkH337yjTaX39VnYQq4+ef5cjLPn1UJ6Hr3XUXEBEhc/ZV01XR5+hEf+rUkW1iN2xQ\nnYQqY906IDwcuPtu1UnoZnpp8Thd9JOTk9GqVSsEBQVh5syZt33NpEmTEBQUhPbt22P37t3lvhf3\nCNEntniMgxus6deAAfKXcnGx2hxOFf3S0lJMnDgRycnJOHjwIBYtWoT09PQbXpOUlITDhw8jIyMD\n8+bNwzPPPFPu+9Wv70wacpfoaDmgXq8HPZMoKZEtsVn09alpUyA4WNqlKjlV9FNTUxEYGIgWLVrA\nx8cHw4cPx6pVq254TWJiIuLj4wEAXbp0wblz51BQUODMZcnD7r0XaNJEvwc9k9i6FbjnHln6T/qk\nh23LvZ354tzcXAQEBFz9tb+/P3bs2HHH1+Tk5MDX1/eW95syZcrV/x8REYGIiAhn4pELlX1Yu3ZV\nnYTKw4OH9C86Ghg8GJg92/H1SCkpKUhxYgtcp4q+rZKptZt27Srv664v+qQvMTHAmDHA66+rTkLl\nWb0a+Pxz1SmoIu3aSU8/PR0IDXXsPW4eEE+dOtWur3eqvePn54fs7Oyrv87Ozob/TfeWN78mJycH\nfn5+zlyWFOjcWc7NzcxUnYRu59AhmVbbsaPqJFQRPexp5VTR79SpEzIyMpCVlYWioiIkJCQg5qZ5\nlzExMfj0008BANu3b0fdunVv29ohfdPbUnK60erV8ufDLUz0T/XUTaeKvre3N+bMmYNHH30UoaGh\nGDZsGEJCQjB37lzMnTsXANC/f3+0bNkSgYGBmDBhAt5//32XBCfPUz1CofJxFa5xREQABw4AJ0+q\nub5Nu7nhrojNZrul90/6cvGizOI5dgyoV091Gipz5ozMsDpxQlZ+kv499pjM2x81yvn3srd26mpF\nLulbjRqyN1JysuokdL21a2X0yIJvHCqnbrLok13Y4tEfTtU0nv79gU2bgMJCz1+bRZ/sEhWlj6Xk\nJIqKgPXr5c+FjKNhQ6B9e2DzZs9fm0Wf7NK0qRxcr3opOYktW2Rpf5MmqpOQvVTdNbPok93Y4tEP\ntnaMq2wjwytXPHtdFn2yW1nR52QrtTSNp80ZWXAwUKsWkJbm2euy6JPd2raVHTcPHlSdxNoOHpQ/\nh7ZtVSchR6nYtpxFn+ymh6XkdG3vfK7CNS4VP0cs+uQQHqyiHvv5xte1K5CTAxw/7rlrsuiTQ3r0\nkPYCj0ZQ46ef5PvP3ceNzctLVuZ6cgDFok8OqVpVDt/+8kvVSawpKQno3RuoVk11EnKWp1s8LPrk\nMLZ41GFrxzz69AG2bQN++cUz1+OGa+Sw06dlo6+CAu774kmFhYCvL3D4MNCokeo05Ar9+skhRY89\nZv/XcsM18pgGDYAOHdQsJbeylBSZpsmCbx6ebPGw6JNT2OLxPLZ2zCc6Wp7TlJS4/1os+uSUslOA\nPL2U3Ko0DVizhkXfbPz9gRYtgK1b3X8tFn1yiqql5Fa1Z4/MnAoJUZ2EXM1TLR4WfXIaWzyes2oV\nV+GaVUyM/Pm6ez4Liz45TeUpQFazciUQG6s6BblDWJjMzPrxR/deh0WfnNa1K5Cd7dml5FaUlSVL\n9rt1U52E3MFTe1qx6JPTvL1lKfmaNaqTmFtZa8fbW3USchcWfTIMtnjcb+VKYNAg1SnInSIigH37\ngJMn3XcNrsgllzh/HvDzA3JzZTYPuRZXP1tHXJyceTxqVOVezxW5pEStWtJrXr9edRJzWrNGNlhj\nwTc/d7d4WPTJZdjicR+2dqyjf39g0ybg0iX3vD+LPrlMTIxsteyJpeRWcvGi7G8UFaU6CXlCw4ZA\nx47Ahg3ueX8WfXKZgACgZUvgm29UJzGXDRuATp2A+vVVJyFPGTwYWL7cPe/Nok8uFRvrvg+rVbG1\nYz2DBslznOJi1783iz651ODBwIoV3IDNVUpK5Id/4EDVSciT3HnXzKJPLnX//UDdukBqquok5vDt\nt0Dz5vIPWYu7Wjws+uRy7uxHWg1bO9Y1eLD8+bv6rplFn1yurOhzrZ1zNI1F38qCg91z18yiTy4X\nFgaUlspycnLc3r1AlSpAmzaqk5Aq7rhrZtEnl7PZOIvHFVaskFE+9863rrKJEa68a2bRJ7co+7CS\n45YuBR57THUKUiksTKZt7t/vuvdk0Se36NpVNgc7fFh1EmNKTwd+/hl48EHVSUglm831LR4WfXIL\nLy9pTXC075hly4AhQ6SnT9bGok+GwambjluyRLbYJSq7a87MdM37seiT20REyHmfubmqkxjLoUNy\niAaPRSRA7poHDnTdXTOLPrlN1apyjCJbPPZZtkzuktjaoTJDhsiDfVfgx4rcKi5OWhVUeUuXsrVD\nN3rkEWnvZGU5/14s+uRWffvKIqO8PNVJjOHIESAnBwgPV52E9MTHR9a+uGIAxaJPblWtmhyu4qpb\nU7Nbtkx+uL28VCchvRk2DPjiC+ffh0Wf3G7oUCAhQXUKY2Brh8rTowdw/Ljzs3hY9MntIiOBH34A\nsrNVJ9G3Y8ekvdOjh+okpEfe3vJA19kWD4s+uV3VqrJQiw90K7ZsmUzN8/FRnYT0auhQ51s8LPrk\nEcOGscVzJ4sWAcOHq05BehYeDuTnAxkZjr+Hw0X/zJkziIyMRHBwMPr06YNz587d9nUtWrRAu3bt\n0KFDBzzIjUQsq2dPaV24YsqZGWVkSPvrkUdUJyE98/KSZz7OjPYdLvozZsxAZGQkDh06hF69emHG\njBm3fZ3NZkNKSgp2796NVJ6hZ1ne3rLgyBWzD8xo0SK5deesHboTZydGOFz0ExMTER8fDwCIj4/H\nypUry32txiOUCGzxlEfTpOiPGKE6CRnBww8Dp0/LTqyO8Hb0wgUFBfD19QUA+Pr6oqCg4Lavs9ls\n6N27N7y8vDBhwgSMHz++3PecMmXK1f8fERGBiIgIR+ORDvXoIfvwHD4MBAaqTqMfe/YAhYXAQw+p\nTkJG8M03KQgISMGzzzo208umVTAMj4yMxIkTJ275/WnTpiE+Ph5nz569+nv169fHmTNnbnltfn4+\nmjZtipMnTyIyMhLvvvsuwm+z3NBms/GOwAImTgSaNgX+8hfVSfTjf/5H2l/Tp6tOQkaxfTswapSM\n9qtUsa92VjjS37BhQ7n/ztfXFydOnECTJk2Qn5+Pxo0b3/Z1TZs2BQA0atQIsbGxSE1NvW3RJ2t4\n/HFg3DjglVd4DCAAXLkCLF4MfPml6iRkJF26AEVFwO7d9n+twz39mJgYzJ8/HwAwf/58DBo06JbX\nXLx4EefPnwcAXLhwAevXr0fbtm0dvSSZQNeuwOXLQFqa6iT6sHUrULs2wB8LsofNJgOozz+3/2sd\nLvovvfQSNmzYgODgYGzevBkvvfQSACAvLw8DBgwAAJw4cQLh4eEICwtDly5dEBUVhT59+jh6STIB\nmw148klgwQLVSfSBD3DJUY8/Lp8fe1XY0/ck9vStIyNDFpnk5Egv26pKSoBmzYBt24D77lOdhozo\n44+BsWPtq51ckUseFxQEtGgBbNyoOola69ZJsWfBJ0eNGWP/17DokxJPPQV89pnqFGp98onMwCDy\nJLZ3SIlTp2SufnY2UKuW6jSed+YMcO+9srNm3bqq05CR2Vs7OdInJRo2BLp3t+75uYsXA/36seCT\n57HokzJWbvGwtUOqsL1DyhQWAv7+wK5dwD33qE7jOQcPysEyx49zgzVyHts7ZBjVq8v+8Z98ojqJ\nZ82fL2sVWPBJBY70Sanvv5fToo4csUYRLCkBmjeX6aqhoarTkBlwpE+GEhYmD3U3bVKdxDM2bpSW\nFgs+qcKiT8qNHQv8+9+qU3jGhx86tqCGyFXY3iHlzp2TFbqHD8uo36zy8oDWreUBrhXXJpB7sL1D\nhlO3LhAdbf5N2D7+WI66Y8EnlTjSJ134+mvgD38A9u0z5z77paVAy5ayGK1jR9VpyEw40idD6t5d\nCuOWLaqGDIv0AAAJT0lEQVSTuEdyMuDry4JP6rHoky7YbDLSnzNHdRL3mDsXmDBBdQoitndIR375\nRR7o7tsH+PmpTuM62dkyNfX4caBmTdVpyGzY3iHDql1bTgOaO1d1EteaN09Ox2LBJz3gSJ90JT0d\n6NlTthyuWlV1GucVFsq+Qt98A9x/v+o0ZEYc6ZOhhYTIatVly1QncY2FC4FOnVjwST9Y9El3Jk4E\n3n1XdQrnaRrw1lvA88+rTkJ0DYs+6U50NJCfLweGG9nmzVL4e/dWnYToGhZ90h1vb+DFF4F//EN1\nEueUjfLNuNiMjIsPckmXLlyQM2S3bDFmP/yHH2TB2bFjwF13qU5DZsYHuWQKNWsCzz4LvPmm6iSO\nmTEDmDSJBZ/0hyN90q1Tp4DgYDlesEkT1Wkq7+hRmbGTmcmDz8n9ONIn02jYEHjiCemNG8kbb8iW\nCyz4pEcc6ZOuHT8OdOggPfJGjVSnubO8PKBNG+DHH42Rl4yPI30ylebN5fD0N95QnaRyZs0CRo5k\nwSf94kifdC83F2jXDjhwQN+9/dxcoG1b820YR/pmb+1k0SdDmDz52gpXvZowAahTxzh3JWQOLPpk\nSidOyPmye/YA/v6q09zq0CGgWzf53/r1VachK2HRJ9N65RV5UPrJJ6qT3GroUHng/PLLqpOQ1bDo\nk2mdPy+rc1etAjp3Vp3mmh07gMGDZZTPPfPJ0zh7h0yrVi3gtddkPxu9jA9KS+WYxxkzWPDJGFj0\nyVDi44FLl4CEBNVJxMcfA9WrA08+qToJUeWwvUOG89//AsOGyRROlatez5yRQ1/WrZMzcIlUYE+f\nLOGZZ4CSEuDDD9VlGDcOqFYNeO89dRmIWPTJEn75RaZwfvop8Mgjnr9+cjLw9NPA3r1yoDuRKnyQ\nS5ZQuzbw/vvA+PEyq8eTzp2T6370EQs+GQ9H+mRoY8cCRUUy4vfECVWaBjz+OFCvnvylQ6QaR/pk\nKe+8A6SlAfPne+Z6//qX7O9v1MNdiDjSJ8M7cACIiJA++wMPuO86qalAVBTw7bdAUJD7rkNkD470\nyXJatwbmzgUGDgSys91zjaNHgdhYYN48FnwyNm/VAYhcYfBg4MgRGYl/9ZVrNz07dQro21f2/hk0\nyHXvS6QCR/pkGi++CPTpA/TqJYXaFX76CYiMBOLiZLsFIqNj0SfTsNlkL/t+/WTu/rFjzr3f0aNA\n9+5ATIzs+UNkBg4X/SVLlqB169bw8vJCWlpaua9LTk5Gq1atEBQUhJkzZzp6OUtJSUlRHUE37P1e\n2GzAtGnA6NHAQw8BmzY5dt3kZPn6iROBqVM9Mx30Tvi5uIbfC8c5XPTbtm2LFStWoHv37uW+prS0\nFBMnTkRycjIOHjyIRYsWIT093dFLWgY/0Nc48r2w2YAXXgAWLJAN2iZMAE6frtzXnjwJjBkji6+W\nLpWirxf8XFzD74XjHC76rVq1QnBwcIWvSU1NRWBgIFq0aAEfHx8MHz4cq1atcvSSRHbp1Uumc3p7\ny4ybF14Adu68dVvmK1dkOuYf/yj79deuLXPxw8PV5CZyJ7f29HNzcxEQEHD11/7+/sjNzXXnJYlu\nUKeObIi2d69sjvbkk0CDBtK66dMH6NoVaNRI2kF16gD798s5vLVqqU5O5B4VTtmMjIzEiRMnbvn9\n6dOnIzo6+o5vbrOzEWrv681s6tSpqiPohju+Fzt23PjrM2dkdD9tmssv5VL8XFzD74VjKiz6GzZs\ncOrN/fz8kH3dapns7Gz4l3OqNVfjEhG5n0vaO+UV7E6dOiEjIwNZWVkoKipCQkICYmJiXHFJIiJy\ngMNFf8WKFQgICMD27dsxYMAA9OvXDwCQl5eHAQMGAAC8vb0xZ84cPProowgNDcWwYcMQEhLimuRE\nRGQ/TbG1a9dq999/vxYYGKjNmDFDdRxljh8/rkVERGihoaFa69attbffflt1JOVKSkq0sLAwLSoq\nSnUUpc6ePasNGTJEa9WqlRYSEqJt27ZNdSRlpk+froWGhmpt2rTRRowYoRUWFqqO5DGjR4/WGjdu\nrLVp0+bq750+fVrr3bu3FhQUpEVGRmpnz5694/soXZHLefzX+Pj4YPbs2Thw4AC2b9+O9957z7Lf\nizJvv/02QkNDLf+A/7nnnkP//v2Rnp6OvXv3WvZuOSsrCx9++CHS0tKwb98+lJaWYvHixapjeczo\n0aORnJx8w+/NmDEDkZGROHToEHr16oUZM2bc8X2UFn3O47+mSZMmCPvtdO27774bISEhyMvLU5xK\nnZycHCQlJWHcuHGWfsj/888/Y8uWLRgzZgwAaZnWqVNHcSo1ateuDR8fH1y8eBElJSW4ePEi/Pz8\nVMfymPDwcNSrV++G30tMTER8fDwAID4+HitXrrzj+ygt+pzHf3tZWVnYvXs3unTpojqKMpMnT8Y/\n/vEPVKli7e2hjh49ikaNGmH06NHo2LEjxo8fj4sXL6qOpUT9+vXx4osvonnz5mjWrBnq1q2L3r17\nq46lVEFBAXx9fQEAvr6+KCgouOPXKP2Jsvpt++38+uuviIuLw9tvv427775bdRwl1qxZg8aNG6ND\nhw6WHuUDQElJCdLS0vDss88iLS0NNWvWrNQtvBllZmbirbfeQlZWFvLy8vDrr79i4cKFqmPphs1m\nq1RNVVr07ZnHbwXFxcUYMmQInnzySQyy8MbtW7duRWJiIu69916MGDECmzdvxsiRI1XHUsLf3x/+\n/v7o3LkzACAuLq7CDQ7NbOfOnejWrRsaNGgAb29vDB48GFu3blUdSylfX9+rC2jz8/PRuHHjO36N\n0qLPefzXaJqGsWPHIjQ0FM8//7zqOEpNnz4d2dnZOHr0KBYvXoyePXvi008/VR1LiSZNmiAgIACH\nDh0CAGzcuBGtW7dWnEqNVq1aYfv27bh06RI0TcPGjRsRGhqqOpZSMTExmP/bAdHz58+v3GDRXdOL\nKispKUkLDg7W7rvvPm369Omq4yizZcsWzWazae3bt9fCwsK0sLAwbe3atapjKZeSkqJFR0erjqHU\n999/r3Xq1Elr166dFhsbq507d051JGVmzpx5dcrmyJEjtaKiItWRPGb48OFa06ZNNR8fH83f31/7\n+OOPtdOnT2u9evWya8qmbg5GJyIi97P21AgiIoth0ScishAWfSIiC2HRJyKyEBZ9IiILYdEnIrKQ\n/wfKbYvCZLFMmgAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 20 }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.special import gammaln\n", "plt.plot(x, gammaln(x))" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "[]" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD9CAYAAABOd5eOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH9xJREFUeJzt3XtY1VW+x/E3jlSO5jXdplCZqYAiUpppmdsULTuko2Tp\nNDqg3Tul9XjyzDRJdhLNsUZPl5ksk9LxMlMZltIRdZM3MsPUMRMlKbxhhZKGSuDv/LFSc7zBvv32\n5fN6nv0Mbth7fWdXHxff3/qtFWFZloWIiAS1WnYXICIinlOYi4iEAIW5iEgIUJiLiIQAhbmISAhQ\nmIuIhIBqhXlaWhoOh4P4+Pgzvjd16lRq1apFaWmp14sTEZHqqVaYp6amkp2dfcbzxcXFLF26lCuv\nvNLrhYmISPVVK8x79OhBo0aNznj+8ccf5/nnn/d6USIiUjNu98zff/99oqKi6NixozfrERERN9R2\n50Xl5eVMnDiRpUuXnnzuXLsCREREuFeZiEiYq8luK27NzAsLCykqKiIhIYFWrVqxa9currvuOvbv\n33/OgvSwGD9+vO01BMpDn4U+C30W53/UlFsz8/j4eEpKSk7+uVWrVnz22Wc0btzYnbcTEREPVWtm\nPnToULp3705BQQHR0dG8+eabp31frRQREXtVa2Y+d+7c837/q6++8koxoc7pdNpdQsDQZ3GKPotT\n9Fm4L8JypzlTkwEiItzq/4iIhLOaZqdu5xcRCQEKcxGREKAwFxEJAQpzEZEQoDAXEQkBCnMRkRCg\nMBcRCQEKcxGREKAwFxEJAQpzEZEA8+WXNX+NwlxEJIAsWAA9e9b8dQpzEZEAMWMGjBkDvzj3p9rc\n2s9cRES86/nn4dVXITcXrrmm5q9XmIuI2Miy4A9/gPffh1WroGVL995HYS4iYpOqKnj4YcjPh48/\nhssuc/+9FOYiIjaoqIDhw2H/fli2DC691LP3U5iLiPhZeTmkpEBkJCxeDJdc4vl7ajWLiIgflZVB\nv36mpfLOO94JcqhmmKelpeFwOIiPjz/53NixY4mNjSUhIYFBgwZRVlbmnYpEREJUSQn06gWJiTBr\nFtT2Ym+kWmGemppKdnb2ac/17duXLVu2sHHjRtq2bUtGRob3qhIRCTE7d8JNN8GAATBtGtTycl+k\nWm/Xo0cPGjVqdNpzSUlJ1Pq5mq5du7Jr1y7vViYiEiI2bYIePcwNQePHQ0SE98fwyt8NM2fOpH//\n/t54KxGRkLJqFSQlwQsvwEMP+W4cjzs2zz33HBdddBHDhg0758+kp6ef/NrpdOJ0Oj0dVkQk4C1a\nBCNHwpw5JtDPx+Vy4XK53B4rwrIsqzo/WFRURHJyMps3bz753KxZs5gxYwbLli3jknNcko2IiKCa\nQ4iIhIzMTHjyScjKguuvr/nra5qdbs/Ms7OzmTJlCrm5uecMchGRcDR1KkyfDi4XxMT4Z8xqzcyH\nDh1Kbm4u3333HQ6Hg2eeeYaMjAwqKipo3LgxAN26deOVV145cwDNzEUkTFgWjBtn2isffQTR0e6/\nV02zs9ptFncpzEUkHFRWwv33w5Yt8OGH0KSJZ+/ntzaLiIgYR47A0KHmf3NyoF49/9eg2/lFRDxQ\nVga33gp16pj2ih1BDgpzERG37dtnjnjr2NEsP7zoIvtqUZiLiLihsNDcnj9okFm54u3b82tKPXMR\nkRrauBH694ennoIHH7S7GkNhLiJSAytWwF13wUsvwZAhdldzitosIiLVNH++CfL58wMryEEzcxGR\napk2DaZMgaVLISHB7mrOpDAXETmP48fNXZ1ZWbB6NVx5pd0VnZ3CXETkHCoqIC0NvvrKBLmnd3X6\nksJcROQsDh2CwYPNzUA5OfDrX9td0fnpAqiIyL/Ztw+cTmjVyhy6HOhBDgpzEZHTbN8ON95ozur8\n61+9e+iyLwVJmSIivrdunQnxZ5+FUaPsrqZmFOYiIsCSJTBiBLzxBiQn211NzanNIiJhb9YsSE01\nyw+DMchBM3MRCWOWBRMnwuuvQ24utGtnd0XuU5iLSFiqqoJHHzXrx1evhhYt7K7IMwpzEQk75eUw\nbJhZS56bCw0a2F2R56rVM09LS8PhcBAfH3/yudLSUpKSkmjbti19+/bl4MGDPitSRMRb9u+HXr1M\ngC9ZEhpBDtUM89TUVLKzs097btKkSSQlJVFQUEDv3r2ZNGmSTwoUEfGWbdugWzdzzNusWfaeDORt\nEVY1j38uKioiOTmZzZs3AxATE0Nubi4Oh4N9+/bhdDr58ssvzxyghidMi4j4wqpVkJJiLnimpdld\nzYXVNDvd7pmXlJTgcDgAcDgclJSUnPNn09PTT37tdDpxOp3uDisiUmP/+Ac8/DC8/Tb062d3NWfn\ncrlwuVxuv97tmXmjRo04cODAye83btyY0tLSMwfQzFxEbGJZMHWq2Yt80SLo1MnuiqrPbzPzE+2V\n5s2bs3fvXpo1a+buW4mIeF1VFYwebVarrFkD0dF2V+Rbbt8Bescdd5CZmQlAZmYmAwcO9FpRIiKe\n+PFHGDQItm6FlStDP8ihmm2WoUOHkpuby3fffYfD4WDChAkMGDCAIUOG8M0333DVVVexYMECGjZs\neOYAarOIiB+VlJhb8uPi4LXXgnfFSk2zs9o9c3cpzEXEX7Ztg9tug+HDYfx4iIiwuyL3+a1nLiIS\nSE4sPczIMJtmhRuFuYgEvQUL4JFHYPZs6NvX7mrsoTAXkaBlWTBpErz6KixdCgkJdldkH4W5iASl\nigp44AHYuBHy8oJ/10NPKcxFJOiUlsLgwVC/Pnz8MdSta3dF9tNJQyISVHbsgO7d4dpr4d13FeQn\nKMxFJGisWgU33WTu7Jw6FX71K7srChxqs4hIUJgzB8aMCe8VK+ejMBeRgGZZ8MwzkJkJy5dDhw52\nVxSYFOYiErCOHoWRI6Gw0KxY+XnXbTkL9cxFJCB9+y306WOWIK5YoSC/EIW5iAScL780x7v16AHz\n50OdOnZXFPgU5iISUJYvh5494Q9/MPus1FJKVYs+JhEJGG+8AUOHwrx5wXFOZyDRBVARsV1lJYwd\nCx9+aO7obNfO7oqCj8JcRGxVVgZ3320C/ZNPoFEjuysKTmqziIhtduyAG26A1q1h8WIFuScU5iJi\nixUr4MYb4dFH4aWXIDLS7oqCm9osIuJ3f/sbPP00zJ0Lt9xidzWhweOZeUZGBu3btyc+Pp5hw4Zx\n7Ngxb9QlIiGostLMxF980WyapSD3Ho/CvKioiBkzZpCfn8/mzZupqqpi3rx53qpNRELIwYNw++3m\n0OW8PGjTxu6KQotHYV6/fn0iIyMpLy+nsrKS8vJyWrZs6a3aRCREbN9uLnTGxJjlhw0b2l1R6PGo\nZ964cWOeeOIJrrjiCurUqUO/fv3o06fPGT+Xnp5+8mun04nT6fRkWBEJIsuWwbBh8OyzcN99dlcT\nuFwuFy6Xy+3XR1iWZbn74sLCQpKTk1m5ciUNGjTgzjvvJCUlhd/+9renBoiIwIMhRCSIvfIKTJhg\n7ujUHK5mapqdHrVZ1q9fT/fu3WnSpAm1a9dm0KBBrFmzxpO3FJEQUFEBDz5olhyuXq0g9wePwjwm\nJoa8vDyOHDmCZVnk5OQQFxfnrdpEJAiVlEDv3rB7t7nQ2bq13RWFB4/CPCEhgeHDh9O5c2c6duwI\nwH1qiomErc8+g+uvh169YOFCqF/f7orCh0c982oNoJ65SFiYM8cctPy3v8GgQXZXE/xqmp26A1RE\nPFJVBePGwbvvmr3I4+Ptrig8KcxFxG2lpWb/8aoqWLcOmjSxu6LwpY22RMQtW7aY/nj79pCdrSC3\nm8JcRGps4UKz3PDpp+GFF6C2fse3nf4RiEi1HT9ubgKaOdPsP96li90VyQkKcxGplkOHYPhw2L/f\n9MebN7e7IvkltVlE5IK2b4du3aBpU7NiRUEeeBTmInJeWVnmRKBHHjFryC++2O6K5GzUZhGRs6qq\ngvR0mDXLBPoNN9hdkZyPwlxEzlBaaratPXoU1q8Hh8PuiuRC1GYRkdNs2ACdO5v14zk5CvJgoTAX\nkZPeegv69oWMDJg6VevHg4n+UYkIFRUwZgwsXQorVkCHDnZXJDWlMBcJc3v2QEqKWXb46afQoIHd\nFYk71GYRCWMff2zu4rz9dnjvPQV5MNPMXCQMWRZMnw4TJ0JmJtx6q90ViacU5iJh5tAhGDUKCgpg\n7Vq4+mq7KxJvUJtFJIxs3myWHTZoAGvWKMhDicJcJExkZsItt8Af/wivvQZ16thdkXiTx2F+8OBB\nUlJSiI2NJS4ujry8PG/UJSJecuSIaatkZJhlh8OH212R+ILHPfPHHnuM/v37889//pPKykp+/PFH\nb9QlIl6wY4dZdhgTY5YdXnqp3RWJr0RYNTn++d+UlZWRmJjIV199de4BanjCtIh4x7vvwgMPwPjx\n8NBDEBFhd0VSEzXNTo9m5jt37qRp06akpqayceNGrrvuOqZNm8avf/3r034uPT395NdOpxOn0+nJ\nsCJyHj/9BE8+acL8gw/MOZ0S+FwuFy6Xy+3XezQzX79+Pd26dWPNmjV06dKF0aNHU79+fSZMmHBq\nAM3MRfxm1y646y5o1Mjss9K4sd0Vibtqmp0eXQCNiooiKiqKLj8fBJiSkkJ+fr4nbykibvq//zN3\ncyYnm/3HFeThxaMwb968OdHR0RQUFACQk5ND+/btvVKYiFRPZaVZbpiaCnPnwrhxUEuLjsOOR20W\ngI0bNzJq1CgqKipo3bo1b775Jg1+scGD2iwivlNcDEOHQt268Pbb0KyZ3RWJt9Q0Oz0O8wsOoDAX\n8YmsLLj3Xnj8cRg7VrPxUOPX1Swi4n/HjpnVKgsXmp0Ou3e3uyIJBApzkSCyY4dZrXLFFZCfr4uc\ncop+MRMJEvPmQbdu8PvfmzXkCnL5Jc3MRQJceTk89hi4XPDRR3DttXZXJIFIM3ORAPbFF+YOzvJy\n01ZRkMu5KMxFApBlwYwZ0LOnWa0ye7Y2yZLzU5tFJMB8/71ZcvjVV5CbC3FxdlckwUAzc5EAsnw5\ndOoEV10Fn3yiIJfq08xcJABUVMCf/mTaKW++CX372l2RBBuFuYjNCgpg2DC4/HL4/HNo2tTuiiQY\nqc0iYhPLgtdfhxtvhLQ0c3u+glzcpZm5iA1KS+G++8ys3OUCbTYqntLMXMTPXC5zkTMqCtatU5CL\nd2hmLuInx47B00+brWpnzoRbb7W7IgklCnMRP9i4EX73O2jd2nyt3rh4m9osIj5UVQWTJkGfPvDE\nE2aDLAW5+IJm5iI+UlgII0ZAZCSsXw9XXml3RRLKNDMX8TLLgtdeg65dYfBgWLZMQS6+p5m5iBft\n2wejRsGePWZfFa1UEX/xysy8qqqKxMREkpOTvfF2IkHpnXfMksNOnSAvT0Eu/uWVmfm0adOIi4vj\n0KFD3ng7kaBy4IA5PGLtWnMu5w032F2RhCOPZ+a7du1i8eLFjBo1qkYnSYuEgg8+gPh4qF/f7Kui\nIBe7eDwzHzNmDFOmTOGHH34458+kp6ef/NrpdOJ0Oj0dVsRWJ2bjq1ebnQ71r7R4yuVy4XK53H59\nhOXBdPqDDz5gyZIlvPzyy7hcLqZOncqiRYtOHyAigqoqi1paNyMhYtEieOABs1IlIwPq1rW7IglF\nERERNep2eDQzX7NmDVlZWSxevJijR4/yww8/MHz4cN56663Tfq6sDBo18mQkEfuVlp7qjf/97+ZI\nN5FA4dHM/Jdyc3P585//fNaZeWGhxdVXe2MUEXtkZcGDD0JKCkycqNm4+J5fZ+ZnG/xsSktRmEtQ\nKi2FRx81Sw3nzoWbb7a7IpGz81onu2fPnmRlZZ31ewcOeGsUEf+wLPjnP81KlcsuM5tjKcglkPnl\nDtDSUn+MIuIdu3fDww/Dtm2wYIE5CUgk0PlljYlm5hIMjh+HV189dRfn558ryCV4+GVmrjCXQLd1\nK9x7rwl0HeMmwcgvM3O1WSRQVVTAhAmmHz50KKxapSCX4KSZuYSttWvNbLxVK8jPh+houysScZ9f\nwvzbb/0xikj1lJXBU0+ZXQ7/8he48044x6pakaDhlzZLcbE/RhE5P8sya8Xj4uDoUfjXv2DIEAW5\nhAa/zMy/+cYfo4ic27ZtZrnhd9+Z9ePdutldkYh3+WVmfvgwlJf7YySR0x05Ak8/bZYY3n67OYtT\nQS6hyC9hHh2tVov435Il0KGDWXb4+ecwZgzU1kGJEqL88q92dLRptbRr54/RJNzt2gWjR8OGDfDy\ny3DrrXZXJOJ7fpmZX301FBb6YyQJZ8eOwfPPm7s327c3FzgV5BIu/DIzj401v+qK+MqSJWav8Xbt\nzPrxNm3srkjEv/wS5nFx8NFH/hhJws2OHaYXvm2bWTPev7/dFYnYwy9tlthY+OILf4wk4eLwYfjv\n/zYHKN98s2mpKMglnPklzK+4Ag4dgu+/98doEsosC+bMgZgY2LMHNm+GsWPhoovsrkzEXn5ps9Sq\nBZ07w7p1cNtt/hhRQtGnn5qWytGj8I9/aL24yC/5ZWYO5tfhvDx/jSah5Ouv4be/hYEDITXVTAoU\n5CKnU5hLwCorM33xa6+Fa64xFzlHjjS/6YnI6Tz6z6K4uJhevXrRvn17OnTowPTp08/5s127mhnV\n8eOejCjhoLISXnnFLDPctw82bYJnnoF69eyuTCRwedQzj4yM5MUXX6RTp04cPnyY6667jqSkJGJj\nY8/42WbNoEUL+Owz6NLFk1ElVFkWfPihuaDZogVkZ5sbgETkwjwK8+bNm9O8eXMA6tWrR2xsLHv2\n7DlrmIO5G2/JEoW5nOnTT2HcONi7F/78Z7PMUFvTilSf11azFBUVsWHDBrp27XrG99LT0wFzSEV2\ntpOnn3Z6a1gJclu3moMiPvkE/vQn0xPXZlgSjlwuFy6Xy+3XR1iWZXlaxOHDh3E6nTz11FMMHDjw\n9AEiIjgxxLFj4HDAl1/CzxN6CVPffAPp6bBoEfzXf8Ejj0CdOnZXJRI4fpmd1eHxuoCffvqJwYMH\nc88995wR5P/u4othwABYsMDTUSVY7d9vdjRMTDR98e3bTY9cQS7iGY/C3LIsRo4cSVxcHKNHj67W\na4YNM3fwSXgpK4Px483WDlVVsGUL/M//QMOGdlcmEho8CvPVq1cze/ZsVqxYQWJiIomJiWRnZ5/3\nNb17m4MqNm/2ZGQJFmVl8OyzZp34zp3mpJ///V+12US8zSs98/MOcJa+z4QJ5gCB117z5chip7Iy\nmD7dPG67zVzkbNvW7qpEgkdNe+a2hHlJibkhpLAQmjTx5ejib78M8f794Y9/VIiLuMPvF0Dd4XCY\nfTZeesmO0cUXDh481U7ZsQPWrIHMTAW5iL/YMjMH8x/8DTdAQQE0buzLCsSX9u41h0K8/jr8x3+Y\ndopO+RHxXFDMzMHM4AYPhkmT7KpAPLFjB9x/vzlr88gRyM83M3EFuYg9bJuZg5nVdewIubnmaDkJ\nfPn5MHkyLF8ODz4I//mf0LSp3VWJhJ6gmZkDXH65uQvw/vu1m2Igsyxzhmu/fnDHHaY9tnOnWZWk\nIBcJDLbOzMHcQNK9O4wYAQ895MtKpKbKy+Htt2HaNIiMNHduDhtm7uQVEd8KiqWJ/66gAG68EXJy\nICHBl9VIdezaBS+/bC5q3nijCfGePbWLoYg/BVWb5YS2bc3sb8gQs8RN/M+yYPVqGDrUXMc4csSc\nDLVwITidCnKRQBcQM/MTxoyBzz83hxLoV3n/KCuD2bPhr381u1o++CCkpUGDBnZXJhLegrLNckJV\nFdx1l+nPzp4Nv/qVLysLb+vXmwB/5x3o29dchO7VSzNwkUAR1GEO5tf7/v2hZUuYNUsHFXjTgQMw\nfz7MmAGlpXDffea0e216JRJ4gj7MwayiGDQI6tY1M3Ttde2+ykqzrDAz89TywrQ0MxvXKfcigSsk\nwhxM//b3vzd3Gr73HkRFeb+2ULZ5swnwOXPgqqvM0s+77oJGjeyuTESqIyhXs5zNxRfD3/8OKSnQ\ntSt4cDRe2Ni2zWx2FR9vWlUXXWQ+t7Vr4YEHFOQioSxgZ+a/9NFHpjUwZAhMnKi2yy/t2GGO4Vuw\nwBzJdued5nPq1k1tFJFgFjJtln/3/ffm0N9PP4UXXoDk5PBceXH8uPkMFi0yj5IS89vLkCFw000K\ncJFQEbJhfkJ2Njz+uDkMePJkuO46r711wPrxR3N37KJF8OGHpl2SnGwe3bppCadIKPJ7zzw7O5uY\nmBjatGnD5MmTPX27C7r1Vti0CX7zGxgwwKzOcLnMHYyh4qefzOEOEyaY2+gdDnNyT4cOsGoVfPGF\n+YvsppsU5CJieDQzr6qqol27duTk5NCyZUu6dOnC3LlziY2NPTWAl2fmv3TsmFm6+PzzZj36qFHw\nu9/BZZf5ZDifOXrUbC27dq35i+njj+Hqq6FPH3MAdo8eZpmmiIQPv7ZZ1q5dyzPPPEN2djYAk34+\naWLcuHFuF+QOyzIB+MYbkJVltmgdMMBs19qypU+HrrHjx6GoCD77zIT32rXmN4127UzL5Oab4ZZb\ntLWsSLiraXZ6dH/l7t27iY6OPvnnqKgoPvnkE0/e0i0REaYd0bMnHDpkVr8sXGgOE3Y4TED26GGW\nOF59tX9aE5ZlLk7u2AFbtsDGjeaxebPZ9yQx0YR3RgZ06aKZt4h4xqMwj6jmcpL09PSTXzudTpxO\npyfDntell5rVHSkpZq+XTZvMrP2998z5lN9+C7Gxpv/cqhVER5sbklq2hIYNoX59qFfv3CtlKivh\n8GHzKC01gb1vn3ns3Wtm3YWF5lGnDrRubU5RSkgwK07i46FJE5/93xeRIOVyuXB5cEONR22WvLw8\n0tPTT7ZZMjIyqFWrFk8++eSpAfzQZqmJH34wM+V//Qu+/trs3V1cDHv2mB0Ef/jB7A9Tt+6pZX4n\ngv3IEXNxsl4985dGgwbmtKTmzc3D4TB3W7ZubR7aeVBE3OXXnnllZSXt2rVj2bJltGjRguuvv96v\nF0B9pbLSLAe0rFOrZCzLzLQvuSQ817eLiH/5tWdeu3ZtXnrpJfr160dVVRUjR448LciDVe3amlWL\nSHAJupuGRETCQchstCUiItWnMBcRCQEKcxGREKAwFxEJAQpzEZEQoDAXEQkBCnMRkRCgMBcRCQEK\ncxGREKAwFxEJAQpzEZEQoDAXEQkBCnMRkRCgMBcRCQEKcxGREKAwFxEJAQpzEZEQoDAXEQkBCnMR\nkRDgdpiPHTuW2NhYEhISGDRoEGVlZd6sKyS5XC67SwgY+ixO0Wdxij4L97kd5n379mXLli1s3LiR\ntm3bkpGR4c26QpL+RT1Fn8Up+ixO0WfhPrfDPCkpiVq1zMu7du3Krl27vFaUiIjUjFd65jNnzqR/\n//7eeCsREXFDhGVZ1rm+mZSUxL59+854fuLEiSQnJwPw3HPPkZ+fzzvvvHP2ASIivFSqiEh4OU88\nn+G8YX4hs2bNYsaMGSxbtoxLLrnE3bcREREP1Xb3hdnZ2UyZMoXc3FwFuYiIzdyembdp04aKigoa\nN24MQLdu3XjllVe8WpyIiFSP2xdAt2/fztdff82GDRvYsGHDGUGenZ1NTEwMbdq0YfLkyR4XGqyK\ni4vp1asX7du3p0OHDkyfPt3ukmxXVVVFYmLiyesu4ergwYOkpKQQGxtLXFwceXl5dpdkm4yMDNq3\nb098fDzDhg3j2LFjdpfkN2lpaTgcDuLj408+V1paSlJSEm3btqVv374cPHjwgu/jkztAq6qqeOSR\nR8jOzuaLL75g7ty5bN261RdDBbzIyEhefPFFtmzZQl5eHi+//HLYfhYnTJs2jbi4uLC/OP7YY4/R\nv39/tm7dyqZNm4iNjbW7JFsUFRUxY8YM8vPz2bx5M1VVVcybN8/usvwmNTWV7Ozs056bNGkSSUlJ\nFBQU0Lt3byZNmnTB9/FJmK9bt45rrrmGq666isjISO6++27ef/99XwwV8Jo3b06nTp0AqFevHrGx\nsezZs8fmquyza9cuFi9ezKhRo2p0pT7UlJWVsXLlStLS0gCoXbs2DRo0sLkqe9SvX5/IyEjKy8up\nrKykvLycli1b2l2W3/To0YNGjRqd9lxWVhYjRowAYMSIESxcuPCC7+OTMN+9ezfR0dEn/xwVFcXu\n3bt9MVRQKSoqYsOGDXTt2tXuUmwzZswYpkyZcvKGs3C1c+dOmjZtSmpqKtdeey333nsv5eXldpdl\ni8aNG/PEE09wxRVX0KJFCxo2bEifPn3sLstWJSUlOBwOABwOByUlJRd8jU/+iwr3X5/P5vDhw6Sk\npDBt2jTq1atndzm2+OCDD2jWrBmJiYlhPSsHqKysJD8/n4ceeoj8/Hzq1q1brV+lQ1FhYSF/+ctf\nKCoqYs+ePRw+fJg5c+bYXVbAiIiIqFam+iTMW7ZsSXFx8ck/FxcXExUV5YuhgsJPP/3E4MGDueee\nexg4cKDd5dhmzZo1ZGVl0apVK4YOHcry5csZPny43WXZIioqiqioKLp06QJASkoK+fn5Nldlj/Xr\n19O9e3eaNGlC7dq1GTRoEGvWrLG7LFs5HI6TN2zu3buXZs2aXfA1Pgnzzp07s337doqKiqioqGD+\n/Pnccccdvhgq4FmWxciRI4mLi2P06NF2l2OriRMnUlxczM6dO5k3bx633HILb731lt1l2aJ58+ZE\nR0dTUFAAQE5ODu3bt7e5KnvExMSQl5fHkSNHsCyLnJwc4uLi7C7LVnfccQeZmZkAZGZmVm8SaPnI\n4sWLrbZt21qtW7e2Jk6c6KthAt7KlSutiIgIKyEhwerUqZPVqVMna8mSJXaXZTuXy2UlJyfbXYat\nPv/8c6tz585Wx44drd/85jfWwYMH7S7JNpMnT7bi4uKsDh06WMOHD7cqKirsLslv7r77buvyyy+3\nIiMjraioKGvmzJnW999/b/Xu3dtq06aNlZSUZB04cOCC7+PR7fwiIhIYwntJgYhIiFCYi4iEAIW5\niEgIUJiLiIQAhbmISAhQmIuIhID/B1c7ZbX+EHwEAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Some interesting properties of UFuncs\n", "\n", "UFuncs have some methods built-in, which allow for some very interesting, flexible, and fast operations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.arange(5)\n", "y = np.arange(1, 6)\n", "np.add(x, y)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 22, "text": [ "array([1, 3, 5, 7, 9])" ] } ], "prompt_number": 22 }, { "cell_type": "code", "collapsed": false, "input": [ "np.add.accumulate(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "array([ 0, 1, 3, 6, 10])" ] } ], "prompt_number": 23 }, { "cell_type": "code", "collapsed": false, "input": [ "np.multiply.accumulate(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 24, "text": [ "array([0, 0, 0, 0, 0])" ] } ], "prompt_number": 24 }, { "cell_type": "code", "collapsed": false, "input": [ "np.multiply.accumulate(y)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 25, "text": [ "array([ 1, 2, 6, 24, 120])" ] } ], "prompt_number": 25 }, { "cell_type": "code", "collapsed": false, "input": [ "np.add.identity" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "0" ] } ], "prompt_number": 26 }, { "cell_type": "code", "collapsed": false, "input": [ "np.multiply.identity" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 27, "text": [ "1" ] } ], "prompt_number": 27 }, { "cell_type": "code", "collapsed": false, "input": [ "np.add.outer(x, y)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 28, "text": [ "array([[1, 2, 3, 4, 5],\n", " [2, 3, 4, 5, 6],\n", " [3, 4, 5, 6, 7],\n", " [4, 5, 6, 7, 8],\n", " [5, 6, 7, 8, 9]])" ] } ], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [ "# make a times-table\n", "x = np.arange(1, 13)\n", "print np.multiply.outer(x, x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 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]]\n" ] } ], "prompt_number": 29 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Ufunc mini-exercises\n", "\n", "Each of the following functions take an array as input, and return an array as output. They are implemented using loops, which is not very efficient.\n", "\n", "1. For each function, implement a fast version which uses ufuncs to calculate the result more efficiently. Double-check that you get the same result for several different arrays.\n", "\n", "2. use the %timeit magic to time the execution of the two implementations for a large array (say, 1000 elements). " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# 1. computing the element-wise sine + cosine\n", "\n", "from math import sin, cos\n", "def slow_sincos(x):\n", " \"\"\"x is a 1-dimensional array\"\"\"\n", " y = np.zeros_like(x)\n", " for i in range(len(x)):\n", " y[i] = sin(x[i]) + cos(x[i])\n", " return y\n", "\n", "x = np.random.random(5)\n", "print slow_sincos(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 1.41379464 1.24311597 1.39711558 1.1674103 1.41168674]\n" ] } ], "prompt_number": 30 }, { "cell_type": "code", "collapsed": false, "input": [ "# write a fast_sincos function\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 31 }, { "cell_type": "code", "collapsed": false, "input": [ "# 2. computing the difference between adjacent squares\n", "\n", "def slow_sqdiff(x):\n", " \"\"\"x is a 1-dimensional array\"\"\"\n", " y = np.zeros(len(x) - 1)\n", " for i in range(len(y)):\n", " y[i] = x[i + 1] ** 2 - x[i] ** 2\n", " return y\n", "\n", "x = np.random.random(5)\n", "print slow_sqdiff(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[-0.11965629 -0.25739054 0.01840098 0.2762583 ]\n" ] } ], "prompt_number": 32 }, { "cell_type": "code", "collapsed": false, "input": [ "# write a fast_sqdiff function\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 33 }, { "cell_type": "code", "collapsed": false, "input": [ "# 3. computing the outer-product of each consecutive pair\n", "\n", "def slow_pairprod(x):\n", " \"\"\"x is a 1-dimensional array\"\"\"\n", " if len(x) % 2 != 0:\n", " raise ValueError(\"length of x must be even\")\n", " N = len(x) / 2\n", " y = np.zeros((N, N))\n", " for i in range(N):\n", " for j in range(N):\n", " y[i, j] = x[2 * i] * x[2 * j + 1]\n", " return y\n", "\n", "x = np.arange(1, 9)\n", "print slow_pairprod(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 2. 4. 6. 8.]\n", " [ 6. 12. 18. 24.]\n", " [ 10. 20. 30. 40.]\n", " [ 14. 28. 42. 56.]]\n" ] } ], "prompt_number": 34 }, { "cell_type": "code", "collapsed": false, "input": [ "# write a fast_pairprod function\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 35 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategy 2. Using Numpy Aggregates\n", "\n", "Aggregates are functions over arrays which return smaller arrays.\n", "\n", "Numpy has several built-in" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# 10 x 10 array drawn from a standard normal\n", "x = np.random.randn(10, 10)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 36 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.mean()" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.00323045925586\n" ] } ], "prompt_number": 37 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.mean(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.00323045925586\n" ] } ], "prompt_number": 38 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.std()" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "1.0801379167\n" ] } ], "prompt_number": 39 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.var()" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "1.1666979191\n" ] } ], "prompt_number": 40 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.sum()" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.323045925586\n" ] } ], "prompt_number": 41 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.prod()" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "-5.72610473992e-23\n" ] } ], "prompt_number": 42 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.median(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.00389134929792\n" ] } ], "prompt_number": 43 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.percentile(x, 50)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.00389134929792\n" ] } ], "prompt_number": 44 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.percentile(x, (25, 75))" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[-0.78582974895964608, 0.7374780554465159]\n" ] } ], "prompt_number": 45 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Aggregates along certain dimensions" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.rand(3, 5)\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0.20461757 0.63148584 0.87816143 0.60442118 0.79228768]\n", " [ 0.40654761 0.97038926 0.96468488 0.27046867 0.17623957]\n", " [ 0.08734478 0.38353736 0.16706956 0.17744179 0.14886722]]\n" ] } ], "prompt_number": 46 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.sum(0) # sum along rows" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.69850996 1.98541246 2.00991587 1.05233165 1.11739447]\n" ] } ], "prompt_number": 47 }, { "cell_type": "code", "collapsed": false, "input": [ "print x.sum(1) # sum along columns" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 3.11097371 2.78832999 0.96426071]\n" ] } ], "prompt_number": 48 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.median(x, 1)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.63148584 0.40654761 0.16706956]\n" ] } ], "prompt_number": 49 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.mean(x, 1)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.62219474 0.557666 0.19285214]\n" ] } ], "prompt_number": 50 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Binary ufuncs as aggregates\n", "\n", "Any binary ufunc (a ufunc taking two arguments) can be turned into an aggregate using the ``reduce()`` method:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.sum(x, 1)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 51, "text": [ "array([ 3.11097371, 2.78832999, 0.96426071])" ] } ], "prompt_number": 51 }, { "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": 52, "text": [ "array([ 3.11097371, 2.78832999, 0.96426071])" ] } ], "prompt_number": 52 }, { "cell_type": "code", "collapsed": false, "input": [ "np.prod(x, 1)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 53, "text": [ "array([ 0.05433798, 0.01814108, 0.00014784])" ] } ], "prompt_number": 53 }, { "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": 54, "text": [ "array([ 0.05433798, 0.01814108, 0.00014784])" ] } ], "prompt_number": 54 }, { "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": 55, "text": [ "array([ 0.77051726, 9.11086414, 51.60324064])" ] } ], "prompt_number": 55 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "A caution: for ``reduce`` methods, the default axis is 0:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.add.reduce(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 56, "text": [ "array([ 0.69850996, 1.98541246, 2.00991587, 1.05233165, 1.11739447])" ] } ], "prompt_number": 56 }, { "cell_type": "code", "collapsed": false, "input": [ "np.sum(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 57, "text": [ "6.8635644110229004" ] } ], "prompt_number": 57 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### A quick efficiency note:\n", "### Beware the built-in Python aggregates!\n", "\n", "Python has a ``min``, ``max``, and ``sum`` aggregate built-in. These are much more general than the versions in NumPy:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.random(10000)\n", "\n", "%timeit np.sum(x)\n", "%timeit sum(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "100000 loops, best of 3: 16.4 \u00b5s per loop\n", "100 loops, best of 3: 4 ms per loop" ] }, { "output_type": "stream", "stream": "stdout", "text": [ "\n" ] } ], "prompt_number": 58 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "**Dynamic type-checking is slow**.\n", "\n", "Make sure to use Numpy's ``sum``, ``min``, and ``max``." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Aggregate Mini-exercises\n", "\n", "Take the following functions, and convert them into an efficient form using aggregates. Each function expects a 1-dimensional array as input. Double-check that your function returns the same result as the original" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def slow_cubesum(x):\n", " \"\"\"x is a 1D array\"\"\"\n", " result = 0\n", " for i in range(len(x)):\n", " result += x[i] ** 3\n", " return result\n", "\n", "x = np.random.random(100)\n", "print slow_cubesum(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "28.3126613072\n" ] } ], "prompt_number": 59 }, { "cell_type": "code", "collapsed": false, "input": [ "# implement fast_cubesum\n", "\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 60 }, { "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", "print slow_rms(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.298246500479\n" ] } ], "prompt_number": 61 }, { "cell_type": "code", "collapsed": false, "input": [ "# implement fast_rms\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 62 }, { "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", "print slow_sillyfunc(100)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "-33.0\n" ] } ], "prompt_number": 63 }, { "cell_type": "code", "collapsed": false, "input": [ "# Implement fast_sillyfunc using ufuncs & aggragates\n", "\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 64 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategy 3: Using Numpy Broadcasting\n", "\n", "We've taken a look at broadcasting previously. But it's important enough that we'll review it quickly here:\n", "\n", "\n", "\n", "\n", "([image source](http://www.astroml.org/book_figures/appendix/fig_broadcast_visual.html))\n", "\n", "**Broadcasting rules:**\n", "\n", "1. If the two arrays differ in their number of dimensions, the shape of the array with fewer dimensions is *padded* with ones on its leading (left) side.\n", "\n", "2. If the shape of the two arrays does not match in any dimension, the array with shape equal to 1 in that dimension is *stretched* to match the other shape.\n", "\n", "3. If in any dimension the sizes disagree and neither is equal to 1, an error is raised." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Some Broadcasting examples..." ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.arange(10)\n", "print x ** 2" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0 1 4 9 16 25 36 49 64 81]\n" ] } ], "prompt_number": 65 }, { "cell_type": "code", "collapsed": false, "input": [ "Y = x * x[:, np.newaxis]\n", "print Y" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 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]]\n" ] } ], "prompt_number": 66 }, { "cell_type": "code", "collapsed": false, "input": [ "print Y + 10 * x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 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]]\n" ] } ], "prompt_number": 67 }, { "cell_type": "code", "collapsed": false, "input": [ "print Y + 10 * x[:, np.newaxis]" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 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]]\n" ] } ], "prompt_number": 68 }, { "cell_type": "code", "collapsed": false, "input": [ "Y = np.random.random((2, 3, 4))\n", "x = 10 * np.arange(3)\n", "\n", "print Y + x[:, np.newaxis]" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[[ 0.86229625 0.64617272 0.81395431 0.57150281]\n", " [ 10.72496019 10.95945758 10.0420753 10.37791557]\n", " [ 20.70117792 20.47498435 20.92576569 20.26448722]]\n", "\n", " [[ 0.98847073 0.84611754 0.62249531 0.97384757]\n", " [ 10.29122012 10.81186259 10.9161937 10.14179008]\n", " [ 20.65886012 20.99372243 20.88986833 20.12829632]]]\n" ] } ], "prompt_number": 69 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Quick Broadcasting Exercise\n", "\n", "Now, assume you have $N$ points in $D$ dimensions, represented by an array of shape ``[N, D]``.\n", "\n", "1. Compute the mean of the distribution of points efficiently using the built-in ``np.mean`` aggregate (that is, find the ``D``-dimensional point which is the mean of the rest of the points)\n", "2. Compute the mean of the distribution of points efficiently using the ``np.add`` ufunc.\n", "3. Compute the standard error of the mean $\\sigma_{mean} = \\sigma N^{-1/2}$, where $\\sigma$ is the standard-deviation, using the ``np.std`` aggregate.\n", "4. Compute this again using the ``np.add`` ufunc.\n", "5. Construct the matrix ``M``, the centered and normalized version of the ``X`` array: $$ M_{ij} = (X_{ij} - \\mu_j) / \\sigma_j $$ This is one version of *whitening* the array.\n", " " ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = np.random.random((1000, 5)) # 1000 points in 5 dimensions" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 70 }, { "cell_type": "code", "collapsed": false, "input": [ "# 1. Compute the mean of the 1000 points in X\n", "\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 71 }, { "cell_type": "code", "collapsed": false, "input": [ "# 2. Compute the mean using np.add\n", "\n" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 72 }, { "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": 73 }, { "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": 74 }, { "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": 75 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Strategy 4: Fancy Indexing and Masking\n", "\n", "The last strategy we will cover is fancy indexing and masking.\n", "\n", "For example, imagine you have an array of data where negative values indicate some kind of error." ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.array([1, 2, 3, -999, 2, 4, -999])" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 76 }, { "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", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[1 2 3 0 2 4 0]\n" ] } ], "prompt_number": 77 }, { "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", "print mask" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[False False False True False False True]\n" ] } ], "prompt_number": 78 }, { "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", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[1 2 3 0 2 4 0]\n" ] } ], "prompt_number": 79 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Typically this is done directly:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.array([1, 2, 3, -999, 2, 4, -999])\n", "x[x < 0] = 0\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[1 2 3 0 2 4 0]\n" ] } ], "prompt_number": 80 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Useful masking functions" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.random(5)\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.97053187 0.72813759 0.95240005 0.04301341 0.77007516]\n" ] } ], "prompt_number": 81 }, { "cell_type": "code", "collapsed": false, "input": [ "x[x > 0.5] = np.nan\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ nan nan nan 0.04301341 nan]\n" ] } ], "prompt_number": 82 }, { "cell_type": "code", "collapsed": false, "input": [ "x[np.isnan(x)] = np.inf\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ inf inf inf 0.04301341 inf]\n" ] } ], "prompt_number": 83 }, { "cell_type": "code", "collapsed": false, "input": [ "np.nan == np.nan" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 84, "text": [ "False" ] } ], "prompt_number": 84 }, { "cell_type": "code", "collapsed": false, "input": [ "x[np.isinf(x)] = 0\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0. 0. 0. 0.04301341 0. ]\n" ] } ], "prompt_number": 85 }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.array([1, 0, -np.inf, np.inf, np.nan])\n", "print \"input \", x\n", "print \"x < 0 \", (x < 0)\n", "print \"x > 0 \", (x > 0)\n", "print \"isinf \", np.isinf(x)\n", "print \"isnan \", np.isnan(x)\n", "print \"isposinf\", np.isposinf(x)\n", "print \"isneginf\", np.isneginf(x)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "input [ 1. 0. -inf inf nan]\n", "x < 0 [False False True False False]\n", "x > 0 [ True False False True False]\n", "isinf [False False True True False]\n", "isnan [False False False False True]\n", "isposinf [False False False True False]\n", "isneginf [False False True False False]\n" ] } ], "prompt_number": 86 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Boolean Operations on Masks\n", "\n", "Use bitwise operators (and make sure to use parentheses!)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.arange(16).reshape((4, 4))\n", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0 1 2 3]\n", " [ 4 5 6 7]\n", " [ 8 9 10 11]\n", " [12 13 14 15]]\n" ] } ], "prompt_number": 87 }, { "cell_type": "code", "collapsed": false, "input": [ "print (x < 5)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ True True True True]\n", " [ True False False False]\n", " [False False False False]\n", " [False False False False]]\n" ] } ], "prompt_number": 88 }, { "cell_type": "code", "collapsed": false, "input": [ "print ~(x < 5)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[False False False False]\n", " [False True True True]\n", " [ True True True True]\n", " [ True True True True]]\n" ] } ], "prompt_number": 89 }, { "cell_type": "code", "collapsed": false, "input": [ "print (x < 10) & (x % 2 == 0)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ True False True False]\n", " [ True False True False]\n", " [ True False False False]\n", " [False False False False]]\n" ] } ], "prompt_number": 90 }, { "cell_type": "code", "collapsed": false, "input": [ "print (x > 3) & (x < 8)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[False False False False]\n", " [ True True True True]\n", " [False False False False]\n", " [False False False False]]\n" ] } ], "prompt_number": 91 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Counting elements with a mask\n", "\n", "Sum over a mask to find the number of ``True`` elements:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.random.random(100)\n", "print \"array is length\", len(x), \"and has\"\n", "print (x > 0.5).sum(), \"elements are greater than 0.5\"" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "array is length 100 and has\n", "50 elements are greater than 0.5\n" ] } ], "prompt_number": 92 }, { "cell_type": "code", "collapsed": false, "input": [ "# clip is a useful function:\n", "x = np.clip(x, 0.3, 0.6)\n", "\n", "print np.sum(x < 0.3)\n", "print np.sum(x > 0.6)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0\n", "0\n" ] } ], "prompt_number": 93 }, { "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": 94 }, { "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", "print x" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0.82518495 0.49805524 0.00840461]\n", " [ 0.85186368 0.98583616 0.77922605]\n", " [ 0.52795826 0.20786051 0.49627606]]\n" ] } ], "prompt_number": 95 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.where(x < 0.3)" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(array([0, 2]), array([2, 1]))\n" ] } ], "prompt_number": 96 }, { "cell_type": "code", "collapsed": false, "input": [ "print x[x < 0.3]" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.00840461 0.20786051]\n" ] } ], "prompt_number": 97 }, { "cell_type": "code", "collapsed": false, "input": [ "print x[np.where(x < 0.3)]" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.00840461 0.20786051]\n" ] } ], "prompt_number": 98 }, { "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", "print X" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0 1 2 3]\n", " [ 4 5 6 7]\n", " [ 8 9 10 11]\n", " [12 13 14 15]]\n" ] } ], "prompt_number": 99 }, { "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": 100, "text": [ "array([1, 4])" ] } ], "prompt_number": 100 }, { "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": 101, "text": [ "array([ 0, 5, 10, 15])" ] } ], "prompt_number": 101 }, { "cell_type": "code", "collapsed": false, "input": [ "X.diagonal()" ], "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() = 100" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "ename": "SyntaxError", "evalue": "can't assign to function call (, line 1)", "output_type": "pyerr", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m1\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": 103 }, { "cell_type": "code", "collapsed": false, "input": [ "X[range(4), range(4)] = 100" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 104 }, { "cell_type": "code", "collapsed": false, "input": [ "print X" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[100 1 2 3]\n", " [ 4 100 6 7]\n", " [ 8 9 100 11]\n", " [ 12 13 14 100]]\n" ] } ], "prompt_number": 105 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Randomizing the rows" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = np.arange(24).reshape((6, 4))\n", "print X" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 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]]\n" ] } ], "prompt_number": 106 }, { "cell_type": "code", "collapsed": false, "input": [ "i = np.arange(6)\n", "np.random.shuffle(i)\n", "print i" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[0 5 3 1 4 2]\n" ] } ], "prompt_number": 107 }, { "cell_type": "code", "collapsed": false, "input": [ "print X[i] # X[i, :] is identical" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0 1 2 3]\n", " [20 21 22 23]\n", " [12 13 14 15]\n", " [ 4 5 6 7]\n", " [16 17 18 19]\n", " [ 8 9 10 11]]\n" ] } ], "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", "print X[i2]" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[[ 0 1 2 3]\n", " [20 21 22 23]]\n", "\n", " [[12 13 14 15]\n", " [ 4 5 6 7]]\n", "\n", " [[16 17 18 19]\n", " [ 8 9 10 11]]]\n" ] } ], "prompt_number": 109 }, { "cell_type": "code", "collapsed": false, "input": [ "print X[i2].shape" ], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(3, 2, 4)\n" ] } ], "prompt_number": 110 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary: Speeding up NumPy\n", "\n", "It's all about **moving loops into compiled code:**\n", "\n", "1. Use Numpy **ufuncs** to your advantage (eliminate loops!)\n", "\n", "2. Use Numpy **aggregates** to your advantage (eliminate loops!)\n", "\n", "3. Use Numpy **broadcasting** to your advantage (eliminate loops!)\n", "\n", "4. Use Numpy **slicing and masking** to your advantage (eliminate loops!)\n", "\n", "5. Use a tool like *SWIG*, *cython* or *f2py* to interface to compiled code." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Homework: asteroid data\n", "\n", "In the [github repository](https://github.com/jakevdp/2013_fall_ASTR599/), there is a file containing measurements of 5000 asteroid orbits, at ``notebooks/data/asteroids5000.csv``.\n", "\n", "These are compiled from a query at [http://ssd.jpl.nasa.gov/sbdb_query.cgi](http://ssd.jpl.nasa.gov/sbdb_query.cgi)\n", "\n", "### Part 1: loading and exploring the data\n", "\n", "1. Use ``np.genfromtxt`` to load the data from the file. This is like ``loadtxt``, but can handle missing data.\n", "\n", " - Remember to set the appropriate ``delimiter`` keyword.\n", "\n", "2. ``genfromtxt`` sets all missing values to ``np.nan``. Use the operations we discussed here to answer these questions:\n", "\n", " - How many values are missing in this data?\n", " - How many complete *rows* are there? i.e. how many objects have **no** missing values?\n", "\n", "3. Create a new array containing only the rows with no missing values.\n", "\n", "4. Compute the maximum, minimum, mean, and standard deviation of the values in each column." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Part 2: Plotting the data\n", "\n", "1. Use the bash ``head`` command to display the first line of the data file: this lists the names of the columns in the dataset. (remember that bash commands in the notebook are indicated by ``!``, and that ``head -n`` displays the first ``n`` lines of a file)\n", "\n", "2. Invoke the ``matplotlib inline`` magic to make figures appear inline in the notebook\n", "\n", "3. Use ``plt.scatter`` to plot the semi-major axis versus the sine of the inclination angle (note that the inclination angle is listed in degrees -- you'll have to convert it to radians to compute the sine). What do you notice about the distribution? What do you think could explain this?\n", "\n", "4. Use ``plt.scatter`` to plot a color-magnitude diagram of the asteroids (H vs B-V). You should see two distinct \"families\" of asteroids indicated in this plot. Over-plot a line that divides these.\n", "\n", "5. Repeat the orbital parameter plot from above, but plot the two \"families\" from #4 in different colors. Note that this magnitude is undefined for many of the asteroids. Do you see any correlation between color and orbit?\n", "\n", "6. Compare what you found to plots in [Parker et al. 2008](http://adsabs.harvard.edu/abs/2008Icar..198..138P). Note that we're not using the same data here, but we're looking at a similar collection of objects." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [] } ], "metadata": {} } ] }