Lecture-numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy -  multidimensional data arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `numpy` package (module) is used in almost all numerical computation using Python. It is a package that provide high-performance vector, matrix and higher-dimensional data structures for Python. It is implemented in C and Fortran so when calculations are vectorized (formulated with vectors and matrices), performance is very good. \n",
    "\n",
    "To use `numpy` you need to import the module, using for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the `numpy` package the terminology used for vectors, matrices and higher-dimensional data sets is *array*. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating `numpy` arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a number of ways to initialize new numpy arrays, for example from\n",
    "\n",
    "* a Python list or tuples\n",
    "* using functions that are dedicated to generating numpy arrays, such as `arange`, `linspace`, etc.\n",
    "* reading data from files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### From lists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, to create new vector and matrix arrays from Python lists we can use the `numpy.array` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 3, 4])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a vector: the argument to the array function is a Python list\n",
    "v = np.array([1,2,3,4])\n",
    "\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2],\n",
       "       [3, 4]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a matrix: the argument to the array function is a nested Python list\n",
    "M = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `v` and `M` objects are both of the type `ndarray` that the `numpy` module provides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(numpy.ndarray, numpy.ndarray)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(v), type(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference between the `v` and `M` arrays is only their shapes. We can get information about the shape of an array by using the `ndarray.shape` property."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4,)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of elements in the array is available through the `ndarray.size` property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equivalently, we could use the function `numpy.shape` and `numpy.size`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.size(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far the `numpy.ndarray` looks awefully much like a Python list (or nested list). Why not simply use Python lists for computations instead of creating a new array type? \n",
    "\n",
    "There are several reasons:\n",
    "\n",
    "* Python lists are very general. They can contain any kind of object. They are dynamically typed. They do not support mathematical functions such as matrix and dot multiplications, etc. Implementing such functions for Python lists would not be very efficient because of the dynamic typing.\n",
    "* Numpy arrays are **statically typed** and **homogeneous**. The type of the elements is determined when the array is created.\n",
    "* Numpy arrays are memory efficient.\n",
    "* Because of the static typing, fast implementation of mathematical functions such as multiplication and addition of `numpy` arrays can be implemented in a compiled language (C and Fortran is used).\n",
    "\n",
    "Using the `dtype` (data type) property of an `ndarray`, we can see what type the data of an array has:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int64')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get an error if we try to assign a value of the wrong type to an element in a numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "invalid literal for long() with base 10: 'hello'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-a09d72434238>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mM\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"hello\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: invalid literal for long() with base 10: 'hello'"
     ]
    }
   ],
   "source": [
    "M[0,0] = \"hello\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want, we can explicitly define the type of the array data when we create it, using the `dtype` keyword argument: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.+0.j,  2.+0.j],\n",
       "       [ 3.+0.j,  4.+0.j]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = np.array([[1, 2], [3, 4]], dtype=complex)\n",
    "\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Common data types that can be used with `dtype` are: `int`, `float`, `complex`, `bool`, `object`, etc.\n",
    "\n",
    "We can also explicitly define the bit size of the data types, for example: `int64`, `int16`, `float128`, `complex128`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using array-generating functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For larger arrays it is inpractical to initialize the data manually, using explicit python lists. Instead we can use one of the many functions in `numpy` that generate arrays of different forms. Some of the more common are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### arange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a range\n",
    "\n",
    "x = np.arange(0, 10, 1) # arguments: start, stop, step\n",
    "\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -1.00000000e+00,  -9.00000000e-01,  -8.00000000e-01,\n",
       "        -7.00000000e-01,  -6.00000000e-01,  -5.00000000e-01,\n",
       "        -4.00000000e-01,  -3.00000000e-01,  -2.00000000e-01,\n",
       "        -1.00000000e-01,  -2.22044605e-16,   1.00000000e-01,\n",
       "         2.00000000e-01,   3.00000000e-01,   4.00000000e-01,\n",
       "         5.00000000e-01,   6.00000000e-01,   7.00000000e-01,\n",
       "         8.00000000e-01,   9.00000000e-01])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(-1, 1, 0.1)\n",
    "\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### linspace and logspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0.        ,   0.41666667,   0.83333333,   1.25      ,\n",
       "         1.66666667,   2.08333333,   2.5       ,   2.91666667,\n",
       "         3.33333333,   3.75      ,   4.16666667,   4.58333333,\n",
       "         5.        ,   5.41666667,   5.83333333,   6.25      ,\n",
       "         6.66666667,   7.08333333,   7.5       ,   7.91666667,\n",
       "         8.33333333,   8.75      ,   9.16666667,   9.58333333,  10.        ])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using linspace, both end points ARE included\n",
    "np.linspace(0, 10, 25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.00000000e+00,   3.03773178e+00,   9.22781435e+00,\n",
       "         2.80316249e+01,   8.51525577e+01,   2.58670631e+02,\n",
       "         7.85771994e+02,   2.38696456e+03,   7.25095809e+03,\n",
       "         2.20264658e+04])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.logspace(0, 10, 10, base=np.e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### mgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, y = np.mgrid[0:5, 0:5] # similar to meshgrid in MATLAB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 0, 0],\n",
       "       [1, 1, 1, 1, 1],\n",
       "       [2, 2, 2, 2, 2],\n",
       "       [3, 3, 3, 3, 3],\n",
       "       [4, 4, 4, 4, 4]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4],\n",
       "       [0, 1, 2, 3, 4],\n",
       "       [0, 1, 2, 3, 4],\n",
       "       [0, 1, 2, 3, 4],\n",
       "       [0, 1, 2, 3, 4]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### random data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.49406293,  0.21425745,  0.61331193,  0.32784842,  0.99591338],\n",
       "       [ 0.24089539,  0.44805359,  0.09345322,  0.56656295,  0.44791492],\n",
       "       [ 0.13505864,  0.88682905,  0.00324599,  0.74313305,  0.50196652],\n",
       "       [ 0.09585222,  0.26228338,  0.02996983,  0.32233701,  0.30819306],\n",
       "       [ 0.77678822,  0.77188318,  0.81466163,  0.01803531,  0.96661312]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# uniform random numbers in [0,1]\n",
    "random.rand(5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.24588781, -0.5974331 , -0.40903861,  1.2327517 ,  0.51265757],\n",
       "       [-0.9178861 ,  0.75621978, -0.11977446,  0.67184885,  1.21928622],\n",
       "       [ 0.29441652, -0.43413505,  0.97153403,  0.14295138, -0.30659839],\n",
       "       [-0.37641243, -1.200776  , -0.97160933,  0.56522347, -1.31024933],\n",
       "       [-0.05640364,  0.28085847,  0.94445629,  0.53184511, -1.74751539]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# standard normal distributed random numbers\n",
    "random.randn(5,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### diag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 0, 0],\n",
       "       [0, 2, 0],\n",
       "       [0, 0, 3]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a diagonal matrix\n",
    "np.diag([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 0, 0],\n",
       "       [0, 0, 2, 0],\n",
       "       [0, 0, 0, 3],\n",
       "       [0, 0, 0, 0]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# diagonal with offset from the main diagonal\n",
    "np.diag([1,2,3], k=1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### zeros and ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.]])"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.zeros((3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  1.,  1.],\n",
       "       [ 1.,  1.,  1.],\n",
       "       [ 1.,  1.,  1.]])"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ones((3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## File I/O"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comma-separated values (CSV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A very common file format for data files is comma-separated values (CSV), or related formats such as TSV (tab-separated values). To read data from such files into Numpy arrays we can use the `numpy.genfromtxt` function. For example, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "echo '\n",
    "1800  1  1    -6.1    -6.1    -6.1 1\n",
    "1800  1  2   -15.4   -15.4   -15.4 1\n",
    "1800  1  3   -15.0   -15.0   -15.0 1\n",
    "1800  1  4   -19.3   -19.3   -19.3 1\n",
    "1800  1  5   -16.8   -16.8   -16.8 1\n",
    "1800  1  6   -11.4   -11.4   -11.4 1\n",
    "1800  1  7    -7.6    -7.6    -7.6 1\n",
    "1800  1  8    -7.1    -7.1    -7.1 1\n",
    "1800  1  9   -10.1   -10.1   -10.1 1\n",
    "1800  1 10    -9.5    -9.5    -9.5 1\n",
    "' > test.dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = np.genfromtxt('test.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 7)"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1800.  1800.  1800.  1800.  1800.  1800.  1800.  1800.  1800.  1800.]\n",
      "[  1.80000000e+03   1.00000000e+00   1.00000000e+00  -6.10000000e+00\n",
      "  -6.10000000e+00  -6.10000000e+00   1.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "print data[:,0]\n",
    "print data[0,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `numpy.savetxt` we can store a Numpy array to a file in CSV format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.56592074,  0.31925451,  0.07071436],\n",
       "       [ 0.40445531,  0.33905924,  0.06890646],\n",
       "       [ 0.43876792,  0.11784734,  0.80457052]])"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = random.rand(3,3)\n",
    "\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.savetxt(\"random-matrix.csv\", M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.659207388768656699e-01 3.192545105089277246e-01 7.071436311328971325e-02\r\n",
      "4.044553055396516994e-01 3.390592443226384534e-01 6.890646403468136949e-02\r\n",
      "4.387679159929198658e-01 1.178473431220541379e-01 8.045705234897031488e-01\r\n"
     ]
    }
   ],
   "source": [
    "!cat random-matrix.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.56592 0.31925 0.07071\r\n",
      "0.40446 0.33906 0.06891\r\n",
      "0.43877 0.11785 0.80457\r\n"
     ]
    }
   ],
   "source": [
    "np.savetxt(\"random-matrix.csv\", M, fmt='%.5f') # fmt specifies the format\n",
    "\n",
    "!cat random-matrix.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.56592  0.31925  0.07071]\n",
      " [ 0.40446  0.33906  0.06891]\n",
      " [ 0.43877  0.11785  0.80457]]\n"
     ]
    }
   ],
   "source": [
    "A = np.loadtxt(\"random-matrix.csv\")\n",
    "print A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "loadtxt can read directly from gzipped (.gz) and bzip2 (.bz2) files.\n",
    "\n",
    "More complicate files can be read using more options. For instance, to read a file\n",
    "skipping all comment lines (starting with a \\#) and the first 5 rows, reading only in\n",
    "columns 0,1,2 using a comma as delimiter we give the command:\n",
    "\n",
    "\n",
    "array = np.loadtxt(filename, dtype=dtype, comments=‘#’, delimiter=‘,’, skiprows=5, usecols=(0, 1, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numpy's native file format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Useful when storing and reading back numpy array data. Use the functions `numpy.save` and `numpy.load`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "random-matrix.npy: data\r\n"
     ]
    }
   ],
   "source": [
    "np.save(\"random-matrix.npy\", M)\n",
    "!file random-matrix.npy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.56592074,  0.31925451,  0.07071436],\n",
       "       [ 0.40445531,  0.33905924,  0.06890646],\n",
       "       [ 0.43876792,  0.11784734,  0.80457052]])"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.load(\"random-matrix.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More properties of the numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.itemsize # bytes per element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "72"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.nbytes # number of bytes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.ndim # number of dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manipulating arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can index elements in an array using square brackets and indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# v is a vector, and has only one dimension, taking one index\n",
    "v[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33905924432263845"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# M is a matrix, or a 2 dimensional array, taking two indices \n",
    "M[1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we omit an index of a multidimensional array it returns the whole row (or, in general, a N-1 dimensional array) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.56592074,  0.31925451,  0.07071436],\n",
       "       [ 0.40445531,  0.33905924,  0.06890646],\n",
       "       [ 0.43876792,  0.11784734,  0.80457052]])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.40445531,  0.33905924,  0.06890646])"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same thing can be achieved with using `:` instead of an index: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.40445531,  0.33905924,  0.06890646])"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M[1,:] # row 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.31925451,  0.33905924,  0.11784734])"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M[:,1] # column 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can assign new values to elements in an array using indexing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M[0,0] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.31925451,  0.07071436],\n",
       "       [ 0.40445531,  0.33905924,  0.06890646],\n",
       "       [ 0.43876792,  0.11784734,  0.80457052]])"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# also works for rows and columns\n",
    "M[1,:] = 0\n",
    "M[:,2] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.31925451, -1.        ],\n",
       "       [ 0.        ,  0.        , -1.        ],\n",
       "       [ 0.43876792,  0.11784734, -1.        ]])"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Index slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Index slicing is the technical name for the syntax `M[lower:upper:step]` to extract part of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 3, 4, 5])"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([1,2,3,4,5])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 3])"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Array slices are *mutable*: if they are assigned a new value the original array from which the slice was extracted is modified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1, -2, -3,  4,  5])"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1:3] = [-2,-3]\n",
    "\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can omit any of the three parameters in `M[lower:upper:step]`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1, -2, -3,  4,  5])"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[::] # lower, upper, step all take the default values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1, -3,  5])"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[::2] # step is 2, lower and upper defaults to the beginning and end of the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1, -2, -3])"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[:3] # first three elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4, 5])"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[3:] # elements from index 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Negative indices counts from the end of the array (positive index from the begining):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = np.array([1,2,3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[-1] # the last element in the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 4, 5])"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[-3:] # the last three elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Index slicing works exactly the same way for multidimensional arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4],\n",
       "       [10, 11, 12, 13, 14],\n",
       "       [20, 21, 22, 23, 24],\n",
       "       [30, 31, 32, 33, 34],\n",
       "       [40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[n+m*10 for n in range(5)] for m in range(5)])\n",
    "\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 12, 13],\n",
       "       [21, 22, 23],\n",
       "       [31, 32, 33]])"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a block from the original array\n",
    "A[1:4, 1:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  2,  4],\n",
       "       [20, 22, 24],\n",
       "       [40, 42, 44]])"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# strides\n",
    "A[::2, ::2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fancy indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fancy indexing is the name for when an array or list is used in-place of an index: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 11, 12, 13, 14],\n",
       "       [20, 21, 22, 23, 24],\n",
       "       [30, 31, 32, 33, 34]])"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "row_indices = [1, 2, 3]\n",
    "A[row_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11, 22, 34])"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col_indices = [1, 2, -1] # remember, index -1 means the last element\n",
    "A[row_indices, col_indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use index masks: If the index mask is an Numpy array of data type `bool`, then an element is selected (True) or not (False) depending on the value of the index mask at the position of each element: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4])"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = np.array([n for n in range(5)])\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 2])"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "row_mask = np.array([True, False, True, False, False])\n",
    "B[row_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 2])"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# same thing\n",
    "row_mask = np.array([1,0,1,0,0], dtype=bool)\n",
    "B[row_mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This feature is very useful to conditionally select elements from an array, using for example comparison operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ,\n",
       "        5.5,  6. ,  6.5,  7. ,  7.5,  8. ,  8.5,  9. ,  9.5])"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(0, 10, 0.5)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False, False, False, False, False, False, False, False, False,\n",
       "       False, False,  True,  True,  True,  True, False, False, False,\n",
       "       False, False], dtype=bool)"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mask = (5 < x) * (x < 7.5)\n",
    "\n",
    "mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.5,  6. ,  6.5,  7. ])"
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions for extracting data from arrays and creating arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index mask can be converted to position index using the `where` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([11, 12, 13, 14]),)"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indices = np.where(mask)\n",
    "\n",
    "indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.5,  6. ,  6.5,  7. ])"
      ]
     },
     "execution_count": 140,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[indices] # this indexing is equivalent to the fancy indexing x[mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### diag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the diag function we can also extract the diagonal and subdiagonals of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0, 11, 22, 33, 44])"
      ]
     },
     "execution_count": 141,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diag(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 21, 32, 43])"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diag(A, -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### take"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `take` function is similar to fancy indexing described above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3, -2, -1,  0,  1,  2])"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v2 = np.arange(-3,3)\n",
    "v2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2,  0,  2])"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "row_indices = [1, 3, 5]\n",
    "v2[row_indices] # fancy indexing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2,  0,  2])"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v2.take(row_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But `take` also works on lists and other objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2,  0,  2])"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.take([-3, -2, -1,  0,  1,  2], row_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### choose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constructs an array by picking elements from several arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5, -2,  5, -2])"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "which = [1, 0, 1, 0]\n",
    "choices = [[-2,-2,-2,-2], [5,5,5,5]]\n",
    "\n",
    "np.choose(which, choices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vectorizing code is the key to writing efficient numerical calculation with Python/Numpy. That means that as much as possible of a program should be formulated in terms of matrix and vector operations, like matrix-matrix multiplication."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scalar-array operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the usual arithmetic operators to multiply, add, subtract, and divide arrays with scalar numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v1 = np.arange(0, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 2, 4, 6, 8])"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1 * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 3, 4, 5, 6])"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1 + 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A*2  [[ 0  2  4  6  8]\n",
      " [20 22 24 26 28]\n",
      " [40 42 44 46 48]\n",
      " [60 62 64 66 68]\n",
      " [80 82 84 86 88]]\n",
      "A+2  [[ 2  3  4  5  6]\n",
      " [12 13 14 15 16]\n",
      " [22 23 24 25 26]\n",
      " [32 33 34 35 36]\n",
      " [42 43 44 45 46]]\n"
     ]
    }
   ],
   "source": [
    "print \"A*2 \",A * 2\n",
    "print \"A+2 \",A + 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Element-wise array-array operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we add, subtract, multiply and divide arrays with each other, the default behaviour is **element-wise** operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   0,    1,    4,    9,   16],\n",
       "       [ 100,  121,  144,  169,  196],\n",
       "       [ 400,  441,  484,  529,  576],\n",
       "       [ 900,  961, 1024, 1089, 1156],\n",
       "       [1600, 1681, 1764, 1849, 1936]])"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A * A # element-wise multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  4,  9, 16])"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1 * v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we multiply arrays with compatible shapes, we get an element-wise multiplication of each row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5, 5), (5,))"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.shape, v1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0,   1,   4,   9,  16],\n",
       "       [  0,  11,  24,  39,  56],\n",
       "       [  0,  21,  44,  69,  96],\n",
       "       [  0,  31,  64,  99, 136],\n",
       "       [  0,  41,  84, 129, 176]])"
      ]
     },
     "execution_count": 162,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A * v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about matrix mutiplication? There are two ways. We can either use the `dot` function, which applies a matrix-matrix, matrix-vector, or inner vector multiplication to its two arguments: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 300,  310,  320,  330,  340],\n",
       "       [1300, 1360, 1420, 1480, 1540],\n",
       "       [2300, 2410, 2520, 2630, 2740],\n",
       "       [3300, 3460, 3620, 3780, 3940],\n",
       "       [4300, 4510, 4720, 4930, 5140]])"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(A, A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 30, 130, 230, 330, 430])"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(A, v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(v1, v1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can cast the array objects to the type `matrix`. This changes the behavior of the standard arithmetic operators `+, -, *` to use matrix algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M = np.matrix(A)\n",
    "v = np.matrix(v1).T # make it a column vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0],\n",
       "        [1],\n",
       "        [2],\n",
       "        [3],\n",
       "        [4]])"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 300,  310,  320,  330,  340],\n",
       "        [1300, 1360, 1420, 1480, 1540],\n",
       "        [2300, 2410, 2520, 2630, 2740],\n",
       "        [3300, 3460, 3620, 3780, 3940],\n",
       "        [4300, 4510, 4720, 4930, 5140]])"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M * M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 30],\n",
       "        [130],\n",
       "        [230],\n",
       "        [330],\n",
       "        [430]])"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M * v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[30]])"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# inner product\n",
    "v.T * v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 30],\n",
       "        [131],\n",
       "        [232],\n",
       "        [333],\n",
       "        [434]])"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# with matrix objects, standard matrix algebra applies\n",
    "v + M*v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we try to add, subtract or multiply objects with incomplatible shapes we get an error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v = np.matrix([1,2,3,4,5,6]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5, 5), (6, 1))"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.shape, v.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "shapes (5,5) and (6,1) not aligned: 5 (dim 1) != 6 (dim 0)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-174-995fb48ad0cc>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mM\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/home/dario/anaconda/lib/python2.7/site-packages/numpy/matrixlib/defmatrix.pyc\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    341\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    342\u001b[0m             \u001b[0;31m# This promotes 1-D vectors to row vectors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 343\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0masmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    344\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misscalar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'__rmul__'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    345\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: shapes (5,5) and (6,1) not aligned: 5 (dim 1) != 6 (dim 0)"
     ]
    }
   ],
   "source": [
    "M * v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See also the related functions: `inner`, `outer`, `cross`, `kron`, `tensordot`. Try for example `help(kron)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array/Matrix transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we have used the `.T` to transpose the matrix object `v`. We could also have used the `transpose` function to accomplish the same thing. \n",
    "\n",
    "Other mathematical functions that transform matrix objects are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.+1.j,  0.+2.j],\n",
       "        [ 0.+3.j,  0.+4.j]])"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = np.matrix([[1j, 2j], [3j, 4j]])\n",
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.-1.j,  0.-2.j],\n",
       "        [ 0.-3.j,  0.-4.j]])"
      ]
     },
     "execution_count": 176,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.conjugate(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hermitian conjugate: transpose + conjugate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.-1.j,  0.-3.j],\n",
       "        [ 0.-2.j,  0.-4.j]])"
      ]
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can extract the real and imaginary parts of complex-valued arrays using `real` and `imag`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.,  0.],\n",
       "        [ 0.,  0.]])"
      ]
     },
     "execution_count": 179,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.real(C) # same as: C.real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.,  2.],\n",
       "        [ 3.,  4.]])"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.imag(C) # same as: C.imag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or the complex argument and absolute value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.78539816,  1.10714872],\n",
       "       [ 1.24904577,  1.32581766]])"
      ]
     },
     "execution_count": 181,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.angle(C+1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.,  2.],\n",
       "        [ 3.,  4.]])"
      ]
     },
     "execution_count": 182,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix computations\n",
    "\n",
    "Let's explore the linear algebra (linalg) methods:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.+2.j ,  0.-1.j ],\n",
       "        [ 0.-1.5j,  0.+0.5j]])"
      ]
     },
     "execution_count": 183,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(C) # equivalent to C.I "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[  1.00000000e+00+0.j,   0.00000000e+00+0.j],\n",
       "        [  1.11022302e-16+0.j,   1.00000000e+00+0.j]])"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.I * C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determinant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.0000000000000004+0j)"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.50000000000000011+0j)"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(C.I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often it is useful to store datasets in Numpy arrays. Numpy provides a number of functions to calculate statistics of datasets in arrays. \n",
    "\n",
    "For example, let's calculate some properties from the dataset defined before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 7)"
      ]
     },
     "execution_count": 190,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-11.83"
      ]
     },
     "execution_count": 192,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "np.mean(data[:,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### standard deviations and variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.3011742582694783, 18.5001)"
      ]
     },
     "execution_count": 193,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(data[:,3]), np.var(data[:,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### min and max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-19.300000000000001"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[:,3].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.0999999999999996"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[:,3].max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### sum, prod, and trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 195,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = np.arange(0, 10)\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "45"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# sum up all elements\n",
    "np.sum(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3628800"
      ]
     },
     "execution_count": 197,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# product of all elements\n",
    "np.prod(d+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  3,  6, 10, 15, 21, 28, 36, 45])"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# cummulative sum\n",
    "np.cumsum(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([      1,       2,       6,      24,     120,     720,    5040,\n",
       "         40320,  362880, 3628800])"
      ]
     },
     "execution_count": 199,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# cummulative product\n",
    "np.cumprod(d+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "110"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# same as: diag(A).sum()\n",
    "np.trace(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If an array as *NaN* elements, it is possible to use the same operations with the nan in their name:\n",
    "for instance, nanmax, nansum, etc. These operations mask all the NaN's in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nan\n",
      "33.3\n"
     ]
    }
   ],
   "source": [
    "a = np.array([0,2,33.3,np.NaN,2.])\n",
    "print np.max(a)\n",
    "print np.nanmax(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another important feature, is finding the index of the maximum of an array.\n",
    "This is done with the function argmax (or nanargmax if the array contains NaN's)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n"
     ]
    }
   ],
   "source": [
    "print np.argmax(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computations on subsets of arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute with subsets of the data in an array using indexing, fancy indexing, and the other methods of extracting data from an array (described above).\n",
    "\n",
    "For example, let's consider a temperature dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1800  1  1    -6.1    -6.1    -6.1 1\r\n",
      "1800  1  2   -15.4   -15.4   -15.4 1\r\n",
      "1800  1  3   -15.0   -15.0   -15.0 1\r\n"
     ]
    }
   ],
   "source": [
    "!head -n 3 stockholm_td_adj.dat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataformat is: year, month, day, daily average temperature, low, high, location.\n",
    "\n",
    "If we are interested in the average temperature only in a particular month, say February, then we can create a index mask and use it to select only the data for that month using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,\n",
       "        12.])"
      ]
     },
     "execution_count": 297,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = np.loadtxt('stockholm_td_adj.dat.bz2')\n",
    "np.unique(data[:,1]) # the month column takes values from 1 to 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mask_feb = data[:,1] == 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.2121095707365961"
      ]
     },
     "execution_count": 225,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the temperature data is in column 3\n",
    "np.mean(data[mask_feb,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these tools we have very powerful data processing capabilities at our disposal. For example, to extract the average monthly average temperatures for each month of the year only takes a few lines of code: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFUhJREFUeJzt3X+w7HV93/HnCxij8ktM9Z5EmnvBJFAhSO50jAlJczRh\npLUBx0lSMWlUMhmbBLU1k4FoO/fa6USxkVRrbYsig4zoKNUKlShSPE6oraJw+S3BUH4p9/gDm4Bt\nqXLf/WO/95tzr+fs2d2zu9/dvc/HzA67393z2fe5nHNe+/l+Pt/PJ1WFJEkAh3VdgCRpdhgKkqSW\noSBJahkKkqSWoSBJahkKkqTWVEIhyfFJbkhyZ5Lbk7y+OX5ckuuS3JPk00mOnUY9kqT1ZRrXKSRZ\nApaqak+So4AvA+cArwG+XVVvT3IBcFxVXTjxgiRJ65pKT6Gq9lbVnub+48DdwPH0guHy5mWXAy+b\nRj2SpPVNpadwwBsmO4AV4FTgoao6bs1zj1bVM6dakCSpNdWB5ubU0VXAG5oew8GJ5JobktShI6b1\nRkmOoBcIV1TVJ5rDq0m2VdVqM+7wjQ2+1rCQpBFUVYZ5/TR7Cu8H7qqqd645djXw6ub+q4BPHPxF\n+1XV3N527drVeQ2Hav3zXLv1d3+b9/pHMZWeQpIzgN8Abk9yC73TRG8CLgI+kuQ84AHg16dRjyRp\nfVMJhar6b8DhGzz9y9OoQZK0Oa9onoLl5eWuS9iSea5/nmsH6+/avNc/iqlPSR1FkpqHOiVpliSh\nZnigWZI04wwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLL\nUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJC2YGlpB0nGclta2tH1tyO5R7O0FUmAcf1sBn/O\nNU7u0SxJ2hJDQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhSkGeYV05o2r2iWtmDSVzR7\nxbS2wiuaJUlbYihIklqGgiSpZShIklpTCYUklyZZTXLbmmO7kjyc5ObmdtY0apEkbWxaPYXLgJes\nc/ziqtrZ3D41pVokSRuYSihU1Y3Ad9Z5aqipUpKkyep6TOH8JHuSvC/JsR3XIkmHvCM6fO/3AP+y\nqirJvwIuBn57oxfv3r27vb+8vMzy8vKk65OkubKyssLKysqW2pjaFc1JtgPXVNVpwzzXPO8VzZpJ\nXtGsWTbrVzSHNWMISZbWPPdy4I4p1iJJWsdUTh8luRJYBn44yYPALuBFSU4H9gH3A6+dRi06tCwt\n7WB19YGxtbdt23b27r1/bO1Js8YF8bTQxnv6BQ4+BePpI82yWT99JEmacYaCJKllKEiSWoaCJKll\nKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKk1Uigk\n+S/jLkSS1L2Rdl5L8iNV9cgE6tno/dx5TSNx57X+7WuxjbLz2kB7NCd5CnAyvZ/Oe6YZCJKk6dk0\nFJK8FPgPwF8CAU5I8tqq+rNJFydJmq5NTx8l+QrwD6vqq83j5wKfrKqTp1Df/ho8faSRePqof/ta\nbKOcPhpkoPmx/YHQuA94bKjKJElzYZAxhS8luRb4CL2PLL8G3JTk5QBV9bEJ1idJmqJBTh9d1ufp\nqqrzxlvSujV4+kgj8fRR//a12EY5fTTSlNRpMxQ0KkOhf/tLSztYXX1gLK1v27advXvvH0tbGo+J\nhEKSE4DXATtYc7qpqs4eocaRGAoalaHQbfvq1qSuU/jPwKXANcC+UQqTJM2HQULh/1bVuyZeiSSp\nc4OcPnol8BPAdcAT+49X1c2TLe2AGjx9pJF4+qjb9tWtSZ0++ingHwMv5m9OH1XzWJK0QAYJhV8D\nTqyq/zfpYiRJ3RrkiuY7gGdMuhBJUvcG6Sk8A/hKkps4cExhalNSJUnTMUgo7Jp4FZKkmTDQFc1J\ntgM/UVXXJ3k6cHhVTW1RPGcfaVTOPuq2fXVrIqukJvkd4CrgPzaHnkPvgrZhCrs0yWqS29YcOy7J\ndUnuSfLpJMcO06YkafwGGWj+feAM4K8Bqupe4NlDvs9lwEsOOnYhcH1VnQTcAPzRkG1KksZskFB4\nYu101CRHMGR/s6puBL5z0OFzgMub+5cDLxumTUnS+A0SCp9L8ibgaUnOBD5Kbx2krXp2Va0CVNVe\nhu99SJLGbJBQuBD4JnA78Frg2qp68wRqcYRKkjo2yJTU11XVO4H37j+Q5A3Nsa1YTbKtqlaTLAHf\n6Pfi3bt3t/eXl5dZXl7e4ttL0mJZWVlhZWVlS20MsiDezVW186Bjt1TVTw/1RskO4Jqq+qnm8UXA\no1V1UZILgOOq6sINvtYpqRqJU1K7bV/dGuuCeEnOBV4JnJDk6jVPHQ08OmRhVwLLwA8neZDeBXFv\nAz6a5DzgAeDXh2lTkjR+G/YUmgvWTgDeSm9cYb/HgNuq6vuTL6+txZ6CRmJPodv21S33aNZcmuQ+\nwYZCt+2rW4aC5tIk/zAZCt22r25NZJkLSdKhw1CQJLVGCoUku8dchyRpBozaU/jyWKuQJM0EB5rV\nOQea1297EdpXt8Z68dqaRt+1zuG/Ar5UVZ8Y5s0kSbNtkNNHTwVOB+5tbqcBxwO/neTfTLA2SdKU\nDbIg3mnAGVX1JECSfw/8OfDz9FZOlSQtiEF6CscBR615fCTwzCYknphIVZKkTgzSU3g7sCfJChDg\n7wF/nORI4PoJ1iZJmrKBZh8l+RHgBc3Dm6rq6xOt6gff39lHC8zZR+u3vQjtq1uTmn10DXAlcHVV\nfXfU4iRJs2+QMYU/AX4BuCvJVUl+NclTJ1yXJKkDA1+8luRw4MXA7wBnVdUxkyzsoPf29NEC8/TR\n+m0vQvvq1kROHzUNPw34FeAfATuBy4cvT5I06wYZU/gIvUHmTwHvBj5XVfsmXZgkafoG6SlcCpy7\n/+I1SdLiGnRK6qnA8+gteQFAVX1ggnUd/P6OKSwwxxTWb3sR2le3JjUldRewTC8UrgX+PnAjMLVQ\nkCRNxyBTUn8V+CVgb1W9Bng+cOxEq5IkdWKQUPg/zcDy95McA3wD+NuTLUuS1IVBBpq/lOQZwHvp\n7bj2OPDfJ1qVJKkTQ+28lmQHcExV3TapgjZ4XweaF5gDzeu3vQjtq1ujDDQPtUdzVd0/7UCQNL+W\nlnaQZCy3paUdXX87hwT3aFbn7Cms37btb96++pt4T0GStNg2DYUk70hyyjSKkSR1a5Cewt3AJUm+\nkOSfJPEaBUlaUJuGQlW9r6rOAH4L2AHcluTKJC+adHGSpOkaaEyh2Uvh5Ob2LeBW4I1JPjzB2iRJ\nU7bp7KMkf0pvL4X/ClxaVV9c89w9VXXSZEt09tGic/bR+m3b/ubtq79JbbJzG/DPN9if+QXDvJkk\nabZt2FNIsrPfF1bVzROpaP1a7CksMHsK67dt+5u3r/7G3VN4R5/nit5+zZKkBdL5Fc1J7gf+CtgH\nfK+qfuCUlD2FxWZPYf22bX/z9tXfpMYUSPJz9Kajtq8f485r+4DlqvrOmNqTJI1okJ3XrgCeC+wB\n9u/TXIxv57XgchuSNBMG6Sn8XeB5Ezx/U8BnkjwJXFJV753Q+0iSNjFIKNwBLAGPTKiGM6rqkSTP\nohcOd1fVjRN6L0lSHxuGQpJr6H2KPxq4K8kXgSf2P19VZ4+jgKp6pPnvN5N8nN61Dz8QCrt3727v\nLy8vs7y8PI63l6SFsbKywsrKypba6Hedwi/2+8Kq+tyW3rn3Hk8HDquqx5McCVwHvKWqrjvodc4+\nWmDOPlq/bdvfvH31N9bZR/v/6Ce5qKouOOiNLgK2HArANuDjSaqp5YMHB4IkaXoGWfvo5qraedCx\n26rqtIlWduD72VNYYPYU1m/b9jdvX/2NtaeQ5HeB3wNOTLJ2X+ajgc+PVqIkaZb1G1M4FjgOeCtw\n4ZqnHquqR6dQ29pa7CksMHsK67dt+5u3r/5G6SkMtMxFs5/CNg68ovnBoSsckaGw2AyF9du2/c3b\nV38TWeYiyfnAbmCV3pIU0Pu/PLUxBUnSdAxy8do/BU6qqm9PuhhJUrcGWXPoIXqrmOoQtbS0gyRj\nuS0t7ej625HUxyBTUi8FTgI+yYFXNF882dIOqMExhQ7N83lnxxQWu331N6mlsx9sbk9pbpKkBTXw\nJjtJjgKoqscnWtH6721PoUPz/GnSnsJit6/+RukpbDqmkOTUJLcAdwJ3JvlyklNGLVKSNLsGGWi+\nBHhjVW2vqu3AHwDueSBJC2iQUDiyqj67/0FVrQBHTqwiSVJnBhlovi/JvwCuaB7/JnDf5EqSJHVl\nkJ7CecCzgI81t2c1xyRJC2bg2UddcvZRt+Z5hoqzjxa7ffU37qWzr+73hePajlOSNDv6jSn8LL0l\nLj4EfAEYKm0kSfOnXygsAWcC5wKvpLfMxYeq6s5pFCZJmr4NB5qr6smq+lRVvQp4IfBVYKVZSluS\ntID6TklN8kPAS+n1FnYA7wI+PvmyJEld6DfQ/AHgVOBa4C1VdcfUqpIkdaLfHs37gO82D9e+KEBV\n1TETrm1tLU5J7dA8T1t0Supit6/+xjoltaoGubBNkrRA/MMvaW65K+D4eUWzNjXPpxg8fWT7W2l/\n3k1kPwVJ0qHDUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVKr81BI\nclaSryT5iyQXdF2PJB3KOg2FJIcB7wZeApwCnJvk5C5rkqRDWdc9hRcA91bVA1X1PeDDwDkd1yRJ\nh6yuQ+E5wENrHj/cHJMkdWDDnddmze7du9v7y8vLLC8vj63tpaUdrK4+MJa2tm3bzt699y9U+9u2\nbWd1dagl2fu2v96xSbU/zrYn3f60/21sf/P2J/27NW4rKyusrKxsqY1ON9lJ8kJgd1Wd1Ty+kN7+\nzxcd9LqJbrIz7xuBuNGINBnz/rs1j5vs3AT8eJLtSZ4CvAK4uuOaJOmQ1enpo6p6Msn5wHX0AurS\nqrq7y5ok6VDmHs3M/+mdee/iSrNq3n+35vH0kSRphhgKkqSWoSBJahkKC6A3vzpjua03V1vSocOB\nZhwIlrS+ef/ddaBZkrQlhoIkbeBQPDXr6SM8fSRpMXn6SJK0JYaCJKllKEiSWoaCJKllKEiSWoaC\nJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKHBorpkuSetxP4UpcD8FSV1w\nP4UZZU9E0rywpyBJC8qegiRpSwwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktToL\nhSS7kjyc5ObmdlZXtUiSerruKVxcVTub26c6rmViVlZWui5hS+a5/nmuHay/a/Ne/yi6DoWh1uSY\nV/P+gzXP9c9z7WD9XZv3+kfRdSicn2RPkvclObbjWiTpkDfRUEjymSS3rbnd3vz3V4D3ACdW1enA\nXuDiSdYiSdrcTCydnWQ7cE1VnbbB890XKUlzaNils4+YVCGbSbJUVXubhy8H7tjotcN+U5Kk0XQW\nCsDbk5wO7APuB17bYS2SJGbk9JEkaTZ0PfuoryRnJflKkr9IckHX9QwjyfFJbkhyZzPA/vquaxpF\nksOaiwuv7rqWYSU5NslHk9zd/H/4ma5rGkaSf5bkjmZyxgeTPKXrmvpJcmmS1SS3rTl2XJLrktyT\n5NOzPMtwg/rf3vz87Enyn5Ic02WN/axX/5rn/iDJviTP3KydmQ2FJIcB7wZeApwCnJvk5G6rGsr3\ngTdW1SnAzwK/P2f17/cG4K6uixjRO4Frq+rvAM8H7u64noEl+VHgdcDOZgLGEcAruq1qU5fR+31d\n60Lg+qo6CbgB+KOpVzW49eq/DjilmSV5L/NXP0mOB84EHhikkZkNBeAFwL1V9UBVfQ/4MHBOxzUN\nrKr2VtWe5v7j9P4gPafbqobT/DD9A+B9XdcyrOYT3S9U1WUAVfX9qvrrjssa1uHAkUmOAJ4OfL3j\nevqqqhuB7xx0+Bzg8ub+5cDLplrUENarv6qur6p9zcP/ARw/9cIGtMG/P8CfAn84aDuzHArPAR5a\n8/hh5uyP6n5JdgCnA1/otpKh7f9hmseBpxOAbyW5rDn9dUmSp3Vd1KCq6uvAO4AHga8B/6uqru+2\nqpE8u6pWofdBCXh2x/VsxXnAn3VdxDCSnA08VFW3D/o1sxwKCyHJUcBVwBuaHsNcSPJSYLXp7YT5\nW5LkCGAn8O+qaifwv+mdypgLSZ5B71P2duBHgaOSvLLbqsZiHj9gkOTNwPeq6squaxlU8yHoTcCu\ntYc3+7pZDoWvAT+25vHxzbG50XT7rwKuqKpPdF3PkM4Azk5yH/Ah4EVJPtBxTcN4mN4npC81j6+i\nFxLz4peB+6rq0ap6EvgY8HMd1zSK1STboHdtEvCNjusZWpJX0zuNOm+h/FxgB3Brkv9J72/ol5P0\n7a3NcijcBPx4ku3NrItXAPM2A+b9wF1V9c6uCxlWVb2pqn6sqk6k929/Q1X9Vtd1Dao5ZfFQkp9s\nDv0S8zVg/iDwwiRPTRJ69c/DQPnBvcqrgVc3918FzPqHowPqb5b0/0Pg7Kp6orOqBtfWX1V3VNVS\nVZ1YVSfQ+6D001XVN5hnNhSaT0fn0xv9vxP4cFXNwy8FAEnOAH4DeHGSW9wzohOvBz6YZA+92Ud/\n3HE9A6uqL9Lr3dwC3ErvF/2STovaRJIrgc8DP5nkwSSvAd4GnJnkHnrB9rYua+xng/r/LXAU8Jnm\nd/g9nRbZxwb1r1UMcPrIi9ckSa2Z7SlIkqbPUJAktQwFSVLLUJAktQwFSVLLUJAktQwFCWiWFf7A\nmseHJ/nmqEuGN8t2/+6ax7+Y5Jpx1CpNkqEg9XwXODXJDzWPz+TABRmHdRzwewcd86IgzTxDQfob\n1wIvbe6fS2/NJ6DdLObjSW5N8vkkpzbHdzWbm3w2yVeTnN98yVuBE5urYC9qjh29ZtOfK6b1TUnD\nMBSknqK3Z8e5TW/hNA5c6vwtwM1V9XzgzcDaP+on0etZ/AywO8nh9FZk/cuq2llV+3cNPJ3e0hvP\nA56bZB4XuNOCMxSkRlXdQW9VyXOBT3LgOjE/TxMEVfVZ4JnNsugAn2w28fk2sAps2+AtvlhVj1Rv\nbZk9zXtJM+WIrguQZszVwL8GloG/NeDXrF09cx8b/16tfd2TfV4ndcaegtSzv1fwfuAtVXXnQc//\nOfCbAEmWgW9tsmnSY8DR4y5SmjQ/qUg9BVBVXwPevc7zu4H3J7mV3kyljfaW2N/Oo82A9G30tnC8\ndr3XSbPGpbMlSS1PH0mSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKn1/wE1DyL6+dxh\noQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7feef2500390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "months = np.arange(1,13)\n",
    "monthly_mean = [np.mean(data[data[:,1] == month, 3]) for month in months]\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(months, monthly_mean)\n",
    "ax.set_xlabel(\"Month\")\n",
    "ax.set_ylabel(\"Monthly avg. temp.\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculations with higher-dimensional data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When functions such as `min`, `max`, etc. are applied to a multidimensional arrays, it is sometimes useful to apply the calculation to the entire array, and sometimes only on a row or column basis. Using the `axis` argument we can specify how these functions should behave: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.90258753,  0.44483294,  0.59624183],\n",
       "       [ 0.26680406,  0.87446277,  0.61826539],\n",
       "       [ 0.47307416,  0.79684284,  0.5340052 ]])"
      ]
     },
     "execution_count": 227,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = random.rand(3,3)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.90258752648090268"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# global max\n",
    "m.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.90258753,  0.87446277,  0.61826539])"
      ]
     },
     "execution_count": 229,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# max in each column\n",
    "m.max(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.90258753,  0.87446277,  0.79684284])"
      ]
     },
     "execution_count": 230,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# max in each row\n",
    "m.max(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many other functions and methods in the `array` and `matrix` classes accept the same (optional) `axis` keyword argument."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reshaping, resizing and stacking arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The shape of an Numpy array can be modified without copying the underlaying data, which makes it a fast operation even for large arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4],\n",
       "       [10, 11, 12, 13, 14],\n",
       "       [20, 21, 22, 23, 24],\n",
       "       [30, 31, 32, 33, 34],\n",
       "       [40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 231,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n, m = A.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,\n",
       "        32, 33, 34, 40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 233,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = A.reshape((1,n*m))\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5,  5,  5,  5,  5, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,\n",
       "        32, 33, 34, 40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 234,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[0,0:5] = 5 # modify the array\n",
    "\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5,  5,  5,  5,  5],\n",
       "       [10, 11, 12, 13, 14],\n",
       "       [20, 21, 22, 23, 24],\n",
       "       [30, 31, 32, 33, 34],\n",
       "       [40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 235,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A # and the original variable is also changed. B is only a different view of the same data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the function `flatten` to make a higher-dimensional array into a vector. But this function create a copy of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5,  5,  5,  5,  5, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,\n",
       "       32, 33, 34, 40, 41, 42, 43, 44])"
      ]
     },
     "execution_count": 236,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = A.flatten()\n",
    "\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 10, 10, 10, 10, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,\n",
       "       32, 33, 34, 40, 41, 42, 43, 44])"
      ]
     },
     "execution_count": 237,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[0:5] = 10\n",
    "\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5,  5,  5,  5,  5],\n",
       "       [10, 11, 12, 13, 14],\n",
       "       [20, 21, 22, 23, 24],\n",
       "       [30, 31, 32, 33, 34],\n",
       "       [40, 41, 42, 43, 44]])"
      ]
     },
     "execution_count": 238,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A # now A has not changed, because B's data is a copy of A's, not refering to the same data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, a few more examples with some graphics.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "array1 = np.arange(90).reshape((-1, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/redux1.png)\n",
    "\n",
    "In this case, one of the axis length is left as -1. Numpy will compute this dimension for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 10)"
      ]
     },
     "execution_count": 291,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array1.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array has 10 columns and 9 rows. \n",
    "\n",
    "We can select one element with two indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "57"
      ]
     },
     "execution_count": 292,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array1[5,7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/redux2.png)\n",
    "\n",
    "Or we can select an entire row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([50, 51, 52, 53, 54, 55, 56, 57, 58, 59])"
      ]
     },
     "execution_count": 293,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array1[5,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/redux3.png)\n",
    "\n",
    "Part of a row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([51, 52, 53, 54, 55, 56, 57])"
      ]
     },
     "execution_count": 294,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array1[5,1:8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/redux4.png)\n",
    "\n",
    "Or a rectangle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[51, 52, 53, 54, 55, 56, 57],\n",
       "       [61, 62, 63, 64, 65, 66, 67],\n",
       "       [71, 72, 73, 74, 75, 76, 77]])"
      ]
     },
     "execution_count": 295,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array1[5:-1,1:8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/redux5.png)\n",
    "\n",
    "We can even use non-sequential indexing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([56, 78])"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array1[(5, 7), (6, 8)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/redux6.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding a new dimension: newaxis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `newaxis`, we can insert new dimensions in an array, for example converting a vector to a column or row matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v = np.array([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 241,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1],\n",
       "       [2],\n",
       "       [3]])"
      ]
     },
     "execution_count": 243,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make a column matrix of the vector v\n",
    "v[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 1)"
      ]
     },
     "execution_count": 244,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# column matrix\n",
    "v[:,np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 3)"
      ]
     },
     "execution_count": 245,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# row matrix\n",
    "v[np.newaxis,:].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stacking and repeating arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using function `repeat`, `tile`, `vstack`, `hstack`, and `concatenate` we can create larger vectors and matrices from smaller ones:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### tile and repeat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = np.array([[1, 2], [3, 4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4])"
      ]
     },
     "execution_count": 247,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# repeat each element 3 times\n",
    "np.repeat(a, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 1, 2, 1, 2],\n",
       "       [3, 4, 3, 4, 3, 4]])"
      ]
     },
     "execution_count": 248,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# tile the matrix 3 times \n",
    "np.tile(a, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### concatenate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "b = np.array([[5, 6]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2],\n",
       "       [3, 4],\n",
       "       [5, 6]])"
      ]
     },
     "execution_count": 250,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.concatenate((a, b), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 5],\n",
       "       [3, 4, 6]])"
      ]
     },
     "execution_count": 251,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.concatenate((a, b.T), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### hstack and vstack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2],\n",
       "       [3, 4],\n",
       "       [5, 6]])"
      ]
     },
     "execution_count": 252,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.vstack((a,b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 5],\n",
       "       [3, 4, 6]])"
      ]
     },
     "execution_count": 253,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.hstack((a,b.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Copy and \"deep copy\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To achieve high performance, assignments in Python usually do not copy the underlaying objects. This is important for example when objects are passed between functions, to avoid an excessive amount of memory copying when it is not necessary (technical term: pass by reference). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2],\n",
       "       [3, 4]])"
      ]
     },
     "execution_count": 254,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# now B is referring to the same array data as A \n",
    "B = A "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10,  2],\n",
       "       [ 3,  4]])"
      ]
     },
     "execution_count": 256,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# changing B affects A\n",
    "B[0,0] = 10\n",
    "\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10,  2],\n",
       "       [ 3,  4]])"
      ]
     },
     "execution_count": 257,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to avoid this behavior, so that when we get a new completely independent object `B` copied from `A`, then we need to do a so-called \"deep copy\" using the function `copy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B = A.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-5,  2],\n",
       "       [ 3,  4]])"
      ]
     },
     "execution_count": 262,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now, if we modify B, A is not affected\n",
    "B[0,0] = -5\n",
    "\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10,  2],\n",
       "       [ 3,  4]])"
      ]
     },
     "execution_count": 263,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterating over array elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, we want to avoid iterating over the elements of arrays whenever we can (at all costs). The reason is that in a interpreted language like Python (or MATLAB), iterations are really slow compared to vectorized operations. \n",
    "\n",
    "However, sometimes iterations are unavoidable. For such cases, the Python `for` loop is the most convenient way to iterate over an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "v = np.array([1,2,3,4])\n",
    "\n",
    "for element in v:\n",
    "    print(element)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('row', array([1, 2]))\n",
      "1\n",
      "2\n",
      "('row', array([3, 4]))\n",
      "3\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "M = np.array([[1,2], [3,4]])\n",
    "\n",
    "for row in M:\n",
    "    print(\"row\", row)\n",
    "    \n",
    "    for element in row:\n",
    "        print(element)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we need to iterate over each element of an array and modify its elements, it is convenient to use the `enumerate` function to obtain both the element and its index in the `for` loop: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('row_idx', 0, 'row', array([1, 2]))\n",
      "('col_idx', 0, 'element', 1)\n",
      "('col_idx', 1, 'element', 2)\n",
      "('row_idx', 1, 'row', array([3, 4]))\n",
      "('col_idx', 0, 'element', 3)\n",
      "('col_idx', 1, 'element', 4)\n"
     ]
    }
   ],
   "source": [
    "for row_idx, row in enumerate(M):\n",
    "    print(\"row_idx\", row_idx, \"row\", row)\n",
    "    \n",
    "    for col_idx, element in enumerate(row):\n",
    "        print(\"col_idx\", col_idx, \"element\", element)\n",
    "       \n",
    "        # update the matrix M: square each element\n",
    "        M[row_idx, col_idx] = element ** 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  4],\n",
       "       [ 9, 16]])"
      ]
     },
     "execution_count": 271,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# each element in M is now squared\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vectorizing functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned several times by now, to get good performance we should try to avoid looping over elements in our vectors and matrices, and instead use vectorized algorithms. The first step in converting a scalar algorithm to a vectorized algorithm is to make sure that the functions we write work with vector inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Theta(x):\n",
    "    \"\"\"\n",
    "    Scalar implemenation of the Heaviside step function.\n",
    "    \"\"\"\n",
    "    if x >= 0:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-274-d55419725688>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mTheta\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-273-9a0cb13d93d4>\u001b[0m in \u001b[0;36mTheta\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mScalar\u001b[0m \u001b[0mimplemenation\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mHeaviside\u001b[0m \u001b[0mstep\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \"\"\"\n\u001b[0;32m----> 5\u001b[0;31m     \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
     ]
    }
   ],
   "source": [
    "Theta(np.array([-3,-2,-1,0,1,2,3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, that didn't work because we didn't write the `Theta` function so that it can handle a vector input... \n",
    "\n",
    "To get a vectorized version of Theta we can use the Numpy function `vectorize`. In many cases it can automatically vectorize a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Theta_vec = np.vectorize(Theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 1, 1, 1, 1])"
      ]
     },
     "execution_count": 276,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Theta_vec(np.array([-3,-2,-1,0,1,2,3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also implement the function to accept a vector input from the beginning (requires more effort but might give better performance):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Theta(x):\n",
    "    \"\"\"\n",
    "    Vector-aware implemenation of the Heaviside step function.\n",
    "    \"\"\"\n",
    "    return 1 * (x >= 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 1, 1, 1, 1])"
      ]
     },
     "execution_count": 278,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Theta(np.array([-3,-2,-1,0,1,2,3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 279,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# still works for scalars as well\n",
    "Theta(-1.2), Theta(2.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def funct3(val):\n",
    "    if val > 3: # Doesn’t work with array\n",
    "        x = 2\n",
    "    else:\n",
    "        x = 5\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function does not work with arrays.\n",
    "\n",
    "So, we have to redefine it with `vectorize`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 2, 2, 5, 5])"
      ]
     },
     "execution_count": 289,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vfunct3 = np.vectorize(funct3)\n",
    "vfunct3(np.array([3,4,5,1,2])) # Now works!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using arrays in conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using arrays in conditions,for example `if` statements and other boolean expressions, one needs to use `any` or `all`, which requires that any or all elements in the array evalutes to `True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  4],\n",
       "       [ 9, 16]])"
      ]
     },
     "execution_count": 280,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "at least one element in M is larger than 5\n"
     ]
    }
   ],
   "source": [
    "if (M > 5).any():\n",
    "    print(\"at least one element in M is larger than 5\")\n",
    "else:\n",
    "    print(\"no element in M is larger than 5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "all elements in M are not larger than 5\n"
     ]
    }
   ],
   "source": [
    "if (M > 5).all():\n",
    "    print(\"all elements in M are larger than 5\")\n",
    "else:\n",
    "    print(\"all elements in M are not larger than 5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Type casting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since Numpy arrays are *statically typed*, the type of an array does not change once created. But we can explicitly cast an array of some type to another using the `astype` functions (see also the similar `asarray` function). This always create a new array of new type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int64')"
      ]
     },
     "execution_count": 283,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.,   4.],\n",
       "       [  9.,  16.]])"
      ]
     },
     "execution_count": 284,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M2 = M.astype(float)\n",
    "\n",
    "M2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('float64')"
      ]
     },
     "execution_count": 285,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M2.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True],\n",
       "       [ True,  True]], dtype=bool)"
      ]
     },
     "execution_count": 286,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M3 = M.astype(bool)\n",
    "\n",
    "M3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* http://numpy.scipy.org\n",
    "* http://scipy.org/Tentative_NumPy_Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "2.7.12 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]"
        },
        {
         "module": "IPython",
         "version": "5.1.0"
        },
        {
         "module": "OS",
         "version": "Linux 4.4.0 47 generic x86_64 with debian jessie sid"
        },
        {
         "module": "numpy",
         "version": "1.11.2"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>2.7.12 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]</td></tr><tr><td>IPython</td><td>5.1.0</td></tr><tr><td>OS</td><td>Linux 4.4.0 47 generic x86_64 with debian jessie sid</td></tr><tr><td>numpy</td><td>1.11.2</td></tr><tr><td colspan='2'>Tue Dec 20 17:36:46 2016 PST</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 2.7.12 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] \\\\ \\hline\n",
       "IPython & 5.1.0 \\\\ \\hline\n",
       "OS & Linux 4.4.0 47 generic x86\\_64 with debian jessie sid \\\\ \\hline\n",
       "numpy & 1.11.2 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Tue Dec 20 17:36:46 2016 PST} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 2.7.12 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]\n",
       "IPython 5.1.0\n",
       "OS Linux 4.4.0 47 generic x86_64 with debian jessie sid\n",
       "numpy 1.11.2\n",
       "Tue Dec 20 17:36:46 2016 PST"
      ]
     },
     "execution_count": 298,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%reload_ext version_information\n",
    "\n",
    "%version_information numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}