lectures/L02_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arrays: numpy\n",
    "\n",
    "*Davide Gerosa (Milano-Bicocca)*\n",
    "\n",
    "**Sources**: Michael Zingale at Stony Brook University: https://sbu-python-class.github.io"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numpy library provides a class for n-dimensional arrays of data. This\n",
    "notebook is based on the numpy tutorial from the scipy webpage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "Note that the traditional way to import numpy is to rename it np.  This saves on typing and makes your code a little more compact.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy provides a multidimensional array class as well as a _large_ number of functions that operate on arrays.\n",
    "\n",
    "Numpy arrays allow you to write fast (optimized) code that works on arrays of data.  To do this, there are some restrictions on arrays:\n",
    "\n",
    "* all elements are of the same data type (e.g. float)\n",
    "* the size of the array is fixed in memory, and specified when you create the array (e.g., you cannot grow the array like you do with lists)\n",
    "\n",
    "\n",
    "The nice part is that arithmetic operations work on entire arrays&mdash;this means that you can avoid writing loops in python (which tend to be slow).  Instead the \"looping\" is done in the underlying compiled code. \n",
    "\n",
    "The usual criticism against python is that it's slow (because it's an interpreted language, not a compiled one like C). Well, that's not true; numpy achieves C-level speed. \n",
    "\n",
    "Always remember:\n",
    "\n",
    "`Time to paper` = `time you need to write code` + `time the code runs`. \n",
    "\n",
    "No point wasting days debugging code to gain 1ms at run time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array Creation and Properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a lot of ways to create arrays.  Let's look at a few\n",
    "\n",
    "Here we create an array using `arange` and then change its shape to be 3 rows and 5 columns.  Note the _row-major ordering_&mdash;you'll see that the rows are together in the inner `[]` (more on this in a bit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15,)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3  4]\n",
      " [ 5  6  7  8  9]\n",
      " [10 11 12 13 14]]\n"
     ]
    }
   ],
   "source": [
    "a = np.arange(15).reshape(3,5)\n",
    "\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 5)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An array is an object of the `ndarray` class, and it has methods that know how to work on the array data.  Here, `.reshape()` is one of the many methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A numpy array has a lot of meta-data associated with it describing its shape, datatype, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "(3, 5)\n",
      "15\n",
      "int64\n",
      "8\n",
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "print(a.ndim)\n",
    "print(a.shape)\n",
    "print(a.size)\n",
    "print(a.dtype)\n",
    "print(a.itemsize)\n",
    "print(type(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#help(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also create an array from a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.  2.1 3.  4. ]\n",
      "float64\n",
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "b = np.array( [1, 2.1, 3, 4] )\n",
    "print(b)\n",
    "print(b.dtype)\n",
    "print(type(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a multi-dimensional array of a specified size initialized all to 0 easily, using the `zeros()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.zeros((10,8),dtype=int)\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also an analogous `ones()` and `empty()` array routine.  Note that here we can explicitly set the datatype for the array in this function if we wish. \n",
    "\n",
    "Unlike lists in python, all of the elements of a numpy array are of the same datatype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[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.]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = np.eye(10, dtype=np.float64)\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`linspace` creates an array with evenly space numbers.  The `endpoint` optional argument specifies whether the upper range is in the array or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.         0.11111111 0.22222222 0.33333333 0.44444444 0.55555556\n",
      " 0.66666667 0.77777778 0.88888889 1.        ]\n"
     ]
    }
   ],
   "source": [
    "d = np.linspace(0, 1, 10, endpoint=True)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<span class=\"fa fa-flash\"></span> Quick Exercise:\n",
    "\n",
    "Analogous to `linspace()`, there is a `logspace()` function that creates an array with elements equally spaced in log. Use `help(np.logspace)` to see the arguments, and create an array with 10 elements from $10^{-6}$ to $10^3$. Then do it again with `geomspace` and figure out the difference.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.00000000e-06 7.94328235e-06 6.30957344e-05 5.01187234e-04\n",
      " 3.98107171e-03 3.16227766e-02 2.51188643e-01 1.99526231e+00\n",
      " 1.58489319e+01 1.25892541e+02]\n"
     ]
    }
   ],
   "source": [
    "#help(np.logspace)\n",
    "e = np.logspace(-6, 3, num=10, endpoint=False)\n",
    "print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.00000000e-06 7.94328235e-06 6.30957344e-05 5.01187234e-04\n",
      " 3.98107171e-03 3.16227766e-02 2.51188643e-01 1.99526231e+00\n",
      " 1.58489319e+01 1.25892541e+02]\n"
     ]
    }
   ],
   "source": [
    "#help(np.geomspace)\n",
    "f = np.logspace(-6, 3, num=10, endpoint=False)\n",
    "print(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array Operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most operations (`+`, `-`, `*`, `/`) will work on an entire array at once, element-by-element.\n",
    "\n",
    "Note that that the multiplication operator is not a matrix multiply (there is a new operator in python 3.5+, `@`, to do matrix multiplicaiton.\n",
    "\n",
    "Let's create a simple array to start with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "l=[1,2,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 2, 3, 1, 2, 3]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 4, 6])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(l)*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n"
     ]
    }
   ],
   "source": [
    "a = np.arange(12).reshape(3,4)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplication by a scalar multiplies every element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  2,  4,  6],\n",
       "       [ 8, 10, 12, 14],\n",
       "       [16, 18, 20, 22]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a*2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding two arrays adds element-by-element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  2,  4,  6],\n",
       "       [ 8, 10, 12, 14],\n",
       "       [16, 18, 20, 22]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplying two arrays multiplies element-by-element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0,   1,   4,   9],\n",
       "       [ 16,  25,  36,  49],\n",
       "       [ 64,  81, 100, 121]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a*a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<span class=\"fa fa-flash\"></span> Quick Exercise:\n",
    "\n",
    "What do you think `1./a` will do?  Try it and see\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/71/dmgz5bhj609bj3ygp2xzd_1r0000gn/T/ipykernel_7882/2189751180.py:1: RuntimeWarning: divide by zero encountered in divide\n",
      "  1./a\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[       inf, 1.        , 0.5       , 0.33333333],\n",
       "       [0.25      , 0.2       , 0.16666667, 0.14285714],\n",
       "       [0.125     , 0.11111111, 0.1       , 0.09090909]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1./a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can think of our 2-d array as a 3 x 4 matrix (3 rows, 4 columns).  We can take the transpose to get a 4 x 3 matrix, and then we can do a matrix multiplication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do some basic linear algebra operations on arrays.  The operator `@` is a matrix multiplication / dot-product operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  4,  8],\n",
       "       [ 1,  5,  9],\n",
       "       [ 2,  6, 10],\n",
       "       [ 3,  7, 11]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = a.transpose() ## Also a.T\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 14,  38,  62],\n",
       "       [ 38, 126, 214],\n",
       "       [ 62, 214, 366]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a @ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can sum the elements of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "66"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 6, 22, 38])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<span class=\"fa fa-flash\"></span> Quick Exercise:\n",
    "\n",
    "`sum()` takes an optional argument, `axis=N`, where `N` is the axis to sum over.  Sum the elements of `a` across rows to create an array with just the sum along each column of `a`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also easily get the extrema"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 11\n"
     ]
    }
   ],
   "source": [
    "print(a.min(), a.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Universal functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Universal functions work element-by-element.  Let's create a new array scaled by `pi / 12`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.26179939 0.52359878 0.78539816]\n",
      " [1.04719755 1.30899694 1.57079633 1.83259571]\n",
      " [2.0943951  2.35619449 2.61799388 2.87979327]]\n"
     ]
    }
   ],
   "source": [
    "b = a*np.pi/12.0\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00  9.65925826e-01  8.66025404e-01  7.07106781e-01]\n",
      " [ 5.00000000e-01  2.58819045e-01  6.12323400e-17 -2.58819045e-01]\n",
      " [-5.00000000e-01 -7.07106781e-01 -8.66025404e-01 -9.65925826e-01]]\n"
     ]
    }
   ],
   "source": [
    "c = np.cos(b)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = b + c "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.         1.22772521 1.38962418 1.49250494]\n",
      " [1.54719755 1.56781598 1.57079633 1.57377667]\n",
      " [1.5943951  1.64908771 1.75196847 1.91386744]]\n"
     ]
    }
   ],
   "source": [
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<span class=\"fa fa-flash\"></span> Quick Exercise:\n",
    "\n",
    "We will often want to write our own function that operates on an array and returns a new array.  We can do this just like we did with functions previously&mdash;the key is to use the methods from the `np` module to do any operations, since they work on, and return, arrays.\n",
    "\n",
    "Write a simple function that returns $\\sin(2\\pi x)$ for an input array `x`.  Then test it \n",
    "by passing in an array `x` that you create via `linspace()`\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slicing"
   ]
  },
  {
   "attachments": {
    "slicing.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAA84AAADzCAYAAABXGU+2AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAAsTAAALEwEAmpwYAAAAB3RJTUUH4QIMEzMzdAGXUAAAIABJREFUeNrt3Xl4FFXWx/FfJwQCIWwhAQ0KCoqKssoOgoiIbIKAoAODqIjiuAu4gaOIII7jwriNDiOCIzjgvMOooMLIKiCMgsqmKCqLLGGHBBKS8/5x6XSHbJ2kO3TC9/M8/aTSVb3UqeqqPn3rnusxMxMAAAAAAMhRBCEAAAAAAIDEGQAAAAAAEmcAAAAAAEicAQAAAAAgcQYAAAAAgMQZAAAAAAASZwAAAAAASJxLkUWLFhEE4k/8QfyJP4g/8QfxJ/4g/pk8ZmZsLgAAAAAAckaLMwAAAAAAJM4AAAAAAJA4AwAAAABA4gwAAAAAAIkzAAAAAAAkzgAAAAAAkDgDAAAAAEDiXJowADrxJ/4g/sQfxJ/4g/gTfxB/fx4zMzYXAAAAAAA5o8UZAAAAAAASZwAAAAAASJwBAAAAACBxBgAAAACAxBkAAAAAABJnAAAAAABInAEAAAAAIHEuTRgAnfgTfxB/4g/iT/xB/Ik/iH+2xDmQNxvoCgVzueJ+rkWLFqljx45Bfb5wjAXxJ/7En/gTf+JP/Ik/8Sf+xJ/4E//Al/GYmfE7h19APB4REuJP/EH8iT+IP/EH8Sf+IP5eXKoNAAAAAACJMwAAAAAAJM4AAAAAAJA4AwAAAABA4gwAAAAAAIkzAAAAAAAkzgAAAAAAkDgDAAAAAEDiDAAAAAAAiTMAAAAAACBxBgAAAACAxBkAAAAAABJnAAAAAABInAEAAAAAIHEGAAAAAIDEGQAAAAAAEmcAAAAAAEicAQAAAAAgcQYAAAAAgMQZAAAAAACQOAMAAAAAQOIMAAAAAACJMwAAAAAAJM4AAAAAAJA4AwAAAABA4gwAAAAAAIkzAAAAAAAkzgAAAAAAkDijZDpxQnrnHalXL+ncc6XoaCk2VrrwQun3v5f+9S8pPZ04FaeFCyWPx90QWuvXSw89JLVsKcXFSVFRUtWqUqNG0j33SOvWEaNQ2rdP+tOfpCuvlOLjpbJlpRo1pMsvlx5/3G0fnB7TpnEcCpU9e3yxLcgNofkONH261K+fdN55UkyMu9WvL912m7R6NTEKtsLs+3wGgi8jQ/r3v6X+/aXataXy5d2+X6+eNHiw9N//EqNg7vZmZoTB/zjgUYkLyebNUu/e+ScHDRpIs2ZJF11E/ENt61apVStpxw73fwlZpxIX//R06cEHpcmT3ckjrxP8Qw9JEyZIkZHEP5hmzZKGD3fJc24iIqTbb5f+/Gd3Uif+xeOHH6SmTaUjR0rEcajExX/hQvdjUUGF6TqW2P1/9WqXIGzcmPc54O67pRdecMcj4h+cxLmwwnA9S+T+v2uXNGCAtGhR3sv16SO9/bZUqRLH/yKixbmkO3RIuuaawFrU1q2Trr7aPQahs2WL1KmTL2lG6Nx6q/TSS3knzd6T9HPPSXfdRcyC6d133Uk7r6RZctvn9delvn1dyxBC7/hxt228STOCjytZTr9Fi6SOHfNOmr3ngJdflu69l5idbhdcQAyC4dgxqUuX/JNmyV112rMn518SZ+jpp6WffvL9/7vfSV984b4sHTokLVkiXX+9b/62bdJrrxG3UPnPf6QWLdxVAAitTz+Vpk71/X/ttdKcOdIvv7ikYdcul9j5n6TfeMM9DkW3bZtrafb+aNGggfTmm9KPP0opKdJvv7lY9+jhe8zcudKUKcSuOIwaJX39NXEorsT5pZdcchbIDcGxf7/7znP0qPu/SRN3DNq8WUpOln79VZo923XZ8frLX6Tly4ldMAS6v5tJkya5x1SrJn38MbELhhdekL75xk2XLevygR9+cN9/jh1zPyaNGSOVKeOWWbzY/YCNou738FeiQpKcbFa5su/Q9MoruS87apRvuXPPNcvIIP7BtGeP2R135H7aYP8Pvuuv98V31Kjclzt40KxhQ9+y111H/IPhoYd8MW3Xzuzw4dyXffhh37JNmxL/UJszp0Qeh0pc/Dt08MV2yRK+/xS3u+/2xb9+fbN9+3JeLi3N7MorfcsOGkT8i9M//uHi7vGYffQR+3+wXHqpb59+/vnclxs3zrdcs2bEv4jo43yKEtXH4cMP3aUXknTJJdJ33+Xe5yQtzRXJ2LLF/b9ihSukRPyLbuxY98uf95LIyEj3y98jj2T9ZZb9P7hq1nStypKUlOSKguVmwQKpc2c3nZDgexzxL7y6dX1Xuyxf7vr059UyVK2am46JCdvLh0tFH+ft210L29697jOyc2eJOQ6VuPjHx7tjT0SEdPCgVLEi33+Ky7Fj0llnSQcOuO89ixdL7drlvvyqVe5qMElKTHRXzBD/0Fu7Vmrd2l2F9MAD0vPPs/8HS/ny7nMgSbt3u+NRbueEWrXcdIUKvis0iH+hcKl2SbZ4sW+6d++8CzVERblqk15cKhM848b5EoHLL3dJxMMPE5dQ27nT15aWV9IsSW3b+qYPHCB2wfDjj77455U0S+5Lk1f16sQuVNLTpZtucklzRISrMozQ2LPHJc2S6w5SwpPmEuejj3zH8m7d8k6aJal5c9/xKgyT5lIpOdnVWUhJcQ0348cTk2CqUsU3nVff5bJlfdOxscSNxDmMjBzpWhkXLPD9ChRK337rm27dOv/l27f3Ta9axfYKposucv1pV650J2iE3wncq2ZN4lGcMjJcPyuv664jJqEybpzvB9WHH5auuoqYhIp//+ZmzVxNkeefdz/SVaniWoPq1JEGDnRDxXBxX3CtXOmb7tuXeISj0aOlTZvc9GuvuWFSETz+DQL/+Efuy82Y4Ztu04a4FVEZQhBEv/4qvf++NHGiO2m2besuD+3c2RWtCPYQCD//7Js+//z8lz/vPN/0jz+yvYLls8/cF1TGJgxfS5b4pjt2JB7F8UPFtm3SvHmuGNjate7+unWlxx4jPqGweLHrIuL9cvTkk8SkuBLn3btdq/Pu3VmX+eUXd5s50/24PWOGdO65xC4YvvrKN920qfu7fr0rDjZvnvs+FhXlxrLt2dMNReXtLoLQW7FCevVVNz1wYOGGbUPennjCFdxMTnbn1RMn3I9I557rvo9u2yZ98IHvh+vy5d1jUCT0cT41IEW5xn74cOmvf815XvXq7sDRubMbEso/iS2satVc30HJDQdTtWrey+/d67tMsmrV/IeQKWnxD7+V8U3Tx/n06tjRN2TD0qVZf6kl/sHVubO76uZU3bq5ip7nnMP+H2x797p+zdu3u2P7mjW+BK0EHYdKVPxHjCj4CBU1ariuPME4/5/p8a9Xz9cAkJQkvfKK++EoLS3n5atUcV0Xundn/y8OLVtKX37pLhPeuDFs9/kSH//ly6VBg7KOrpOTc85x+/8VVxD/IuJS7WCqXDn3eUlJ0j//6ZLr8893LS/Dh7v7CpvA+nfwL18+/+X9l2FsT5wp/vpXX9I8cGDYJs2lhv+VMJlnmgg3JMYvvxCfUBg61CXNkvTWW7RqFof167P+36ePa/Xfv9+1/Bw65LpTvfCCdPbZbpldu1w9kvR04ldUBw/6pp991rWk5ZY0S64/dK9erm80Quvf/3ZJs+S+55aApLnEat3aDcuZ19jY55zjco0wTppLElqcTw1IUX7xGD9eevzxgj8uMtL1kercWbrjjsBbZCIjfWOoZmTkf6lwRoZ7jFtR32NLS/zDb2V807Q4nx7//a8b3zk1VbrwQte3v1Il4h8qZq4lqE4d19KwZYuvj5v3M/HSS+6ySeIfHC+9JN13n5seMcK1vJXQ41CJir+3orb33P/oo7kvu2OHO8d7K5xPnSr9/vfEvyiio914tV5RUdKdd0qDB7tRRiIj3VjyCxZIzzzja5GrUsW1VIfhZdul4vxrJjVu7MYXjopysQ7jq4xKdPyTkqSbbw78x6Df/95dJVOhAvEvAlqcg8m/wl1BREW5ioN9+xbsAOO/8weSBPsv419lDyiNli51hahSU92wJR9/HNZJc6ng8bgvSgsWuL5XGze6vqDefuVmLsn73/+IVTB8/bU0apSbbtgwrId6KXX27PFVac4raZZci7N/3/5//pP4FZX/95lKlVxr/0svuZEtKlSQypVzP+Ddequ0erX7QU9yLc9/+QvxC5VZs1zSLEk33lhikuYS5/hxqUsXX9J88cXuB7mff3bzkpOlH35w3aNq13bLvPOOuzImDBvNShIS52B+YfzDHwr2mEsucYXEduxwO7S3wEWg/Ie/CGRcNv/Ls2Ni2GYovRYskLp2dft8YqK0cKHrHoHid8kl0qef+pLnjAzpxReJS1EdOeKGeklNdcfzmTOpWhvOunXzTfsXtkLh+A+rM3583kPiVa0qPfWU738u1w6dN97wTd95J/EIlb/8xf1wKrkW/lWrXIty7dquYax8efdj0fDh7nhz4YVu2U8/dQk2SJxPS6LsfzOT/vOf/B8XHS317+8qMa9b58r151fUK6+TgZf3krG8+C9TowbbEKXT+++7L6lHj7oTx9KlvpMGTo+oKGnCBN//3j7nKLwRI1yLgiRNnuyGxEP4qlWrYOdr5M2/pkz//vkv7z80m/dzg+Datk36/HM3fcklef+YgaJ57z3f9LPP5t0YVq2aW8ZryhTiR+JczEmyN1H2v516ID/VxRe71uXt290X+86di/6+Lr7YN51TQZ5Tbdnim65Th+2K0iU93Q3Bc+ONrhXu8sulZcvY18NFkya+aW9fTxTetGm+6Vtuyfk8dWrdi9zuR+ilpvqmueKr6PyLIfm3PufGv6EhJYX4hcL06b7LgAcOJB6htGGDb7pdu/yX79TJN+0dHhIkziFLlE9NknPrvH5qH2f/1uX1613rcjALUlx2mW969er8l/dfpmFDtjNKjz173DBvf/yjO3HfdJPr85aQQGzChX8V3EC+6AKlyfff+6bpNlJ0/l3bvvsu/+V37fJN0+82NP7xD9/09dcTj1Dyz0MC6bNcpoxv+sQJ4kfiHLxE2bw7ZH5Jck68Lc4XXeRal7dtC17rck78S8vPmpX/8u+/75vu0IHtjtLzhbRlS3eJWESEuyTp3XcDG6INhXf22b5j59Kl+S/v35Xl0kuJH0rP/h9IsTv/gmCh+k5wJrn6at90IMW+5s3zTfs3OiA4duxww69J7gfrBg2ISSh5C35JgbUg+5+jGa6QxLnICbNfklykC9gSEtyloRs2uNbluLjQvv8OHXx9lVetkl5+Ofdln3tOWrPGTVepIl15JXs/SsfJ+sorXTeEmBiXnHmrDKP4vriOGiUdO5b7sr/8knWoPi7jK7qcroTK7+qowvwgjJz5n0PHjs17bOZvvnEVnyU3TNLttxO/YMTfOz7wtGlZ+3yeatcu143Ha8gQ4hds3r7NUmCXDqNoevTwTU+cmPcY5snJ0iOP+P7v3p34kTgXIln2T5iDJTpaatOm+NYlMlIaOtT3/733ui+ky5a5wkhHj7rp/v2zJhO33071VZR8GRnucrAdO9z+/MknWSvXIrQefNAV/ZKk5ctdK9rHH0u7d7skIiXF/Yj4/POuv7m3X3ODBtJttxE/lGz33eeucJHcft+jh+uWlZTk9v+jR10L3Pjx7uowb7/akSN9CR+K9l1u0iTf/4MGSfff7yoNHzzoEolff5XefFNq3txdASi5WgucJ4JvxQrfdNu2xCPUHnjA1+Xpww+lRo2kF15wx5wjR1xNhV9/df3OW7TwVfKvWNGdu1H4Q4/ZGfDTs38hlHxWt8QNgH7okCuSsXt3YMufd55reQ7T8WxLXPyDtN8R/0KYOlW6+ebCPTZM17HE7f+TJ0v33BP48mef7S4ZC9PEoVQdf0rgcajExf/pp6UxYwJfvkcPafZsN1wM8Q+Oe+5xx6FAxMa6Wi9hOspCiT7+dOjg6opI7sov/xbREnOoLGHx/89/XONBoH2Wy5Rxx59evYh/EZTuFuecinuVNpUqSf/6V2DFLi6+2LXKhWnSDBTIK68Qg9Pt7rvdL9qnFkbMyVVXuS4ltLahtHj8cde/Nr8q2RERbviwME6aS6yXX3Z1LfKraXHxxa5VlKEJQ+PHH33T/v1vETo9e7pxmRMT81+2Vi23bJgmzSUqtSx1Lc5F/HW9xP7il5LivsDOmeMuVdq7113KGh/vLk3q21f63e98l1aG7eajxZn4B6h8+bz71eaFFufg2rdP+vvf3SWr334rHTjgLp8/+2zX323QIKljR/Z/jkOlM/6//ebGRv30U9c94cABl0zXquVqAQwZ4i6lZP8PnW3b3GXZc+dKP/0kHT4snXWWS5gHD3bfgcqVI/6hEh0tHT/upg8cyHt4VuIfXMeOuYKoc+a4S7K948THx7vq8z17unMw+z+Jc45fDoq4OqXyixMHLhB/4g/iT/xB/Ik/iD/xL7QyJMwAAAAAAJTGxJmEGQAAAABA4pxH0kzCDAAAAAAgcc4hYZZImgEAAAAAJM45Js0kzAAAAAAAEuccEmaJpBkAAAAAQOJMwgwAAAAAIHEONGkmYQYAAAAAhIEIkmYAAAAAAEpK4kzSDAAAAAAgcSZpBgAAAACQOJM0AwAAAABInEmaAQAAAAAkzgAAAAAAkDifJrQ2AwAAAABInEmaAQAAAAAkziTNAAAAAAASZ5JmAAAAAACJM0kzAAAAAACnIXEmaQYAAAAAkDgDAAAAAEDiXHC0NgMAAAAASJxJmgEAAAAAJM4kzQAAAAAAEmcAAAAAABDqxJnWZgAAAAAAiTNJMwAAAACAxJmkGQAAAABA4gwAAAAAAEKdONPaDAAAAAAgcQYAAAAAgMS54GhtBgAAAACQOAMAAAAAQOJccLQ2AwAAAABInAEAAAAAIHEuOFqbAQAAAAAkzgAAAAAAkDgXHK3NAAAAAAASZwAAAAAASJwLjtZmAAAAAACJMwAAAAAAJM4FR2szAAAAAOBMYWa2cOFCy0+WZdzD8l8u0OcLk+dauHChSeLGjRs3bty4cePGjRs3bsV0C2Y+F6r80ePy4AIqxS3OixYtUocOHfhFhfgTfxB/4g/iT/xB/Ik/iL9LgUmcAQAAAADIXcH7OJM0AwAAAABInAEAAAAAAIkzAAAAAAAkzgAAAAAAFFfiTP9mAAAAAACJMwAAAAAAIHEGAAAAAIDEueAWLVpEEIg/8QfxJ/4g/sQfxJ/4g/hn8pgF2GmZ/s0AAAAAgDMQLc4AAAAAAJA4AwAAAABA4gwAAAAAAIkzAAAAAAAkzgAAAAAAkDgDAAAAAEDiDAAAAAAAiXNpwgDoxJ/4g/gTfxB/4g/iT/xB/P15zMzyX8ojBbAYAAAAAAClDS3OAAAAAACQOAMAAAAAQOIMAAAAAACJMwAAAAAAJM4AAAAAAJA4AwAAAABA4gwAAAAAAIkzAAAAAAAkzqXQ559/rj59+qhWrVoqV66catWqpc6dO2v58uXsGafB8OHD5fF4CEQxWbBggXr37q2aNWsqKipKcXFxatasmd5//32CE2Lp6el69dVX1bx5c8XExCg2NlYtWrTQfffdp/T0dAJ0GvzrX/+Sx+PhGFQMtm3blhnr3G4IrYyMDE2ZMkXt27dXlSpVFB0drdq1a+uPf/yjkpOTCVCI5Lff8xkoHnv27NGIESNUq1YtlS1bVrVr11b//v114MABglMMlixZop49eyouLk4xMTFq1KiR7rjjDh05ciT837wFIsDFSpIpU6aYx+MxSdluHo/H/va3vxmKz7vvvpu5PRB6Y8aMyXHf996GDh1qGRkZBCoEMjIyrF+/frnGvkOHDnbs2DECVYwOHjxoZ599duY2QGi9//77eR5/2AahlZKSYt26dcs19o0bN7aDBw8SqBDIb7/33ipVqkSwQiQpKcnq1q2bY9wvvfRS9v0Qe/nll3PNvy699FLbtm1bWL//M7LF+dChQ3rggQdkZvrd736n9evX6/jx41q/fr2uv/56mZkeeughpaam8rNQMfjb3/6mIUOGyJ1TEGrz58/XuHHjVKZMGT3yyCPatGmTkpOT9f3332v06NEqU6aM/v73v+uNN94gWCHwl7/8RbNmzVJiYqJmzZqlvXv36uDBg/r3v/+tWrVqadGiRRo/fjyBKkajR4/Wjh07CEQx8V7V9eyzz8rMcrwhdB577DF9/PHHat26tT799FPt379f+/bt06xZs3TOOedozZo1GjduHIEKTWNVnrcpU6ZIkqZOnUqwQuSJJ57Qjz/+qOuuu06//PKL0tLStGXLFnXr1k3fffednnrqKYIUIitXrtR9990nj8ejsWPHasuWLUpNTdXOnTv14osvasOGDRo4cGB4nwPOxBbnqVOnmiS79tprs807ceKENWvWzCTZf/7zH34aCqGdO3fakCFDTJKVL1+eloZi0qtXL5NkkyZNynH+s88+a5KsUaNGBCsEGjZsaJLsiy++yDbv888/N0lWp04dAlVMli5dah6Px+Lj4zkGFZNWrVqZJFu1ahXBKGbbt2+3yMhIK1u2rG3ZsiXb/C+++MIkWWJiIsEqZmvXrrXy5cvbrbfeSjBCyHt10caNG7Pc/+2335okO/fccwlSiPTv398k2V133ZXj/KFDh5okmz17Ni3O4eSLL76QJA0YMCDbvMjISN12222SpHnz5vHzUAhdcMEFmjp1qurUqaPPP/+cgBTz/j9o0KAc5998882SpE2bNhGsEFi7dq3MTK1bt842r2XLlpKk3bt3E6hikJqaqmHDhsnM9OKLLxKQYnD8+HF9/fXXqlq1qpo2bUpAitmMGTOUnp6u3//+96pTp062+a1bt5aZadu2bQSrGHm3SUxMjCZNmkRAQmjv3r2SlG3/v+CCCyRJu3btIkghsmTJEknSHXfckef3zw8++CBs1yGwxNlMKkWFCr766itJUuPGjXOc36JFC0nSl19+yV4eQidOnNDo0aP17bffZiYMCL09e/bIzHTWWWflmkxIUmxsLME6TT9qNGnShGAUg2eeeUYbNmxQ165dddNNNxGQYjr/Hj9+XG3bttXkyZPVqFEjRUdHKz4+Xv3799eaNWsIUggtXbpUktS3b1+CEUZeffVVrV27VmPHjlW1atUISAidf/75kqSffvopy/0///xzjgk1gicpKUmSVLdu3Rzne7/7fP311yU8cS5lvH3ZatWqleN87/2//PILe3kIbd26VRMnTlTFihUJRhj59NNPJUmdOnUiGMVk586dev3113XjjTeqXLlymjhxIkEJsQ0bNmjChAmKiYnRa6+9RkCKibd/89y5c3Xffffpm2++0fHjx5WUlKRZs2apZcuWVPYPoXXr1kmSGjRooOXLl6tXr16qVq2aYmJi1K5dO/373/8mSMVs//79Gjt2rBITE3NtiUPw3HvvvZKkRx55RFu3blV6erq2bdumRx99NMt8BF/lypWzJNA55QX+eVo4KnMmbrhDhw5Jyr1FzXs/ZelDKy4ujiCE4Wdj7NixiomJ0ZgxYwhIMbjyyiu1cOFCSVLVqlW1YMECtW3blsCEkJlp2LBhSk1N1cSJE09LC0NKSor279+vAwcOaP/+/dmm/f+PiIgI60vXCpM4e3+w6NKliypXrqydO3fqo48+0oMPPqghQ4aoUaNGql+/PjtrkO3Zs0eS+4F0+PDhWYa/W7ZsmZYtW6YnnnhCf/zjHwlWMXn11Vd14MABPfzww4qKiiIgITZ8+HClpaXpscce07nnnpt5f0xMjP785z/rzjvvJEgh0rJlS3388cd6//339eCDD2ab/9JLL0lSeA9LFXBv6FJUMCUiIsIk5TrcTnp6ukmyiIgIevIXI1GY57RKTU21Ll26WPXq1W3p0qUEpJi0bt3aLr/8cktISDBJlpCQYDNnziQwIfTqq6+aJLv88svtxIkTIT0GzZs3z7p3725t2rSxiy++2GrWrGnR0dEBD0sjyfr27VtqYp+YmJhn8Zfx48ebJAokhUhUVJRJsujoaBs0aJCtWbPGUlJSbNeuXTZlyhSLjY01SfbZZ58RrGJw7Ngxq1GjhlWoUMH27dtHQIrBmjVr7JprrsnMBby3iIgI69Wrl33//fcEKUTmzZtnkiwmJsZef/11++233+z48eO2fv16Gz58eOa5MTIyMnxzlTMxca5QoYJJspSUlGzzFi5caMnJyYyjR+J8RklLS7O+ffta9erVbcOGDQTkNHn66aetSpUqJsnmzZtHQEJg+/btVqlSJStTpox9/fXXxXIMmjdvnp177rkFSpb9b6+99toZsW0WLlxomzZtMklWu3ZtdtYQKFu2rEmy+++/P8f4v/DCCybJunXrRrCKwfTp002S3XzzzbZw4UICUgxJc8WKFa169er23nvv2Z49eywtLc127txpDz/8sFWuXNni4+Nt8+bNBCtEHn744RzPcxERETZlyhSTZLGxsSTO4aRmzZomyfbu3Zvj/D179pgkq1WrFns4iXOpd/z4cevdu7dFRUXZmjVrCEiYfJFq27YtwQiB3r17myQbNWpUsR6Djh49aqNHj87WyhHI7UxqATl69KhJsrJly7KzhkBcXJxJsq1bt+Y4/+effzZJFhcXR7CKwXXXXWeS7MMPPyQYxaBHjx4myaZNm5bj/DfffNMk2YABAwhWCM2ePdvatWtnFSpUsCpVqli3bt3siy++sJ07d5okq1u3bti+d8/JLwv583hcde1SoE2bNlq+fLnWrl2rhg0bZpu/Zs0aNWnSRC1atNDKlSvplFBMPCcrt1sp2c9KgqSkJA0cOFALFizQ2LFj9eSTTxKU02z//v2qVq2aYmNjM+sxIPjHmUD7QgfbsmXLdNttt2njxo0BLZ+YmHhGDQ3k3f8rVaqkgwcPssMGWbNmzfTVV18pJSVF0dHR2eYfO3ZM5cuXV7ly5XTs2DECFkJHjhxRfHy8oqKitGfPHpUrV46ghFj16tW1d+9e7d27N8fq5bt371aNGjVUo0YN7dy5k4AVs7lz56pbt27q1q2bPvroo7B8j2dkVe1LL71UkhtPNSfe4TCaNWvGXoxjAi3pAAAc1ElEQVRSa/ny5WrcuLH++9//aty4cSTNxSQhIUEejyfXxMlbFKNMmTIEqxRq27atvvrqK40ePVqRkZH5Lr99+/ZStf7x8fHyeDzZhoLx2rJlS5bzNILLO3b2N998k+f+lpiYSLBCbNGiRTp27JiuuuoqkuZikpKSkuf51Xv/4cOHCdZp4B3VpV27dmH7Hs/IxLlLly6SpKlTp+Y433t/69at2YtRKn300Ufq1KmTDh48qA8++ECPP/44QSkmHTt2lCS9/fbbOc5/7733JEnt27cnWCFwsotSjrdTlwmV8uXLa+LEiVq9enVmIpObadOmyePxFKilPJxdccUVkqQZM2bkOP8f//iHJKlfv37srCHQp08fScp1yDvvdrnmmmsIVoh5K8y3aNGCYBQT7zjBuY0XP3fuXElSo0aNCFYI1K9fXx6PR19++WW2ecePH9cHH3ygcuXKadCgQWH9JaIgnVBLxbX1R44csapVq5oku+mmm2zDhg2WmppqGzZssAEDBpgkq1Klih06dIiOCPRxLnXWr19vMTExVrlyZVuxYgUBKWZffvmlRUREWLly5WzKlCm2d+9eO3HihO3YscMmT55s5cuXtzJlytjKlSsJ1hlwDEpNTbWJEydauXLlcuzfvG3btizvr6RbsmRJZlXV6dOn2/79+y0tLc22bt1qf/7zny0qKsrOO+88O3r0KDtkCKSnp1vLli1Nko0YMcJ++uknS01NtT179tibb75pFSpUsHLlytn69esJVohdffXVJsk++eQTglFMZs2aZZLs7LPPttdff91+/vlnS0tLs6SkJJsxY4bFx8ebJEa2CJHHHnvMJFm9evVs/vz5duTIETtw4IDNnTvX2rdvb5Js5MiR4f094UxMnM3M/vSnP+VaiMXj8diUKVPYw0mcS6U+ffoEXJQIofHKK6/kWiSqbNmy9s477xCkM+wY9O2331qLFi2y7AsXX3xxju+xpJswYUKux5yEhARbt24dO2MIbdmyxc4///wc4x8ZGWlvvfUWQSoGdevWNUmMZHGakrfcbqNHjyZIIZKcnGzNmjXLNfYDBw60tLS0sF6HwIuDSaWqQJgkffDBB3rjjTe0evVqHT58WPHx8WrevLnuuecederUiWsqihnFwYpHTEyMkpOTA76sFaGxbNkyTZo0ScuWLdOhQ4dUo0YNXXXVVRo5cqQaNGhAgM7AY1BGRobeeustPfDAAzp69Kj+8Ic/aPLkyaXyWPnZZ5/phRde0IoVK3T06FHVqlVLPXv21GOPPab4+Hh2xhA7dOiQnnvuOc2aNUs///yzYmNj1bp1a40aNUpt27YlQMWgYsWKOnr0qA4dOqTY2FgCUozmz5+vl19+WStWrND+/ftVtWpVtWrVSnfffbeuvvpqAhRChw8f1oQJEzR79mxt2bJFVatWVYMGDTRixAj17ds37LslndGJMwAA4eann37SsGHDdPfdd6t37955Jvr8uAUAQDgmziTPAACEnJkpLS1NZcuWzeeUTPIMAEBxiCAEWS1atIggEH/iD+JP/E8rj8eTb9JM/EH8iT+IP/EvxnMzLc4AAJRctDoDABB6tDgDAFCCmVmpGecZAIDSkzibuVZnAAAAAABInAEAQLij1RkAABJnAAAAAABInAEAQOHR6gwAQLglzvRzBgAAAACQOAMAgJKEVmcAAMItcabVGQAAAABA4nzmWbRoEUEg/sQfxJ/4l1gltdWZ/Z/4E38Qf+IfzjxmZkV7Bo9rfQYAAOFxcvd4ZJybAQAIGlqcAQAAAAAgcQYAAAAA4HQlzhQJAwAgrFBdGwCAcEucAQAAAAAgcc4Hrc4AAAAAABJnAAAAAABInAuPVmcAAAAAAIlz6ccA6MSf+IP4E38Qf+IP4k/8Qfz9eczMgvuMHtf6DAAATt8J3uORcT4GACAogt/izCXbAAAAAAASZ5JnAADCGWM5AwAQ7okzAAAAAAAkzvmg1RkAAAAAQOJM8gwAAAAAIHEGAAAAAIDEOWRodQYAAAAAkDiTPAMAAAAASJzPiOR50aJF7BXEn/iD+BN/EH/iD+JP/EH8M3nMzIr3FT0uiQYAACE+5XpknHMBACiy4i8OxmXbAAAAAAASZwAAAAAASJwLj1ZnAAAAAACJM8kzAAAAAIDEmeQZAAAAAEDiTPIMAAAAACBxDofkmQQ6u82bXVyqVy95733uXKlpUyk6WkpIkLZvZ3sGYt48t81btSIWAIp4emUoKgAASk/i7E2ew6D1OeABuI8c8SX7O3eyJ+WU/HXrJq1ZI7VoIXXtKlWuHLz44/Tu/6Ewd677PK1eTfxB/Ik/iD/xB/En/mHFY+H4c7TH45LocHbkiBQb66Z/+02qWTM0r7N5s3TBBVJcnJSUVHI+Ac2buwTo9del4cM5IhT0R4drr5VatpRWrDhz1nvSJGn0aGnVKunyy9kPAAAAEDbCcxxnLt0u+b77zv0dOJBYIDA//UQMAAAAQOJc4OSZwmEl17FjUmRkQJdnAyTOQNHNnj1bL774Yub/zz//vP71r3+dUTFIT0/XggULMv9PSUnRkiVL2DkAAKU4cfZPoEtS6/P06e69vveeu7T6rrukxERXIKthQ2nmzNzXc8oUqVkzqUIFVwzs5pul3btzf63586Xu3d1l3NHR0iWXSOPHu6TVy9tvND4+e1/sCRPcvP79A1u3BQuknj3deytXTjrvPLd+W7f6luna1bet0tN9287jkQ4cyPv5jxxx76lpU3cZfMWK0tVXSytXZl924kSXmM+b5y5nbtNGKl/ed4lvfvMDfb0+fdx7f+qpnN/zU0+5+WPH5h+/QLZXsJ7DW2Ds1VelbdukG2+UqlVzRdoeeEA6ccItN22ab59LSJBuuSXvLgGBvn5BPgfPPSfVqeP2L8ld5u/dZ2bM4CiNQtu8ebM8Ho+ql8DiinPnzlXTpk0VHR2thIQEbc+nuOKOHTs0dOhQnTj52V65cqVGjRqlqKioMyquw4YNU+fOnbXiZDeXDh06qHPnztq3bx8fCABAUfPSEsTbDh0ODh/2vZ/ffvPdP22au69fP7Natcxq1jTr2dOscWPf8jNnZn++YcPcvMhIs/btzXr1MktIMKte3d0fF5d1+aeecvdHRZm1aWPWvbvZOee4+9q3Nzt+3Lds9+7u/qFDffft2GEWE2NWoYLZr7/mv77jxrnniIhwr9e7t1m9eu6+qlXNvvzSLTd9utlVV/nW9frrzQYMcLejR3N//o0bzc47zz2mVi2zbt3MmjRx/0dHm339ddblJ0xw826/3a3DpZea3XCDuz+Q+YG+3ocfuvsSE83S07O+h4wMF3OPx+zHH/OOX0G219y57v6WLYv+HF26uH3wwgvNevQwq1HD3f/YY2bPP++mmzUz69rVrHJl93+TJmapqUVbh4J8DmbOdPtHxYq+9+zdZ774wlD6HD582CSZJPvN//gZZD/88INJsrhTj59hbu7cuSbJPB6PtW/f3gYPHmyHDx/O8zFjx461iIgIO3jwoJmZ3XDDDZaYmHjGxfWvf/2rSbJnnnnGzMzuuecek2Rz5szhgwcAKFoqWjLftU5/Ep1f4iyZPfRQ1gTk0Ufd/ZddlvW55sxx91eq5EtAzcySk13ycGriPH++u69evawJW1qa2b33unkTJ/ru//57s7JlXYK3erW779Zb3XLjx+e/rgsW+N7f8uVZE8cnn3Tzatc2O3Ys+zbavz+weG7d6hLFN94wO3HCd/+DD7rn6ds358RZMnvkkezPl9/8QF/vxAmXNEtmn3+e9TkWL3b3d+yY97oVdHvllDgX9jm82zgjw92/f79Z3bpuW8bEmL32WtaYxMe7x/zzn0Vbh8J8DurWdfevWsWRmcT5jE6cL7/8cpNkr7/+esCP6dy5syUkJGT+X6tWLbviiivOuLh+9913JslGjRplZmYzZswwSfbOO+/wwQMAnIGJczgk0fklzh06ZH/MwYO+xxw44Lu/Sxd33+OP5/QNJXvifPXV7r7PPsu+fGqqazmsXz/r/SNHuse0bWv2zTeuZbtevazJbm66dnWPnTQp5/nNmrn506YVPnHOzS+/5Nzi7k2M69TJmvgGOr8gr+dN9G67LeuyI0a4+6dOzfs5C7q9ckqcC/scTZtmX977Y0fdutnn3X+/m3fffUVbh8J8DkicSZxJnM3MLDo62iTZAf/PRz6uuOIKK1eunCUnJ5uZWUJCgl1wwQVnXFzXrl1rkuzVV181M7N3333XJNnixYv54AEAiqRMKbjW3Dft3w/6dA9nVbFi9vsqVZKqVHF9fffudYWz0tOlpUvd/EAqUKelSYsXu+k33pDeeiv7Mh6PtGmTlJoqlS3r7hszxvVnXbbMja+cni699JLrq5yX9HTJO7Zat245L3PttdL//uf6vw4aFNw41qrl/u7dK2VkSBGndMvv0sX1Zc5NfvMDeb1bbnF9oWfPll55xcU0Pd39X6mS1K9f8LdXsJ4jp/6NiYnu70UXZZ9Xu7Zv/YPx+oF+DgBkOnbsmCIjI1W5AJ+Niy66SIsXL9acOXM0YMCAzP/Xrl2rRo0anTGxW7BggaKiotSrV6/M/88++2y1atWKHQsAUCQRpWpt/Nuf/YtSFaC4WMgH4PYmcRkZ7m9SkpScnDVpyUtSknT8uJueNcsVWTr1duCAK4SVluZ7XGys9PTTbnrbNqlz59wT4VNfLyXFTZ9zTs7LeO//9dcih2fJ/PnSZ59Jzz4rjRghDR3qm+mN2akJW17ym3/8eP6vV7eu1LGjtH+/9PHH7r7//lfatUsaMMAV1gr29gr2c+S0D5bJ/rvZ5i1bssc62K+f0+cAxXP8CXPTp0+Xx+PRe++9p6SkJN11111KTExUdHS0GjZsqJm5FFc0M02ZMkXNmjVThQoVVL16dd18883anUdxxfnz56t79+6Ki4tTdHS0LrnkEt1222065lfobu7cufJ4PIqPj9fOU4orTpgwQR6PR/0DLK64YMEC9ezZU9WrV1e5cuV03nnn6a677tJWv+KKXbt2lefkMSs9PV0ejyfzdiCf4oqDBw9Wly5ddNHJH8SGDh2qG264QXFxcfmu9/jx47Osd1Hi6t1uZ511liIiIrKsg8fjUZ06dXJ9P2XLls3z/SxZskRXXHGFYmJiVLFiRV1xxRX68MMPsywzefJkPfLII0pMTNTBgwc1ffp0vfjiiyEpksbxB8Sf+OMMi7+Z2cKFC/Ntmg5kmWAvF/TnyppaZ7vEe+HChZm3/Czx70Oa06Xa3bvn/DxxcW7+Dz+4/3fs8D2Ptx9qDpdqp1aq5P5PSnLLVqhQ8Fh4L8OVLLVKlayXyeb2XDt3+t7fyaIzp9r0wANufqdOuV6qnW9M09Ptx9tvt7RKlXLfTmlpmYv/dPvt7r7hw3N+Pu+l2rnNT093fXKrVQvo9dY/9pi7r39/d8ctt7j/TylelW09c9leecbj5L518OKLC7/NT7ncO8syf/+7m3fdddme5wdvf+Xf/S7Pdch3e/p9DnLk9znIfK48LtUuFcefAJYJ9PgTju+/IMt5i1/plEu1p02bZpKsX79+Fh8fbzVr1rSePXta48aNM5efeUpxxYULF9qwYcNMkkVGRlr79u2tV69elpCQYNWrVzdJVsl7/DzpqaeeMkkWFRVlbdq0se7du1tCQoJJsoYNG9pxv0J33bt3N0k21K+44uzZsy0mJsYqVKhgv+ZSXNE/FuPGjTNJFhERYW3atLHevXtbvXr1TJLFxsbalydrW0yfPt2uuuqqzHW9/vrrbcCAAdapUyc7mldxxQDj713vMmXKZK73OeecY5Ksffv2WdY7kLjGxcVlec1Dhw5Z3bp1TZJdcskl1rNnT6tdu3bm+kiyjn41Ify3w6WXXmqtWrXK3A6nvp/PP//cIiIizOPx2GWXXWZdu3a1SpUqmSRbsWJF5nKbN2+2jIyMzPf1g/dcG0b7P8cf4k/8iT/xL5nxL/l9nIMhl0Q6T/n1cQ4gYTAzs5QUV6laMtu1K9fEObPPbUaGr/rxtm2Br+PGja4a8vnn+wqOPfhg/o87ccKsfHm3/Lp1OS/z+ONu/pAhuSbO+TrZBzs5MdFs8mSzr74y27fPvX4OiWy+iXF+8719vuvVC+z1UlLMqlRxFbeTklwl8Ysuyn+9CrO9Tu3jHIzn8JdL4mxmOSfOhXn9gn4O8kmczxSBHuxLutz6OHsTZ0n20EMPWapfUblHH33UJNllpxSVmzNnTmZy/KVfccXk5GQbMGBAtr648+fPN0lWr149+9Gv0F1aWpr169fPJNlEv0J333//vZUtW9Y8Ho+tPllc8dZbbzVJNj6A4ooLFizIfH/L/YorZmRk2JNPPmmSrHbt2nbMr96ENwb7i1ojwk9e633vvfdmW++CxtXM7LHHHjNJduedd1rGyR+Cjx8/btdcc41JsjfeeCPz/pzez8KFC3N9P+3atTNJ9re//S3LfjRr1iy+w3D8If4g/sS/eFJGNlWhv/kFJ3E2M7v4Ynff++9nX37VquzFqvr2dfd5h1YKxLXXuse8+66rihwV5Spt5/NrvJn5ioM9/3zO85s3L3pxsNhYV/U7px8PQpE4F/T1zHzFwAYOzLtY2qkKur1ySnqD8RwBJM42eXL2xLkwr1+Yz4F3eDOKg53xiXOHHIrKHTx4MPMx/kWzunTpYpLs8RyKK+ZUxOrqq682SfZZDoXuUlNTrXLlylb/lOKKI0eONEnWtm1b++abbywyMtLq1auXJdnN/fDZ1STZpFyOF82aNTNJNs3v+BmKxLmg613QuJpZZoL8v//9L8v93srWTzzxRKHfz0UXXWSSbP369XyAAACnRQRX1YcBb1GwJ56Q9u3z3b9xo3TTTdmXHz3a9REdN06aNy/rvK1bpdtvl959V34d9dytSRPpxhul88+Xhg93hZwefDD/9zdypPv79NPS6tU6paOftGqV658dYF+/HGVkuHR1zx7ffYcPS3/4Q2hiXpjXu/VW93fGDNc/ePDgwF6roNsrVM9RFMXx+tWqub/eftZB6DOPkqliDkXlKlWqpCpVqkiS9p4sXpeenq6lJ4srDgyguGJaWpoWnyx098Ybb2jgwIFZboMHD5bH49GmTZuUmpqa+bgxY8aoZs2aWrZsmbp166b09HS99NJLKpdPccX09PTMflvdcqkpce2112b29Q2Vgq53QePq5d0++/zPY37bq+zJwoGF2Q7e4l7Dhg3T119/zYcEAFDsyhCCMPDAAy7p2LDBFaJq0cIVY1q2THr0UVe4yl/z5tILL0j33usqWjduLJ13nrR9u6tunZ4unXWW9xuTdP/9bnriRF+xrDFjpLfflubMcdWwO3fO/f116uQSpjFjpFatpNatpfh4ad066fvvpapVpfffz79Cd14GDXIVm9u0cYW4UlKkFSukq65y8fjyy+DGvDCv17Sp+/Hh669d3GvWDOy1CrK9QvkcRVEcr9+pk4v7sGHSX//qKnm/844rwAZIijxZVC7jZFG5pKQkJZ8srlg7gOKKSUlJOn6y0N2sWbNyXa58+fJKS0vLTPRiY2P19NNP67bbbtO2bdvUuXPnXBPhU18v5WRxxXNyKa7ovf/XEP5QVND13r9/f4Hi6nXDDTdo5syZeuihh/T222+rfv36WrVqlZ555hl5PB5179690Nth0qRJWr9+vZYtW6amTZuqYcOGuvXWWzV8+PB8f8AAACAYaHEOBxUruiGfBg1ywx8tXepag2fNkp580iUpp7r7bmnJEql3b1cl+8MPXRLTt69LOJ580i03ebIbJqhHDzc0k1dCgkuEvYl7enre7/Hxx12Cfe21LsH/8EOX3N95p7R2rUs2i+L5511SVrGiq3K9e7dL9D/4ILDq38X1es2bu7+33FKw1wt0e4X6OYoi1K//6KOuFT8iwl3ZMGCAdMUVHB+Qqwy/quwxMTH5Lu9NhCtUqKCTXZVyvCUnJ2d7vnXr1mVOr127VgcPHgzKOtjJoRM9AY78UBgFXe+CxtXr+uuv18iRI7V27Vo1adJEFSpUUIcOHbRr1y5NmDAhc1iswmyH+Ph4ffHFF5ozZ45uuOEGbdq0Sffee69at26to0eP8mEAAIScx+x0D3gMlBDJyW4M5LJlXeLI8CZAoRw5ckSxsbGSpN9++001T169MX36dA0ePFjdu3fPNsyQJFWvXl179+7VDz/8oHr16unYsWOZid6uXbuUkJCQZfnNmzfrggsuUFxcnJKSkmRmqlq1qg4ePKht27Yp0TumeT42bdqkyy67TOecc46aN2+umTNn6sEHH9Sf/vSnPB+Xnp6u2NhYpaSkaN26dbrkkkuyLTNmzBg9/fTTGjJkiN5+++0sSfT+/fszL38uanJekPUuaFy90tLS1Lt3b/3www/q16+ftm/frrPPPltDhgzJHCarMO8nJ7t27VKvXr305Zdf6pVXXtGIESP4YAEAQooWZyBQ06a58YoHDSJpBsJAdHS06tevLynnMSBPHfvY4/Go88luKdOmTQv4de6//36lpaVp3LhxeuaZZxQVFaXJkydr8+bNeT4uMjJSHTp0kCTNO7U2wEmffPKJJGW+r1Ao6HoXNK5e//d//6ePP/5Yt99+u5555hlNnTpVEyZMyJI0F2U7+KtRo4aGDx8uSfrpp5/4MAAASJyBsPHKK+5vQS/TBhAy3uJVTzzxRJaiVBs3btRNORRXHD16tCIjIzVu3LhsyezWrVt1++23612/Qndz587V3Llz1aRJE9144406//zzNXz4cKWmpurBAIorjjxZXPHpp5/W6lOKK06YMEGrVq1S7dq11b8oxRUDUND1LmhcJdcn2ftarVq1Ur9+/XTLLbfo4Ycf1vTp07Nc3l7Q9/PCCy9oxowZmf9nZGRkXpXQoEEDPggAgNCjsDgQgEWL3FBJnTsTCyDEw1F1z2UYs7i4OJNkP/gNY3b48GG78MILTZJVqVLFunTpYh06dLAyZcrY2LFjrVy5ctmGTXr55ZfN4/GYJGvcuLH16dPHWrRoYZGRkSbJxo4dmzksUv369U2SffLJJ5mP37Vrl1WsWDHX4ZRONW7cOJNkkZGR1q5dO+vTp0/me65ataqtXLkyy/IKwXBUBVnvwsY1NTXVOnbsaNHR0VapUqXM9fDeatSoYRs2bCjw+0lOTrazzjrLJNmFF15ovXr1ynxvDRo0sJSUFD5UAICQI3E+BQOgE3/iD+JfchJnM7PffvvNBg0aZNWqVbMKFSpY69at7f/+7//MzKxly5bZEjwzs6VLl1rv3r2tevXqFhUVZYmJiXbllVfa4sWLM5d5/vnnTZL16NEj2+OfffZZk2SXXXaZnThxIt91nj9/vvXo0cPi4uIsKirKateubXfeeaf9+uuv2U/MIUqcc1vvG264Ict6FzauXbt2tYSEBPv2228z7zt27JitX7/errvuOpNkw4YNy/X9lClTJtf3s2PHDrvjjjssMTEx832PGDHC9u7dy4GD4w/xB/En/sWC4mAAAKBIvvnmGzVq1Ei33nqr3nrrrWzz582bp2uvvVaDBw/WO++8Q8AAACUOfZwBAECRePs3z58/P9uY1Pv379dzzz0nSerTpw/BAgCUSLQ4AwCAIrv55ps1depUlSlTRs2bN1fNmjW1b98+rVq1SsnJybrjjjv02muvESgAAIkzAAA4M2VkZGjq1Kl65513tHbtWh06dEhVqlRRixYtdOedd6pnz54ECQBA4gwAAAAAQGlEH2cAAAAAAEicAQAAAAAgcQYAAAAAgMQ51BYtWkQQiD/xB/En/iD+xB/En/iD+GeiOBgAAAAAAHmgxRkAAAAAABJnAAAAAABInAEAAAAACLr/B0JLFr33HW+QAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Slicing works very similarly to how we saw with strings.  Remember, python uses 0-based indexing\n",
    "\n",
    "![slicing.png](attachment:slicing.png)\n",
    "\n",
    "Let's create this array from the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.arange(9)\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now look at accessing a single element vs. a range (using slicing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Giving a single (0-based) index just references a single value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Giving a range uses the range of the edges to return the values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 3]\n"
     ]
    }
   ],
   "source": [
    "print(a[2:4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `:` can be used to specify all of the elements in that dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2],\n",
       "       [ 3,  4,  5],\n",
       "       [ 6,  7,  8],\n",
       "       [ 9, 10, 11],\n",
       "       [12, 13, 14]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(15).reshape(5,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(15).reshape(5,3)[1:2,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multidimensional Arrays\n",
    "\n",
    "Multidimensional arrays are stored in a contiguous space in memory -- this means that the columns / rows need to be unraveled (flattened) so that it can be thought of as a single one-dimensional array.  Different programming languages do this via different conventions:\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/4/4d/Row_and_column_major_order.svg/340px-Row_and_column_major_order.svg.png)\n",
    "\n",
    "Storage order:\n",
    "\n",
    "* Python/C use *row-major* storage: rows are stored one after the other\n",
    "* Fortran/matlab use *column-major* storage: columns are stored one after another\n",
    "\n",
    "The ordering matters when \n",
    "\n",
    "* Passing arrays between languages \n",
    "* Looping over arrays -- you want to access elements that are next to one-another in memory\n",
    "  * e.g, in Fortran:\n",
    "  \n",
    "    ```\n",
    "    double precision :: A(M,N)\n",
    "    do j = 1, N\n",
    "       do i = 1, M\n",
    "          A(i,j) = …\n",
    "       enddo\n",
    "    enddo\n",
    "    ```\n",
    "  \n",
    "  * in C\n",
    "  \n",
    "    ```\n",
    "    double A[M][N];\n",
    "    for (i = 0; i < M; i++) {\n",
    "       for (j = 0; j < N; j++) {\n",
    "          A[i][j] = …\n",
    "       }\n",
    "    }  \n",
    "    ```\n",
    "  \n",
    "\n",
    "In python, using numpy, we'll try to avoid explicit loops over elements as much as possible\n",
    "\n",
    "Let's look at multidimensional arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4],\n",
       "       [ 5,  6,  7,  8,  9],\n",
       "       [10, 11, 12, 13, 14]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.arange(15).reshape(3,5)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the output of `a` shows the row-major storage.  The rows are grouped together in the inner `[...]`\n",
    "\n",
    "Giving a single index (0-based) for each dimension just references a single value in the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing slices will access a range of elements.  Think of the start and stop in the slice as referencing the left-edge of the slots in the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1],\n",
       "       [5, 6]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[0:2,0:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Access a specific column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  6, 11])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes we want a one-dimensional view into the array -- here we see the memory layout (row-major) more explicitly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4],\n",
       "       [ 5,  6,  7,  8,  9],\n",
       "       [10, 11, 12, 13, 14]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can also iterate -- this is done over the first axis (rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4]\n",
      "[5 6 7 8 9]\n",
      "[10 11 12 13 14]\n"
     ]
    }
   ],
   "source": [
    "for row in a:\n",
    "    print(row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n"
     ]
    }
   ],
   "source": [
    "for row in a:\n",
    "    for el in row:\n",
    "        print(el)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or element by element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n"
     ]
    }
   ],
   "source": [
    "for e in a.flat:\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally speaking, we want to avoid looping over the elements of an array in python&mdash;this is slow.  Instead we want to write and use functions that operate on the entire array at once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<span class=\"fa fa-flash\"></span> Quick Exercise:\n",
    "\n",
    "Consider the array defined as:\n",
    "\n",
    "```\n",
    "\n",
    " q = np.array([[1, 2, 3, 2, 1],\n",
    "               [2, 4, 4, 4, 2],\n",
    "               [3, 4, 4, 4, 3],\n",
    "               [2, 4, 4, 4, 2],\n",
    "               [1, 2, 3, 2, 1]])\n",
    " \n",
    "```\n",
    "\n",
    "  * Using slice notation, create an array that consists of only the `4`'s in `q` (this will be called a _view_, as we'll see shortly)\n",
    "  * Zero out all of the elements in your view\n",
    "  * How does `q` change?\n",
    "</dov> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Copying Arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simply using \"=\" does not make a copy, but much like with lists, you will just have multiple names pointing to the same ndarray object\n",
    "\n",
    "Therefore, we need to understand if two arrays, `A` and `B` point to:\n",
    "* the same array, including shape and data/memory space\n",
    "* the same data/memory space, but perhaps different shapes (a _view_)\n",
    "* a separate copy of the data (i.e. stored completely separately in memory)\n",
    "\n",
    "All of these are possible:\n",
    "* `B = A`\n",
    "  \n",
    "  this is _assignment_.  No copy is made. `A` and `B` point to the same data in memory and share the same shape, etc.  They are just two different labels for the same object in memory\n",
    "  \n",
    "\n",
    "* `B = A[:]`\n",
    "\n",
    "  this is a _view_ or _shallow copy_.  The shape info for A and B are stored independently, but both point to the same memory location for data\n",
    "  \n",
    "  \n",
    "* `B = A.copy()`\n",
    "\n",
    "  this is a _deep_ copy.  A completely separate object will be created in memory, with a completely separate location in memory.\n",
    "  \n",
    "Let's look at examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "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": {},
   "source": [
    "Here is assignment&mdash;we can just use the `is` operator to test for equality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = a\n",
    "b is a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since `b` and `a` are the same, changes to the shape of one are reflected in the other&mdash;no copy is made."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 2 3 4]\n",
      " [5 6 7 8 9]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(2, 5)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b.shape = (2, 5)\n",
    "print(b)\n",
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b is a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 2 3 4]\n",
      " [5 6 7 8 9]]\n"
     ]
    }
   ],
   "source": [
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A shallow copy creates a new *view* into the array&mdash;the _data_ is the same, but the array _properties_ can be different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n",
      "[ 0  1  2  3  4  5  6  7  8  9 10 11]\n"
     ]
    }
   ],
   "source": [
    "a = np.arange(12)\n",
    "c = a[:]\n",
    "a.shape = (3,4)\n",
    "\n",
    "print(a)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the underlying data is the same memory, changing an element of one is reflected in the other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0 -1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n"
     ]
    }
   ],
   "source": [
    "c[1] = -1\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even slices into an array are just views, still pointing to the same memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3 4 5 6 7]\n"
     ]
    }
   ],
   "source": [
    "d = c[3:8]\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "d[:] = 0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0 -1  2  0]\n",
      " [ 0  0  0  0]\n",
      " [ 8  9 10 11]]\n",
      "[ 0 -1  2  0  0  0  0  0  8  9 10 11]\n",
      "[0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "print(a)\n",
    "print(c)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are lots of ways to inquire if two arrays are the same, views, own their own data, etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "print(c is a)\n",
    "print(c.base is a)\n",
    "#print(c.flags.owndata)\n",
    "#print(a.flags.owndata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a copy of the data of the array that you can deal with independently of the original, you need a _deep copy_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0 -1  2  0]\n",
      " [ 0  0  0  0]\n",
      " [ 8  9 10 11]]\n",
      "[[0 0 0 0]\n",
      " [0 0 0 0]\n",
      " [0 0 0 0]]\n"
     ]
    }
   ],
   "source": [
    "d = a.copy()\n",
    "d[:,:] = 0.0\n",
    "\n",
    "print(a)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boolean indexing (aka masking)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are lots of fun ways to index arrays to access only those elements that meet a certain condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11]])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.arange(12).reshape(3,4)\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we set all the elements in the array that are > 4 to zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False, False, False, False],\n",
       "       [False,  True,  True,  True],\n",
       "       [ True,  True,  True,  True]])"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a > 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5,  6,  7,  8,  9, 10, 11])"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[a > 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3],\n",
       "       [4, 0, 0, 0],\n",
       "       [0, 0, 0, 0]])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[a > 4] = 0\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(12)\n",
    "b = np.arange(12)*123"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   0,  123,  246,  369,  492,  615,  738,  861,  984, 1107, 1230,\n",
       "       1353])"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False, False, False, False, False, False, False,  True,  True,\n",
       "        True,  True,  True])"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a>6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 861,  984, 1107, 1230, 1353])"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b[a>6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and now, all the zeros to -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11]])"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[a == 0] = -1\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have 2 tests, we need to use `logical_and()` or `logical_or()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 0,  0,  0,  0],\n",
       "       [ 0,  0, 10, 11]])"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.arange(12).reshape(3,4)\n",
    "a[np.logical_and(a > 3, a <= 9)] = 0.0\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our test that we index the array with returns a boolean array of the same shape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False, False, False, False],\n",
       "       [False, False, False, False],\n",
       "       [False, False, False, False]])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a > 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have three or more tests, look into `reduce()`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Avoiding loops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, you want to avoid loops over elements on an array.\n",
    "\n",
    "Here, let's create 1-d x and y coordinates and then try to fill some larger array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3200,)\n",
      "(6400,)\n"
     ]
    }
   ],
   "source": [
    "M = 3200\n",
    "N = 6400\n",
    "xmin = ymin = 0.0\n",
    "xmax = ymax = 1.0\n",
    "\n",
    "x = np.linspace(xmin, xmax, M, endpoint=False)\n",
    "y = np.linspace(ymin, ymax, N, endpoint=False)\n",
    "\n",
    "print(x.shape)\n",
    "print(y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll time our code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time elapsed: 11.915959119796753 s\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "\n",
    "g = np.zeros((M, N))\n",
    "\n",
    "for i in range(M):\n",
    "    for j in range(N):\n",
    "        g[i,j] = np.sin(2.0*np.pi*x[i]*y[j])\n",
    "        \n",
    "t1 = time.time()\n",
    "print(\"time elapsed: {} s\".format(t1-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's instead do this using all array syntax.  First will extend our 1-d coordinate arrays to be 2-d.  Numpy has a function for this (`meshgrid()`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. ... 0. 0. 0.]\n",
      "[0.         0.26179939 0.52359878 0.78539816 1.04719755 1.30899694\n",
      " 1.57079633 1.83259571 2.0943951  2.35619449 2.61799388 2.87979327\n",
      " 3.14159265 3.40339204 3.66519143 3.92699082 4.1887902  4.45058959\n",
      " 4.71238898 4.97418837 5.23598776 5.49778714 5.75958653 6.02138592\n",
      " 6.28318531]\n",
      "[0.0000000e+00 1.5625000e-04 3.1250000e-04 ... 9.9953125e-01 9.9968750e-01\n",
      " 9.9984375e-01]\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0.]\n"
     ]
    }
   ],
   "source": [
    "x2d, y2d = np.meshgrid(x, y, indexing=\"ij\")\n",
    "\n",
    "print(x2d[:,0])\n",
    "print(x2d[0,:])\n",
    "\n",
    "print(y2d[:,0])\n",
    "print(y2d[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time elapsed: 0.011608123779296875 s\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "g2 = np.sin(2.0*np.pi*x2d*y2d)\n",
    "t1 = time.time()\n",
    "print(\"time elapsed: {} s\".format(t1-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "print(np.max(np.abs(g2-g)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to construct a derivative, \n",
    "$$\n",
    "\\frac{d f}{dx}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 2*np.pi, 25)\n",
    "f = np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.26179939, 0.52359878, 0.78539816, 1.04719755,\n",
       "       1.30899694, 1.57079633, 1.83259571, 2.0943951 , 2.35619449,\n",
       "       2.61799388, 2.87979327, 3.14159265, 3.40339204, 3.66519143,\n",
       "       3.92699082, 4.1887902 , 4.45058959, 4.71238898, 4.97418837,\n",
       "       5.23598776, 5.49778714, 5.75958653, 6.02138592, 6.28318531])"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2617993877991494"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5.49778714, 5.75958653])"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[-4:-2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to do this without loops&mdash;we'll use views into arrays offset from one another.  Recall from calculus that a derivative is approximately:\n",
    "$$\n",
    "\\frac{df}{dx} = \\frac{f(x+h) - f(x)}{h}\n",
    "$$\n",
    "Here, we'll take $h$ to be a single adjacent element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  2.58819045e-01,  5.00000000e-01,  7.07106781e-01,\n",
       "        8.66025404e-01,  9.65925826e-01,  1.00000000e+00,  9.65925826e-01,\n",
       "        8.66025404e-01,  7.07106781e-01,  5.00000000e-01,  2.58819045e-01,\n",
       "        1.22464680e-16, -2.58819045e-01, -5.00000000e-01, -7.07106781e-01,\n",
       "       -8.66025404e-01, -9.65925826e-01, -1.00000000e+00, -9.65925826e-01,\n",
       "       -8.66025404e-01, -7.07106781e-01, -5.00000000e-01, -2.58819045e-01,\n",
       "       -2.44929360e-16])"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.58819045e-01,  5.00000000e-01,  7.07106781e-01,  8.66025404e-01,\n",
       "        9.65925826e-01,  1.00000000e+00,  9.65925826e-01,  8.66025404e-01,\n",
       "        7.07106781e-01,  5.00000000e-01,  2.58819045e-01,  1.22464680e-16,\n",
       "       -2.58819045e-01, -5.00000000e-01, -7.07106781e-01, -8.66025404e-01,\n",
       "       -9.65925826e-01, -1.00000000e+00, -9.65925826e-01, -8.66025404e-01,\n",
       "       -7.07106781e-01, -5.00000000e-01, -2.58819045e-01, -2.44929360e-16])"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  2.58819045e-01,  5.00000000e-01,  7.07106781e-01,\n",
       "        8.66025404e-01,  9.65925826e-01,  1.00000000e+00,  9.65925826e-01,\n",
       "        8.66025404e-01,  7.07106781e-01,  5.00000000e-01,  2.58819045e-01,\n",
       "        1.22464680e-16, -2.58819045e-01, -5.00000000e-01, -7.07106781e-01,\n",
       "       -8.66025404e-01, -9.65925826e-01, -1.00000000e+00, -9.65925826e-01,\n",
       "       -8.66025404e-01, -7.07106781e-01, -5.00000000e-01, -2.58819045e-01])"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = x[1] - x[0]\n",
    "dfdx = (f[1:] - f[:-1])/dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.98861593,  0.92124339,  0.79108963,  0.60702442,  0.38159151,\n",
       "        0.13015376, -0.13015376, -0.38159151, -0.60702442, -0.79108963,\n",
       "       -0.92124339, -0.98861593, -0.98861593, -0.92124339, -0.79108963,\n",
       "       -0.60702442, -0.38159151, -0.13015376,  0.13015376,  0.38159151,\n",
       "        0.60702442,  0.79108963,  0.92124339,  0.98861593])"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfdx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the exams, prepare 2 of these exercises. Or only one if you like Q6 (actually, if you go for Q6 feel free to skip a question from another lecture as well). Pick the ones you like the most."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some  of these come from / are inspired from https://github.com/rougier/numpy-100 and http://www.scipy-lectures.org/intro/numpy/exercises.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q1: Fun with arrays\n",
    "\n",
    "**A.**  Create the array: \n",
    "```\n",
    "[[1,  6, 11],\n",
    " [2,  7, 12],\n",
    " [3,  8, 13],\n",
    " [4,  9, 14],\n",
    " [5, 10, 15]]\n",
    "```\n",
    "with out explicitly typing it in.\n",
    "\n",
    "Now create a new array containing only its 2nd and 4th rows.\n",
    "\n",
    "\n",
    "**B.** Create a 2d array with `1` on the border and `0` on the inside, e.g., like:\n",
    "```\n",
    "1 1 1 1 1\n",
    "1 0 0 0 1\n",
    "1 0 0 0 1\n",
    "1 1 1 1 1\n",
    "```\n",
    "\n",
    "Do this using array slice notation to let it work for an arbitrary-sized array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q2: Histograms\n",
    "\n",
    "Here we will read in columns of numbers from a file and create a histogram, using NumPy routines.  Make sure you have the data file\n",
    "\"`sample.txt`\" in the same directory as this notebook. You download it from  https://github.com/raw/sbu-python-summer/python-tutorial/master/day-3/sample.txt (and use python to download a file!)\n",
    "\n",
    "  * Use `np.loadtxt()` to read this file in.  \n",
    "\n",
    "  * Next, use `np.histogram()` to create a histogram array.  The output returns both the count and an array of edges.\n",
    "  \n",
    "  * Finally, loop over the bins and print out the bin center (averaging the left and right edges of the bin) and the count for that bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q3: Are you faster than numpy?\n",
    "\n",
    "Numpy of course has a standard deviation function, `np.std()`, but here we'll write our own that works on a 1-d array (vector).  The standard\n",
    "deviation is a measure of the \"width\" of the distribution of numbers\n",
    "in the vector.\n",
    "\n",
    "Given an array, $a$, and an average $\\bar{a}$, the standard deviation\n",
    "is:\n",
    "$$\n",
    "\\sigma = \\left [ \\frac{1}{N} \\sum_{i=1}^N (a_i - \\bar{a})^2 \\right ]^{1/2}\n",
    "$$\n",
    "\n",
    "Write a function to calculate the standard deviation for an input array, `a`:\n",
    "\n",
    "  * First compute the average of the elements in `a` to define $\\bar{a}$\n",
    "  * Next compute the sum over the squares of $a - \\bar{a}$\n",
    "  * Then divide the sum by the number of elements in the array\n",
    "  * Finally take the square root (you can use `np.sqrt()`)\n",
    "  \n",
    "Test your function on a random array, and compare to the built-in `np.std()`. Check the runtime as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q4: Bootstraps\n",
    "\n",
    "Bootstrapping is a prominent statistical technique to estimate errors when you're stuck and don't know what else to o: https://erikjanvankesteren.nl/blog/2018-03-01-baron-bootstrap.html\n",
    "\n",
    "Generate a 1D array (say from a Gaussian, but then add some outliers for fun). Compute the bootstrapped 95% confidence intervals for the mean of that array (i.e., resample the elements of an array with replacement N times, compute the mean of each sample, and then compute percentiles over the means). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q5: Einstein summation\n",
    "\n",
    "einsum is a powerful (but often painful) numpy thing:\n",
    "- https://numpy.org/doc/stable/reference/generated/numpy.einsum.html\n",
    "- https://stackoverflow.com/questions/26089893/understanding-numpys-einsum\n",
    "\n",
    "Take 2 vectors A and B. Write the einsum equivalent of inner, outer, sum, and mul function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q6: Conway's Game of Life\n",
    "\n",
    "**Exercise**: Code up Conway's Game of Life using numpy \n",
    "\n",
    "The Game of Life is a cellular automaton devised by mathematician John Horton Conway in 1970. It is a zero-player game, meaning that its evolution is determined by its initial state, requiring no further input. One interacts with the Game of Life by creating an initial configuration and observing how it evolves. It is Turing complete and can simulate a universal constructor or any other Turing machine.\n",
    "\n",
    "https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life\n",
    "\n",
    "The Game of Life is *really* (really, really) cool. There are just four extremely simple rules, and these result in an immense richness of behaviour and complexity.\n",
    "\n",
    "https://www.youtube.com/watch?v=C2vgICfQawE&t=221s&ab_channel=RationalAnimations\n",
    "\n",
    "https://www.youtube.com/watch?v=jvSp6VHt_Pc&ab_channel=TheDevDoctor\n",
    "\n",
    "Here some web apps to play:\n",
    "\n",
    "https://conwaylife.com/\n",
    "\n",
    "https://playgameoflife.com/\n",
    "\n",
    "Some computational hints:\n",
    "\n",
    "https://blog.datawrapper.de/game-of-life/\n",
    "\n",
    "For instance, here is a Game-of-Life structure that sends a message at fixed intervals (that little spaceship leaving toward the bottom right)\n",
    "\n",
    "![](https://blog.datawrapper.de/wp-content/uploads/2021/06/game-of-life-loop-cropped.gif)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}