<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A Workshop Introduction to NumPy\n",
    "\n",
    "The Python language is an excellent tool for general-purpose programming, with a highly readable syntax, rich and powerful data types (strings, lists, sets, dictionaries, arbitrary length integers, etc) and a very comprehensive standard library. It was not, however, designed specifically for mathematical and scientific computing. Neither the language nor its standard library have facilities for the efficient representation of multidimensional datasets, tools for linear algebra and general matrix manipulations (essential building blocks of virtually all scientific computing).\n",
    "\n",
    "For example, Python lists are very flexible containers that can be nested arbitrarily deep and can hold any Python object in them, but they are poorly suited to represent common mathematical constructs like vectors and matrices. It is for this reason that NumPy exists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Workshop Aims\n",
    "\n",
    "The aim of this workshop is to enable you use NumPy to:\n",
    "\n",
    "  - manipulate numerical arrays, and\n",
    "  - perform efficient array calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Table of Contents\n",
    "\n",
    "* [Getting Started with NumPy](#getting-started)\n",
    "* [Motivation: what are arrays good for?](#motivation)\n",
    "* [The Array Object](#array-object)\n",
    "* [Application: Arithmetic and Broadcasting](#arithmetic_and_broadcasting)\n",
    "* [Application: Statistics and Masked Arrays](#statistics)\n",
    "* [Application: Efficiency](#efficiency)\n",
    "* [Further Reading](#further_reading)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting Started with NumPy <a class=\"anchor\" id=\"getting-started\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will be able to create NumPy arrays with the aid of the NumPy documentation.\n",
    "\n",
    "NumPy is the fundamental package for scientific computing with Python. Its primary purpose is to provide a powerful N-dimensional array object; the focus for this workshop.\n",
    "\n",
    "To begin with let's import NumPy, check where it is being imported from and check the version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NumPy is generally imported as 'np'.\n",
    "import numpy as np\n",
    "print(np)\n",
    "print(np.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Documentation\n",
    "Here is a link to the NumPy documentation for v1.11:\n",
    "https://docs.scipy.org/doc/numpy-1.11.0/reference/\n",
    "\n",
    "There are many online forums with tips and suggestions for solving problems with NumPy, such as http://stackoverflow.com/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explore: creating arrays\n",
    "\n",
    "NumPy provides many different ways to create arrays. These are listed in the documentation at: https://docs.scipy.org/doc/numpy-1.11.0/user/basics.creation.html#arrays-creation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "##### Part 1\n",
    "\n",
    "Try out some of these ways of creating NumPy arrays. See if you can produce:\n",
    "\n",
    "* a NumPy array from a list of numbers,\n",
    "* a 3-dimensional NumPy array filled with a constant value -- either `0` or `1`,\n",
    "* a NumPy array filled with a constant value -- **not** `0` or `1`. (**Hint:** this can be achieved using the last array you created, or you could use `np.empty` and find a way of filling the array with a constant value),\n",
    "* a NumPy array of 8 elements with a range of values starting from `0` and a spacing of `3` between each element, and\n",
    "* a NumPy array of 10 elements that are logarithmically spaced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part 2\n",
    "\n",
    "How could you change the shape of the 8-element array you created previously to have shape `(2, 2, 2)`? **Hint:** this can be done without creating a new array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Motivation: what are arrays good for? <a class=\"anchor\" id=\"motivation\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will understand key benefits of using NumPy arrays for numerical processing in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extensive Features\n",
    "\n",
    "\n",
    "NumPy provides routines for fast operations on arrays, including mathematical, logical, shape manipulation, sorting, selecting, I/O, discrete Fourier transforms, basic linear algebra, basic statistical operations, random simulation and much more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fast Calculations\n",
    "\n",
    "It is a lot faster than Python alone for numerical computing tasks.\n",
    "\n",
    "Element-by-element operations are the “default mode” when an ndarray is involved, but the element-by-element operation is speedily executed by pre-compiled C code. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clear Syntax\n",
    "\n",
    "In NumPy,\n",
    "\n",
    "```python\n",
    "c = a * b\n",
    "```\n",
    "\n",
    "calculates the element-wise product of `a` and `b`, at near-C speeds, but with the code simplicity we expect from the Python language.\n",
    "\n",
    "This demonstrates a core feature of NumPy, called **vectorization**. This removes the need for loops iterating through elements of arrays, which can make code easier to read, as well as performing fast calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interfacing to other Libraries\n",
    "\n",
    "Many scientific Python libraries use NumPy as their core array representation. From plotting libraries such as matplotlib, to parallel processing libraries such as Dask, to data interoperability libraries such as Iris, NumPy arrays are at the core of how these libraries operate and communicate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## The Array Object <a class=\"anchor\" id=\"array-object\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will be able to manipulate NumPy array objects through indexing and explain key features and properties of NumPy arrays.\n",
    "\n",
    "The multidimensional array object is at the core of all of NumPy's functionality.  Let's explore this object some more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array properties\n",
    "\n",
    "Let's create a NumPy array and take a look at some of its properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.ones((3, 2, 4))\n",
    "\n",
    "print(\"Array shape:\", arr.shape)\n",
    "print(\"Array element dtype:\", arr.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "##### Part 1\n",
    "\n",
    "Use the code cells below to explore features of this array.\n",
    "\n",
    "* What is the `type` of the array created above? (Hint: you can find the type of an object in Python using `type(<object>)`, where `<object>` is replaced with the name of the variable containing the object.)\n",
    "* Look this type up in the NumPy documentation (see https://docs.scipy.org/doc/numpy-1.11.0/reference/generated/numpy.ndarray.html#numpy.ndarray)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part 2\n",
    "\n",
    "Consider the following NumPy array properties:\n",
    "\n",
    "* `ndim`,\n",
    "* `nbytes`,\n",
    "* `size`,\n",
    "* `T`.\n",
    "\n",
    "For each of these properties:\n",
    "\n",
    "* Use information from the array documentation to find out more about each array property.\n",
    "* Apply each property to the array `arr` defined above. Can you explain the results you get in each case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexing\n",
    "\n",
    "You can index NumPy arrays in the same way as other Python objects, by using square brackets `[]`. This means we can index to retrieve a single element, multiple consecutive elements, or a more complex sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.array([1, 2, 3, 4, 5, 6])\n",
    "print(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"arr[2] = {}\".format(arr[2]))\n",
    "print(\"arr[2:5] = {}\".format(arr[2:5]))\n",
    "print(\"arr[::2] = {}\".format(arr[::2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also index multidimensional arrays using an enhanced indexing syntax, which allows for multi-element indexing. (Sometimes this is referred to as \"extended slicing\".) Remember that Python uses zero-based indexing!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lst_2d = [[1, 2, 3], [4, 5, 6]]\n",
    "arr_2d = arr.reshape(2, 3)\n",
    "\n",
    "print(\"2D list:\")\n",
    "print(lst_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"2D array:\")\n",
    "print(arr_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Single array element:\")\n",
    "print(arr_2d[1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Single row:\")\n",
    "print(arr_2d[1])\n",
    "print(\"First two columns:\")\n",
    "print(arr_2d[:, :2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you only provide one index to slice a multidimensional array, then the slice will be expanded to \"`:`\" for all of the remaining dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Second row: {} is equivalent to {}'.format(arr_2d[1], arr_2d[1, :]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is known as **ellipsis**. Ellipsis can be specified explicitly using \"`...`\", which automatically expands to \"`:`\" for each dimension unspecified in the slice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.empty((4, 6, 3))\n",
    "print('Original shape: ', arr1.shape)\n",
    "print(arr1[...].shape)\n",
    "print(arr1[..., 0:2].shape)\n",
    "print(arr1[2:4, ..., ::2].shape)\n",
    "print(arr1[2:4, :, ..., ::-1].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Boolean Indexing\n",
    "\n",
    "NumPy provides syntax to index conditionally, based on the data in the array.\n",
    "\n",
    "You can pass in an array of True and False values (a boolean array), or, more commonly, a condition that returns a boolean array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(arr_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bools = arr_2d % 2 == 0\n",
    "print(bools)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(arr_2d[bools])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "##### Part 1\n",
    "\n",
    "Why do these indexing examples give the stated results?\n",
    "\n",
    " * result of `arr_2d[1, 0]` is `4`\n",
    " * result of `arr_2d[0]` is `[1, 2, 3]`\n",
    " * result of `arr_2d[1, 1:]` is `[5, 6]`\n",
    " * result of `arr_2d[0:, ::2]` is `[[1, 3], [4, 6]]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part 2\n",
    "\n",
    "How would you index `arr_2d` to retrieve:\n",
    "\n",
    " * the third value: resulting in `3`\n",
    " * the second row: resulting in `[4 5 6]` \n",
    " * the first column: resulting in `[1 4]`\n",
    " * the first column, retaining the outside dimension: resulting in `[[1] [4]]`\n",
    " * only values greater than or equal to 3: resulting in `[3 4 5 6]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arrays are not lists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: why do the following examples produce different results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lst_2d[0:2][1])\n",
    "print(arr_2d[0:2, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result we just received points to an important piece of learning, which is that in most cases NumPy arrays behave very differently to Python lists. Let's explore the differences (and some similarities) between the two."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### dtype\n",
    "\n",
    "A NumPy array has a fixed data type, called `dtype`.  This is the type of all the elements of the array.  \n",
    "\n",
    "This is in contrast to Python lists, which can hold elements of different types.\n",
    "\n",
    "#### Exercise\n",
    "\n",
    "* What happens in Python when you add an integer to a float?\n",
    "* What happens when you put an integer into a NumPy float array?\n",
    "* What happens when you do numerical calculations between arrays of different types?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating 2D coordinate arrays\n",
    "\n",
    "A common requirement of NumPy arrays is to generate arrays that represent the coordinates of our data.\n",
    "\n",
    "When orthogonal 1d coordinate arrays already exist, NumPy's meshgrid function is very useful:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 9, 3)\n",
    "y = np.linspace(4, 8, 3)\n",
    "x2d, y2d = np.meshgrid(x, y)\n",
    "print(x2d)\n",
    "print(y2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The Array Object: Summary of key  points\n",
    " * properties : `shape`, `dtype`.\n",
    " * arrays are homogeneous; all elements have the same type: `dtype`.\n",
    " * indexing arrays to produce further arrays: views on the original arrays.\n",
    " * multi-dimensional indexing (slicing) and boolean indexing.\n",
    " * combinations to form 2D arrays: `meshgrid`\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Arithmetic and Broadcasting<a class=\"anchor\" id=\"app-calc\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section you will be able to explain how broadcasting allows mathematical operations between NumPy arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elementwise Arithmetic <a class=\"anchor\" id=\"arithmetic_and_broadcasting\"></a>\n",
    "\n",
    "You can use NumPy to perform arithmetic operations between two arrays in an element-by-element fashion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.arange(4)\n",
    "arr2 = np.arange(4)\n",
    "\n",
    "print('{} + {} = {}'.format(arr1, arr2, arr1 + arr2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise \n",
    "\n",
    "Explore arithmetic operations between arrays:\n",
    "\n",
    "* Create a number of arrays of different shapes and dtypes. Make sure some of them include values of `0`.\n",
    "* Make use of all the basic Python mathematical operators (i.e. `+`, `-`, `*`, `/`, `//`, `%`).\n",
    "\n",
    "What operations work? What operations do not work? Can you explain why operations do or do not work?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It makes intrinsic sense that you should be able to add a constant to all values in an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.arange(4)\n",
    "const = 5\n",
    "\n",
    "print(\"Original array: {}\".format(arr))\n",
    "print(\"\")\n",
    "print(\"Array + const: {}\".format(arr + const))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting\n",
    "\n",
    "There are times when you need to perform calculations between NumPy arrays of different sizes.\n",
    "\n",
    "For example, suppose you have maximum temperatures from each of three recording stations, recorded on two separate days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "daily_records = np.array([[12, 14, 11], [11, 12, 15]])\n",
    "\n",
    "print('raw data:')\n",
    "print(daily_records)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each station is known to overstate the maximum recorded temperature by a different known constant value. You wish to subtract the appropriate offset from each station's values.\n",
    "\n",
    "You can do that like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "offset = np.array([2, 1, 4])\n",
    "\n",
    "corrected_records = daily_records - offset\n",
    "\n",
    "print('corrected values:')\n",
    "print(corrected_records)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy allows you to do this easily using a powerful piece of functionality called **broadcasting**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Broadcasting is a way of treating the arrays ***as if*** they had the same dimensions, and thus have elements all corresponding.  It is then easy to perform the calculation, element-wise.  \n",
    "It does this by matching dimensions in one array to the other where possible, and using repeated values where there is no corresponding dimension in the other array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rules of Broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Broadcasting applies these three rules:\n",
    "\n",
    "1.    If the two arrays differ in their number of dimensions, the shape of the array with fewer dimensions is padded with ones on its leading (left) side.\n",
    "\n",
    "1.    If the shape of the two arrays does not match in any dimension, either array with shape equal to 1 in a given dimension is stretched to match the other shape.\n",
    "\n",
    "1.    If in any dimension the sizes disagree and neither has shape equal to 1, an error is raised.\n",
    "\n",
    "Note that all of this happens without ever actually creating the expanded arrays in memory! This broadcasting behavior is in practice enormously powerful, especially given that when NumPy broadcasts to create new dimensions or to 'stretch' existing ones, it doesn't actually duplicate the data. In the example above the operation is carried out as if the scalar 1.5 was a 1D array with 1.5 in all of its entries, but no actual array is ever created. This can save lots of memory in cases when the arrays in question are large. As such this can have significant performance implications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Illustration of broadcasting](../images/fig_broadcast_visual_1.png)\n",
    "\n",
    "([image source](http://www.astroml.org/book_figures/appendix/fig_broadcast_visual.html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "For the following cases:\n",
    "\n",
    "* What will be the result of adding `arr1` to `arr2`? \n",
    "* What will be the shape of the resulting array?\n",
    "* What rules of broadcasting are being used?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((2, 3))\n",
    "arr2 = np.ones((2, 1))\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((2, 3))\n",
    "arr2 = np.ones(3)\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((1, 3))\n",
    "arr2 = np.ones((2, 1))\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((1, 3))\n",
    "arr2 = np.ones((1, 2))\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Reshaping arrays to aid broadcasting\n",
    "\n",
    "NumPy allows you to change the shape of an array, so long as the total number of elements in the array does not change. For example, we could reshape a flat array with 12 elements to a 2D array with shape `(2, 6)`, or `(3, 2, 2)`, or even `(3, 4, 1)`. We could not, however, reshape it to have shape `(2, 5)`, because the total number of elements would not be kept constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose you want to apply a correction for each day.\n",
    "\n",
    "For example, you might try :\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "days_adjust = np.array([1.5, 3.7])\n",
    "adjusted = daily_records - days_adjust"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but that results in a `ValueError`. Clearly, this doesn't work !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "With reference to the above rules of broadcasting:\n",
    " 1. describe why the above attempt to subtract `days_adjust` fails\n",
    " 1. work out how you can modify the 'days_adjust' array to get the desired result.  \n",
    "    Hint: imagine how the 'days_adjust' values *should* look, when expanded to match the dimensions of 'daily_records'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise \n",
    "\n",
    "Sometimes an operation will produce a result, but not the one you wanted.\n",
    "\n",
    "For example, suppose we have :\n",
    "```python\n",
    "A = np.array([[1, 2, 3],\n",
    "              [4, 5, 6],\n",
    "              [7, 8, 9]])\n",
    "```\n",
    "\n",
    "and we wish to add 0, 100 and 400 ***to the rows***.\n",
    "That is,\n",
    "\n",
    "```python\n",
    "B = np.array([0, 100, 400])\n",
    "# and...\n",
    "desired_result = np.array([[  1,   2,   3],\n",
    "                           [104, 105, 406],\n",
    "                           [407, 408, 409]])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Questions:\n",
    "\n",
    " 1. What will be the result of adding A to B directly?\n",
    " 1. How can you perform the correct calculation on `A` and `B` to get `desired_result`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Arithmetic and Broadcasting: Summary of key  points\n",
    "\n",
    " * arithmetic operations are performed in an element-by-element fashion,\n",
    " * operations can be performed between arrays of different shapes,\n",
    " * the arrays' dimensions are aligned according to fixed rules; where one input lacks a given dimension, values are repeated,\n",
    " * reshaping can be used to get arrays to combine as required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Statistics and Masked Arrays <a class=\"anchor\" id=\"statistics\"></a>\n",
    "\n",
    "**Learning outcome:** By the end of this section, you will be able to apply statistical operations and masked arrays to real-world problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics\n",
    "\n",
    "Numpy arrays support many common statistical calculations. For a list of common operations, see: https://docs.scipy.org/doc/numpy/reference/routines.statistics.html.\n",
    "\n",
    "The simplest operations consist of calculating a single statistical value from an array of numbers -- such as a mean value, a variance or a minimum.\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(12).reshape((3, 4))\n",
    "mean = np.mean(a)\n",
    "print(a)\n",
    "print(mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Used without any further arguments, statistical functions simply reduce the whole array to a single value.  In practice, however, we very often want to calculate statistics over only *some* of the dimensions.\n",
    "\n",
    "The most common requirement is to calculate a statistic along a single array dimension, while leaving all the other dimensions intact.   This is referred to as \"collapsing\" or \"reducing\" the chosen dimension.\n",
    "\n",
    "This is done by adding an \"axis\" keyword specifying which dimension, such as `np.min(data, axis=1)`. \n",
    "\n",
    "For example, recall that the above \"daily_records\" data varies over 2 timepoints and 3 locations, mapped to array dimensions 0 and 1 respectively:\n",
    "\n",
    "```python\n",
    "print(daily_records)\n",
    "[[12 14 11]\n",
    " [11 12 15]]\n",
    "```\n",
    "\n",
    "For this data, by adding the `axis` keyword, we can calculate either...\n",
    "  * using `axis=0`: a statistic **over both days, for each station**, or\n",
    "  * using `axis=1`: a statistic **over all three stations, for each day**.\n",
    "\n",
    "For most statistical functions (but not all), the \"axis\" keyword will also accept multiple dimensions, e.g. `axis=(0, 2)`.  \n",
    "Remember also that the *default* statistical operation, as seen above, is to collapses over _all_ dimensions, reducing the array to a single value.  This is equivalent to specifying `axis=None`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise \n",
    "\n",
    "* What other similar statistical operations exist (see above link)?\n",
    "* A mean value can also be calculated with `<array>.mean()`.  Is that the same thing?\n",
    "* produce the two kinds of single-axis average mentioned above -- that is:\n",
    "   1. over all days at each station, \n",
    "   2. over all stations for each day.\n",
    "* Create a sample 3-D array (for example, this could represent air temperature at given [time, X, Y] positions):\n",
    "   1. how can you form a mean over all X and Y at each timestep?  \n",
    "   1. what shape does the result have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masked Arrays\n",
    "\n",
    "Real-world measurements processes often result in certain datapoint values being uncertain or simply \"missing\".  This is usually indicated by additional data quality information, stored alongside the data values.\n",
    "\n",
    "In these cases we often need to make calculations that count only the ***valid*** datapoints.  NumPy provides a special \"masked array\" type for this type of calculation. Here's a link to the documentation for NumPy masked arrays: https://docs.scipy.org/doc/numpy-1.11.0/reference/maskedarray.generic.html#maskedarray-generic-constructing.\n",
    "\n",
    "For example, we might know that in our previous `daily_records` data, the value for station-2 on day-1 is invalid.\n",
    "To represent this missing datapoint, we can make a masked version of the data :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "daily_records = np.array([[12, 14, 11], [11, 12, 15]])\n",
    "masked_data = np.ma.masked_array(daily_records)\n",
    "masked_data[0, 1] = np.ma.masked\n",
    "print('masked data:')\n",
    "print(masked_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The statistics of the masked data version are different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('unmasked average = ', np.mean(daily_records))\n",
    "print('masked average = ', np.ma.mean(masked_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `np.ma.masked_array()` function seen above is a simple creation method for masked data.  \n",
    "The sub-module `np.ma` contains all the NumPy masked-data routines.\n",
    "\n",
    "Instead of masking selected points, as above, a mask is often specified as a \"mask array\":  This is a boolean array of the same size and shape as the data, where `False` means a good datapoint and `True` means a masked or missing point.  Such a 'boolean mask' can be passed in the `np.ma.masked_array` creation method, and can be extracted with the `np.ma.getmaskarray()` function.\n",
    "\n",
    "Note that most file storage formats represent missing data in a _different_ way, using a distinct \"missing data\" value appearing in the data.  \n",
    "There is special support for converting between this type of representation and NumPy masked arrays :  Every masked array has a 'fill_value' property and a 'filled()' method to convert between the two."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "  * Create a masked array from the numbers 0-11, where all the values less than 5 are masked.\n",
    "  * Create a masked array of positive values, with a value of `-1.0` to represent missing points.\n",
    "  * Look up the masked array creation documentation. What routines exist to produce masked arrays like the ones you've just created more efficiently?\n",
    "  * Use `np.ma.filled()` to create a 'plain' (i.e. unmasked) array from a masked array.\n",
    "  * How can you create a plain array from a masked array, but using a _different_ fill-value for masked points?\n",
    "  * Try performing a mathematical operation between two masked arrays. What happens to the 'fill_value' properties when you do so?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Statistics and Masked Arrays: Summary of key points\n",
    " * most statistical functions are available in two different forms, as in `array.mean()` and also `np.mean(array)`,\n",
    "   the choice being mostly a question of style.\n",
    " * statistical operations operate over, and remove (or \"collapse\") the array dimensions that they are applied to.\n",
    " * an \"axis\" keyword specifies operation over dimensions : this can be one; multiple; or all.\n",
    "   * NOTE: not all operations permit operation over specifically selected dimensions\n",
    " * Statistical operations are not really part of NumPy itself, but are defined by the higher-level Scipy project.\n",
    " * Missing datapoints can be represented using \"masked arrays\"\n",
    "   * these are useful for calculation, but usually require converting to another form for data storage\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Efficiency <a class=\"anchor\" id=\"efficiency\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will be able to apply concepts that allow you to perform fast and efficient calculations on NumPy arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Views on Arrays\n",
    "\n",
    "NumPy attempts to not make copies of arrays. Many NumPy operations will produce a reference to an existing array, known as a \"view\", instead of making a whole new array. For example, indexing and reshaping provide a view of the same memory wherever possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.arange(8)\n",
    "arr_view = arr.reshape(2, 4)\n",
    "\n",
    "# Print the \"view\" array from reshape.\n",
    "print('Before\\n{}'.format(arr_view))\n",
    "\n",
    "# Update the first element of the original array.\n",
    "arr[0] = 1000\n",
    "\n",
    "# Print the \"view\" array from reshape again,\n",
    "# noticing the first value has changed.\n",
    "print('After\\n{}'.format(arr_view))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this means is that if one array (`arr`) is modified, the other (`arr_view`) will also be updated :  the same memory is being shared. This is a valuable tool which enables the system memory overhead to be managed, which is particularly useful when handling lots of large arrays. The lack of copying allows for very efficient vectorized operations.\n",
    "\n",
    "Remember, this behaviour is automatic in most of NumPy, so it requires some consideration in your code, it can lead to some bugs that are hard to track down. For example, if you are changing some elements of an array that you are using elsewhere, you may want to explicitly copy that array before making changes. If in doubt, you can always copy the data to a different block of memory with the `copy()` method.\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.arange(8)\n",
    "arr_view = arr.reshape(2, 4).copy()\n",
    "\n",
    "# Print the \"view\" array from reshape.\n",
    "print('Before\\n{}'.format(arr_view))\n",
    "\n",
    "# Update the first element of the original array.\n",
    "arr[0] = 1000\n",
    "\n",
    "# Print the \"view\" array from reshape again,\n",
    "# noticing the first value has changed.\n",
    "print('After\\n{}'.format(arr_view))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loops and Vectorised Operations\n",
    "\n",
    "We will now explore calculation performance and consider efficiency in terms of processing time.\n",
    "\n",
    "Firstly let's look at a simple processing time tool that is provided in notebooks; ``%%timeit`` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit \n",
    "x = range(500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat that, specifying only 100 loops and fastest time of 5 runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 100 -r 5\n",
    "x = range(500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us an easy way to evaluate performance for implementations ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rands = np.random.random(1000000).reshape(100, 100, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 10 -r 5\n",
    "overPointEightLoop = 0\n",
    "for i in range(100):\n",
    "    for j in range(100):\n",
    "        for k in range(100):\n",
    "            if rands[i, j, k] > 0.8:\n",
    "                overPointEightLoop +=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 10 -r 5\n",
    "overPointEightWhere = rands[rands > 0.8].size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Efficiency: Summary of key  points\n",
    "\n",
    " * NumPy enables vectorised calculations that are fast:\n",
    "   * Python loops are much slower.\n",
    " * NumPy can only use the system memory available to hold arrays:\n",
    "   * very large arrays can result in a `MemoryError`.\n",
    " * Many NumPy operations return views on exisiting array:\n",
    "   * a view shares memory with the original array,\n",
    "   * changing array contents in place will affect all related views,\n",
    "   * sometimes it is useful to explicitly `copy` arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: trapezoidal integration\n",
    "\n",
    "In this exercise, you are tasked with implementing the simple trapezoid rule\n",
    "formula for numerical integration. If we want to compute the definite integral\n",
    "\n",
    "$$\n",
    "     \\int_{a}^{b}f(x)dx\n",
    "$$\n",
    "\n",
    "we can partition the integration interval $[a,b]$ into smaller subintervals. We then\n",
    "approximate the area under the curve for each subinterval by calculating the area of\n",
    "the trapezoid created by linearly interpolating between the two function values\n",
    "at each end of the subinterval:\n",
    "\n",
    "![Illustration of the trapezoidal rule](../images/trapezoidal_rule.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a pre-computed $y$ array (where $y = f(x)$ at discrete samples) the trapezoidal rule equation is:\n",
    "\n",
    "$$\n",
    "   \\int_{a}^{b}f(x)dx\\approx\\frac{1}{2}\\sum_{i=1}^{n}\\left(x_{i}-x_{i-1}\\right)\\left(y_{i}+y_{i-1}\\right).\n",
    "$$\n",
    "\n",
    "In pure python, this can be written as:\n",
    "\n",
    "    def trapz_slow(x, y):\n",
    "        area = 0.\n",
    "        for i in range(1, len(x)):\n",
    "            area += (x[i] - x[i-1]) * (y[i] + y[i-1])\n",
    "        return area / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 1\n",
    "\n",
    "Create two arrays $x$ and $y$, where $x$ is a linearly spaced array in the interval $[0, 3]$ of length 11, and $y$ represents the function $f(x) = x^2$ sampled at $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 2\n",
    "\n",
    "Use indexing (not a for loop) to find the 10 values representing $y_{i}+y_{i-1}$ for $i$ between 1 and 11.\n",
    "\n",
    "Hint: What indexing would be needed to get all but the last element of the 1d array **``y``**. Similarly what indexing would be needed to get all but the first element of a 1d array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 3\n",
    "\n",
    "Write a function `trapz(x, y)`, that applies the trapezoid formula to pre-computed values, where `x` and `y` are 1-d arrays. The function should not use a for loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 4\n",
    "\n",
    "Verify that your function is correct by using the arrays created in #1 as input to ``trapz``. Your answer should be a close approximation of $\\int_0^3 x^2$ which is $9$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 5 (extension)\n",
    "\n",
    "``numpy`` and ``scipy.integrate`` provide many common integration schemes.  Find the documentation for NumPy's own version of the trapezoidal integration scheme and check its result with your own:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 6 (extension)\n",
    "\n",
    "Write a function  `trapzf(f, a, b, npts=100)` that accepts a function `f`, the endpoints `a` and `b` and the number of samples to take `npts`.  Sample the function uniformly at these\n",
    "points and return the value of the integral.\n",
    "\n",
    "Use the trapzf function to identify the minimum number of sampling points needed to approximate the integral $\\int_0^3 x^2$ with an absolute error of $<=0.0001$. (A loop is necessary here.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Reading <a class=\"anchor\" id=\"further_reading\"></a>\n",
    "\n",
    "Here are a selection of topics we chose not to include in this workshop, or only touched upon briefly.  You may come across some of them or find these links and terms useful in the future.\n",
    "\n",
    "### Concepts:\n",
    "\n",
    "  * [NaN](https://docs.scipy.org/doc/numpy/user/misc.html)s and NaN functions\n",
    "  * [Random number arrays](https://docs.scipy.org/doc/numpy-1.11.0/reference/routines.random.html)\n",
    "  * matrix multiplication (for example using [dot](https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html))\n",
    "  * C and Fortran storage order -- see \"C order\" and \"Fortran order\" in [the glossary](https://docs.scipy.org/doc/numpy/glossary.html)\n",
    "  * [file I/O](https://docs.scipy.org/doc/numpy/reference/routines.io.html)\n",
    "  * [Universal Functions (ufunc)](https://docs.scipy.org/doc/numpy/reference/ufuncs.html)\n",
    "  * [indexing with arrays of indices](https://docs.scipy.org/doc/numpy/user/basics.indexing.html#index-arrays)\n",
    "\n",
    "### Functions:\n",
    "\n",
    "  * [empty](https://docs.scipy.org/doc/numpy/reference/generated/numpy.empty.html)\n",
    "  * [stack](https://docs.scipy.org/doc/numpy/reference/generated/numpy.stack.html), [concatenate](https://docs.scipy.org/doc/numpy/reference/generated/numpy.concatenate.html#numpy.concatenate) and [split](https://docs.scipy.org/doc/numpy/reference/generated/numpy.split.html#numpy.split)\n",
    "  * [transpose](https://docs.scipy.org/doc/numpy/reference/generated/numpy.transpose.html) ([`array.T`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.T.html#numpy.ndarray.T))\n",
    "  * [roll](https://docs.scipy.org/doc/numpy/reference/generated/numpy.roll.html) and [rollaxis](https://docs.scipy.org/doc/numpy/reference/generated/numpy.rollaxis.html#numpy.rollaxis)"
   ]
  }
 ],
 "metadata": {
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
<commit_msg>Shorten NumPy course (#139)<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A Workshop Introduction to NumPy\n",
    "\n",
    "The Python language is an excellent tool for general-purpose programming, with a highly readable syntax, rich and powerful data types (strings, lists, sets, dictionaries, arbitrary length integers, etc) and a very comprehensive standard library. It was not, however, designed specifically for mathematical and scientific computing. Neither the language nor its standard library have facilities for the efficient representation of multidimensional datasets, tools for linear algebra and general matrix manipulations (essential building blocks of virtually all scientific computing).\n",
    "\n",
    "For example, Python lists are very flexible containers that can be nested arbitrarily deep and can hold any Python object in them, but they are poorly suited to represent common mathematical constructs like vectors and matrices. It is for this reason that NumPy exists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Workshop Aims\n",
    "\n",
    "The aim of this workshop is to enable you use NumPy to:\n",
    "\n",
    "  - manipulate numerical arrays, and\n",
    "  - perform efficient array calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Table of Contents\n",
    "\n",
    "* [Getting Started with NumPy](#getting-started)\n",
    "* [Motivation: what are arrays good for?](#motivation)\n",
    "* [The Array Object](#array-object)\n",
    "* [Application: Arithmetic and Broadcasting](#arithmetic_and_broadcasting)\n",
    "* [Application: Statistics and Masked Arrays](#statistics)\n",
    "* [Application: Efficiency](#efficiency)\n",
    "* [Further Reading](#further_reading)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting Started with NumPy <a class=\"anchor\" id=\"getting-started\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will be able to create NumPy arrays with the aid of the NumPy documentation.\n",
    "\n",
    "NumPy is the fundamental package for scientific computing with Python. Its primary purpose is to provide a powerful N-dimensional array object; the focus for this workshop.\n",
    "\n",
    "To begin with let's import NumPy, check where it is being imported from and check the version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NumPy is generally imported as 'np'.\n",
    "import numpy as np\n",
    "print(np)\n",
    "print(np.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Documentation\n",
    "Here is a link to the NumPy documentation for v1.11:\n",
    "https://docs.scipy.org/doc/numpy-1.11.0/reference/\n",
    "\n",
    "There are many online forums with tips and suggestions for solving problems with NumPy, such as http://stackoverflow.com/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explore: creating arrays\n",
    "\n",
    "NumPy provides many different ways to create arrays. These are listed in the documentation at: https://docs.scipy.org/doc/numpy-1.11.0/user/basics.creation.html#arrays-creation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "##### Part 1\n",
    "\n",
    "Try out some of these ways of creating NumPy arrays. See if you can produce:\n",
    "\n",
    "* a NumPy array from a list of numbers,\n",
    "* a 3-dimensional NumPy array filled with a constant value -- either `0` or `1`,\n",
    "* a NumPy array filled with a constant value -- **not** `0` or `1`. (**Hint:** this can be achieved using the last array you created, or you could use `np.empty` and find a way of filling the array with a constant value),\n",
    "* a NumPy array of 8 elements with a range of values starting from `0` and a spacing of `3` between each element, and\n",
    "* a NumPy array of 10 elements that are logarithmically spaced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part 2\n",
    "\n",
    "How could you change the shape of the 8-element array you created previously to have shape `(2, 2, 2)`? **Hint:** this can be done without creating a new array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## Motivation: what are arrays good for? <a class=\"anchor\" id=\"motivation\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will understand key benefits of using NumPy arrays for numerical processing in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extensive Features\n",
    "\n",
    "\n",
    "NumPy provides routines for fast operations on arrays, including mathematical, logical, shape manipulation, sorting, selecting, I/O, discrete Fourier transforms, basic linear algebra, basic statistical operations, random simulation and much more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fast Calculations\n",
    "\n",
    "It is a lot faster than Python alone for numerical computing tasks.\n",
    "\n",
    "Element-by-element operations are the “default mode” when an ndarray is involved, but the element-by-element operation is speedily executed by pre-compiled C code. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clear Syntax\n",
    "\n",
    "In NumPy,\n",
    "\n",
    "```python\n",
    "c = a * b\n",
    "```\n",
    "\n",
    "calculates the element-wise product of `a` and `b`, at near-C speeds, but with the code simplicity we expect from the Python language.\n",
    "\n",
    "This demonstrates a core feature of NumPy, called **vectorization**. This removes the need for loops iterating through elements of arrays, which can make code easier to read, as well as performing fast calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interfacing to other Libraries\n",
    "\n",
    "Many scientific Python libraries use NumPy as their core array representation. From plotting libraries such as matplotlib, to parallel processing libraries such as Dask, to data interoperability libraries such as Iris, NumPy arrays are at the core of how these libraries operate and communicate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {},
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "## The Array Object <a class=\"anchor\" id=\"array-object\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will be able to manipulate NumPy array objects through indexing and explain key features and properties of NumPy arrays.\n",
    "\n",
    "The multidimensional array object is at the core of all of NumPy's functionality.  Let's explore this object some more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array properties\n",
    "\n",
    "Let's create a NumPy array and take a look at some of its properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.ones((3, 2, 4))\n",
    "\n",
    "print(\"Array shape:\", arr.shape)\n",
    "print(\"Array element dtype:\", arr.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "##### Part 1\n",
    "\n",
    "Use the code cells below to explore features of this array.\n",
    "\n",
    "* What is the `type` of the array created above? (Hint: you can find the type of an object in Python using `type(<object>)`, where `<object>` is replaced with the name of the variable containing the object.)\n",
    "* Look this type up in the NumPy documentation (see https://docs.scipy.org/doc/numpy-1.11.0/reference/generated/numpy.ndarray.html#numpy.ndarray)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part 2\n",
    "\n",
    "Consider the following NumPy array properties:\n",
    "\n",
    "* `ndim`,\n",
    "* `nbytes`,\n",
    "* `size`,\n",
    "* `T`.\n",
    "\n",
    "For each of these properties:\n",
    "\n",
    "* Use information from the array documentation to find out more about each array property.\n",
    "* Apply each property to the array `arr` defined above. Can you explain the results you get in each case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexing\n",
    "\n",
    "You can index NumPy arrays in the same way as other Python objects, by using square brackets `[]`. This means we can index to retrieve a single element, multiple consecutive elements, or a more complex sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.array([1, 2, 3, 4, 5, 6])\n",
    "print(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"arr[2] = {}\".format(arr[2]))\n",
    "print(\"arr[2:5] = {}\".format(arr[2:5]))\n",
    "print(\"arr[::2] = {}\".format(arr[::2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also index multidimensional arrays using an enhanced indexing syntax, which allows for multi-element indexing. (Sometimes this is referred to as \"extended slicing\".) Remember that Python uses zero-based indexing!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lst_2d = [[1, 2, 3], [4, 5, 6]]\n",
    "arr_2d = arr.reshape(2, 3)\n",
    "\n",
    "print(\"2D list:\")\n",
    "print(lst_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"2D array:\")\n",
    "print(arr_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Single array element:\")\n",
    "print(arr_2d[1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Single row:\")\n",
    "print(arr_2d[1])\n",
    "print(\"First two columns:\")\n",
    "print(arr_2d[:, :2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you only provide one index to slice a multidimensional array, then the slice will be expanded to \"`:`\" for all of the remaining dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Second row: {} is equivalent to {}'.format(arr_2d[1], arr_2d[1, :]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is known as **ellipsis**. Ellipsis can be specified explicitly using \"`...`\", which automatically expands to \"`:`\" for each dimension unspecified in the slice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.empty((4, 6, 3))\n",
    "print('Original shape: ', arr1.shape)\n",
    "print(arr1[...].shape)\n",
    "print(arr1[..., 0:2].shape)\n",
    "print(arr1[2:4, ..., ::2].shape)\n",
    "print(arr1[2:4, :, ..., ::-1].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Boolean Indexing\n",
    "\n",
    "NumPy provides syntax to index conditionally, based on the data in the array.\n",
    "\n",
    "You can pass in an array of True and False values (a boolean array), or, more commonly, a condition that returns a boolean array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(arr_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bools = arr_2d % 2 == 0\n",
    "print(bools)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(arr_2d[bools])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "##### Part 1\n",
    "\n",
    "Why do these indexing examples give the stated results?\n",
    "\n",
    " * result of `arr_2d[1, 0]` is `4`\n",
    " * result of `arr_2d[0]` is `[1, 2, 3]`\n",
    " * result of `arr_2d[1, 1:]` is `[5, 6]`\n",
    " * result of `arr_2d[0:, ::2]` is `[[1, 3], [4, 6]]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part 2\n",
    "\n",
    "How would you index `arr_2d` to retrieve:\n",
    "\n",
    " * the third value: resulting in `3`\n",
    " * the second row: resulting in `[4 5 6]` \n",
    " * the first column: resulting in `[1 4]`\n",
    " * the first column, retaining the outside dimension: resulting in `[[1] [4]]`\n",
    " * only values greater than or equal to 3: resulting in `[3 4 5 6]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arrays are not lists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: why do the following examples produce different results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lst_2d[0:2][1])\n",
    "print(arr_2d[0:2, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result we just received points to an important piece of learning, which is that in most cases NumPy arrays behave very differently to Python lists. Let's explore the differences (and some similarities) between the two."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### dtype\n",
    "\n",
    "A NumPy array has a fixed data type, called `dtype`.  This is the type of all the elements of the array.  \n",
    "\n",
    "This is in contrast to Python lists, which can hold elements of different types.\n",
    "\n",
    "#### Exercise\n",
    "\n",
    "* What happens in Python when you add an integer to a float?\n",
    "* What happens when you put an integer into a NumPy float array?\n",
    "* What happens when you do numerical calculations between arrays of different types?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating 2D coordinate arrays\n",
    "\n",
    "A common requirement of NumPy arrays is to generate arrays that represent the coordinates of our data.\n",
    "\n",
    "When orthogonal 1d coordinate arrays already exist, NumPy's meshgrid function is very useful:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 9, 3)\n",
    "y = np.linspace(4, 8, 3)\n",
    "x2d, y2d = np.meshgrid(x, y)\n",
    "print(x2d)\n",
    "print(y2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The Array Object: Summary of key  points\n",
    " * properties : `shape`, `dtype`.\n",
    " * arrays are homogeneous; all elements have the same type: `dtype`.\n",
    " * indexing arrays to produce further arrays: views on the original arrays.\n",
    " * multi-dimensional indexing (slicing) and boolean indexing.\n",
    " * combinations to form 2D arrays: `meshgrid`\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Arithmetic and Broadcasting<a class=\"anchor\" id=\"app-calc\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section you will be able to explain how broadcasting allows mathematical operations between NumPy arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elementwise Arithmetic <a class=\"anchor\" id=\"arithmetic_and_broadcasting\"></a>\n",
    "\n",
    "You can use NumPy to perform arithmetic operations between two arrays in an element-by-element fashion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.arange(4)\n",
    "arr2 = np.arange(4)\n",
    "\n",
    "print('{} + {} = {}'.format(arr1, arr2, arr1 + arr2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise \n",
    "\n",
    "Explore arithmetic operations between arrays:\n",
    "\n",
    "* Create a number of arrays of different shapes and dtypes. Make sure some of them include values of `0`.\n",
    "* Make use of all the basic Python mathematical operators (i.e. `+`, `-`, `*`, `/`, `//`, `%`).\n",
    "\n",
    "What operations work? What operations do not work? Can you explain why operations do or do not work?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It makes intrinsic sense that you should be able to add a constant to all values in an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.arange(4)\n",
    "const = 5\n",
    "\n",
    "print(\"Original array: {}\".format(arr))\n",
    "print(\"\")\n",
    "print(\"Array + const: {}\".format(arr + const))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting\n",
    "\n",
    "There are times when you need to perform calculations between NumPy arrays of different sizes. NumPy allows you to do this easily using a powerful piece of functionality called **broadcasting**.\n",
    "\n",
    "Broadcasting is a way of treating the arrays as if they had the same dimensions and thus have elements all corresponding.  It is then easy to perform the calculation, element-wise. It does this by matching dimensions in one array to the other where possible, and using repeated values where there is no corresponding dimension in the other array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rules of Broadcasting\n",
    "\n",
    "Broadcasting applies these three rules:\n",
    "\n",
    "1.    If the two arrays differ in their number of dimensions, the shape of the array with fewer dimensions is padded with ones on its leading (left) side.\n",
    "\n",
    "1.    If the shape of the two arrays does not match in any dimension, either array with shape equal to 1 in a given dimension is stretched to match the other shape.\n",
    "\n",
    "1.    If in any dimension the sizes disagree and neither has shape equal to 1, an error is raised.\n",
    "\n",
    "Note that all of this happens without ever actually creating the expanded arrays in memory! This broadcasting behavior is in practice enormously powerful, especially given that when NumPy broadcasts to create new dimensions or to 'stretch' existing ones, it doesn't actually duplicate the data. In the example above the operation is carried out as if the scalar 1.5 was a 1D array with 1.5 in all of its entries, but no actual array is ever created. This can save lots of memory in cases when the arrays in question are large. As such this can have significant performance implications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Illustration of broadcasting](../images/fig_broadcast_visual_1.png)\n",
    "\n",
    "([image source](http://www.astroml.org/book_figures/appendix/fig_broadcast_visual.html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "For the following cases:\n",
    "\n",
    "* What will be the result of adding `arr1` to `arr2`? \n",
    "* What will be the shape of the resulting array?\n",
    "* What rules of broadcasting are being used?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((2, 3))\n",
    "arr2 = np.ones((2, 1))\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((2, 3))\n",
    "arr2 = np.ones(3)\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((1, 3))\n",
    "arr2 = np.ones((2, 1))\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.ones((1, 3))\n",
    "arr2 = np.ones((1, 2))\n",
    "\n",
    "# (arr1 + arr2).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Reshaping arrays to aid broadcasting\n",
    "\n",
    "NumPy allows you to change the shape of an array, so long as the total number of elements in the array does not change. For example, we could reshape a flat array with 12 elements to a 2D array with shape `(2, 6)`, or `(3, 2, 2)`, or even `(3, 4, 1)`. We could not, however, reshape it to have shape `(2, 5)`, because the total number of elements would not be kept constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise, continued\n",
    "\n",
    "For the failing example above, what reshape operation could you apply to `arr2` so that it can be broadcast with `arr1`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Arithmetic and Broadcasting: Summary of key  points\n",
    "\n",
    " * arithmetic operations are performed in an element-by-element fashion,\n",
    " * operations can be performed between arrays of different shapes,\n",
    " * the arrays' dimensions are aligned according to fixed rules; where one input lacks a given dimension, values are repeated,\n",
    " * reshaping can be used to get arrays to combine as required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Statistics and Masked Arrays <a class=\"anchor\" id=\"statistics\"></a>\n",
    "\n",
    "**Learning outcome:** By the end of this section, you will be able to apply statistical operations and masked arrays to real-world problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics\n",
    "\n",
    "Numpy arrays support many common statistical calculations. For a list of common operations, see: https://docs.scipy.org/doc/numpy/reference/routines.statistics.html.\n",
    "\n",
    "The simplest operations consist of calculating a single statistical value from an array of numbers -- such as a mean value, a variance or a minimum.\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(12).reshape((3, 4))\n",
    "mean = np.mean(a)\n",
    "print(a)\n",
    "print(mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Used without any further arguments, statistical functions simply reduce the whole array to a single value.  In practice, however, we very often want to calculate statistics over only *some* of the dimensions. The most common requirement is to calculate a statistic along a single array dimension, while leaving all the other dimensions intact.   This is referred to as \"collapsing\" or \"reducing\" the chosen dimension.\n",
    "\n",
    "This is done by adding an \"axis\" keyword specifying the dimension to collapse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean(a, axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise \n",
    "\n",
    "* What other similar statistical operations exist (see above link)?\n",
    "* A mean value can also be calculated with `<array>.mean()`.  Is that the same thing?\n",
    "* Create a 3D array (that could be considered to describe `[time, x, y]`) and find the mean over all `x` and `y` at each timestep.\n",
    "* What shape does the result have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masked Arrays\n",
    "\n",
    "Real-world measurements processes often result in certain datapoint values being uncertain or simply \"missing\".  This is usually indicated by additional data quality information, stored alongside the data values.\n",
    "\n",
    "In these cases we often need to make calculations that count only the valid datapoints.  NumPy provides a special \"masked array\" type for this type of calculation. Here's a link to the documentation for NumPy masked arrays: https://docs.scipy.org/doc/numpy-1.11.0/reference/maskedarray.generic.html#maskedarray-generic-constructing.\n",
    "\n",
    "To construct a masked array we need some data and a mask. The data can be any kind of NumPy array, but the mask array must contain a boolean-type values only (either `True` and `False` or `1` and `0`). Let's make each of these and convert them together into a NumPy masked array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.arange(4)\n",
    "mask = np.array([0, 0, 1, 0])\n",
    "print('Data: {}'.format(data))\n",
    "print('Mask: {}'.format(mask))\n",
    "masked_data = np.ma.masked_array(data, mask)\n",
    "print('Masked data: {}'.format(masked_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mask is applied where the values in the mask array are **`True`**. Masked arrays are printed with a double-dash `--` denoting the locations in the array where the mask has been applied.\n",
    "\n",
    "The statistics of masked data are different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Unmasked average: {}'.format(np.mean(data)))\n",
    "print('Masked average: {}'.format(np.ma.mean(masked_data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that most file formats represent missing data in a _different_ way, using a distinct \"missing data\" value appearing in the data. There is special support for converting between this type of representation and NumPy masked arrays. Every masked array has a `fill_value` property and a `filled()` method to fill the masked points with the fill value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "  * Create a masked array from the numbers 0-11, where all the values less than 5 are masked.\n",
    "  * Create a masked array of positive values, with a value of `-1.0` to represent missing points.\n",
    "  * Look up the masked array creation documentation. What routines exist to produce masked arrays like the ones you've just created more efficiently?\n",
    "  * Use `np.ma.filled()` to create a 'plain' (i.e. unmasked) array from a masked array.\n",
    "  * How can you create a plain array from a masked array, but using a _different_ fill-value for masked points?\n",
    "  * Try performing a mathematical operation between two masked arrays. What happens to the 'fill_value' properties when you do so?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Statistics and Masked Arrays: Summary of key points\n",
    " * most statistical functions are available in two different forms, as in `array.mean()` and also `np.mean(array)`,\n",
    "   the choice being mostly a question of style.\n",
    " * statistical operations operate over, and remove (or \"collapse\") the array dimensions that they are applied to.\n",
    " * an \"axis\" keyword specifies operation over dimensions : this can be one; multiple; or all.\n",
    "   * NOTE: not all operations permit operation over specifically selected dimensions\n",
    " * Statistical operations are not really part of NumPy itself, but are defined by the higher-level Scipy project.\n",
    " * Missing datapoints can be represented using \"masked arrays\"\n",
    "   * these are useful for calculation, but usually require converting to another form for data storage\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Efficiency <a class=\"anchor\" id=\"efficiency\"></a>\n",
    "\n",
    "**Learning outcome:** by the end of this section, you will be able to apply concepts that allow you to perform fast and efficient calculations on NumPy arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Views on Arrays\n",
    "\n",
    "NumPy attempts to not make copies of arrays. Many NumPy operations will produce a reference to an existing array, known as a \"view\", instead of making a whole new array. For example, indexing and reshaping provide a view of the same memory wherever possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.arange(8)\n",
    "arr_view = arr.reshape(2, 4)\n",
    "\n",
    "# Print the \"view\" array from reshape.\n",
    "print('Before\\n{}'.format(arr_view))\n",
    "\n",
    "# Update the first element of the original array.\n",
    "arr[0] = 1000\n",
    "\n",
    "# Print the \"view\" array from reshape again,\n",
    "# noticing the first value has changed.\n",
    "print('After\\n{}'.format(arr_view))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this means is that if one array (`arr`) is modified, the other (`arr_view`) will also be updated :  the same memory is being shared. This is a valuable tool which enables the system memory overhead to be managed, which is particularly useful when handling lots of large arrays. The lack of copying allows for very efficient vectorized operations.\n",
    "\n",
    "Remember, this behaviour is automatic in most of NumPy, so it requires some consideration in your code, it can lead to some bugs that are hard to track down. For example, if you are changing some elements of an array that you are using elsewhere, you may want to explicitly copy that array before making changes. If in doubt, you can always copy the data to a different block of memory with the `copy()` method.\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.arange(8)\n",
    "arr_view = arr.reshape(2, 4).copy()\n",
    "\n",
    "# Print the \"view\" array from reshape.\n",
    "print('Before\\n{}'.format(arr_view))\n",
    "\n",
    "# Update the first element of the original array.\n",
    "arr[0] = 1000\n",
    "\n",
    "# Print the \"view\" array from reshape again,\n",
    "# noticing the first value has changed.\n",
    "print('After\\n{}'.format(arr_view))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loops and Vectorised Operations\n",
    "\n",
    "We will now explore calculation performance and consider efficiency in terms of processing time.\n",
    "\n",
    "Firstly let's look at a simple processing time tool that is provided in notebooks; ``%%timeit`` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit \n",
    "x = range(500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat that, specifying only 100 loops and fastest time of 5 runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 100 -r 5\n",
    "x = range(500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us an easy way to evaluate performance for implementations ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rands = np.random.random(1000000).reshape(100, 100, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 10 -r 5\n",
    "overPointEightLoop = 0\n",
    "for i in range(100):\n",
    "    for j in range(100):\n",
    "        for k in range(100):\n",
    "            if rands[i, j, k] > 0.8:\n",
    "                overPointEightLoop +=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 10 -r 5\n",
    "overPointEightWhere = rands[rands > 0.8].size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Efficiency: Summary of key  points\n",
    "\n",
    " * NumPy enables vectorised calculations that are fast:\n",
    "   * Python loops are much slower.\n",
    " * NumPy can only use the system memory available to hold arrays:\n",
    "   * very large arrays can result in a `MemoryError`.\n",
    " * Many NumPy operations return views on exisiting array:\n",
    "   * a view shares memory with the original array,\n",
    "   * changing array contents in place will affect all related views,\n",
    "   * sometimes it is useful to explicitly `copy` arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: trapezoidal integration\n",
    "\n",
    "In this exercise, you are tasked with implementing the simple trapezoid rule\n",
    "formula for numerical integration. If we want to compute the definite integral\n",
    "\n",
    "$$\n",
    "     \\int_{a}^{b}f(x)dx\n",
    "$$\n",
    "\n",
    "we can partition the integration interval $[a,b]$ into smaller subintervals. We then\n",
    "approximate the area under the curve for each subinterval by calculating the area of\n",
    "the trapezoid created by linearly interpolating between the two function values\n",
    "at each end of the subinterval:\n",
    "\n",
    "![Illustration of the trapezoidal rule](../images/trapezoidal_rule.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a pre-computed $y$ array (where $y = f(x)$ at discrete samples) the trapezoidal rule equation is:\n",
    "\n",
    "$$\n",
    "   \\int_{a}^{b}f(x)dx\\approx\\frac{1}{2}\\sum_{i=1}^{n}\\left(x_{i}-x_{i-1}\\right)\\left(y_{i}+y_{i-1}\\right).\n",
    "$$\n",
    "\n",
    "In pure python, this can be written as:\n",
    "\n",
    "    def trapz_slow(x, y):\n",
    "        area = 0.\n",
    "        for i in range(1, len(x)):\n",
    "            area += (x[i] - x[i-1]) * (y[i] + y[i-1])\n",
    "        return area / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 1\n",
    "\n",
    "Create two arrays $x$ and $y$, where $x$ is a linearly spaced array in the interval $[0, 3]$ of length 11, and $y$ represents the function $f(x) = x^2$ sampled at $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 2\n",
    "\n",
    "Use indexing (not a for loop) to find the 10 values representing $y_{i}+y_{i-1}$ for $i$ between 1 and 11.\n",
    "\n",
    "Hint: What indexing would be needed to get all but the last element of the 1d array **``y``**. Similarly what indexing would be needed to get all but the first element of a 1d array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 3\n",
    "\n",
    "Write a function `trapz(x, y)`, that applies the trapezoid formula to pre-computed values, where `x` and `y` are 1-d arrays. The function should not use a for loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 4\n",
    "\n",
    "Verify that your function is correct by using the arrays created in #1 as input to ``trapz``. Your answer should be a close approximation of $\\int_0^3 x^2$ which is $9$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 5 (extension)\n",
    "\n",
    "``numpy`` and ``scipy.integrate`` provide many common integration schemes.  Find the documentation for NumPy's own version of the trapezoidal integration scheme and check its result with your own:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part 6 (extension)\n",
    "\n",
    "Write a function  `trapzf(f, a, b, npts=100)` that accepts a function `f`, the endpoints `a` and `b` and the number of samples to take `npts`.  Sample the function uniformly at these\n",
    "points and return the value of the integral.\n",
    "\n",
    "Use the trapzf function to identify the minimum number of sampling points needed to approximate the integral $\\int_0^3 x^2$ with an absolute error of $<=0.0001$. (A loop is necessary here.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Reading <a class=\"anchor\" id=\"further_reading\"></a>\n",
    "\n",
    "Here are a selection of topics we chose not to include in this workshop, or only touched upon briefly.  You may come across some of them or find these links and terms useful in the future.\n",
    "\n",
    "### Concepts:\n",
    "\n",
    "  * [NaN](https://docs.scipy.org/doc/numpy/user/misc.html)s and NaN functions\n",
    "  * [Random number arrays](https://docs.scipy.org/doc/numpy-1.11.0/reference/routines.random.html)\n",
    "  * matrix multiplication (for example using [dot](https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html))\n",
    "  * C and Fortran storage order -- see \"C order\" and \"Fortran order\" in [the glossary](https://docs.scipy.org/doc/numpy/glossary.html)\n",
    "  * [file I/O](https://docs.scipy.org/doc/numpy/reference/routines.io.html)\n",
    "  * [Universal Functions (ufunc)](https://docs.scipy.org/doc/numpy/reference/ufuncs.html)\n",
    "  * [indexing with arrays of indices](https://docs.scipy.org/doc/numpy/user/basics.indexing.html#index-arrays)\n",
    "\n",
    "### Functions:\n",
    "\n",
    "  * [empty](https://docs.scipy.org/doc/numpy/reference/generated/numpy.empty.html)\n",
    "  * [stack](https://docs.scipy.org/doc/numpy/reference/generated/numpy.stack.html), [concatenate](https://docs.scipy.org/doc/numpy/reference/generated/numpy.concatenate.html#numpy.concatenate) and [split](https://docs.scipy.org/doc/numpy/reference/generated/numpy.split.html#numpy.split)\n",
    "  * [transpose](https://docs.scipy.org/doc/numpy/reference/generated/numpy.transpose.html) ([`array.T`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.T.html#numpy.ndarray.T))\n",
    "  * [roll](https://docs.scipy.org/doc/numpy/reference/generated/numpy.roll.html) and [rollaxis](https://docs.scipy.org/doc/numpy/reference/generated/numpy.rollaxis.html#numpy.rollaxis)"
   ]
  }
 ],
 "metadata": {
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

{'fraction_non_alphanumeric': 0.17656669773562148, 'fraction_numerical': 0.009309720927682438, 'mean_word_length': 2.8734284743459058, 'pattern_counts': {'":': 1111, '<': 43, '<?xml version=': 0, '>': 45, 'https://': 40, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3736636', 'n_tokens_mistral': 31084, 'n_tokens_neox': 29859, 'n_words': 9459}