NumPy Lectures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# NumPy\n",
    "\n",
    "[NumPy](http://www.numpy.org/), short for Numerical Python, is the core Python package for numerical computing. The main features of NumPy are:\n",
    "\n",
    "* $n$-dimensional array object `ndarray`\n",
    "* Vectorized operations and functions which broadcast across arrays for fast computation\n",
    "\n",
    "To get started with NumPy, let's adopt the standard convention and import it using the name `np`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## NumPy Arrays\n",
    "\n",
    "The fundamental object provided by the NumPy package is the `ndarray`. We can think of a 1D (1-dimensional) `ndarray` as a list, a 2D (2-dimensional) `ndarray` as a matrix, a 3D (3-dimensional) `ndarray` as a 3-tensor (or a \"cube\" of numbers), and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Creating Arrays\n",
    "\n",
    "The function `np.array` creates a NumPy array from a Python sequence such as a list, a tuple or a list of lists. For example, create a 1D NumPy array from a Python list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4 5]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1, 2, 3, 4, 5])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note that when we print a NumPy array it looks a lot like a Python list except that the entries are separated by spaces whereas entries in a Python list are separated by commas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 3, 4, 5]\n"
     ]
    }
   ],
   "source": [
    "print([1, 2, 3, 4, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note also that a NumPy array is displayed slightly differently when output by a cell (as opposed to being explicitly printed to output by the `print` function):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 3, 4, 5])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Use the built-in function `type` to verify the type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Nested sequences, like a list of lists with the same length, are converted into a 2D NumPy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "M = np.array([[1, 2, 3], [4, 5, 6]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create an $n$-dimensional NumPy array from nested Python lists. For example, the following is a 3D NumPy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 1  2]\n",
      "  [ 3  4]]\n",
      "\n",
      " [[ 5  6]\n",
      "  [ 7  8]]\n",
      "\n",
      " [[ 9 10]\n",
      "  [11 12]]]\n"
     ]
    }
   ],
   "source": [
    "N = np.array([ [[1, 2], [3, 4]] , [[5, 6], [7, 8]] , [[9, 10], [11, 12]] ])\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are several NumPy functions for [creating arrays](https://docs.scipy.org/doc/numpy/user/quickstart.html#array-creation):\n",
    "\n",
    "| Function | Description |\n",
    "| ---: | :--- |\n",
    "| `np.array(a)` | Create $n$-dimensional NumPy array from sequence `a` |\n",
    "| `np.linspace(a,b,N)` | Create 1D NumPy array with `N` equally spaced values from `a` to `b` (inclusively)|\n",
    "| `np.arange(a,b,step)` | Create 1D NumPy array with values from `a` to `b` (exclusively) incremented by `step`|\n",
    "| `np.zeros(N)` | Create 1D NumPy array of zeros of length $N$ |\n",
    "| `np.zeros((n,m))` | Create 2D NumPy array of zeros with $n$ rows and $m$ columns |\n",
    "| `np.ones(N)` | Create 1D NumPy array of ones of length $N$ |\n",
    "| `np.ones((n,m))` | Create 2D NumPy array of ones with $n$ rows and $m$ columns |\n",
    "| `np.eye(N)` | Create 2D NumPy array with $N$ rows and $N$ columns with ones on the diagonal (ie. the identity matrix of size $N$) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create a 1D NumPy array with 11 equally spaced values from 0 to 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(0, 1, 11)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Create a 1D NumPy array with values from 0 to 20 (exclusively) incremented by 2.5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.   2.5  5.   7.5 10.  12.5 15.  17.5]\n"
     ]
    }
   ],
   "source": [
    "y = np.arange(0, 20, 2.5)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "These are the functions that we'll use most often when creating NumPy arrays. The function `np.linspace` works best when we know the *number of elements* we want in the array, and `np.arange` works best when we know *step size* between values in the array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create a 1D NumPy array of zeros of length 5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "z = np.zeros(5)\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Create a 2D NumPy array of zeros with 2 rows and 5 columns (note that the argument is a tuple):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "M = np.zeros((2, 5))\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create a 1D NumPy array of ones of length 7:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "w = np.ones(7)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Create a 2D NumPy array of ones with 3 rows and 2 columns (note that the argument is a tuple):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]]\n"
     ]
    }
   ],
   "source": [
    "N = np.ones((3, 2))\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create the identity matrix of size 10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "I = np.eye(10)\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "When dealing with NumPy array, often a dimension are also called **axis**. The only difference is that while dimensions are numbered starting at 1, axis are numbered starting at 0. So, the first dimension correponds to axis 0, the second dimension to axis 1, and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Dimension, Shape and Size\n",
    "\n",
    "We can think of a 1D NumPy array as a list of numbers, a 2D NumPy array as a matrix, a 3D NumPy array as a cube of numbers, and so on. Given a NumPy array, we can find out how many dimensions it has by accessing its `.ndim` attribute. The result is a number telling us how many dimensions it has.\n",
    "\n",
    "For example, create a 2D NumPy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2]\n",
      " [3 4]\n",
      " [5 6]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2], [3, 4], [5, 6]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.ndim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The result tells us that `A` has 2 dimensions. The first dimension corresponds to the vertical direction counting the rows and the second dimension corresponds to the horizontal direction counting the columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can find out how many rows and columns `A` has by accessing its `.shape` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 2)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The result is a tuple `(3, 2)` of length 2 which means that `A` is a 2D array with 3 rows and 2 columns.\n",
    "\n",
    "We can also find out how many entries `A` has in total by accessing its `.size` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is the expected result since we know that `A` has 3 rows and 2 columns and therefore 2 $\\times$ 3 = 6 total entries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The attributes `.ndim`, `.shape` and `.size` return the same results of the corresponding functions `np.ndim`, `np.shape` and `np.size` when used with the same argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[9 3 1 7]\n"
     ]
    }
   ],
   "source": [
    "r = np.array([9, 3, 1, 7])\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ndim(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4,)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.size(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note that `np.shape(r)` is a tuple with a single entry `(4,)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Array Datatype\n",
    "\n",
    "NumPy arrays are *homogeneous*: all entries in the array are the same datatype. There are different kinds of datatypes provided by NumPy for different applications but we'll mostly be working with the default integer type `np.int64` and the default float type `np.float64`. These are very similar to the built-in Python datatypes `int` and `float` but with some differences that we won't go into. Check out the NumPy documentation on [numeric datatypes](https://docs.scipy.org/doc/numpy/user/basics.types.html) for more information.\n",
    "\n",
    "The most important point for now is to know how to determine if a NumPy array contains integers elements or float elements. We can access the datatype of a NumPy array by its `.dtype` attribute. For example, create a 2D NumPy array from a list of lists of integers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3], [4, 5, 6]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We expect the datatype of `A` to be integers and we verify:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int64')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Most of the other NumPy functions which create arrays use the `np.float64` datatype by default. For example, using `np.linspace`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.   0.25 0.5  0.75 1.  ]\n"
     ]
    }
   ],
   "source": [
    "u = np.linspace(0, 1, 5)\n",
    "print(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('float64')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Notice that numbers are printed with a decimal point when the datatype of the NumPy array is any kind of float."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can explicitly convert or *cast* a NumPy array from one data type to another using the `.astype` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 1]\n"
     ]
    }
   ],
   "source": [
    "u_int = u.astype(np.int64)\n",
    "print(u_int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int64')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_int.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can explicitly define the type of the array data when we create it, using the `dtype` keyword argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[1, 2], [3, 4]], dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 2.]\n",
      " [3. 4.]]\n"
     ]
    }
   ],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Array Operations\n",
    "\n",
    "[Arithmetic operators](https://docs.scipy.org/doc/numpy/user/quickstart.html#basic-operations) including addition `+`, subtraction `-`, multiplication `*`, division `/` and exponentiation `**` are applied to arrays *element-wise*. For addition and substraction, these are the familiar vector operations we see in linear algebra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "v = np.array([1,2,3])\n",
    "w = np.array([1,0,-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 2 2]\n"
     ]
    }
   ],
   "source": [
    "print(v + w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2 4]\n"
     ]
    }
   ],
   "source": [
    "print(v - w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the same way, array multiplication and division are performed element by element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1  0 -3]\n"
     ]
    }
   ],
   "source": [
    "print(v * w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.          0.         -0.33333333]\n"
     ]
    }
   ],
   "source": [
    "print(w / v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Notice that the datatype of both `v` and `w` is `numpy.int64` however division `w / v` returns an array with datatype `numpy.float64`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The exponent operator `**` also acts element by element in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 4 9]\n"
     ]
    }
   ],
   "source": [
    "print(v ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's see these operations for 2D arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[3,1],[2,-1]])\n",
    "B = np.array([[2,-2],[5,1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5 -1]\n",
      " [ 7  0]]\n"
     ]
    }
   ],
   "source": [
    "print(A + B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1  3]\n",
      " [-3 -2]]\n"
     ]
    }
   ],
   "source": [
    "print(A - B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.5 -0.5]\n",
      " [ 0.4 -1. ]]\n"
     ]
    }
   ],
   "source": [
    "print(A / B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 6 -2]\n",
      " [10 -1]]\n"
     ]
    }
   ],
   "source": [
    "print(A * B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[9 1]\n",
      " [4 1]]\n"
     ]
    }
   ],
   "source": [
    "print(A ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Notice that array multiplication and exponentiation are performed element-wise.\n",
    "\n",
    "In Python 3.5+, the symbol `@` computes matrix multiplication for NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, -5],\n",
       "       [-1, -5]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A @ B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "[Matrix powers](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.matrix_power.html) are performed by the function `numpy.linalg.matrix_power`. It's a long function name and so it's convenient to import it with a shorter name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from numpy.linalg import matrix_power as mpow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Compute $A^3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[37,  9],\n",
       "       [18,  1]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpow(A,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Equivalently, use the `@` operator to compute $A^3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[37,  9],\n",
       "       [18,  1]])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A @ A @ A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Conditional operations between arrays of the same size yield boolean array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ True False  True]\n",
      " [False  True False]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([[1, 2, 3], [4, 5, 6]])\n",
    "b = np.array([[0., 4., 1.], [7., 2., 12.]])\n",
    "print(a > b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Broadcasting\n",
    "\n",
    "We know from linear algebra that we can only add matrices of the same size. [Broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) is a set of NumPy rules which relaxes this constraint and allows us to combine a smaller array with a bigger when it makes sense.\n",
    "\n",
    "For example, suppose we want to create a 1D NumPy array of $y$ values for $x=0.0,0.25,0.5,0.75,1.0$ for the function $y = x^2 + 1$. From what we've seen so far, it makes sense to create `x`, then `x**2` and then add an array of ones `[1. 1. 1. 1. 1.]`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.     1.0625 1.25   1.5625 2.    ]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([0, 0.25, 0.5, 0.75, 1.0])\n",
    "y = x ** 2 + np.array([1, 1, 1, 1, 1])\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "An example of broadcasting in NumPy is the following equivalent operation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.     1.0625 1.25   1.5625 2.    ]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([0, 0.25, 0.5, 0.75, 1.0])\n",
    "y = x ** 2 + 1\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The number 1 is a scalar and we are adding it to a 1D NumPy array of length 5. The broadcasting rule in this case is to broadcast the scalar value 1 across the larger array. The result is a simpler syntax for a very comman operation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's try another example. What happens when we try to add a 1D NumPy array of length 4 to a 2D NumPy array of size 3 by 4?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4]\n"
     ]
    }
   ],
   "source": [
    "u = np.array([1, 2, 3, 4])\n",
    "print(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 1 1 1]\n",
      " [2 2 2 2]\n",
      " [3 3 3 3]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 3 4 5]\n",
      " [3 4 5 6]\n",
      " [4 5 6 7]]\n"
     ]
    }
   ],
   "source": [
    "result = A + u\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The 1D NumPy array is broadcast across the 2D array because the length of the first dimension in each array are equal!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Broadcasting can also be used with comparison operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ True  True  True False  True False False False]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1, 2, 3, -1, 2, 0, -3, -1])\n",
    "print(a > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[False False False  True False False False  True]\n"
     ]
    }
   ],
   "source": [
    "print(a == -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Broadcasting rules are more complex than in these examples, but we will not discuss them in detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Accessing Arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Accessing the entries in an array is called *indexing* and accessing rows and columns (or subarrays) is called *slicing*. See the NumPy documentation for more information about [indexing and slicing](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Slicing\n",
    "\n",
    "Create a 1D NumPy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.  0.5 1.  1.5 2.  2.5 3.  3.5 4.  4.5 5. ]\n"
     ]
    }
   ],
   "source": [
    "v = np.linspace(0, 5, 11)\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Access the entries in a 1D array using the square brackets notation just like a Python list. For example, access the entry at index 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Notice that NumPy array indices start at 0 just like Python sequences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create a 2D array of integers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[6 5 3 1 1]\n",
      " [1 0 4 0 1]\n",
      " [5 9 2 2 9]]\n"
     ]
    }
   ],
   "source": [
    "B = np.array([[6, 5, 3, 1, 1], [1, 0, 4, 0, 1], [5, 9, 2, 2, 9]])\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Access the entries in a 2D array using the square brackets with 2 indices. In particular, access the entry at row index 1 and column index 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[1, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Access the top left entry in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Negative indices work for NumPy arrays as they do for Python sequences. Access the bottom right entry in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[-1, -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Access the row at index 2 using the colon `:` syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 9, 2, 2, 9])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[2, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Access the column at index 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 2])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[:, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Select the subarray of rows at index 1 and 2, and columns at index 2, 3 and 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[4 0 1]\n",
      " [2 2 9]]\n"
     ]
    }
   ],
   "source": [
    "subB = B[1:3, 2:5]\n",
    "print(subB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Slices of NumPy arrays are again NumPy arrays but possibly of a different dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subB.ndim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 3)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subB.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(subB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The variable `subB` is assigned to a 2D NumPy array of shape 2 by 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's do the same for the column at index 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3 4 2]\n"
     ]
    }
   ],
   "source": [
    "colB = B[:,2]\n",
    "print(colB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "colB.ndim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "colB.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(colB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The variable `colB` is assigned to a 1D NumPy array of length 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One important distinction from Python's lists is that array slices are *views* on the original array. This means that the data is not copied, and any modification to the view will be reflected in the source array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[6 5 3 1 1]\n",
      " [1 0 4 0 1]\n",
      " [5 9 2 2 9]]\n"
     ]
    }
   ],
   "source": [
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "colB[0] = 17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 6  5 17  1  1]\n",
      " [ 1  0  4  0  1]\n",
      " [ 5  9  2  2  9]]\n"
     ]
    }
   ],
   "source": [
    "print(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "That is, **array slices are mutable**: if they are assigned a new value the original array from which the slice was extracted is modified:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If we omit an index of a $n$-dimensional array it returns, in general, a $(n-1)$-dimensional array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 0 4 0 1]\n"
     ]
    }
   ],
   "source": [
    "print(B[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For 2D arrays, it is particularly useful remembering that *rows* correspond to the first dimension, and *columns* correspond to the second dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "> So far an 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",
    "> * 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. Numpy arrays are **memory efficient**.\n",
    ">\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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fancy Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Fancy indexing is the NumPy name for when an integer array or list is used in-place of an index:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([n for n in range(5)])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3 1 2 3]\n"
     ]
    }
   ],
   "source": [
    "row_indices = [3, 1, 2, 3]\n",
    "print(a[row_indices])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "You can also use negative indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 4 3 2]\n"
     ]
    }
   ],
   "source": [
    "row_indices = [-3, -1, -2, -3]\n",
    "print(a[row_indices])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Unlike slicing, selecting data from an array by **fancy indexing always creates a copy of the data** in the new array.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Boolean Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Instead of using integers or arrays/lists of integers to access the element of a NumPy array, we can also use arrays of boolean, called **index masks**. An index masks is a list/array of boolean dat type; when used to index a NumPy array, an element is selected or not depending on the value of the index mask at the position of each element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([n for n in range(5)])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2]\n"
     ]
    }
   ],
   "source": [
    "index_mask = np.array([True, False, True, False, False])\n",
    "print(a[index_mask])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2]\n"
     ]
    }
   ],
   "source": [
    "# same thing\n",
    "index_mask = np.array([1, 0, 1, 0, 0], dtype=bool)\n",
    "print(a[index_mask])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Boolean indexing is very useful to conditionally select elements from an array, using for example comparison operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Alice' 'Charlie' 'David']\n"
     ]
    }
   ],
   "source": [
    "names = np.array(['Alice', 'Bob', 'Bob', 'Charlie', 'Bob', 'David'])\n",
    "index_mask = names != 'Bob'\n",
    "print(names[index_mask])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "It is also possible to mix boolean array using logical conditions. However, the Python keywords `not`, `and` and `or` do not work with boolean arrays. Use `~`, `&` and `|` instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Charlie' 'David']\n"
     ]
    }
   ],
   "source": [
    "print(names[~(names == 'Bob') & (names != 'Alice')])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Mixing up boolean indexing and broadcasting allows us many common-sense applications. For example, to set all negative values in an array to 0 we need only to write:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1  2  3 -2 -3  4  0]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([-1, 2, 3, -2, -3, 4, 0])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2 3 0 0 4 0]\n"
     ]
    }
   ],
   "source": [
    "a[a < 0] = 0\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Unlike slicing, selecting data from an array by **boolean indexing always creates a copy of the data** in the new array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Copies versus Views\n",
    "\n",
    "If we want to obtain a copy of the data using slicing too, 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 `np.copy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3, 4, 5])\n",
    "b = np.copy(a[:3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0]\n"
     ]
    }
   ],
   "source": [
    "b[:] = 0\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4 5]\n"
     ]
    }
   ],
   "source": [
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Array Functions\n",
    "\n",
    "There are *many* [array functions](https://docs.scipy.org/doc/numpy/reference/routines.html) we can use to compute with NumPy arrays. The following is a partial list and we'll look closer at mathematical functions in the next section.\n",
    "\n",
    "| Array Function | Description |\n",
    "| ---: | :--- |\n",
    "| `np.sum` | Sum of all the elements in the array. |\n",
    "| `np.prod` | Product of all the elements in the array.|\n",
    "| `np.mean` | Arithmetic mean.|\n",
    "| `np.std`, `np.var` | Standard deviation and variance. |\n",
    "| `np.min`, `np.max` | Minimum and maximum. |\n",
    "| `np.argmin`, `np.argmax` | Indices of the minimum and maximum elements. |\n",
    "| `np.cumsum` | Cumulative sum of elements, starting from 0. |\n",
    "| `np.cumprod` | Cumulative product of elements, starting from 1. |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create a 1D NumPy array with random values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 8 -2  4  7 -3]\n"
     ]
    }
   ],
   "source": [
    "arr = np.array([8, -2, 4, 7, -3])\n",
    "print(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Compute the mean of the values in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.8"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Verify the mean once more:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.8\n"
     ]
    }
   ],
   "source": [
    "m = np.sum(arr) / arr.size\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Find the index of the maximum element in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "arr = np.array([8, -2, 4, 7, -3])\n",
    "max_i = np.argmax(arr)\n",
    "print(max_i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Verify the maximum value in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr[max_i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Array functions apply to 2D arrays as well (and $n$-dimensional arrays in general) with the added feature that we can choose to apply array functions to the entire array, down the columns or across the rows (or any axis).\n",
    "\n",
    "Create a 2D NumPy array with random values and compute the sum of all the entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 4 2]\n",
      " [2 1 1]\n",
      " [3 2 0]\n",
      " [0 6 2]]\n"
     ]
    }
   ],
   "source": [
    "M = np.array([[2, 4, 2], [2, 1, 1], [3, 2, 0], [0, 6, 2]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The function `numpy.sum` also takes a keyword argument `axis` which determines along which dimension to compute the sum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7, 13,  5])"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(M, axis=0) # Sum of the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8, 4, 5, 8])"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(M, axis=1) # Sum of the rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Mathematical Functions\n",
    "\n",
    "[Mathematical functions](http://docs.scipy.org/doc/numpy/reference/routines.math.html) in NumPy are called [**universal functions**](https://docs.scipy.org/doc/numpy/user/quickstart.html#universal-functions) and are *vectorized*. Vectorized functions operate *element-wise* on arrays producing arrays as output and are built to compute values across arrays *very* quickly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The following table contains a list of the most important unary ufuncs.\n",
    "\n",
    "|Function| Description |\n",
    "|:-------|:---------|\n",
    "|`np.abs`|Compute the absolute value element-wise for integer, floating-point, or complex values|\n",
    "|`np.sqrt`|Compute the square root of each element|\n",
    "|`np.exp`|Compute the exponent $e^x$ of each element|\n",
    "|`np.log`, `np.log10`, `np.log2`, `np.log1p`|Natural logarithm (base e), log base 10, log base 2, and log(1 + x), respectively|\n",
    "|`np.sign`|Compute the sign of each element: 1 (positive), 0 (zero), or –1 (negative)|\n",
    "|`np.ceil`|Compute the ceiling of each element|\n",
    "|`np.floor`|Compute the floor of each element|\n",
    "|`np.modf`|Return fractional and integral parts of array as a separate array|\n",
    "|`np.isnan`|Return boolean array indicating whether each value is `NaN` (Not a Number)|\n",
    "|`np.cos`, `np.cosh`, `np.sin`, `np.sinh`, `np.tan`, `np.tanh`|Regular and hyperbolic trigonometric functions|\n",
    "|`np.arccos`, `np.arccosh`, `np.arcsin`, `np.arcsinh`, `np.arctan`, `np.arctanh`|Inverse trigonometric functions|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The following table contains a list of the most important binary ufuncs.\n",
    "\n",
    "|Function| Description |\n",
    "|:-------|:---------|\n",
    "|`np.add`|Element-wise addition|\n",
    "|`np.subtract`|Element-wise subtraction|\n",
    "|`np.multiply`|Element-wise multiplication|\n",
    "|`np.divide`|Element-wise division|\n",
    "|`np.mod`|Element-wise modulus|\n",
    "|`np.power`|Raise elements in first array to powers indicated in second array |\n",
    "|`np.maximum`, `np.fmax`|Element-wise maximum; `np.fmax` ignores `NaN`|\n",
    "|`np.minimum`, `np.fmin`|Element-wise minimum; `np.fmin` ignores `NaN`|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Compute the values $\\sin(2 \\pi x)$ for $x = 0, 0.25, 0.5, \\dots, 1.75$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.   0.25 0.5  0.75 1.  ]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(0, 1.25, 0.25)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.0000000e+00,  1.0000000e+00,  1.2246468e-16, -1.0000000e+00,\n",
       "       -2.4492936e-16])"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sin(2 * np.pi * x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We expect the array `[0. 1. 0. -1. 0.]` however there is (as always with floating point numbers) some rounding errors in the result. In numerical computing, we can interpret a number such as $10^{-16}$ as $0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Compute the values $\\log_{10}(x)$ for $x = 1,10,100,1000,10000$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[    1    10   100  1000 10000]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1, 10, 100, 1000, 10000])\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1., 2., 3., 4.])"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log10(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that we can also evaluate mathematical functions with scalar values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sin(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "NumPy also provides familiar mathematical constants such as $\\pi$ and $e$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589793"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.718281828459045"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For example, verify the limit\n",
    "\n",
    "$$\n",
    "\\lim_{x \\to \\infty} \\arctan(x) = \\frac{\\pi}{2}\n",
    "$$\n",
    "\n",
    "by evaluating $\\arctan(x)$ for some (arbitrary) large value $x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5706963267952299"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arctan(10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5707963267948966"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.pi/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Ternary Expressions\n",
    "\n",
    "The `np.where`function is a vectorized version of the ternary expression `x if <condition> else y`. Suppose we have the following arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])\n",
    "y = np.array([0.1, 0.2, 0.3, 0.4, 0.5])\n",
    "condition = [False, True, False, False, True]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We want to take a value from `x` whenever the corresponding condition in `condition` is `True`, and otherwise take the value from `y`. With `np.where` we can write:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.1 2.  0.3 0.4 5. ]\n"
     ]
    }
   ],
   "source": [
    "result = np.where(condition, x, y)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sorting\n",
    "\n",
    "Like Python's lists, NumPy array can be sorted in-place:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2, 1, 2, 4, 7, 88]\n"
     ]
    }
   ],
   "source": [
    "a = [4, 2, 7, 88, -2, 1]\n",
    "a.sort()\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2, 1, 2, 4, 7, 88]\n"
     ]
    }
   ],
   "source": [
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Using the function `np.sort` the input array is not modified and a new, sorted array is created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2  1  2  4  7 88]\n",
      "[4, 2, 7, 88, -2, 1]\n"
     ]
    }
   ],
   "source": [
    "a = [4, 2, 7, 88, -2, 1]\n",
    "print(np.sort(a))\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Special Methods for Boolean Arrays\n",
    "\n",
    "There are two additional methods for boolean arrays:\n",
    "* `any` to check if at least one value in the array is `True`;\n",
    "* `all` to test if all values in the array are `True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "a = np.array([False, False, True, False])\n",
    "print(a.any())\n",
    "print(a.all())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Set Functions\n",
    "\n",
    "NumPy has some basic set operations working on 1D arrays only:\n",
    "* `np.unique` returns the sorted unqiue values in an array;\n",
    "* `np.in1d` tests membership of the values in one array in another, returning a boolean array.\n",
    "\n",
    "\n",
    "|Function| Description |\n",
    "|:-------|:---------|\n",
    "|`np.unique`|Compute the sorted, unique elements in an array.|\n",
    "|`np.intersect1d`|Compute the sorted, common elements in two arrays.|\n",
    "|`np.union1d`|Compute the sorted union elements in two arrays.|\n",
    "|`np.in1d`|Test membership of the values in one array in another, returning a boolean array.|\n",
    "|`np.setdiff1d`|Set difference, return elements in and array that are not in the other.|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Alice' 'Bob' 'Charlie' 'David']\n"
     ]
    }
   ],
   "source": [
    "names = np.array(['Alice', 'Bob', 'Bob', 'Charlie', 'Bob', 'David'])\n",
    "print(np.unique(names))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ True False False  True  True False  True]\n"
     ]
    }
   ],
   "source": [
    "values = np.array([6, 0, 0, 3, 2, 5, 6])\n",
    "print(np.in1d(values, [2, 3, 6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Reshaping, resizing and stacking\n",
    "\n",
    "The shape of an Numpy array can be modified **without copying the underlaying data**, i.e., returing a view, which makes it a fast operation even for large arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "M = np.array([[1, 2, 3], [4, 5, 6]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3 4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "m, n = M.shape\n",
    "B = M.reshape((1, n * m))\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": 118,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4 5 6]\n"
     ]
    }
   ],
   "source": [
    "B = M.flatten()\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Adding a new dimension\n",
    "\n",
    "With `np.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": 119,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,)\n"
     ]
    }
   ],
   "source": [
    "v = np.array([1, 2, 3])\n",
    "print(v.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 1)\n"
     ]
    }
   ],
   "source": [
    "# make a column matrix of the vector v\n",
    "w = v[:, np.newaxis]\n",
    "print(w.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 3)\n"
     ]
    }
   ],
   "source": [
    "# make a row matrix of the vector v\n",
    "w = v[np.newaxis, :]\n",
    "print(w.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Stacking\n",
    "\n",
    "We can build bigger arrays out of smaller arrays by [stacking](https://docs.scipy.org/doc/numpy/user/quickstart.html#stacking-together-different-arrays) along different dimensions using the functions `numpy.hstack` and `numpy.vstack`.\n",
    "\n",
    "Stack 3 different 1D NumPy arrays of length 3 vertically forming a 3 by 3 matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([1, 1, 1])\n",
    "y = np.array([2, 2, 2])\n",
    "z = np.array([3, 3, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 1 1]\n",
      " [2 2 2]\n",
      " [3 3 3]]\n"
     ]
    }
   ],
   "source": [
    "vstacked = np.vstack((x,y,z))\n",
    "print(vstacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Stack 1D NumPy arrays horizontally to create another 1D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1 1 2 2 2 3 3 3]\n"
     ]
    }
   ],
   "source": [
    "hstacked = np.hstack((x,y,z))\n",
    "print(hstacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Use `np.hstack` and `np.vstack` to build the matrix $T$ where\n",
    "\n",
    "$$\n",
    "T = \n",
    "\\begin{bmatrix}\n",
    "1 & 1 & 2 & 2 \\\\ \n",
    "1 & 1 & 2 & 2 \\\\ \n",
    "3 & 3 & 4 & 4 \\\\ \n",
    "3 & 3 & 4 & 4\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 2. 2.]\n",
      " [1. 1. 2. 2.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.ones((2, 2))\n",
    "B = 2 * np.ones((2, 2))\n",
    "C = 3 * np.ones((2, 2))\n",
    "D = 4 * np.ones((2, 2))\n",
    "A_B = np.hstack((A, B))\n",
    "print(A_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3. 3. 4. 4.]\n",
      " [3. 3. 4. 4.]]\n"
     ]
    }
   ],
   "source": [
    "C_D = np.hstack((C, D))\n",
    "print(C_D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 2. 2.]\n",
      " [1. 1. 2. 2.]\n",
      " [3. 3. 4. 4.]\n",
      " [3. 3. 4. 4.]]\n"
     ]
    }
   ],
   "source": [
    "T = np.vstack((A_B, C_D))\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## File I/O with Arrays\n",
    "\n",
    "The `np.save` and `np.load` are the two functions for efficiently saving and loading array data from and to disk. Arrays are saved by default in an uncompressed raw binary format with file extension `.npy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4 5 6 7 8 9]\n"
     ]
    }
   ],
   "source": [
    "a = np.arange(10)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If your save filename does not end with `.npy`, the extensions will be automatically added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "np.save('test_a', a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "To correctly load a saved file, the filename **must include the `.npy` extensions**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4 5 6 7 8 9]\n"
     ]
    }
   ],
   "source": [
    "b = np.load('test_a.npy')\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Random Number Generators\n",
    "\n",
    "The subpackage `numpy.random` contains functions to generate NumPy arrays of [random numbers](https://docs.scipy.org/doc/numpy/reference/routines.random.html) sampled from different distributions. The following is a partial list of distributions:\n",
    "\n",
    "| Function | Description |\n",
    "| :--- | :--- |\n",
    "| `numpy.random.rand(d1,...,dn)` | Create a NumPy array (with shape `(d1,...,dn)`) with entries sampled uniformly from `[0,1)` |\n",
    "| `numpy.random.randn(d1,...,dn)` | Create a NumPy array (with shape `(d1,...,dn)`) with entries sampled from the standard normal distribution |\n",
    "| `numpy.random.randint(a,b,size)` | Create a NumPy array (with shape `size`) with integer entries from `low` (inclusive) to `high` (exclusive) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Sample a random number from the [uniform distribution](https://en.wikipedia.org/wiki/Uniform_distribution_%28continuous%29):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6584711063945924"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.rand()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Sample 3 random numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.1135212 , 0.5790922 , 0.53603493])"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.rand(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Create 2D NumPy array of random samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.10372371, 0.99003381, 0.90302009, 0.3193851 ],\n",
       "       [0.08293488, 0.80402934, 0.74208849, 0.61407595]])"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.rand(2,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Random samples from the [standard normal distribution](https://en.wikipedia.org/wiki/Normal_distribution):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.14613974413530106"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.5368617 ,  0.5607    ,  0.46728486])"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randn(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.46546945],\n",
       "       [-0.54656086],\n",
       "       [-0.29600639]])"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randn(3,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Random integers sampled uniformly from various intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randint(-10,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 1, 1, 0, 1, 1, 0],\n",
       "       [0, 0, 0, 0, 1, 1, 0, 0],\n",
       "       [0, 1, 1, 1, 1, 0, 0, 0],\n",
       "       [0, 1, 0, 0, 0, 1, 0, 1]])"
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randint(0,2,(4,8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4,  5],\n",
       "       [ 6,  3],\n",
       "       [-2, -6],\n",
       "       [-2,  8],\n",
       "       [ 1,  1]])"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randint(-9,10,(5,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Iterating over arrays\n",
    "\n",
    "In general, **avoid iterating over the elements of arrays whenever possible (at all costs)**. The reason is that in a interpreted language like Python iterations are really slow compared to vectorized operations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercises (solved)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Brute Force Optimization\n",
    "\n",
    "Find the [absolute maximum and minimum](https://en.wikipedia.org/wiki/Maxima_and_minima) values of the function\n",
    "\n",
    "$$\n",
    "f(x) = x\\sin(x)+\\cos(4x)\n",
    "$$\n",
    "\n",
    "on the interval $[0,2\\pi]$. We know that the maximum and minimum values must occur at either the endpoints $x=0,2\\pi$ or at critical points where $f'(x)=0$. However, the derivative is given by\n",
    "\n",
    "$$\n",
    "f'(x) = \\sin(x) + x \\cos(x) - 4\\sin(4x)\n",
    "$$\n",
    "\n",
    "and the equation $f'(x) = 0$ is impossible to solve explicitly.\n",
    "\n",
    "Instead, create a 1D NumPy array of $x$ values from $0$ to $2\\pi$ of length $N$ (for some arbitrarily large value $N$) and use the functions `numpy.min` and `numpy.max` to find maximum and minimum $y$ values, and the functions `numpy.argmin` and `numpy.argmax` to find the indices of the corresponding $x$ values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Absolute maximum value is y = 2.5992726072887007 at x = 1.628136126702901\n",
      "Absolute minimum value is y = -5.129752039182 at x = 5.34187001663503\n"
     ]
    }
   ],
   "source": [
    "N = 10000\n",
    "x = np.linspace(0, 2*np.pi, N)\n",
    "y = x * np.sin(x) + np.cos(4*x)\n",
    "y_max = np.max(y)\n",
    "y_min = np.min(y)\n",
    "x_max = x[np.argmax(y)]\n",
    "x_min = x[np.argmin(y)]\n",
    "print('Absolute maximum value is y =', y_max, 'at x =', x_max)\n",
    "print('Absolute minimum value is y =', y_min, 'at x =', x_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Riemann Sums\n",
    "\n",
    "Write a function called `exp_int` which takes input parameters $b$ and $N$ and returns the (left) [Riemann sum](https://en.wikipedia.org/wiki/Riemann_sum)\n",
    "\n",
    "$$\n",
    "\\int_0^b e^{-x^2} dx \\approx \\sum_{k=0}^{N-1} e^{-x_k^2} \\Delta x\n",
    "$$\n",
    "\n",
    "for $\\Delta x = b/N$ and the partition $x_k=k \\, \\Delta x$, $k=0,\\dots,N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def exp_int(b,N):\n",
    "    \"Compute left Riemann sum of exp(-x^2) from 0 to b with N subintervals.\"\n",
    "    x = np.linspace(0,b,N+1)\n",
    "    x_left_endpoints = x[:-1]\n",
    "    Delta_x = b/N\n",
    "    I = Delta_x * np.sum(np.exp(-x_left_endpoints**2))\n",
    "    return I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The infinite integral satisfies the [beautiful identity](https://en.wikipedia.org/wiki/Gaussian_integral)\n",
    "\n",
    "$$\n",
    "\\int_0^{\\infty} e^{-x^2} dx = \\frac{\\sqrt{\\pi}}{2}\n",
    "$$\n",
    "\n",
    "Compute the integral with large values of $b$ and $N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.886726925452758"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp_int(100,100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Compare to the true value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8862269254527579"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.pi**0.5/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Infinite Products\n",
    "\n",
    "The cosine function has the following infinite product representation\n",
    "\n",
    "$$\n",
    "\\cos x = \\prod_{k = 1}^{\\infty} \\left(1 - \\frac{4 x^2}{\\pi^2 (2k - 1)^2} \\right)\n",
    "$$\n",
    "\n",
    "Write a function called `cos_product` which takes input parameters $x$ and $N$ and returns the $N$th partial product\n",
    "\n",
    "$$\n",
    "\\prod_{k = 1}^{N} \\left(1 - \\frac{4 x^2}{\\pi^2 (2k - 1)^2} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def cos_product(x,N):\n",
    "    \"Compute the product \\prod_{k=1}^N (1 - 4x^2/(pi^2 (2k - 1)^2).\"\n",
    "    k = np.arange(1,N+1)\n",
    "    terms = 1 - 4*x**2 / (np.pi**2 * (2*k - 1)**2)\n",
    "    return np.prod(terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Verify our function using values for which we know the result. For example, $\\cos(0)=1$, $\\cos(\\pi)=-1$ and $\\cos(\\pi/4) = \\frac{1}{\\sqrt{2}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos_product(0,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0001000050002433"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos_product(np.pi,10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7071067856245614"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos_product(np.pi/4,10000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7071067811865475"
      ]
     },
     "execution_count": 148,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/2**0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercises (unsolved)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "*  The natural log satisfies the following definite integral\n",
    "\n",
    "    $$\n",
    "    \\int_1^e \\frac{\\ln x \\ dx}{(1 + \\ln x)^2} = \\frac{e}{2} - 1\n",
    "    $$\n",
    "\n",
    "    Write a function called `log_integral` which takes input parameters $c$ and $N$ and returns the value of the (right) Riemann sum\n",
    "\n",
    "    $$\n",
    "    \\int_1^c \\frac{\\ln x \\ dx}{(1 + \\ln x)^2} \\approx \\sum_{k=1}^N \\frac{\\ln x_k \\ \\Delta x}{(1 + \\ln x_k)^2} \\ , \\ \\ \\Delta x = \\frac{c - 1}{N}\n",
    "    $$\n",
    "    \n",
    "    for the partition $x_k = 1 + k \\Delta x$, for $k = 0, \\dots , N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Write a function called `k_sum` which takes input parameters `k` and `N` and returns the partial sum\n",
    "\n",
    "    $$\n",
    "    \\sum_{n=1}^{N} \\frac{1}{n^k}\n",
    "    $$\n",
    "    \n",
    "    Verify your function by comparing to the infinite series identity\n",
    "\n",
    "    $$\n",
    "    \\sum_{n=1}^{\\infty} \\frac{(-1)^{n+1}}{n^2} = \\frac{\\pi^2}{12}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Write a function called `dot` which takes 3 inputs `M`, `i` and `j` where `M` is a square NumPy array and the function returns the dot product of the $i$th row and the $j$th column of $M$."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}