Numpy_scipy_with_exercise_answers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DyznzL8W2WDt"
   },
   "source": [
    "NumPy (SciPy) Workshop - 20 July 2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xJ1k-vir2pTm"
   },
   "source": [
    "# Import NumPy <br>\n",
    "NumPy is an external library; traditionally imported as 'np'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "AYgwYUQw2cZw",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Import\n",
    "import numpy as np\n",
    "import sys #Also import sys to check python version\n",
    "\n",
    "print (\"Numpy version:\", np.__version__)\n",
    "print (\"Python version:\", sys.version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is NumPy (and SciPy)?\n",
    "NumPy is the fundamental package for scientific computing with Python. It contains among other things: <br> <br>\n",
    "\n",
    "-a powerful N-dimensional *array* object <br>\n",
    "-fast, parallelized math functions <br>\n",
    "-tools for integrating C/C++ and Fortran code <br>\n",
    "-useful linear algebra, Fourier transform, and random number capabilities <br> <br>\n",
    "\n",
    "NumPy was originally forked (derived) from SciPy, so that one could avoid importing the large SciPy package just to get an array object; therefore, it contains many redundancies with SciPy, but is also incomplete with regards to some functionality. <br> <br>\n",
    "\n",
    "Broadly speaking, NumPy is the basis for the scientific computing and data science stacks in Python, including SciPy and Pandas. <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "CEe2L6ta3cDl"
   },
   "source": [
    "# Arrays <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3RSNQSstmAvg"
   },
   "source": [
    "The fundamental 'object' in NumPy is the *array* <br>\n",
    "Compared to default Python, *lists* $\\approx$ *arrays*, except every element of an *array* **must be of the same type** <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "a-GpWzzN4T7c",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Arrays can be created from lists...\n",
    "l1 = [1, 2, 3, 4]\n",
    "a1 = np.array(l1)\n",
    "a1.reshape(4,1)\n",
    "print('l1 is of ' + str(type(l1))[1:-1])\n",
    "print('a1 is of ' + str(type(a1))[1:-1])\n",
    "\n",
    "# Or lists can be created from arrays...\n",
    "a2 = np.array([1, 2, 3, 4])\n",
    "l2 = list(a2)\n",
    "a2.reshape(1,4)\n",
    "print('l2 is of ' + str(type(l2))[1:-1])\n",
    "print('a2 is of ' + str(type(a2))[1:-1])\n",
    "\n",
    "print('It is ' + str(np.all(a1 == a2)) + #np.all is described later\n",
    "      ' that a1 is identical to a2') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "a5Ix2_pW-Ws0"
   },
   "source": [
    "But, again, all elements of an *array* must be of the same type! <br>\n",
    "\n",
    "So, if you feed an array something of mixed data type, it will force everything into one type and normal math operations may not function\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "hmJBY5Qo-xJ7",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "l3 = ['test', 2, 3, 4] # a list with mixed type\n",
    "a3 = np.array(l3) # a numpy array from the list of mixed type\n",
    "print(str(a3[1]*2) + ' does not equal ' + str(l3[1]*2))\n",
    "print('This is because in a NumPy array, \"2\" is being treated as ' +\n",
    "      str(type(a3[1]))) \n",
    "print('In the list, \"2\" remains ' +\n",
    "      str(type(l3[1])))\n",
    "a3 #full output, unsuppressed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7Jfbdmx3Agn9"
   },
   "source": [
    "Note the 'dtype' output in the full description of the 'a3' array above <br>\n",
    "'dtype' = data type <br>\n",
    "'dtype' is both a standalone class object in numpy (see [Data Type Objects](https://docs.scipy.org/doc/numpy-1.14.0/reference/arrays.dtypes.html) and [np.dtype](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.dtype.html) ) and an attribute of an *array* ('a4.dtype', for ex) describing the **datatype of the elements in that array**<br>\n",
    "There are various defaults for assigning datatypes to arrays without returning errors (as in the above, a mixed *string* and *int* list was converted to unicode; two floats of different accuracy will be converted to the one of higher accuracy, etc.)<br>\n",
    "See [Data Types](https://docs.scipy.org/doc/numpy-1.13.0/user/basics.types.html) for a concise list of datatypes <br> <br>\n",
    "\n",
    "If you want to be explicit about your *array*  datatype, you can call *dtype* when forming an *array*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "i_y7fMezC27v",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "a4 = np.array([1, 2, 3, 4], dtype=float) #note: could be 'float64'\n",
    "print('It is ' + str(type(a4) == type(a1)) +\n",
    "      ' that both a4 and a1 are arrays of ' + str(type(a4)))\n",
    "print('The elements in the original array a1 were ' + str((a1[0]).dtype))\n",
    "print('The elements in a4 were \"cast\" into ' + str((a4[0]).dtype))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Oru6qQPhKCl9"
   },
   "source": [
    "Note: though at times one can simply recast the dtype with (for example)\n",
    "\n",
    "```\n",
    "a4.dtype = 'np.int'\n",
    "```\n",
    "The safer and more consistent way to recast dtype is with the *array* method *array*.astype\n",
    "```\n",
    "a4.astype(np.int)\n",
    "```\n",
    "as done in the code below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "s1tXN67QKAaf",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "a4 = np.array([1, 2, 3, 4], dtype=float)\n",
    "print('a4 is ' + str(a4))\n",
    "print('elements of a4 are of type ' + str(a4.dtype))\n",
    "print('Though the elements can be cast back into the original form with dtype,')\n",
    "print('in this case they\\'re re-written (as a copy) onto the original with astype')\n",
    "a4 = a4.astype(np.int)\n",
    "print('a4 is now ' + str(a4))\n",
    "print('elements of a4 are now of type ' + str(a4.dtype))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "LKsKm-00Low-"
   },
   "source": [
    "The above example with ```array.astype(...)``` utilizes the concept of 'copying' an array<br>\n",
    "It is worth noting that,  like Python lists, **arrays are not immutable** <br>\n",
    "When an new array is made from another array, for example, it points directly to that original array <br>\n",
    "And when that new array is modified, that old array is modified as well <br>\n",
    "You can use [*array*.copy()](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.copy.html) to be safer <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "HQ_91SilNhX1",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "a5 = a4\n",
    "print('The array a4 is')\n",
    "print(a4)\n",
    "a5[-1] = int(5)\n",
    "print('After modifying the last element of a5, both a5 and a4 are')\n",
    "print(a4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "HQ_91SilNhX1",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "a4 = np.array([1, 2, 3, 4], dtype=int)\n",
    "a5 = a4.copy()\n",
    "a5[-1] = int(5)\n",
    "print('After modifying a5, a copy of a4, its values are ')\n",
    "print(a5)\n",
    "print('BUT a4 is the same')\n",
    "print(a4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JzJRdtSJOpmL"
   },
   "source": [
    "Okay-- <br>\n",
    "Everything we have done so far with one-dimensional *arrays* is basically identical to Python *lists*...<br>\n",
    "One of the well-known benefits of NumPy, however, is more computationally efficient operations (it better exploits *C*). This extends even to 1-D arrays vs lists--but it is most obvious with large arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example highlighting the difference between Python's built-in ```range()``` and ```np.arange()```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "jmHwWfE_Ou6S",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#modified from https://stackoverflow.com/questions/10698858/built-in-range-or-numpy-arange-which-is-more-efficient\n",
    "\n",
    "time_size = int(1E6)\n",
    "\n",
    "# Loop index example\n",
    "print('range() loop-index example timed')\n",
    "%timeit for x in range(time_size): x * 2\n",
    "#out: 1 loop, best of 3: 287 ms per loop\n",
    "print('np.arange() loop-index example timed')\n",
    "%timeit for x in np.arange(time_size): x * 2\n",
    "#out: 1 loop, best of 3: 162 ms per loop\n",
    "\n",
    "# Math example\n",
    "print('range() math example timed')\n",
    "%timeit list(range(time_size)) * 2\n",
    "#out: 10 loops, best of 3: 49.7 ms per loop\n",
    "print('np.arange() math example timed')\n",
    "%timeit np.arange(time_size) * 2\n",
    "#out: 1000 loops, best of 3: 1.87 ms per loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1zXxRoaLLPye"
   },
   "source": [
    "So, [np.arange](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.arange.html) is **not always faster** Python's built-in [range](https://docs.python.org/3/library/functions.html#func-range)<br>, even when used just for building a loop index--use it whenever you have already imported NumPy! <br>\n",
    "When used for a straight-forward mathematical operation, it is **way faster** <br> \n",
    "Note: in the above example,\n",
    "\n",
    "```\n",
    "%timeit\n",
    "```\n",
    "\n",
    "is 'iPython magic' -- these functions are beyond the scope of this workshop, but as they have default and modifiable integration into Jupyter notebooks, they are worth reading about. <br>\n",
    "<br>\n",
    "For more information, see [Built-in Magic Commands](http://ipython.readthedocs.io/en/stable/interactive/magics.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HhM2bfZuP1DA"
   },
   "source": [
    "## Exercise 1 <br>\n",
    "Using [np.linspace](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.linspace.html), create an *array* with *dtype* ```int``` of length 1 million (1E6) evenly spaced between '0' and '1E6,'' with '0' omitted. That is, $[1, 2, 3, ... 1*10^6]$ <br>\n",
    "\n",
    "Bonus: use *%timeit* to compare how long it takes to create this *array* VS. how long it takes to create a *list* with the same elements from Python's built-in *range* function.\n",
    "\n",
    "Hint: You may need to index for the *range* function to omit '0'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "dtOT6VjmQQk6",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "million = np.linspace(1,time_size,num=time_size, dtype=int)\n",
    "print((million))\n",
    "\n",
    "%timeit np.linspace(1,time_size,num=time_size, dtype=int)\n",
    "#out: 100 loops, best of 3: 4.33 ms per loop\n",
    "\n",
    "size_thousand = 1001\n",
    "%timeit list(range(time_size))[1:]\n",
    "#out: 10 loops, best of 3: 41.6 ms per loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2CFYEAjtv7-x"
   },
   "source": [
    "## NumPy Speed-up, In General\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Nxo8TYJmwMFG"
   },
   "source": [
    "More generically, NumPy will save computational (and coding) time if it replaces coding *loops* with *vectorized* (or any higher-dimensional) math.<br> <br>\n",
    "Here is an example comparing NumPy element-wise multiplication vs. loop-based indexing for element-wise multiplication "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "kWLTo8cCxMBS"
   },
   "outputs": [],
   "source": [
    "#modified from https://stackoverflow.com/questions/47755442/what-is-vectorization\n",
    "a40 = np.array([1., 2., 3.])\n",
    "print('The array')\n",
    "print(a40)\n",
    "# By default all arithmetic operators (like *) in NumPy are element-wise\n",
    "a4040 = a40 * a40\n",
    "print('The array a40 * a40 (squared) is ')\n",
    "print(a4040)\n",
    "print('%timeit of the operation is')\n",
    "%timeit a40 * a40\n",
    "# out: 597 nsec (with cached)\n",
    "# Compare this to python math with loops to complete element wise\n",
    "l40 = a40.tolist()\n",
    "sq = [0., 0., 0.] # pre-allocate empty list (for speed)\n",
    "def square(lst,empty):\n",
    "  for i in np.arange(len(l40)): #use np.arange, even (still slower)!\n",
    "    empty[i] = l40[i] * l40[i]\n",
    "  return empty\n",
    "sq = square(l40,sq)\n",
    "print('The array l40[i] * l40[i] for i in len(l40) (squared) is ')\n",
    "print(sq)\n",
    "%timeit square(l40,sq)\n",
    "#out: 2.5 microsec (with cached)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NZ8QgAxv6AM9"
   },
   "source": [
    "So in general, if you want to speed up your code with NumPy, **avoid loops at all cost**. This includes list-comprehensions (which are loops)--these are included in this workshop for brevity, not for speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FsnrfyCT4PoM"
   },
   "source": [
    "Note: **Do not expect a non-element-wise operator from the default arithmetic operators.** <br>\n",
    "A common confusion with non-element-wise operations is regarding ```*``` vs. the *dot product*, which is actually given by [np.dot(...)](https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html); <br>\n",
    "However, since row and column dimensions are important for linear opertors, you should first understand NumPy array dimensionality (that is, *shape*) to use *np.dot(...)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3U1a21ojRNAM"
   },
   "source": [
    "# Multidimensional Arrays (and Matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WOue93dsqXaV"
   },
   "source": [
    "You may have noticed earlier, that the *type* of an *array* is\n",
    "\n",
    "```\n",
    "numpy.ndarray\n",
    "```\n",
    "What is an 'ndarray'?<br>\n",
    "This refers to NumPy's exceptional ability to build and perform math on *arrays* of arbitrary (*N*) dimensionality. <br>\n",
    "In more concrete terms: <br>\n",
    "a vector is a 1-D array (as we have been working with) <br>\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "G_9HLGnWriXY"
   },
   "outputs": [],
   "source": [
    "a6 = np.array([1, 2, 3, 4])\n",
    "print(a6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "r7iCy1b7vX51"
   },
   "source": [
    "Note that every *array* has a *shape* attribute (tuple) that is tied to it. <br>\n",
    "By default, the dimensionality of a given array is as short as possible (1's omitted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "DSOouJcXvgl3"
   },
   "outputs": [],
   "source": [
    "print('The shape of a6 is')\n",
    "print(a6.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_Yp8ClDPtRTW"
   },
   "source": [
    "In relation to 2-D matrix math (linear algebra): this 1-D array can be either a column vector or row vector <br>\n",
    "And the most generic  way to enforce the direction of the 1-D vector in 2-D space is with [*np.reshape(...)*](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.reshape.html) or with [*array.reshape(...)*](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.ndarray.reshape.html#numpy.ndarray.reshape); the latter having slightly easier syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "chHqbH36t7-G"
   },
   "outputs": [],
   "source": [
    "print('Again, the shape of a6 is')\n",
    "print(a6.shape)\n",
    "a6 = a6.reshape(4,1) # array.reshape automatically acts on a6\n",
    "print('After reshaping to (4,1), a6 prints like this')\n",
    "print(a6)\n",
    "a6 = np.reshape(a6,(1,4)) #np.reshape is wordy, must be directed at a6\n",
    "print('After reshaping to (1,4), a6 prints like this')\n",
    "print(a6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "R4XE1p4f3h8f"
   },
   "source": [
    "# Note: Broadcasting, a perplexing default for element-wise operators\n",
    "Now that we've introduced the concept of *array* dimensionality and *shape*, it is worth mentioning that NumPy's ```*``` (and other element-wise operators) follow a strange rule called [broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) when given two arrays of differing *shapes*. <br>\n",
    "Of course, element-wise operations can only work if the shape of the *arrays* is the same in each dimension. If they are not, it will throw a \n",
    "\n",
    "```\n",
    "ValueError: frames are not aligned\n",
    "```\n",
    "\n",
    "\n",
    "However, **if one of the *array* dimensions is 1 (and the other is greater than 1), the array will be 'stretched' (*tiled*) to match the greater dimension** and the output will be of larger dimensionality than you might have anticipated!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "hk_X8Kqv3hXx"
   },
   "outputs": [],
   "source": [
    "a6_v = a6.reshape(4,1)\n",
    "a6_h = a6_v.T # .T short-hand is an 'attribute' of an array class, described below\n",
    "print('Because of broadcasting, the output is')\n",
    "print(a6_h * a6_v)\n",
    "print('instead of (as you might have been trying to do)')\n",
    "print(a6_h * a6_v.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aXDm_RLmfGLf"
   },
   "source": [
    "## Aside: class attributes are useful; example Matrix Class\n",
    "Note, in the above example, [*array.T*](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ndarray.T.html#numpy.ndarray.T) is an attribute of arrays utilized to perform a transpose (swapping of dimensions) as a short-hand. <br>\n",
    "This will only work *after* the second dimension has been defined by *array.reshape* <br>\n",
    "These attributes are powerful shorthand (we've already used *dtype* and *shape*) for items of class type *ndarray* <br>\n",
    "Another example attribute of *ndarray*s is [*array.real*](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ndarray.real.html#numpy.ndarray.real) , which saves time when working with complex numbers in taking only the 'real' part of array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "DdH1DgYDhr0q"
   },
   "outputs": [],
   "source": [
    "a7 = np.array([3, 1, 4+0.1j, 1, 5-2j, 9, 2, 6+2j, 5])\n",
    "print(a7)\n",
    "print(a7[0].dtype)\n",
    "print(a7.real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_0dOySNuhsAj"
   },
   "source": [
    "Similarly, in relation to matrix math (linear algebra), there is a special *array* class *matrix* which may be called to form a 2-D matrix with [np.mat](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.mat.html) <br>\n",
    "This [*matrix*](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.matrix.html#numpy.matrix) class will carry its own shorthand attributes that are useful <br>\n",
    "\n",
    "*matrix.T* = transpose\n",
    "\n",
    "*matrix.H* = Hermitian transpose (transpose with complex conjugate) A.H=transpose(A*)\n",
    "\n",
    "*matrix.I* = matrix inverse\n",
    "\n",
    "*matrix.A* = matrix as a basic ndarray\n",
    "\n",
    "*matrix.A1* = matrix as a one-dimensional ndarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "VvzCw9gkhslW"
   },
   "outputs": [],
   "source": [
    "#From Christina Maimone's (NU's Research Computing Services) '17 Numpy/Scipy \n",
    "A = np.mat([[3, 1, 4+0.1j], [1, 5-2j, 9], [2, 6+2j, 5]]) #row-wise, bracketing\n",
    "\n",
    "print(type(A))\n",
    "\n",
    "A_tr = A.T\n",
    "A_H = A.H\n",
    "A_I = A.I\n",
    "A_A = A.A\n",
    "A_A1 = A.A1\n",
    "\n",
    "print('\\n Transpose of A=')\n",
    "print(A_tr)\n",
    "print('\\n Hermitian of A=')\n",
    "print(A_H)\n",
    "print('\\n Inverse of A=')\n",
    "print(A_I)\n",
    "print('\\n Standard ndarray=')\n",
    "print(A_A)\n",
    "print('\\n One dimestional array of A=')\n",
    "print(A_A1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix Properties\n",
    "\n",
    "One can also compute matrix properties such as eigenvalues, left and right eigenvectors, determinant and trace using scipy and numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "import scipy.linalg as scipy_la\n",
    "import numpy.linalg as numpy_la\n",
    "\n",
    "A1 = np.mat([[1,2,3],[1,0,2],[0,1,3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the eigenvectors and eigenvalues, one must use the linear algebra module of scipy, scipy.linalg. This module must be imported separately.\n",
    "\n",
    "\n",
    "The right eigenvectors, $v_r$, satisfy the equation, $ A v_r = \\lambda v_r$, and the left eigenvectors, $v_l$, satisfy the equation, $v_l A  = \\lambda v_l.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[w,vl]=scipy_la.eig(A1,left=True,right=False)\n",
    "print(\"The eigenvalues are \\n\", w,\"\\n\")\n",
    "print(\"The left eigenvectors are \\n\", vl,\"\\n\")\n",
    "\n",
    "\n",
    "[w,vr]=scipy_la.eig(A1,left=False,right=True)\n",
    "print(\"The eigenvalues are \\n\", w,\"\\n\")\n",
    "print(\"The right eigenvectors are \\n\", vr,\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute trace and determinant, one can use numpy or scipy. For our purposes, we will be using numpy.\n",
    "\n",
    "Numpy contains the function, *trace*, to compute the trace. \n",
    "\n",
    "However to compute the determinant, one must port the linear algebra module of numpy separately.\n",
    "\n",
    "The same operations are true if one were to use scipy instead of numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The trace is\",np.trace(A1),\"\\n\")\n",
    "print(\"The determinant is\", numpy_la.det(A1),\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sparse Matrices\n",
    "\n",
    "When working in research, the matrices often become very large, sometimes on the order of a million by a million. The matrices, at times, contain a larger number of zeros than non-zero numbers. These are called **sparse** matrices.\n",
    "\n",
    "Since large 2-D array require a large amount of memory to store, scientists have come up with a short cut to store and manipulate these arrays when they are sparse. \n",
    "\n",
    "Sparse matrices are very important when working with big data problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On way of storing sparse matrices is the **COO** (*COO*rdinate) format. It stores the 3 arrays: row, col and data. The row array stores the row index of non-zero entries. The col array stores the column index of non-zero entries. The data array stores the value of non-zero entries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_dense = np.mat([[2,4,0,0], [0,0,0,1], [3,0,0,0]])\n",
    "A_coo_sparse = scipy.sparse.coo_matrix(A_dense)\n",
    "print(\"This is the original matrix: \\n \", A_dense,\"\\n\")\n",
    "print(\"COO matrix row indices: \",A_coo_sparse.row,\"\\n\")\n",
    "print(\"COO matrix col indices:\",A_coo_sparse.col,\"\\n\")\n",
    "print(\"This is the data stored in the COO matrix data array: \", A_coo_sparse.data,\"\\n\")\n",
    "\n",
    "print(\"To convert it back to the original matrix, simply use the todense() function \\n \",A_coo_sparse.todense())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "COO format is less commonly used since it does not accomodate fast matrix arithmetic. The more commonly used formats are CSR (Compressed Sparse Row format) and CSC (Compressed Sparse Column format).\n",
    "\n",
    "CSR like COO stores three arrays: indptr, indices and data. The *indptr* array is a tally of the non-zero entries in the rows. The second array stores the column index of non-zero entries. The data array stores the value of non-zero entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_csr_sparse = scipy.sparse.csr_matrix(A_dense)\n",
    "print(\"This is the original matrix: \\n \", A_dense,\"\\n\")\n",
    "print(\"indptr array:\",A_csr_sparse.indptr,\"\\n\")\n",
    "print(\"indices array: \",A_csr_sparse.indices,\"\\n\")\n",
    "print(\"data array: \", A_csr_sparse.data,\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1.5\n",
    "\n",
    "Similarly, CSC stores three arrays. What are these array and what do they store? See the documentation for [scipy.sparse.csc_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html). Note: CSC is very similar to CSR.\n",
    "\n",
    "Convert A_dense to CSC format using [scipy.sparse.csc_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html) and print each array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The indices array stores the row index of non-zero entries.\n",
    "#The ith entry of the indptr array stores the number \n",
    "#of non-zero entries in the columns with index less than *i*\n",
    "#The third array stores the value of non-zero entries.\n",
    "\n",
    "A_csc_sparse = scipy.sparse.csc_matrix(A_dense)\n",
    "print(\"This is the original matrix: \\n \", A_dense,\"\\n\")\n",
    "print(\"indptr array:\",A_csc_sparse.indptr,\"\\n\")\n",
    "print(\"row array: \",A_csc_sparse.indices,\"\\n\")\n",
    "print(\"data array: \", A_csc_sparse.data,\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: You can convert from one sparse format to another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "check = (scipy.sparse.coo_matrix(A_csr_sparse) != A_coo_sparse)\n",
    "print(\"It is\",check.nnz == 0,\"that A_coo_sparse and A_csr_sparse store the same data.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-81Cj5EVsrlZ"
   },
   "source": [
    "## 2-D Arrays outside the Matrix Class\n",
    "In a non class-specific sense, a matrix can be thought of as a 2-D *array* where both of the dimensions are strictly greater than 1 <br>\n",
    "This can be formed from the get-go by bracketing comma-separated sets of bracket-ed *arrays* or it can be formed by *array.reshape(...)* on an *array* of different dimensionality <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "Jveht2OtsvuJ"
   },
   "outputs": [],
   "source": [
    "a8 = np.array([ [1, 2], [3, 4] ])\n",
    "print('A 2x2 matrix like')\n",
    "print(a8)\n",
    "print('has shape')\n",
    "print(a8.shape)\n",
    "\n",
    "a6 = np.array([1, 2, 3, 4]) # just to remind\n",
    "\n",
    "a8 = a6.reshape(2,2)\n",
    "print('\\n A 2x2 matrix formed by reshaping')\n",
    "print(a8)\n",
    "print('has shape')\n",
    "print(a8.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IKJVMV5JSU83"
   },
   "source": [
    "Note that *array.reshape* only works if the dimensions multiply to the overall size of the array <br> <br>\n",
    "That is, in the example above, $\\text{dimension 1 } (\\text{ex. } 2) * \\text{dimension 2 } (\\text{ex. } 2) = \\text{overall size } (\\text{ex. }4)$ <br>\n",
    "You can query the overall size of any *array* with the method [*array.size*](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ndarray.size.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "i0y7VK5tdCc_"
   },
   "outputs": [],
   "source": [
    "print('The overall size of a8 is')\n",
    "print(a8.size)\n",
    "print('The relationship between dimension in np.shape and np.size is ' +\n",
    "      str(a8.shape[0] * a8.shape[1] == a8.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A73UYfH4RfCh"
   },
   "source": [
    "# Note: 'Length' vs 'size'\n",
    "\n",
    "The built-in Python function ```len``` **will only return the length of the first dimension of a NumPy *array*** <br>\n",
    "When working with any *array*, then, it is safer to use *np.shape()* or the method *array.shape*, and index this tuple to return the length in a specific dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "a1G6L8wmdgCz"
   },
   "outputs": [],
   "source": [
    "a9 = a8.reshape(1,4)\n",
    "print('It is ' + str(len(a9) == a9.shape[1]) + \n",
    "      ' that the len function reads the 2nd dimension')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3F0Dw6g-sv-d"
   },
   "source": [
    "In general, dimensionality can be stepped up to any arbitrary N-dimensions (as long as there are enough elements to satisfy to the overall shape) <br>\n",
    "For ex. a 'hypercube' or 'hyperplane' in *N*-dimensional 'hyperspace' will be of *N*-dimensions <br> <br>\n",
    "Perhaps the two most useful of all commands for creating *N-D arrays* is [np.repeat(...)](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.repeat.html#numpy.repeat) and [np.tile(...)](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.tile.html#numpy.tile)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L1wsJX84l2OG"
   },
   "source": [
    "## Exercise 2\n",
    "\n",
    "Create an 4-dimensional array with ```np.repeat(...)``` and ```np.tile(...)``` <br>\n",
    "The array should be filled with *1*s and be of shape ```(4, 4, 4, 4)``` <br>\n",
    "Bonus: write a function to generalize this to return an array of 1s of shape (N, N, N, etc.) for any N-dimenionsal array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "X1FUmTwuu9Td"
   },
   "outputs": [],
   "source": [
    "#Exercise 2:\n",
    "\n",
    "a_4d = np.tile(np.repeat(1,4),(4,4,4,1))\n",
    "print(a_4d.shape)\n",
    "\n",
    "def ones(N):\n",
    "    np_tile_size = np.repeat(N, N-1)\n",
    "    np_tile_size = np.append(np_tile_size, 1)\n",
    "    array = np.tile(np.repeat(1,N),(np_tile_size))\n",
    "    return array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EXJrAnAVndKs"
   },
   "source": [
    "In fact, there are already several NumPy functions for creating common *N*-dimensional arrays; [np.ones(...)](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ones.html) does what we just did in exercise 2 (and more) <br><br>\n",
    "for example, [np.zeros](https://docs.scipy.org/doc/numpy/reference/generated/numpy.zeros.html) is useful for creating zero-filled *N-D arrays* which you might use to preallocate that *array* before you (*have* to) assign to it in a loop; <br><br>\n",
    "another example, [np.identity](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.identity.html) returns the square identity matrix $I$ which is useful in matrix math (for ex. it can be used to verify the inverse of a matrix as $ I = AA^{-1}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "WOHwjYSroEoc"
   },
   "outputs": [],
   "source": [
    "A_test=A*A.I\n",
    "print('It holds ' +\n",
    "      str(np.all(np.all(np.isclose(A_test, np.identity(3)),axis=0))) + # Boolean operation np.all is described below\n",
    "     ' that np.identity is euqal to A*A.inverse')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FBevWwxt-OxM"
   },
   "source": [
    "Note: The function [np.isclose](https://docs.scipy.org/doc/numpy/reference/generated/numpy.isclose.html) is useful for comparing whether two arrays or values are equal within a certain tolerance <br>\n",
    "For example, because of differences in floating point precision\n",
    "\n",
    "```\n",
    "array1 == array2\n",
    "```\n",
    "May return ```False``` even if the arrays are filled with floating point numbers while ```np.isclose(...)``` would return ```True```<br> <br>\n",
    "The above code could have been simplified to a single functional call of [np.allclose(...)](https://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html#numpy.allclose), but bringing up ```np.all()``` gives us a reason to discuss Booleans in Numpy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "83WcP-CM-C_I"
   },
   "source": [
    "# NumPy Booleans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7dYbUvBF_APK"
   },
   "source": [
    "But there is another reason you want to be careful using ```array1 == array2```<br>\n",
    "To understand this, let's analyze\n",
    "\n",
    "```\n",
    "print(np.all(np.all(np.isclose(A_test, np.identity(3)),axis=0)))\n",
    "```\n",
    "We discussed ```np.isclose(...)``` but not ```np.all(...)``` <br> <br>\n",
    "Just like there are NumPy functions for vector math on arrays, there are also NumPy for boolean interpretation of arrays. <br> <br>\n",
    "The main take-home, though, is that **NumPy booleans are not Python booleans (and vice-versa)**--therefore **use [Numpy boolean functions](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.logic.html) instead of built-in boolean functions** (like ```or``` or ```in```) to stay out of trouble <br> <br>\n",
    "For example, [np.all(...)](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.all.html#numpy.all) in the example above tests whether all array elements along a given axis ('axis=0' in above ex.) evaluates to True. What if we did 'axis=1'?:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "zZ6jxlgbA2gY"
   },
   "outputs": [],
   "source": [
    "A_test=A*A.I\n",
    "# First evaluate np.isclose(A_test, np.identiy(3))\n",
    "print('np.isclose(...) is element-wise, returning a Matrix of the same shape')\n",
    "is_close = np.isclose(A_test,np.identity(3))\n",
    "print(is_close) # Note: this element-wise\n",
    "# Next, evaluate np.all to evaluate whether all is True along 1 dimension (axis)\n",
    "all_ax1 = np.all(is_close, axis=1)\n",
    "print('Evaluating on the other axis...')\n",
    "print(all_ax1)\n",
    "print('The above returns a Boolean array, rather than a scalar')\n",
    "# The last statement evaluates the 1D vector with np.all, doesn't need 'axis='\n",
    "all_all_ax1 = np.all(all_ax1)\n",
    "print('To get a single Boolean value, like...')\n",
    "print(all_all_ax1)\n",
    "print('we had to apply np.all() again')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A note on conditional Logic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "siRz0se0ZnC4"
   },
   "source": [
    "So you see how Booleans evaluations (conditionals) can be *vectorized* in NumPy--and, indeed, NumPy expects it! If you try conditional Python logic on an array, like\n",
    "\n",
    "```\n",
    "if True in all_ax1\n",
    "```\n",
    "It will return an  error message <br>\n",
    "\n",
    "```The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()```<br>\n",
    "\n",
    "We already used *np.all(...)*, and [np.any(...)](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.all.html#numpy.all) is the vector equivalent of\n",
    "```\n",
    "True in all_ax1\n",
    "```\n",
    "That is, if it occurs once (across that axis), it will evaluate as True <br><br>\n",
    "Perhaps the most useful task you can acheive with a Boolean *N-D array* is as a 'mask' (acting like indices) to filter an *array* of the same dimensionality.<br>\n",
    "For this, [np.where(...)](https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html) is invaluable; it is a vectorized (element-wise) version of the satement:  ```x``` if ```condition is True``` else ```y``` with the format\n",
    "\n",
    "```\n",
    "np.where(condition, x, y)\n",
    "```\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "zLJd46UEizrt"
   },
   "outputs": [],
   "source": [
    "arr=np.random.randn(4,4)\n",
    "print(arr)\n",
    "mask_int = np.where(arr > 0, 1, 0)\n",
    "# Unfortunately, this will not work as mask as it is not a Boolean type in NumPy\n",
    "# arr[mask_int] # indexes arr quite strangely\n",
    "mask_bool = mask_int.astype(np.bool)\n",
    "print('The positive elements in the array are')\n",
    "print(arr[mask_bool])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SXT8q6k7d6Tr"
   },
   "source": [
    "## Exercise 3\n",
    "Given the two following 2-D *arrays*, use *np.not_equal(...)* to evaluate where they are equal and return a 0 if they are, or else return a 1 (basically the opposite of the Boolean array *np.isclose()* normally returns). <br>\n",
    "Use the output as a 'mask' to select the un-matching elements from either array.\n",
    "\n",
    "Bonus: use *np.isclose* as a condition (the first argument) in *np.where(...)* to do the same.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "edvkibPLeq8y"
   },
   "outputs": [],
   "source": [
    "#Exercise 3\n",
    "a10 = np.random.randn(4, 4)\n",
    "np.fill_diagonal(a10, 1.)\n",
    "a11 = np.random.randn(4, 4)\n",
    "np.fill_diagonal(a11, 1.)\n",
    "print('a10 is')\n",
    "print(a10)\n",
    "print('a11 is')\n",
    "print(a11)\n",
    "\n",
    "a_mask = np.not_equal(a10,a11)\n",
    "\n",
    "print(a10[a_mask])\n",
    "\n",
    "print(a11[a_mask])\n",
    "\n",
    "#Bonus\n",
    "\n",
    "a_mask_bonus = np.where(np.isclose(a10,a11), False, True)\n",
    "\n",
    "print(a_mask_bonus)\n",
    "# print(a10[a_mask_bonus])\n",
    "# print(a11[a_mask_bonus])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HowphD4HpFWF"
   },
   "source": [
    "# 'Random' Sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-jHo163qibAh"
   },
   "source": [
    "The NumPy routine for handling/generating random events, numbers, etc. is [numpy.random](https://docs.scipy.org/doc/numpy-1.14.0/reference/routines.random.html)\n",
    "\n",
    "Because there are sub-functions for the random routine, e.g. ([random.random_sample(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.random_sample.html#numpy.random.random_sample)) it is often useful to import *np.random* as 'random' (or as *) to make calls to this routine shorter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "sJ5hllVgk01K"
   },
   "outputs": [],
   "source": [
    "import numpy.random as random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomness is extremely useful in numerical evaluations, simulations, and algorithms in general."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OrLfk9pCk3cO"
   },
   "source": [
    "## Sampling Uniform Probability Distributions\n",
    "\n",
    "For example, if we want to randomly sample (with uniform probability) from any given array (discrete values), we would use [numpy.random.choice(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.choice.html#numpy.random.choice)\n",
    "\n",
    "Note: this needs a 1-D array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "yk052lelk4Zs"
   },
   "outputs": [],
   "source": [
    "random_sample = np.array(random.choice(a6.reshape(4,),size=20))\n",
    "print(random_sample)\n",
    "  \n",
    "# or return as an array in a list comprehension\n",
    "print(np.array([random.choice(a6.reshape(4,)) for i in np.arange(20)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VsjiGbMOmxN9"
   },
   "source": [
    "What you'll notice by comparing with another workshop attendee or by running the above cell multiple times is that this sequence of samples appears truly random for each instance. <br>\n",
    "However, at times (especially when comparing or testing code), it is more useful when randomness can be repeated. <br>\n",
    "Because *numpy.random* works from a [pseudo-random number generator](https://en.wikipedia.org/wiki/Pseudorandom_number_generator) (PRNG, specifically the [Mersenne Twister](https://en.wikipedia.org/wiki/Mersenne_Twister) algorithm), it can return the same 'random sample' repeatedly by fixing the 'seed' for the PRNG. This is done with [numpy.random.seed(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.seed.html#numpy.random.seed) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "gZtVyItZpJri"
   },
   "outputs": [],
   "source": [
    "print('A random sample can be recreated with fidelity by fixing the seed')\n",
    "for i in np.arange(5):\n",
    "  random.seed(seed = 51)\n",
    "  random_sample = np.array(random.choice(a6.reshape(4,),size=20))\n",
    "  print(random_sample)\n",
    "#   random_sample = np.array(random.choice(a6.reshape(4,),size=20))\n",
    "#   print(random_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RmWRYIABrIuh"
   },
   "source": [
    "Another important feature of random sampling from a pre-existing array is the option to *replace* while sampling. <br> <br>\n",
    "By default, this option is true, so technically the *np.random.choice(...)* sampling performs [bootstrapping](https://en.wikipedia.org/wiki/Bootstrapping_(statistics), which is a powerful method in statistics. <br> <br>\n",
    "When *replace=False*, the length of the sample size must be less than the length of the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "qmCwCYONtGyE"
   },
   "outputs": [],
   "source": [
    "random_sample = np.array(random.choice(a6.reshape(4,),size=4, replace=False))\n",
    "print(random_sample)\n",
    "\n",
    "# Uncommenting the below code will give an error due to the size of the sample\n",
    "# random_sample = np.array(random.choice(a6.reshape(4,),size=20, replace=False))\n",
    "# print(random_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rHwjs5zjtyUP"
   },
   "source": [
    "Another example: <br>\n",
    "Say you don't want a set of discrete values you want to sample from, but rather, you only know a *continuous* range of numbers or distribution limits you'd like to sample from: *numpy.random* is still useful for randomly sampling (with uniform probability). <br> <br>\n",
    "The most generic method for this is the aptly-named [np.random.random_sample(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.random_sample.html#numpy.random.random_sample) <br>\n",
    "This ''returns random floats in the half-open interval [0, 1)'' (1 not inclusive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "mk80WmZLu4Xd"
   },
   "outputs": [],
   "source": [
    "print(random.random_sample(size=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YaH5LB95vISi"
   },
   "source": [
    "Note: this function is redundant with* np.random.random(...)*, *np.random.ranf(...)*, and *np.random.sample(...)* <br>\n",
    "<br>\n",
    "Okay...but let's say you want your continuous (uniform) sampling from '2' to '10' aka [2, 10)<br>\n",
    "Just multiply the random.random_sample by '(b -a)' and add the offset 'a'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "VOZuX2fPv6HX"
   },
   "outputs": [],
   "source": [
    "b = 10.\n",
    "a = 2.\n",
    "sample_2_10 = (10.-2.)*random.random_sample(size=10) + 2.\n",
    "print(sample_2_10)\n",
    "# and if you wanted these to be rounded to integers...\n",
    "print(np.array([int(sample_2_10[i]) for i in np.arange(np.size(sample_2_10))]))\n",
    "# or you could just use np.floor to round down (or np.ceil to round up)\n",
    "print(np.floor(sample_2_10).astype(int))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also [np.random.randint(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.randint.html#numpy.random.randint) and [np.random.random_integers(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.random_integers.html#numpy.random.random_integers) for similar behavior to that last line of code. <br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UM3vOtjSzKPK"
   },
   "source": [
    "# Non-uniform, continuous probability distributions\n",
    "Often a uniform probability distribution is not useful. <br>\n",
    "Thankfully, *numpy.random* has 35 common distributions built-in from which you can sample. <br>\n",
    "For example, one of the most common distributions, the 'Normal' (Gaussian) distribution, has its own random sampling function ([np.random.randn(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.randn.html#numpy.random.randn)) that is redundant with [np.random.normal(...)](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.normal.html#numpy.random.normal)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "m4dB1bSj2HVE"
   },
   "outputs": [],
   "source": [
    "# for plotting\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "#for sampling from the normal (gaussian) distribution\n",
    "random.seed(6) # note: same as random.seed(seed = 6)\n",
    "mu, sigma = 0.1, 1.\n",
    "n_dist = random.normal(mu, sigma, 500)\n",
    "\n",
    "#for sampling from cauchy (lorentzian) distribution\n",
    "random.seed(6)\n",
    "c_dist = random.standard_cauchy(120)\n",
    "\n",
    "plt.hist(n_dist)\n",
    "plt.hist(c_dist,45)\n",
    "plt.show() #not strictly necessary with %pylab inline, but removes some print-out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ajupZHFDSBi_"
   },
   "source": [
    "Great...but what if you wan't to use a non-standard cauchy distribution? (one not centered on zero and with gamma variance greater than 1) <br>\n",
    "<br>\n",
    "First, note that in these calls of *random.normal*, we have been using shorthand\n",
    "\n",
    "```\n",
    "random.normal(mu, sigma, 500)\n",
    "```\n",
    "which is really\n",
    "\n",
    "```\n",
    "random.normal(loc = mu, scale = sigma, size = 500)\n",
    "```\n",
    "the ```loc``` and ```scale``` keywords are, in general, related to the statistical variables that determine the center (ex. 'mean' for Normal) and variance (ex. 'beta' for Cauchy). <br> <br>\n",
    "But NumPy's built-in [cauchy](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.standard_cauchy.html) doesn't have ```loc``` or ```mu``` keywords! <br>\n",
    "So, back to the question: to sample from a non-standard cauchy distribution, we'll need *SciPy*, namely [scipy.stats.cauchy(...)](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.cauchy.html) which has the ```rvs()``` method for random-variable sampling.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "0VKwhMsfaZUm"
   },
   "outputs": [],
   "source": [
    "import scipy.stats\n",
    "from scipy.stats import cauchy\n",
    "\n",
    "random.seed(6)\n",
    "c_dist = random.standard_cauchy(120)\n",
    "\n",
    "random.seed(6)\n",
    "c_dist2 = scipy.stats.cauchy.rvs(loc=100, scale=2.5, size=120)\n",
    "plt.hist(c_dist,5)\n",
    "plt.hist(c_dist2,45)\n",
    "plt.show()\n",
    "print(type(scipy.stats.cauchy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3AttA1w_bnEW"
   },
   "source": [
    "There are a multitude of other statistical applications that can be evaluated with ```scipy.stats```; for example, on any specific continuous distribution (in this case of the *rv_continuous* class), the mean, variance, skewness, and kurtosis can be readily generated.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "S5AQ0_AHcalP"
   },
   "outputs": [],
   "source": [
    "mean, var, skew, kurt = scipy.stats.norm.stats(moments='mvsk')\n",
    "print('The mean of the normal distribution is ' + str(mean))\n",
    "print('The variance of the normal distribution is ' + str(var))\n",
    "print('The skew of the normal distribution is ' + str(skew))\n",
    "print('The kurtosis of the normal distribution is ' + str(kurt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-7UBr0DMgW-S"
   },
   "source": [
    "The additional method ```stats``` is specific to the ```rv_continuous``` class. If there is a distribution you can't find in ```numpy.random``` or ```scipy.stats```, you can actually generate a custom one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "ii59S_NWhA3L"
   },
   "outputs": [],
   "source": [
    "from scipy.stats import rv_continuous #import rv_continuous class\n",
    "\n",
    "class gaussian_gen(rv_continuous):\n",
    "  \"Gaussian distribution\"                           #comment describing function\n",
    "  def _pdf(self, x):                       #rewrites the pdf method to be called\n",
    "    return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi) #The custom function\n",
    "\n",
    "\n",
    "custom_gauss = gaussian_gen() #instance created for custom rv_continuous class\n",
    "n_dist2 = custom_gauss.rvs(size = 500) #instance method .rvs\n",
    "plt.hist(n_dist)\n",
    "plt.hist(n_dist2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "QDNU34BZkBn9"
   },
   "source": [
    "# Creating + Fitting 'Noisy' Data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Zi5lqOMVLSzA"
   },
   "source": [
    "NumPy also has many built-in math functions. Trigonometric expressions like [np.sin](https://docs.scipy.org/doc/numpy/reference/generated/numpy.sin.html) can be called on numpy *array* for instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "ii6jSfQyLy7S"
   },
   "outputs": [],
   "source": [
    "t = np.linspace(0, 2*np.pi, 101) #array of 101 pts evenly spaced from 0 to 2pi\n",
    "y = np.sin(t) # apply trig sine function to array\n",
    "plt.plot(t, y) # plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qN9BsarcMpL1"
   },
   "source": [
    "Now, let's say we want to add noise to this data to simulate real data that one might collect from an instrument.\n",
    "\n",
    "We can do this with the np.random routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "e25r3tuqNU6b"
   },
   "outputs": [],
   "source": [
    "random.seed(51)\n",
    "noisy_y = random.normal(loc = np.sin(t), scale = 0.3) # Gaussian noise add to y\n",
    "plt.plot(t, noisy_y, 'o') # plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "J9R-SMogO1yp"
   },
   "source": [
    "Now the data generated from a continuous function acting on discrete points looks realistically noisy. <br> <br>\n",
    "Imagine this were from an experiment, you might want to 'fit' it with some form of regression (going from discrete points back to a continuous function).<br>\n",
    "NumPy has some black-box functions for regression, the most common of which are based on least-squares fitting (minimizing the square of the residual between the data and fit). <br>\n",
    "A particularly easy one to use (if we naively assume the data is representative of a polynomial) is [np.polyfit(...)](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html)--used to perform a least-squares polynomial fit--and [np.poly1d](https://docs.scipy.org/doc/numpy/reference/generated/numpy.poly1d.html)--used to calculate a new *y* from the fit.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "hQrC07DFR0LW"
   },
   "outputs": [],
   "source": [
    "fit_obj = np.polyfit(t, noisy_y, 3) #perform least-squares fit\n",
    "fit_fn = np.poly1d(fit_obj) #generate function parameters\n",
    "fit_y = fit_fn(t) # apply function to x to generate fit y\n",
    "plt.plot(t, noisy_y, 'o') #plot as 'o' dots\n",
    "plt.plot(t, fit_y, 'r') #plot as red 'r' line\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks pretty good! <br>\n",
    "But, of course, knowing that this was generated from a sine function, outside of these bounds the fit would be rather poor. <br> <br>\n",
    "\n",
    "Note: If you *do* know your data is periodic, a better *basis* for your fit might be Fourier basis <br>\n",
    "Or, on that same note: if the data is conditioned properly (windowed to zero, constant time interval, etc.), the fast Fourier Transform (of discrete data) is easily implementable in NumPy and could give the frequency of sine or cosine functions. <br> The relevant functions for this are in the Numpy [fft](https://docs.scipy.org/doc/numpy/reference/routines.fft.html) module. (There is also a nearly-identical SciPy module)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear (least-squares) fitting, in general\n",
    "Focusing on the math behind np.poly1d, you realize it is actually solving a linear least-squares problem for the coefficients of the polynomial (the data $y$ is linear in the coefficients even though it is non-linear in $x$). <br> <br>\n",
    "This means we could solve the above fitting more explicitly with a general least-squares problem by considering the features $x$ of the model that are linearly related to $y$ <br>\n",
    "The general module for linear algebra in NumPy is [np.linalg](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.linalg as linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Returning to the problem at hand, in generic notation, linear fitting means you have some data $y$, model parameters $x$, and $A$ a vector holding the data as it relates to each model parameter. In other words: <br>\n",
    "$$y = Ax$$\n",
    "and you want to solve for\n",
    "$$x = A^{-1}y$$\n",
    "<br>\n",
    "If $A$ is invertible, then, it's easy enough to solve this problem with [np.dot](https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html) and [np.linalg.inv](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.inv.html) <br> <br>\n",
    "But $A$ is not invertible (in this case, it is not square), so we minimize the squared difference between $y$ and $Ax$ ('ordinary' least-squares fitting, a.k.a. minimize the L2 norm). This has the general solution:\n",
    "$$ x = (A^TA)^{-1}A^Ty$$\n",
    "First, lets define A more concretely in terms of our sine-wave/3rd order polynomial example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  The model we are trying to fit is y = mt^3 + mt^2 + mt + b\n",
    "#the relationship b/w y and m in each is [t^3,   t^2,   t,  1]\n",
    "A = np.vstack([t**3, t**2, t, np.ones(np.shape(t)[0])]).T\n",
    "print('The shape of A is ' + str(np.shape(A)))\n",
    "print('because there are 101 points (t) and 4 linear parameters')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n",
    "Solve for x with the general least squares solution\n",
    "$$ x = (A^TA)^{-1}A^Ty$$\n",
    "using [np.dot](https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html) and [np.linalg.inv](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.inv.html) <br><br>\n",
    "Assign the coefficients to variables ```m3```, ```m2```, ```m```, and ```c``` and use the next cell to plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 4\n",
    "m3, m2, m, b = np.dot(linalg.inv(np.dot(A.T, A)), np.dot(A.T, y))\n",
    "print(m3)\n",
    "print(m2)\n",
    "print(m)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to make our fit look continuous, we need finer sampling\n",
    "tf = np.linspace(0, 2*np.pi, 501) #finer sampling points\n",
    "fit_y_genlstsq = np.polyval([m3, m2, m, b], tf) #evaluate fit\n",
    "plt.plot(t, noisy_y, 'o') #plot as 'o' dots\n",
    "plt.plot(tf, fit_y_genlstsq, 'r') #plot as red 'r' line\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The black-box way to solve this is just with [np.linalg.lstsq(...)](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.lstsq.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Within precision tolerance, it is ' + \n",
    "      str(np.allclose(np.linalg.lstsq(A, y, rcond=None)[0],\n",
    "                      #index 0th element for coeffecients\n",
    "                      [m3, m2, m, b])) + \n",
    "      ' that numpy.linalg.lstsq returns the same coefficients.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For non-linear fitting (parameters are not linearly related to the data), you would use [scipy.optimize.curve_fit(...)](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy and beyond!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most scipy modules can be accessed by importing the 'whole' scipy library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import scipy #we did this before with 'as spy'\n",
    "#help(scipy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, in the help documents, you'll see (near the top) there are 'Subpackages' and\n",
    "```\n",
    "Using any of these subpackages requires an explicit import.\n",
    "```\n",
    "We utilized this for ```import scipy.linalg as scipy_la``` before <br><br>\n",
    "There is a **lot** of functionality built into Scipy. To briefly describe some of the most important components, defer to a table from (USCB)[https://engineering.ucsb.edu/~shell/che210d/numpy.pdf] <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Module | Code for.. \n",
    " --- | --- \n",
    "scipy.constants| Many mathematical and physical constants. \n",
    "scipy.integrate| Special functions for mathematical physics, such as iry, elliptic, bessel, gamma, beta, hypergeometric, parabolic cylinder, mathieu, spheroidal wave, struve, and kelvin functions. \n",
    "scipy.integrate | Functions for performing numerical integration using trapezoidal, Simpson's, Romberg, and other methods. Also provides methods for integration of ordinary differential equations.\n",
    "scipy.optimize | Standard minimization / maximization routines that operate on generic user-defined objective functions. Algorithms include: Nelder-Mead Simplex, Powell's, conjugate gradient, BFGS, least-squares, constrained optimizers, simulated annealing, brute force, Brent's method, Newton's method, bisection method, Broyden, Anderson, and line search.\n",
    "*scipy.linalg* | Much broader base of linear algebra routines than NumPy. Offers more control for using special, faster routines for specific cases (e.g., tridiagonal matrices). Methods include: inverse, determinant, solving a linear system of equations, computing norms and pseudo/generalized inverses, eigenvalue/eigenvector decomposition, singular value decomposition, LU decomposition, Cholesky decomposition, QR decomposition, Schur decomposition, and various other mathematical operations on matrices.\n",
    "*scipy.sparse* | Routines for working with large, sparse matrices.\n",
    "scipy.interpolate |  Routines and classes for interpolation objects that can be used with discrete numeric data. Linear and spline interpolation available for one- and two-dimensional data sets\n",
    "scipy.fftpack | Fast Fourier transform routines and processing\n",
    "scipy.signal | Signal processing routines, such as convolution, correlation, finite fourier transforms, B-spline smoothing, filtering, etc.\n",
    "*scipy.stats* | Huge library of various statistical distributions and statistical functions for operating on sets of data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already used the modules in this workshop italicized in the table above. <br>\n",
    "Best of luck to you making use of the others!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "default_view": {},
   "name": "Numpy-scipy.ipynb",
   "provenance": [],
   "version": "0.3.2",
   "views": {}
  },
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}