Notebook10_NumpyScipy_part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d95c49d43afa9cf5",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Part 1: Intro to Numpy/Scipy\n",
    "\n",
    "[Numpy](http://www.numpy.org/) is a Python module that provides fast primitives for multidimensional arrays. It's well-suited to implementing numerical linear algebra algorithms, and for those can be much faster than Python's native list and dictionary types when you only need to store and operate on numerical data.\n",
    "\n",
    "Some of the material from this lesson is copied from the following, and more comprehensive, tutorial: [link](http://www.scipy-lectures.org/intro/numpy/index.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-8d1656231e9af5be",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Quick demo.** The recommended importing idiom is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-846ef2892d749d83",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-5a7dbee2c730c482",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Creating a simple numpy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1,2,3,4])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why bother with Numpy? A motivating example\n",
    "\n",
    "We already have lists and dictionary types, which are pretty easy to use and very flexible. So why bother with this special type?\n",
    "\n",
    "**Exercise 0** (ungraded). One reason to consider Numpy is that it \"can be much faster,\" as noted above. But how much faster is that? Run the experiment below to see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "302 ms ± 5.75 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "L = range(n)\n",
    "%timeit [i**2 for i in L]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(10) # Moral equivalent to `range`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "901 µs ± 234 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "A = np.arange(n)\n",
    "%timeit A**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-5783a1fb911dd611",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Creating multidimensional arrays\n",
    "\n",
    "Beyond simple arrays, Numpy supports multidimensional arrays. To do more than one dimension, call `numpy.array()` but nest each new dimension within a list. It's easiest to see by example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n"
     ]
    }
   ],
   "source": [
    "# Create a two-dimensional array of size 3 rows x 4 columns:\n",
    "B = np.array([[0, 1, 2, 3],\n",
    "              [4, 5, 6, 7],\n",
    "              [8, 9, 10, 11]])\n",
    "\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "(3, 4)\n",
      "3\n"
     ]
    }
   ],
   "source": [
    "print(B.ndim) # Gives the number of dimensions (2 dimensions)\n",
    "print(B.shape) # gives the nXm \"shape\" of the matrix  (height, length) (3, 4)\n",
    "print(len (B)) # gives the number of rows or lists needed to make matrix (3 rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0  1  2  3]\n",
      "  [ 4  5  6  7]\n",
      "  [ 8  9 10 11]]\n",
      "\n",
      " [[12 13 14 15]\n",
      "  [16 17 18 19]\n",
      "  [20 21 22 23]]]\n",
      "3\n",
      "(2, 3, 4)\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "C1 = [[0, 1, 2, 3],\n",
    "      [4, 5, 6, 7],\n",
    "      [8, 9, 10, 11]]\n",
    "\n",
    "C2 = [[12, 13, 14, 15],\n",
    "      [16, 17, 18, 19],\n",
    "      [20, 21, 22, 23]]\n",
    "\n",
    "C = np.array([C1, C2])  # create a two dimensional matrix\n",
    "\n",
    "print(C)\n",
    "print(C.ndim)         # gives number of dimentions of the matrix (3 dimensions)\n",
    "print(C.shape)        # gives the actual dimensions (width, height, length)  (2, 3, 4)\n",
    "print(len (C))        # gives the number of items needed to make the shape (2 matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are routines for creating various kinds of structured matrices as well, which are similar to those found in [MATLAB](http://www.mathworks.com/products/matlab/) and [Octave](https://www.gnu.org/software/octave/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.zeros((3, 4)))  # zero matrix, give the dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.ones((3, 4)))   # matrix of only ones, give the dimentions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.  0.]\n",
      " [ 0.  1.  0.]\n",
      " [ 0.  0.  1.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.eye(3))      # identity nXn matrix (Gausian reduced row echelon) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 0 0]\n",
      " [0 2 0]\n",
      " [0 0 3]]\n"
     ]
    }
   ],
   "source": [
    "print(np.diag([1, 2, 3]))      # nXn matrix, give the diagonal elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also create empty (uninitialized) arrays. What does the following produce?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.empty((3, 4)) # An \"empty\" 3 x 4 matrix\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1** (ungraded). The following code creates an identity matrix in two different ways, which are found to be equal according to the assertion. But in fact there is a subtle difference between the `I` and `I_u` matrices created below; can you spot it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> I = eye(n):\n",
      "[[ 1.  0.  0.]\n",
      " [ 0.  1.  0.]\n",
      " [ 0.  0.  1.]]\n",
      "\n",
      "==> u:\n",
      " [1, 1, 1]\n",
      "==> I_u = diag (u):\n",
      " [[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n"
     ]
    }
   ],
   "source": [
    "n = 3\n",
    "I = np.eye(n)\n",
    "\n",
    "print(\"==> I = eye(n):\")  # uses float values\n",
    "print(I)\n",
    "\n",
    "u = [1] * n\n",
    "I_u = np.diag(u)\n",
    "\n",
    "print(\"\\n==> u:\\n\", u)\n",
    "print(\"==> I_u = diag (u):\\n\", I_u)    # uses integer values\n",
    "\n",
    "assert np.all(I_u == I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false
   },
   "source": [
    "**Answer.** Give this some thought before you read the answer that follows!\n",
    "\n",
    "The difference is in the element types. The `eye()` function returns an identity matrix and uses a floating-point type as the element type. By contrast, `diag()`, which expects a list of initializer values upon input, derives the element type from that input. In this case, `u` contains values that will be stored as integers; therefore, `diag()` constructs its output assuming integer elements.\n",
    "\n",
    "Try running `print(I_u.dtype)` and `print(I.dtype)` to confirm that these element types differ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing and slicing\n",
    "\n",
    "The usual 0-based slicing and indexing notation you know and love from lists is also supported for Numpy arrays. In the multidimensional case,  including their natural multidimensional analogues with index ranges separated by commas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0  1  2  3]\n",
      "  [ 4  5  6  7]\n",
      "  [ 8  9 10 11]]\n",
      "\n",
      " [[12 13 14 15]\n",
      "  [16 17 18 19]\n",
      "  [20 21 22 23]]]\n"
     ]
    }
   ],
   "source": [
    "# Recall: C\n",
    "print (C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What part of C will the following slice extract? Run the code to find out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 8  9 10 11]\n"
     ]
    }
   ],
   "source": [
    "print (C[0, 2, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What will the following slice return? Run the code to find out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[15 14 13 12]\n"
     ]
    }
   ],
   "source": [
    "print (C[1, 0, ::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** (5 points). Consider the following $6 \\times 6$ matrix, which has 4 different subsets highlighted.\n",
    "\n",
    "<img src=\"https://github.com/cse6040/labs-fa17/raw/master/lab10-numpy/slicing-exercise.png\" alt=\"Exercise: Extract these slices\" width=\"240\">\n",
    "\n",
    "For each subset illustrated above, write an indexing or slicing expression that extracts the subset. Store the result of each slice into `Z_green`, `Z_red`, `Z_orange`, and `Z_cyan`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "11e5e133ff5caee1a7562795f3d26d2e",
     "grade": false,
     "grade_id": "create_Z",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "Z= np.array([[0,1,2,3,4,5],[10,11,12,13,14,15],[20,21,22,23,24,25],[30,31,32,33,34,35],[40,41,42,43,44,45],[50,51,52,53,54,55]])\n",
    "\n",
    "# Construct `Z_green`, `Z_red`, `Z_orange`, and `Z_cyan`:\n",
    "# YOUR CODE HERE\n",
    "Z_red = Z[:, 2]\n",
    "Z_orange = Z[0, 3:5]\n",
    "Z_cyan = Z[4:, 4:]\n",
    "Z_green = Z[2::2, 0::2]  # start on :: skip this many  (Ex. 0::2 start at 0 and skip every other element)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5ba2d769c7f7534e17eab32ec1e3d2c8",
     "grade": true,
     "grade_id": "check_Z",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> Z:\n",
      " [[ 0  1  2  3  4  5]\n",
      " [10 11 12 13 14 15]\n",
      " [20 21 22 23 24 25]\n",
      " [30 31 32 33 34 35]\n",
      " [40 41 42 43 44 45]\n",
      " [50 51 52 53 54 55]]\n",
      "\n",
      "==> Orange slice:\n",
      " [3 4]\n",
      "\n",
      "==> Red slice:\n",
      " [ 2 12 22 32 42 52]\n",
      "\n",
      "==> Cyan slice:\n",
      " [[44 45]\n",
      " [54 55]]\n",
      "\n",
      "==> Green slice:\n",
      " [[20 22 24]\n",
      " [40 42 44]]\n",
      "\n",
      "(Passed!)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `check_Z`\n",
    "\n",
    "print(\"==> Z:\\n\", Z)\n",
    "assert (Z == np.array([np.arange(0, 6),\n",
    "                       np.arange(10, 16),\n",
    "                       np.arange(20, 26),\n",
    "                       np.arange(30, 36),\n",
    "                       np.arange(40, 46),\n",
    "                       np.arange(50, 56)])).all()\n",
    "\n",
    "print(\"\\n==> Orange slice:\\n\", Z_orange)\n",
    "assert (Z_orange == np.array ([3, 4])).all()\n",
    "\n",
    "print(\"\\n==> Red slice:\\n\", Z_red)\n",
    "assert (Z_red == np.array ([2, 12, 22, 32, 42, 52])).all()\n",
    "\n",
    "print(\"\\n==> Cyan slice:\\n\", Z_cyan)\n",
    "assert (Z_cyan == np.array ([[44, 45], [54, 55]])).all()\n",
    "\n",
    "print(\"\\n==> Green slice:\\n\", Z_green)\n",
    "assert (Z_green == np.array ([[20, 22, 24], [40, 42, 44]])).all()\n",
    "\n",
    "print(\"\\n(Passed!)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slices are views\n",
    "\n",
    "To help save memory, when you slice a Numpy array, you are actually creating a _view_ into that array. That means modifications through the view will modify the original array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> Recall C: (2, 3, 4)\n",
      "[[[ 0  1  2  3]\n",
      "  [ 4  5  6  7]\n",
      "  [ 8  9 10 11]]\n",
      "\n",
      " [[12 13 14 15]\n",
      "  [16 17 18 19]\n",
      "  [20 21 22 23]]]\n"
     ]
    }
   ],
   "source": [
    "print(\"==> Recall C: %s\" % str(C.shape))\n",
    "print(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> C_view: (2, 2)\n",
      "[[13 15]\n",
      " [21 23]]\n"
     ]
    }
   ],
   "source": [
    "C_view = C[1, 0::2, 1::2] # Question: What does this produce?\n",
    "print (\"==> C_view: %s\" % str (C_view.shape))\n",
    "print (C_view)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-23 -21]\n",
      " [-15 -13]]\n"
     ]
    }
   ],
   "source": [
    "C_view[:, :] = -C_view[::-1, ::-1] # Question: What does this do?\n",
    "print (C_view)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[  0   1   2   3]\n",
      "  [  4   5   6   7]\n",
      "  [  8   9  10  11]]\n",
      "\n",
      " [[ 12 -23  14 -21]\n",
      "  [ 16  17  18  19]\n",
      "  [ 20 -15  22 -13]]]\n"
     ]
    }
   ],
   "source": [
    "print (C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can force a copy using the `.copy()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> C_view:\n",
      "[[-23 -21]\n",
      " [-15 -13]]\n",
      "\n",
      "==> C_copy:\n",
      "[[13 15]\n",
      " [21 23]]\n"
     ]
    }
   ],
   "source": [
    "C_copy = C[1, 0::2, 1::2].copy ()\n",
    "C_copy[:, :] = -C_copy[::-1, ::-1]\n",
    "\n",
    "print (\"==> C_view:\")\n",
    "print (C_view)\n",
    "\n",
    "print (\"\\n==> C_copy:\")\n",
    "print (C_copy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to check whether two Numpy array variables point to the same object, you can use the `numpy.may_share_memory()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C and C_view share memory: True\n",
      "C and C_copy share memory: False\n"
     ]
    }
   ],
   "source": [
    "print (\"C and C_view share memory: %s\" % np.may_share_memory (C, C_view))\n",
    "print (\"C and C_copy share memory: %s\" % np.may_share_memory (C, C_copy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indirect addressing\n",
    "\n",
    "Two other common ways to index a Numpy array are to use a boolean mask or to use a set of integer indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10  3  8  0 19 10 11  9 10  6  0 12  7 14 17]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(3)\n",
    "x = np.random.randint(0, 20, 15) # 15 random ints in [0, 20)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  9 10  7]\n"
     ]
    }
   ],
   "source": [
    "# Pull out an arbitrary subset of elements\n",
    "inds = np.array([3, 7, 8, 12])\n",
    "print(x[inds])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before looking at how to use a boolean mask for indexing, let's create one.\n",
    "\n",
    "**Exercise 3** (1 point). Given the input array, `x[:]`, above, create an array, `mask_mult_3[:]` such that `mask_mult_3[i]` is true only if `x[i]` is a positive multiple of 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8e07e2c94176962381e0f8b39cf70602",
     "grade": false,
     "grade_id": "mask_mult_3",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False,  True, False,  True, False, False, False,  True, False,\n",
       "        True,  True,  True, False, False, False], dtype=bool)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# YOUR CODE HERE\n",
    "mask_mult_3 = x%3 == 0\n",
    "mask_mult_3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "9a226b25ac1f6e3e9dedc154362f2342",
     "grade": true,
     "grade_id": "mask_mult_3_test",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x: [10  3  8  0 19 10 11  9 10  6  0 12  7 14 17]\n",
      "mask_mult_3: [False  True False  True False False False  True False  True  True  True\n",
      " False False False]\n",
      "==> x[mask_mult_3]: [ 3  0  9  6  0 12]\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `mask_mult_3_test`\n",
    "\n",
    "print (\"x:\", x)\n",
    "print (\"mask_mult_3:\", mask_mult_3)\n",
    "print (\"==> x[mask_mult_3]:\", x[mask_mult_3])\n",
    "\n",
    "inv_mask_mult_3 = np.invert (mask_mult_3)\n",
    "assert ((x[mask_mult_3] % 3) == np.zeros (sum (mask_mult_3))).all ()\n",
    "assert (((x[inv_mask_mult_3] % 3) != np.zeros (sum (inv_mask_mult_3))) | (x[inv_mask_mult_3] == 0)).all ()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4** (3 points). Complete the prime number sieve algorithm, which is illustrated below.\n",
    "\n",
    "<img src=\"https://github.com/cse6040/labs-fa17/raw/master/lab10-numpy/prime-sieve.png\" alt=\"Exercise: Extract these slices\" width=\"480\">\n",
    "\n",
    "That is, given a positive integer $n$, the algorithm iterates from $i \\in \\{2, 3, 4, \\ldots, \\left\\lfloor\\sqrt{n}\\right\\rfloor\\}$, repeatedly \"crossing out\" values that are strict multiples of $i$. \"Crossing out\" means maintaining an array of, say, booleans, and setting values that are multiples of $i$ to `False`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "2eeb476bfd7cfcf266c8af9efe5cc8de",
     "grade": false,
     "grade_id": "prime_sieve",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> Primes through 20:\n",
      " [ 2  3  5  7 11 13 17 19 23 29 31 37]\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt\n",
    "import numpy as np\n",
    "\n",
    "def sieve(n):\n",
    "    \"\"\"\n",
    "    Returns the prime number 'sieve' shown above.\n",
    "    \n",
    "    That is, this function returns an array `X[0:n+1]`\n",
    "    such that `X[i]` is true if and only if `i` is prime.\n",
    "    \"\"\"\n",
    "    is_prime = np.empty(n+1, dtype=bool) # an empty array to be seeded with bool values T/F.  n+1 to include value 0\n",
    "\n",
    "    # Initial values\n",
    "    is_prime[0:2] = False # {0, 1} are _not_ considered prime\n",
    "    is_prime[2:] = True # All other values might be prime\n",
    "\n",
    "    # Implement the sieving loop\n",
    "    # YOUR CODE HERE\n",
    "    for i in range(2, int(sqrt(n))):\n",
    "            is_prime[2*i::i] = False  # the \"sieve\", the corresponding indices that are not prime to false.  When i = 2,\n",
    "                                        # [4::2] starts at 4 and every second value (evens) are marked false for not prime\n",
    "                                        # [6::3] starts at 6 and every third value (multi of 3s) are false etc. \n",
    "            \n",
    "    return is_prime                     # returns a, array of booleans where non-primes are set to false and primes are true\n",
    "\n",
    "# Prints your primes\n",
    "print(\"==> Primes through 20:\\n\", np.nonzero(sieve(40))[0])  # returns the index values where it is true (primes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2]\n"
     ]
    }
   ],
   "source": [
    "# m = np.empty(20+1, dtype=bool)\n",
    "# m[0:2]= False\n",
    "# m[2:] = True\n",
    "# m[4::2] = False\n",
    "# m[6::3] = False\n",
    "# m[8::4] = False\n",
    "# m\n",
    "# # m = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]\n",
    "# # m[4::2] = False\n",
    "# # np.nonzero(m[4::2])#[0]\n",
    "print(np.nonzero([True, False, True, False])[0])  # prints index values where it is true (or if the value is nonzero)\n",
    "# sieve(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "99a64fc0190abc8cdba6386398befdb5",
     "grade": true,
     "grade_id": "prime_sieve_test",
     "locked": true,
     "points": 3,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# Test cell: `prime_sieve_test`\n",
    "\n",
    "is_prime = sieve(20)\n",
    "assert len (is_prime) == 21\n",
    "assert (is_prime == np.array([False, False, True, True, False, True, False, True, False, False, False, True, False, True, False, False, False, True, False, True, False])).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Fin!** If you've reached this point and everything executed without error, you can submit this part and move on to the next one."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": [],
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}