notebooks/05_NumpyIntro.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:71fb8d99b5d571a5b56a7bfa094d1651546f3596bf6933228c3deb8afa60578e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "skip"
      }
     },
     "source": [
      "<style>\n",
      "div.cell, div.text_cell_render{\n",
      "  max-width:750px;\n",
      "  margin-left:auto;\n",
      "  margin-right:auto;\n",
      "}\n",
      "\n",
      ".rendered_html\n",
      "{\n",
      "  font-size: 140%;\n",
      "  }\n",
      "\n",
      ".rendered_html li\n",
      "{\n",
      "  line-height: 1.8;\n",
      "  }\n",
      "\n",
      ".rendered_html h1, h2 {\n",
      "  text-align:center;\n",
      "  font-familly:\"Charis SIL\", serif;\n",
      "}\n",
      "\n",
      ".input_prompt, .CodeMirror-lines, .output_area\n",
      "{\n",
      "  font-family: Consolas;\n",
      "  font-size: 110%;\n",
      "  }\n",
      "</style>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<small><i>This notebook was put together by [Jake Vanderplas](http://www.vanderplas.com) for UW's [Astro 599](http://www.astro.washington.edu/users/vanderplas/Astr599/) course. Source and license info is on [GitHub](https://github.com/jakevdp/2013_fall_ASTR599/).</i></small>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# An Introduction to Numerical Computing with Python\n",
      "While the Python language is an excellent tool for general-purpose programming, with a highly readable syntax, rich and powerful data types (strings, lists, sets, dictionaries, arbitrary length integers, etc) and a very comprehensive standard library, it was not designed specifically for mathematical and scientific computing.  Neither the language nor its standard library have facilities for the efficient representation of multidimensional datasets, tools for linear algebra and general matrix manipulations (an essential building block of virtually all technical computing), nor any data visualization facilities.\n",
      "\n",
      "In particular, Python lists are very flexible containers that can be nested arbitrarily deep and which can hold any Python object in them, but they are poorly suited to represent efficiently common mathematical constructs like vectors and matrices.  In contrast, much of our modern heritage of scientific computing has been built on top of libraries written in the Fortran language, which has native support for vectors and matrices as well as a library of mathematical functions that can efficiently operate on entire arrays at once."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Basics of Numpy arrays\n",
      "\n",
      "We now turn our attention to the Numpy library, which forms the base layer for the entire 'scipy ecosystem'.  Once you have installed numpy, you can import it as"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### The Numpy array structure\n",
      "<center>\n",
      "<img src=\"files/images/dtype.png\" width=90%>\n",
      "</center>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "As mentioned above, the main object provided by numpy is a powerful array.  We'll start by exploring how the numpy array differs from Python lists.  We start by creating a simple list and an array with the same contents of the list:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "lst = [10, 20, 30, 40]\n",
      "arr = np.array([10, 20, 30, 40])"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Element indexing\n",
      "\n",
      "Elements of a one-dimensional array are accessed with the same syntax as a list:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lst[0]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "10"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr[0]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "10"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr[-1]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "40"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr[2:]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "array([30, 40])"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Differences between arrays and lists\n",
      "\n",
      "The first difference to note between lists and arrays is that arrays are *homogeneous*; i.e. all elements of an array must be of the same type.  In contrast, lists can contain elements of arbitrary type. For example, we can change the last element in our list above to be a string:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lst[-1] = 'a string inside a list'\n",
      "lst"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[10, 20, 30, 'a string inside a list']"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "but the same can not be done with an array, as we get an error message:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr[-1] = 'a string inside an array'"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "invalid literal for long() with base 10: 'a string inside an array'",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-8-29c0bfa5fa8a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0marr\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'a string inside an array'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;31mValueError\u001b[0m: invalid literal for long() with base 10: 'a string inside an array'"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Array Attributes\n",
      "\n",
      "The information about the type of an array is contained in its *dtype* attribute:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr.dtype"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "dtype('int64')"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Once an array has been created, its dtype is fixed and it can only store elements of the same type.  For this example where the dtype is integer, if we store a floating point number it will be automatically converted into an integer:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr[-1] = 1.234\n",
      "arr"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "array([10, 20, 30,  1])"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Creating Arrays\n",
      "\n",
      "Above we created an array from an existing list; now let us now see other ways in which we can create arrays, which we'll illustrate next.  A common need is to have an array initialized with a constant value, and very often this value is 0 or 1 (suitable as starting value for additive and multiplicative loops respectively); `zeros` creates arrays of all zeros, with any desired dtype:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.zeros(5, dtype=float)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "array([ 0.,  0.,  0.,  0.,  0.])"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.zeros(3, dtype=int)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([0, 0, 0])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.zeros(3, dtype=complex)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "array([ 0.+0.j,  0.+0.j,  0.+0.j])"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "and similarly for `ones`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print '5 ones:', np.ones(5)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "5 ones: [ 1.  1.  1.  1.  1.]\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "If we want an array initialized with an arbitrary value, we can create an empty array and then use the fill method to put the value we want into the array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = np.empty(4)\n",
      "a.fill(5.5)\n",
      "a"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "array([ 5.5,  5.5,  5.5,  5.5])"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Filling arrays with sequences\n",
      "\n",
      "Numpy also offers the `arange` function, which works like the builtin `range` but returns an array instead of a list:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.arange(5)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "array([0, 1, 2, 3, 4])"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "and the `linspace` and `logspace` functions to create linearly and logarithmically-spaced grids respectively, with a fixed number of points and including both ends of the specified interval:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"A linear grid between 0 and 1:\"\n",
      "print np.linspace(0, 1, 4)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "A linear grid between 0 and 1:\n",
        "[ 0.          0.33333333  0.66666667  1.        ]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"A logarithmic grid between 10**1 and 10**3:\"\n",
      "print np.logspace(1, 3, 4)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "A logarithmic grid between 10**1 and 10**3:\n",
        "[   10.            46.41588834   215.443469    1000.        ]\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Creating random arrays\n",
      "\n",
      "Finally, it is often useful to create arrays with random numbers that follow a specific distribution.  The `np.random` module contains a number of functions that can be used to this effect, for example this will produce an array of 5 random samples taken from a standard normal distribution (0 mean and variance 1):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.randn(5)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "array([ 0.87809648,  0.39991   ,  1.31517435, -0.3666163 ,  0.61513923])"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "whereas this will also give 5 samples, but from a normal distribution with a mean of 10 and a variance of 3:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "norm10 = np.random.normal(10, 3, 5)\n",
      "norm10"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "array([  7.1999085 ,   4.46144008,   9.62660515,  15.59534389,  10.6923198 ])"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Indexing with other arrays\n",
      "Above we saw how to index arrays with single numbers and slices, just like Python lists.  But arrays allow for a more sophisticated kind of indexing which is very powerful: you can index an array with another array, and in particular with an array of boolean values.  This is particluarly useful to extract information from an array that matches a certain condition.\n",
      "\n",
      "Consider for example that in the array `norm10` we want to replace all values above 9 with the value 0.  We can do so by first finding the *mask* that indicates where this condition is true or false:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mask = norm10 > 9\n",
      "mask"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "array([False, False,  True,  True,  True], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Now that we have this mask, we can use it to either read those values or to reset them to 0:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Values above 9:', norm10[mask]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Values above 9: [  9.62660515  15.59534389  10.6923198 ]\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Resetting all values above 9 to 0...'\n",
      "norm10[mask] = 0\n",
      "print norm10"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Resetting all values above 9 to 0...\n",
        "[ 7.1999085   4.46144008  0.          0.          0.        ]\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Arrays with more than one dimension\n",
      "Up until now all our examples have used one-dimensional arrays.  But Numpy can create arrays of aribtrary dimensions, and all the methods illustrated in the previous section work with more than one dimension.  For example, a list of lists can be used to initialize a two dimensional array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lst2 = [[1, 2], [3, 4]]\n",
      "arr2 = np.array([[1, 2], [3, 4]])\n",
      "arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "array([[1, 2],\n",
        "       [3, 4]])"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "With two-dimensional arrays we start seeing the power of numpy: while a nested list can be indexed using repeatedly the `[ ]` operator, multidimensional arrays support a much more natural indexing syntax with a single `[ ]` and a set of indices separated by commas:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print lst2[0][1]\n",
      "print arr2[0,1]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2\n",
        "2\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Arrays with more than one dimension\n",
      "\n",
      "Most of the array creation functions listed above can be used with more than one dimension, for example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.zeros((2,3))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "array([[ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.]])"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.normal(10, 3, (2, 4))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "array([[ 11.21418405,   9.56604053,  11.01177765,  13.46153439],\n",
        "       [ 10.57163512,  13.23242428,   9.99134628,  10.87468603]])"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "In fact, the shape of an array can be changed at any time, as long as the total number of elements is unchanged.  For example, if we want a 2x4 array with numbers increasing from 0, the easiest way to create it is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr = np.arange(8).reshape(2, 4)\n",
      "print arr"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[0 1 2 3]\n",
        " [4 5 6 7]]\n"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Views, not Copies\n",
      "\n",
      "Note that reshaping (like most numpy operations) provides a **view** of *the same memory*:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr = np.arange(8)\n",
      "arr2 = arr.reshape(2, 4)\n",
      "\n",
      "arr[0] = 1000\n",
      "print arr\n",
      "print arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[1000    1    2    3    4    5    6    7]\n",
        "[[1000    1    2    3]\n",
        " [   4    5    6    7]]\n"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "This lack of copying allows for very efficient vectorized operations."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Slices\n",
      "\n",
      "With multidimensional arrays, you can also use slices, and you can mix and match slices and single indices in the different dimensions (using the same array as above):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Slicing in the second row:', arr2[1, 2:4]\n",
      "print 'All rows, third column   :', arr2[:, 2]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Slicing in the second row: [6 7]\n",
        "All rows, third column   : [2 6]\n"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "If you only provide one index, then you will get an array with one less dimension containing that row:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'First row:  ', arr2[0]\n",
      "print 'Second row: ', arr2[1]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "First row:   [1000    1    2    3]\n",
        "Second row:  [4 5 6 7]\n"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Array Properties and Methods\n",
      "\n",
      "Now that we have seen how to create arrays with more than one dimension, it's a good idea to look at some of the most useful properties and methods that arrays have.  The following provide basic information about the size, shape and data in the array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Data type                :', arr.dtype\n",
      "print 'Total number of elements :', arr.size\n",
      "print 'Number of dimensions     :', arr.ndim\n",
      "print 'Shape (dimensionality)   :', arr.shape\n",
      "print 'Memory used (in bytes)   :', arr.nbytes"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Data type                : int64\n",
        "Total number of elements : 8\n",
        "Number of dimensions     : 1\n",
        "Shape (dimensionality)   : (8,)\n",
        "Memory used (in bytes)   : 64\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Arrays also have many useful methods, some especially useful ones are:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Minimum and maximum             :', arr.min(), arr.max()\n",
      "print 'Sum and product of all elements :', arr.sum(), arr.prod()\n",
      "print 'Mean and standard deviation     :', arr.mean(), arr.std()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Minimum and maximum             : 1 1000\n",
        "Sum and product of all elements : 1028 5040000\n",
        "Mean and standard deviation     : 128.5 329.401350938\n"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "For these methods, the above operations area all computed on all the elements of the array.  But for a multidimensional array, it's possible to do the computation along a single dimension, by passing the `axis` parameter; for example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'For the following array:\\n', arr2\n",
      "print 'The sum of elements along the rows is    :', arr2.sum(axis=1)\n",
      "print 'The sum of elements along the columns is :', arr2.sum(axis=0)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "For the following array:\n",
        "[[1000    1    2    3]\n",
        " [   4    5    6    7]]\n",
        "The sum of elements along the rows is    : [1006   22]\n",
        "The sum of elements along the columns is : [1004    6    8   10]\n"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "As you can see in this example, the value of the `axis` parameter is the dimension which will be *consumed* once the operation has been carried out.  This is why to sum along the rows we use `axis=0`.  \n",
      "\n",
      "This can be easily illustrated with an example that has more dimensions; we create an array with 4 dimensions and shape `(3,4,5,6)` and sum along the axis number 2 (i.e. the *third* axis, since in Python all counts are 0-based).  That consumes the dimension whose length was 5, leaving us with a new array that has shape `(3,4,6)`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.zeros((3,4,5,6)).sum(2).shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 35,
       "text": [
        "(3, 4, 6)"
       ]
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Another widely used property of arrays is the `.T` attribute, which allows you to access the transpose of the array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Array:\\n', arr2\n",
      "print 'Transpose:\\n', arr2.T"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Array:\n",
        "[[1000    1    2    3]\n",
        " [   4    5    6    7]]\n",
        "Transpose:\n",
        "[[1000    4]\n",
        " [   1    5]\n",
        " [   2    6]\n",
        " [   3    7]]\n"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## More array properties\n",
      "\n",
      "We don't have time here to look at all the methods and properties of arrays, here's a complete list.  Simply try exploring some of these IPython to learn more, or read their description in the full Numpy documentation:\n",
      "\n",
      "    arr.T             arr.copy          arr.getfield      arr.put           arr.squeeze\n",
      "    arr.all           arr.ctypes        arr.imag          arr.ravel         arr.std\n",
      "    arr.any           arr.cumprod       arr.item          arr.real          arr.strides\n",
      "    arr.argmax        arr.cumsum        arr.itemset       arr.repeat        arr.sum\n",
      "    arr.argmin        arr.data          arr.itemsize      arr.reshape       arr.swapaxes\n",
      "    arr.argsort       arr.diagonal      arr.max           arr.resize        arr.take\n",
      "    arr.astype        arr.dot           arr.mean          arr.round         arr.tofile\n",
      "    arr.base          arr.dtype         arr.min           arr.searchsorted  arr.tolist\n",
      "    arr.byteswap      arr.dump          arr.nbytes        arr.setasflat     arr.tostring\n",
      "    arr.choose        arr.dumps         arr.ndim          arr.setfield      arr.trace\n",
      "    arr.clip          arr.fill          arr.newbyteorder  arr.setflags      arr.transpose\n",
      "    arr.compress      arr.flags         arr.nonzero       arr.shape         arr.var\n",
      "    arr.conj          arr.flat          arr.prod          arr.size          arr.view\n",
      "    arr.conjugate     arr.flatten       arr.ptp           arr.sort          "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Operating with arrays\n",
      "Arrays support all regular arithmetic operators, and the numpy library also contains a complete collection of basic mathematical functions that operate on arrays.  It is important to remember that in general, all operations with arrays are applied *element-wise*, i.e., are applied to all the elements of the array at the same time.  Consider for example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr1 = np.arange(4)\n",
      "arr2 = np.arange(10, 14)\n",
      "print arr1, '+', arr2, '=', arr1+arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[0 1 2 3] + [10 11 12 13] = [10 12 14 16]\n"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Importantly, you must remember that even the multiplication operator is by default applied element-wise, it is *not* the matrix multiplication from linear algebra (as is the case in Matlab, for example):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print arr1, '*', arr2, '=', arr1*arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[0 1 2 3] * [10 11 12 13] = [ 0 11 24 39]\n"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "We may also multiply an array by a scalar:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "1.5 * arr1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 40,
       "text": [
        "array([ 0. ,  1.5,  3. ,  4.5])"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "This is a first example of **broadcasting**"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Broadcasting\n",
      "\n",
      "While this means that in principle arrays must always match in their dimensionality in order for an operation to be valid, numpy will *broadcast* dimensions when possible.  Here is an example of broadcasting a scalar to a 1D array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.arange(3)\n",
      "print np.arange(3) + 5"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[0 1 2]\n",
        "[5 6 7]\n"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "We can also broadcast a 1D array to a 2D array, in this case adding a vector to all rows of a matrix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.ones((3, 3)) + np.arange(3)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 42,
       "text": [
        "array([[ 1.,  2.,  3.],\n",
        "       [ 1.,  2.,  3.],\n",
        "       [ 1.,  2.,  3.]])"
       ]
      }
     ],
     "prompt_number": 42
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "We can also broadcast in two directions at a time:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.arange(3).reshape((3, 1)) + np.arange(3)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 43,
       "text": [
        "array([[0, 1, 2],\n",
        "       [1, 2, 3],\n",
        "       [2, 3, 4]])"
       ]
      }
     ],
     "prompt_number": 43
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Rules of Broadcasting\n",
      "\n",
      "Broadcasting rules can do the following:\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.\n",
      "\n",
      "Note that all of this happens without ever actually creating the stretched arrays in memory! This broadcasting behavior is in practice enormously powerful, especially because when numpy broadcasts to create new dimensions or to 'stretch' existing ones, it doesn't actually replicate the data.  In the example above the operation is carried *as if* the 1.5 was a 1-d array with 1.5 in all of its entries, but no actual array was ever created.  This can save lots of memory in cases when the arrays in question are large and can have significant performance implications."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Broadcasting Examples:\n",
      "\n",
      "So when we do\n",
      "\n",
      "    np.arange(3) + 5\n",
      "    \n",
      "The scalar 5 is\n",
      "\n",
      "* first 'promoted' to a 1-dimensional array of length 1\n",
      "* then, this array is 'stretched' to length 3 to match the first array.\n",
      "\n",
      "After these two operations are complete, the addition can proceed as now both operands are one-dimensional arrays of length 3."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "When we do\n",
      "\n",
      "    np.ones((3, 3)) + np.arange(3)\n",
      "    \n",
      "The second array is\n",
      "\n",
      "- first 'promoted' to a 2-dimensional array of shape (1, 3)\n",
      "- then 'stretched' to length 3 to match the first array"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "When we do\n",
      "\n",
      "    np.arange(3).reshape((3, 1)) + np.arange(3)\n",
      "    \n",
      "- the second array is 'promoted' to a 2-dimensional array of shape (1, 3)\n",
      "- the second array is 'stretched' to shape (3, 3)\n",
      "- the first array is 'stretched' to shape (3, 3)\n",
      "\n",
      "Then the operation proceeds as if on two 3 $\\times$ 3 arrays"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "The general rule is: when operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing dimensions, and works its way forward, creating dimensions of length 1 as needed. Two dimensions are considered compatible when\n",
      "\n",
      "* they are equal to begin with, or\n",
      "* one of them is 1; in this case numpy will do the 'stretching' to make them equal.\n",
      "\n",
      "If these conditions are not met, a `ValueError: frames are not aligned` exception is thrown, indicating that the arrays have incompatible shapes. The size of the resulting array is the maximum size along each dimension of the input arrays."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Visualizing Broadcasting\n",
      "\n",
      "<img src=\"http://www.astroml.org/_images/fig_broadcast_visual_1.png\">\n",
      "\n",
      "([image source](http://www.astroml.org/book_figures/appendix/fig_broadcast_visual.html))"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Questions:\n",
      "\n",
      "Will the following broadcasting operations work?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr1 = np.ones((2, 3))\n",
      "arr2 = np.ones((2, 1))\n",
      "\n",
      "arr1 + arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 44,
       "text": [
        "array([[ 2.,  2.,  2.],\n",
        "       [ 2.,  2.,  2.]])"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr1 = np.ones((2, 3))\n",
      "arr2 = np.ones(2)\n",
      "\n",
      "arr1 + arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "operands could not be broadcast together with shapes (2,3) (2,) ",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-45-133a817ed55f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0marr2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0marr1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0marr2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (2,3) (2,) "
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr1 = np.ones((2, 3))\n",
      "arr2 = np.ones((2, 1))\n",
      "\n",
      "arr1 + arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "array([[ 2.,  2.,  2.],\n",
        "       [ 2.,  2.,  2.]])"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Quick Exercise:\n",
      "\n",
      "Use ``np.arange`` and ``reshape``\n",
      "\n",
      "    A = [[1 2 3 4]\n",
      "         [5 6 7 8]]\n",
      "\n",
      "Use ``np.arange`` to create the array\n",
      "\n",
      "    B = [1 2]\n",
      "    \n",
      "Use broadcasting to add ``B`` to each **column** of ``A`` to create the final array\n",
      "\n",
      "    A + B = [[2  3  4  5]\n",
      "             [7  8  9 10]\n",
      "\n",
      "Hint: what shape does ``B`` have to be changed to?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Answers:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.arange(1, 9).reshape((2, 4))\n",
      "B = np.arange(1, 3)\n",
      "A + B.reshape((2, 1))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 47,
       "text": [
        "array([[ 2,  3,  4,  5],\n",
        "       [ 7,  8,  9, 10]])"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Another way to change the shape of ``B`` is to use the ``newaxis`` keyword:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print B.shape\n",
      "print B[:, np.newaxis].shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(2,)\n",
        "(2, 1)\n"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A + B[:, np.newaxis]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "array([[ 2,  3,  4,  5],\n",
        "       [ 7,  8,  9, 10]])"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Element-wise Functions\n",
      "\n",
      "As we mentioned before, Numpy ships with a full complement of mathematical functions that work on entire arrays, including logarithms, exponentials, trigonometric and hyperbolic trigonometric functions, etc.  Furthermore, scipy ships a rich special function library in the `scipy.special` module that includes Bessel, Airy, Fresnel, Laguerre and other classical special functions.  For example, sampling the sine function at 100 points between $0$ and $2\\pi$ is as simple as:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "x = np.linspace(0, 2*np.pi, 100)\n",
      "y = np.sin(x)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Linear algebra in numpy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Numpy ships with a basic linear algebra library, and all arrays have a `dot` method whose behavior is that of the scalar dot product when its arguments are vectors (one-dimensional arrays) and the traditional matrix multiplication when one or both of its arguments are two-dimensional arrays:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "v1 = np.array([2, 3, 4])\n",
      "v2 = np.array([1, 0, 1])\n",
      "\n",
      "print v1, '.', v2, '=', np.dot(v1, v2)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[2 3 4] . [1 0 1] = 6\n"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Here is a regular matrix-vector multiplication, note that the array `v1` should be viewed as a *column* vector in traditional linear algebra notation; numpy makes no distinction between row and column vectors and simply verifies that the dimensions match the required rules of matrix multiplication, in this case we have a $2 \\times 3$ matrix multiplied by a 3-vector, which produces a 2-vector:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.arange(6).reshape(2, 3)\n",
      "print A, 'x', v1, '=', np.dot(A, v1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[0 1 2]\n",
        " [3 4 5]] x [2 3 4] = [11 38]\n"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "For matrix-matrix multiplication, the same dimension-matching rules must be satisfied, e.g. consider the difference between $A \\times A^T$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.dot(A, A.T)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 5 14]\n",
        " [14 50]]\n"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "and $A^T \\times A$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.dot(A.T, A)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 9 12 15]\n",
        " [12 17 22]\n",
        " [15 22 29]]\n"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Furthermore, the `numpy.linalg` module includes additional functionality such as determinants, matrix norms, Cholesky, eigenvalue and singular value decompositions, etc.  For even more linear algebra tools, `scipy.linalg` contains the majority of the tools in the classic LAPACK libraries as well as functions to operate on sparse matrices.  We refer the reader to the Numpy and Scipy documentations for additional details on these."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Reading and writing arrays to disk\n",
      "Numpy lets you read and write arrays into files in a number of ways.  In order to use these tools well, it is critical to understand the difference between a *text* and a *binary* file containing numerical data.  In a text file, the number $\\pi$ could be written as \"3.141592653589793\", for example: a string of digits that a human can read, with in this case 15 decimal digits.  In contrast, that same number written to a binary file would be encoded as 8 characters (bytes) that are not readable by a human but which contain the exact same data that the variable `pi` had in the computer's memory.  \n",
      "\n",
      "The tradeoffs between the two modes are thus:\n",
      "\n",
      "* Text mode: occupies more space, precision can be lost (if not all digits are written to disk), but is readable and editable by hand with a text editor.  Can *only* be used for one- and two-dimensional arrays.\n",
      "\n",
      "* Binary mode: compact and exact representation of the data in memory, can't be read or edited by hand.  Arrays of any size and dimensionality can be saved and read without loss of information."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "#### Text data\n",
      "\n",
      "First, let's see how to read and write arrays in text mode.  The `np.savetxt` function saves an array to a text file, with options to control the precision, separators and even adding a header:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr = np.arange(10).reshape(2, 5)\n",
      "np.savetxt('test.out', arr, fmt='%.2e', header=\"My dataset\")\n",
      "!cat test.out"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "# My dataset\r\n",
        "0.00e+00 1.00e+00 2.00e+00 3.00e+00 4.00e+00\r\n",
        "5.00e+00 6.00e+00 7.00e+00 8.00e+00 9.00e+00\r\n"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "And this same type of file can then be read with the matching `np.loadtxt` function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr2 = np.loadtxt('test.out')\n",
      "print arr2"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.  1.  2.  3.  4.]\n",
        " [ 5.  6.  7.  8.  9.]]\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "#### Binary Data\n",
      "\n",
      "For binary data, Numpy provides the `np.save` and `np.savez` routines.  The first saves a single array to a file with `.npy` extension, while the latter can be used to save a *group* of arrays into a single file with `.npz` extension.  The files created with these routines can then be read with the `np.load` function.\n",
      "\n",
      "Let us first see how to use the simpler `np.save` function to save a single array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.save('test.npy', arr2)\n",
      "# Now we read this back\n",
      "arr2n = np.load('test.npy')\n",
      "# Let's see if any element is non-zero in the difference.\n",
      "# A value of True would be a problem.\n",
      "print 'Any differences?', np.any(arr2-arr2n)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Any differences? False\n"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Now let us see how the `np.savez` function works.  You give it a filename and either a sequence of arrays or a set of keywords.  In the first mode, the function will auotmatically name the saved arrays in the archive as `arr_0`, `arr_1`, etc:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.savez('test.npz', arr, arr2)\n",
      "arrays = np.load('test.npz')\n",
      "arrays.files"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "['arr_1', 'arr_0']"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### .npz: multiple binary outputs in one file\n",
      "\n",
      "Alternatively, we can explicitly choose how to name the arrays we save:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.savez('test.npz', array1=arr, array2=arr2)\n",
      "arrays = np.load('test.npz')\n",
      "arrays.files"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 59,
       "text": [
        "['array2', 'array1']"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "The object returned by `np.load` from an `.npz` file works like a dictionary, though you can also access its constituent files by attribute using its special `.f` field; this is best illustrated with an example with the `arrays` object from above:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'First row of first array:', arrays['array1'][0]\n",
      "# This is an equivalent way to get the same field\n",
      "print 'First row of first array:', arrays.f.array1[0]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "First row of first array: [0 1 2 3 4]\n",
        "First row of first array: [0 1 2 3 4]\n"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "This `.npz` format is a very convenient way to package compactly and without loss of information, into a single file, a group of related arrays that pertain to a specific problem.  At some point, however, the complexity of your dataset may be such that the optimal approach is to use one of the standard formats in scientific data processing that have been designed to handle complex datasets, such as NetCDF or HDF5.  \n",
      "\n",
      "Fortunately, there are tools for manipulating these formats in Python, and for storing data in other ways such as databases.  A complete discussion of the possibilities is beyond the scope of this discussion, but of particular interest for scientific users we at least mention the following:\n",
      "\n",
      "**Matlab files**\n",
      "\n",
      "The `scipy.io` module contains routines to read and write Matlab files in `.mat` format and files in the NetCDF format that is widely used in certain scientific disciplines.\n",
      "\n",
      "** HDF5 files**\n",
      "\n",
      "For manipulating files in the HDF5 format, there are two excellent options in Python: The PyTables project offers a high-level, object oriented approach to manipulating HDF5 datasets, while the h5py project offers a more direct mapping to the standard HDF5 library interface.  Both are excellent tools; if you need to work with HDF5 datasets you should read some of their documentation and examples and decide which approach is a better match for your needs."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Breakout: trapezoidal integration\n",
      "\n",
      "**Illustrates**: basic array slicing, functions as first class objects.\n",
      "\n",
      "In this exercise, you are tasked with implementing the simple trapezoid rule\n",
      "formula for numerical integration. If we want to compute the definite integral\n",
      "\n",
      "$$\n",
      "     \\int_{a}^{b}f(x)dx\n",
      "$$\n",
      "\n",
      "we can partition the integration interval $[a,b]$ into smaller subintervals,\n",
      "and approximate the area under the curve for each subinterval by the area of\n",
      "the trapezoid created by linearly interpolating between the two function values\n",
      "at each end of the subinterval:\n",
      "\n",
      "<img src=\"http://upload.wikimedia.org/wikipedia/commons/thumb/d/dd/Trapezoidal_rule_illustration.png/316px-Trapezoidal_rule_illustration.png\" /img>\n",
      "\n",
      "The blue line represents the function $f(x)$ and the red line\n",
      "is the linear interpolation.  By subdividing the interval $[a,b]$, the area under $f(x)$ can thus be approximated as the sum of the areas of all\n",
      "the resulting trapezoids. \n",
      "\n",
      "If we denote by $x_{i}$ ($i=0,\\ldots,n,$ with $x_{0}=a$ and\n",
      "$x_{n}=b$) the abscissas where the function is sampled, then\n",
      "\n",
      "$$\n",
      "   \\int_{a}^{b}f(x)dx\\approx\\frac{1}{2}\\sum_{i=1}^{n}\\left(x_{i}-x_{i-1}\\right)\\left(f(x_{i})+f(x_{i-1})\\right).\n",
      "$$\n",
      "\n",
      "The common case of using equally spaced abscissas with spacing $h=(b-a)/n$ reads simply\n",
      "\n",
      "$$\n",
      "   \\int_{a}^{b}f(x)dx\\approx\\frac{h}{2}\\sum_{i=1}^{n}\\left(f(x_{i})+f(x_{i-1})\\right).\n",
      "$$\n",
      "\n",
      "One frequently receives the function values already precomputed, $y_{i}=f(x_{i}),$\n",
      "so the equation above becomes\n",
      "\n",
      "$$\n",
      "   \\int_{a}^{b}f(x)dx\\approx\\frac{1}{2}\\sum_{i=1}^{n}\\left(x_{i}-x_{i-1}\\right)\\left(y_{i}+y_{i-1}\\right).\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Exercises\n",
      "\n",
      "#### Part 1\n",
      "\n",
      "Write a function `trapz(x, y)`, that applies the trapezoid formula to pre-computed values, where `x` and `y` are 1-d arrays.  (Use numpy operations, not loops!)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Part 2 \n",
      "\n",
      "Write a function  `trapzf(f, a, b, npts=100)` that accepts a function `f`, the endpoints `a`\n",
      "and `b` and the number of samples to take `npts`.  Sample the function uniformly at these\n",
      "points and return the value of the integral."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Part 3\n",
      "\n",
      "Verify that both functions above are correct by showing that they produce correct values for a simple integral such as $\\int_0^3 x^2$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "#### Part 4\n",
      "\n",
      "Preview of what's to come: use the documentation features of IPython to explore the submodule ``scipy.integrage``.  Can you find a suitable function to perform this integration above?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 60
    }
   ],
   "metadata": {}
  }
 ]
}