{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson A11 – NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial will focus on the NumPy module. NumPy is fundamental for scientific computing and will provide the basis for our everyday work, as well as offer advanced functionalities for numerical operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refer also to the NumPy [Quickstart](https://docs.scipy.org/doc/numpy/user/quickstart.html) tutorial online or to\n",
    "\n",
    "\"Travis E. Oliphant. __2015__. *Guide to NumPy* (2nd. ed.). CreateSpace Independent Publishing Platform, North Charleston, SC, USA.\"\n",
    "\n",
    "if you prefer books."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To import the module run the cell below. NumPy is a third-party Python module that you might need to install first (e.g. using the conda package manager with `conda install numpy`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.869224Z",
     "start_time": "2021-05-06T06:24:08.856739Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now call any function `x` implemented into NumPy with `np.x()`. Lets generate some lists to get started. You should aready know about lists. Recall that a list in Python is a mutable, ordered collection of objects with no restriction on the type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.891981Z",
     "start_time": "2021-05-06T06:24:08.871528Z"
    }
   },
   "outputs": [],
   "source": [
    "list1 = [1, 3, 5]\n",
    "list2 = [0.0, 2.0, 4.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For numerical operations, this data type is of limited use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arrays basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to turn these lists into *NumPy arrays*. Arrays are a type of collection better suited for certain tasks and are essential for our day-to-day-work. They offer a lot of functionality. To convert these lists to arrays we run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.906625Z",
     "start_time": "2021-05-06T06:24:08.895167Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [1 3 5] <class 'numpy.ndarray'>\n",
      "B =  [0. 2. 4.] <class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "A = np.asarray(list1)\n",
    "B = np.asarray(list2)\n",
    "print('A = ', A, type(A))\n",
    "print('B = ', B, type(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have also created a NumPy array directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.917838Z",
     "start_time": "2021-05-06T06:24:08.909362Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [1 3 5] <class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "A = np.array([1, 3, 5])\n",
    "print('A = ', A, type(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a first glance, the arrays we created look like lists when printed, but they are in fact very different.\n",
    "The type of objects collected in an array is constrained to one specific type. We can discover the data type associated with an array from the `dtype` array-attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.929590Z",
     "start_time": "2021-05-06T06:24:08.920244Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type stored in A:  [1 3 5] int64\n",
      "Type stored in B:  [0. 2. 4.] float64\n"
     ]
    }
   ],
   "source": [
    "print('Type stored in A: ', A, A.dtype)\n",
    "print('Type stored in B: ', B, B.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:** Arrays are mutable, ordered sequences in which all stored objects\n",
    "    need to have the same type.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Advanced:** The Python standard library also offers a type of collection\n",
    "    called array. These arrays are similar to the arrays from the NumPy module,\n",
    "    with respect to the data type restriction, but they are handled\n",
    "    differently. You can use the array module with `import array`. Do not\n",
    "    confuse arrays of the type `numpy.ndarray` with those of the type\n",
    "    `array.array`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of the restriction to only one type, lists and other sequences can only\n",
    "be converted to NumPy arrays straight away if they fulfill this criterion. Otherwise the objects they contain are transformed into a common basic type if possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.947145Z",
     "start_time": "2021-05-06T06:24:08.932010Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['2', 'True', 'x'], dtype='<U21')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flexible_collection = [2, True, \"x\"]\n",
    "np.asarray(flexible_collection)  # Converted to string object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.958663Z",
     "start_time": "2021-05-06T06:24:08.949063Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, True, 'x', None], dtype=object)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flexible_collection = [2, True, \"x\", None]\n",
    "np.asarray(flexible_collection)  # Converted to `object`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T05:45:25.054595Z",
     "start_time": "2021-05-06T05:45:25.039860Z"
    }
   },
   "source": [
    "To quite a large extend arrays can be treated just as other sequences. We can for\n",
    "example iterate over them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:08.971426Z",
     "start_time": "2021-05-06T06:24:08.961843Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "9\n",
      "25\n"
     ]
    }
   ],
   "source": [
    "A = np.array([1, 3, 5])\n",
    "for i in A:\n",
    "    print(i**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can access their elements by indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:09.079062Z",
     "start_time": "2021-05-06T06:24:08.974123Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 5]\n",
      "1\n",
      "5\n",
      "[1 3]\n",
      "[3]\n",
      "[1 5]\n",
      "[5 3 1]\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "index 10 is out of bounds for axis 0 with size 3",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-48fc44548521>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m   \u001b[1;31m# print every second character\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m# print every character in reversed order\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 8\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m    \u001b[1;31m# raises IndexError\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m: index 10 is out of bounds for axis 0 with size 3"
     ]
    }
   ],
   "source": [
    "print(A)        # print the whole string\n",
    "print(A[0])     # print the first character of the string\n",
    "print(A[-1])    # print the last character of the string\n",
    "print(A[:2])    # print the first two characters of the string\n",
    "print(A[1:-1])  # print the second to the last character of the string\n",
    "print(A[::2])   # print every second character\n",
    "print(A[::-1])  # print every character in reversed order\n",
    "print(A[10])    # raises IndexError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The indexing capabilities of arrays go, however, even further. We can for example use another sequence (e.g. another array or list)\n",
    "to specifiy the indices of elements we want to get."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:19.774967Z",
     "start_time": "2021-05-06T06:24:19.769032Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3 7]\n"
     ]
    }
   ],
   "source": [
    "# Using a sequence of indices to index an array\n",
    "A = np.array([1, 3, 5, 7])  # Array of elements\n",
    "B = [1, 3]                  # List of indices\n",
    "print(A[B]) # Get all elements in A with indices given in B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not working for lists for example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:20.516082Z",
     "start_time": "2021-05-06T06:24:20.504970Z"
    }
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "list indices must be integers or slices, not list",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-11-935edb6c284e>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m7\u001b[0m\u001b[1;33m]\u001b[0m  \u001b[1;31m# List of elements\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mB\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m        \u001b[1;31m# List of indices\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m# Can't get elements in A with indices given in B\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: list indices must be integers or slices, not list"
     ]
    }
   ],
   "source": [
    "# Trying to use a sequence of indices to index a list\n",
    "A = [1, 3, 5, 7]  # List of elements\n",
    "B = [1, 3]        # List of indices\n",
    "print(A[B])  # Can't get elements in A with indices given in B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also very useful, we can use a sequence of objects of the boolean type (`True` and `False`)\n",
    "with the same length as the array we want to index to get only specific values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:20.854295Z",
     "start_time": "2021-05-06T06:24:20.848088Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 7]\n"
     ]
    }
   ],
   "source": [
    "# Using a boolean sequence of indices to index an array\n",
    "A = np.asarray([1, 3, 5, 7])\n",
    "B = [True, True, False, True]\n",
    "# A and B must be of same length\n",
    "print(A[B])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:21.036285Z",
     "start_time": "2021-05-06T06:24:21.024870Z"
    }
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "list indices must be integers or slices, not list",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-13-5b4bb0a00302>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m7\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mB\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: list indices must be integers or slices, not list"
     ]
    }
   ],
   "source": [
    "# Trying to use a boolean sequence of indices to index a list\n",
    "A = [1, 3, 5, 7]\n",
    "B = [True, True, False, True]\n",
    "print(A[B])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This extends to testing an array against a condition to get a boolean sequence for indexing (see also [Boolean expressions and np.where() to find them](#sec:Boolean_expr))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:21.662038Z",
     "start_time": "2021-05-06T06:24:21.656285Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3 5 7]\n"
     ]
    }
   ],
   "source": [
    "A = np.asarray([1, 3, 5, 7])\n",
    "print(A[A > 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the most striking advantages of arrays over other sequences, is that we can perform\n",
    "certain operations on all of their elements at the same time. We can add, subtract, multiply or divide scalars and the arrays. The operation will be performed element-wise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:22.060039Z",
     "start_time": "2021-05-06T06:24:22.049941Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A + 5.0 =  [ 6.  8. 10. 12.]\n",
      "A - 5,0 =  [-4. -2.  0.  2.]\n",
      "A * 5.0 =  [ 5. 15. 25. 35.]\n",
      "A / 5.0 =  [0.2 0.6 1.  1.4]\n"
     ]
    }
   ],
   "source": [
    "A = np.asarray([1, 3, 5, 7])\n",
    "print('A + 5.0 = ', A + 5.0)\n",
    "print('A - 5,0 = ', A - 5.0)\n",
    "print('A * 5.0 = ', A * 5.0)\n",
    "print('A / 5.0 = ', A / 5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also add, subtract, multiply or divide arrays (of same length/size). The operation will be performed again element-wise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:22.821734Z",
     "start_time": "2021-05-06T06:24:22.810819Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A + B =  [ 1.   5.   9.  15.5]\n",
      "A - B =  [ 1.   1.   1.  -1.5]\n",
      "A * B =  [ 0.   6.  20.  59.5]\n",
      "A / B =  [       inf 1.5        1.25       0.82352941]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-16-42687086e92f>:6: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  print('A / B = ', A / B)\n"
     ]
    }
   ],
   "source": [
    "A = np.asarray([1, 3, 5, 7])\n",
    "B = np.asarray([0, 2.0, 4.0, 8.5])\n",
    "print('A + B = ', A + B)\n",
    "print('A - B = ', A - B)\n",
    "print('A * B = ', A * B)\n",
    "print('A / B = ', A / B)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The division by zero in the last line is the reason behind the warning message (`RuntimeWarning`). Note that a division by zero outside NumPy would result in a `ZeroDivisionError`. NumPy can handle this and has a type to represent infinity (`np.inf`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:23.139584Z",
     "start_time": "2021-05-06T06:24:23.131913Z"
    }
   },
   "outputs": [
    {
     "ename": "ZeroDivisionError",
     "evalue": "division by zero",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mZeroDivisionError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-17-bc757c3fda29>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;36m1\u001b[0m \u001b[1;33m/\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mZeroDivisionError\u001b[0m: division by zero"
     ]
    }
   ],
   "source": [
    "1 / 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The operations above can be translated to vector operations in linear algebra. Let's define $\\mathbf{a}$ and $\\mathbf{b}$ as the vectors corresponding to the NumPy arrays `A` and `B`.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\mathbf{a} =  \\begin{pmatrix} 1\\\\ 3\\\\ 5\\\\ 7 \\end{pmatrix},~~~~~~~\n",
    "\\mathbf{b} =  \\begin{pmatrix} 0\\\\ 2.0\\\\ 4.0\\\\ 8.57 \\end{pmatrix}\n",
    "\\end{equation*}\n",
    "\n",
    "Let $\\mathbf{e} = \\mathbf{1}$ be a vector of all ones and $\\mathbf{I}$ the identity matrix.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\mathbf{e} =  \\begin{pmatrix} 1\\\\ 1\\\\ 1\\\\ 1 \\end{pmatrix},~~~~~~~\n",
    "\\mathbf{I} =  \\begin{pmatrix} 1& 0& 0& 0\\\\ 0& 1& 0& 0\\\\ 0& 0& 1& 0\\\\ 0& 0& 0& 1 \\end{pmatrix}\n",
    "\\end{equation*}\n",
    "\n",
    "Then, the NumPy operations have the following algebraic expressions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`A + B` $\\rightarrow$ $\\mathbf{a} + \\mathbf{b}$\n",
    "\n",
    "`A - B` $\\rightarrow$ $\\mathbf{a} - \\mathbf{b}$\n",
    "\n",
    "`A * 5` $\\rightarrow$ $5 \\mathbf{a}$\n",
    "\n",
    "`A / 5` $\\rightarrow$ $\\frac{1}{5} \\mathbf{a}$\n",
    "\n",
    "`A + 5` $\\rightarrow$ $\\mathbf{a} + 5\\cdot\\mathbf{e}$\n",
    "\n",
    "`A - 5` $\\rightarrow$ $\\mathbf{a} - 5\\cdot\\mathbf{e}$\n",
    "\n",
    "`A * B` $\\rightarrow$ $\\mathbf{a} \\cdot \\mathbf{I} \\cdot \\mathbf{b}$\n",
    "\n",
    "`A / B` $\\rightarrow$ $\\mathbf{a} \\cdot (\\mathbf{I} \\cdot \\mathbf{b})^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $(\\mathbf{I} \\cdot \\mathbf{b})^{-1}$ being the inverse of $\\mathbf{I} \\cdot \\mathbf{b}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:** Be careful, `A * B` is not the scalar product $\\mathbf{a} \\cdot \\mathbf{b}$. This you get with `np.dot(A, B)` or `A @ B`.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrays come with many handy methods like `sum()`, `mean()`, `std()` etc. Here are some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:** Remember that you can get an overview over all these methods by using\n",
    "    `dir(object)`. You can also explore them by typing `\"A.\"` and then pressing\n",
    "    tab to use the interactive help in a Jupyter notebook.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:24.306038Z",
     "start_time": "2021-05-06T06:24:24.299000Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A.sum() =  16\n",
      "A.max() =  7\n",
      "A.argmax() =  3\n"
     ]
    }
   ],
   "source": [
    "print(\"A.sum() = \", A.sum())\n",
    "# The largest Entry of the Array\n",
    "print(\"A.max() = \", A.max())\n",
    "# The index of the largest entry\n",
    "print(\"A.argmax() = \", A.argmax())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrays can be multidimensional. This is how you represent matrices in NumPy. The dimensionality of an array can be viewed with the `shape` property. We are going to make a 2D-array `C`, which will be the array made from a list of lists. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:24.645082Z",
     "start_time": "2021-05-06T06:24:24.632288Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C = \n",
      "[[1. 3. 5.]\n",
      " [0. 2. 4.]]\n",
      "\n",
      "\n",
      "C.shape (rows, columns) =  (2, 3)\n",
      "\n",
      "\n",
      "C[0] =  [1. 3. 5.]\n",
      "C[1] =  [0. 2. 4.]\n",
      "C[:, 0] =  [1. 0.]\n"
     ]
    }
   ],
   "source": [
    "list1 = [1, 3, 5]\n",
    "list2 = [0.0, 2.0, 4.0]\n",
    "\n",
    "C = np.array([\n",
    "    list1,\n",
    "    list2,\n",
    "    ])\n",
    "\n",
    "print(\"C = \")\n",
    "print(C)\n",
    "print(\"\\n\")\n",
    "print('C.shape (rows, columns) = ', C.shape)\n",
    "print(\"\\n\")\n",
    "print('C[0] = ', C[0])\n",
    "print('C[1] = ', C[1])\n",
    "print('C[:, 0] = ', C[:, 0])  # [rowindex, columnindex]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this is fundamentally different to having a nested list of lists. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:24.973263Z",
     "start_time": "2021-05-06T06:24:24.963422Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nested_list[0][0] (nested index) => 1\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "list indices must be integers or slices, not tuple",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-20-e19173c88ae2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"nested_list[0][0] (nested index) =>\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnested_list\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m \u001b[0mnested_list\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m  \u001b[1;31m# raises TypeError\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: list indices must be integers or slices, not tuple"
     ]
    }
   ],
   "source": [
    "nested_list = [\n",
    "    list1,\n",
    "    list2,\n",
    "    ]\n",
    "\n",
    "print(\"nested_list[0][0] (nested index) =>\", nested_list[0][0])\n",
    "nested_list[:, 0]  # raises TypeError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you depend on fast row- or column-wise indexing, e.g. when working with matrices,\n",
    "NumPy arrays can make your life much easier. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boolean expressions and np.where() to find them"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another essential function that deserves extra mentioning is `np.where()`. In the most simple case, it takes a boolean sequence, i.e. a sequence of `True` and `False`, as an input and returns the indices where the sequence is *true*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:25.665788Z",
     "start_time": "2021-05-06T06:24:25.660079Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function where in module numpy:\n",
      "\n",
      "where(...)\n",
      "    where(condition, [x, y])\n",
      "    \n",
      "    Return elements chosen from `x` or `y` depending on `condition`.\n",
      "    \n",
      "    .. note::\n",
      "        When only `condition` is provided, this function is a shorthand for\n",
      "        ``np.asarray(condition).nonzero()``. Using `nonzero` directly should be\n",
      "        preferred, as it behaves correctly for subclasses. The rest of this\n",
      "        documentation covers only the case where all three arguments are\n",
      "        provided.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    condition : array_like, bool\n",
      "        Where True, yield `x`, otherwise yield `y`.\n",
      "    x, y : array_like\n",
      "        Values from which to choose. `x`, `y` and `condition` need to be\n",
      "        broadcastable to some shape.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    out : ndarray\n",
      "        An array with elements from `x` where `condition` is True, and elements\n",
      "        from `y` elsewhere.\n",
      "    \n",
      "    See Also\n",
      "    --------\n",
      "    choose\n",
      "    nonzero : The function that is called when x and y are omitted\n",
      "    \n",
      "    Notes\n",
      "    -----\n",
      "    If all the arrays are 1-D, `where` is equivalent to::\n",
      "    \n",
      "        [xv if c else yv\n",
      "         for c, xv, yv in zip(condition, x, y)]\n",
      "    \n",
      "    Examples\n",
      "    --------\n",
      "    >>> a = np.arange(10)\n",
      "    >>> a\n",
      "    array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])\n",
      "    >>> np.where(a < 5, a, 10*a)\n",
      "    array([ 0,  1,  2,  3,  4, 50, 60, 70, 80, 90])\n",
      "    \n",
      "    This can be used on multidimensional arrays too:\n",
      "    \n",
      "    >>> np.where([[True, False], [True, True]],\n",
      "    ...          [[1, 2], [3, 4]],\n",
      "    ...          [[9, 8], [7, 6]])\n",
      "    array([[1, 8],\n",
      "           [3, 4]])\n",
      "    \n",
      "    The shapes of x, y, and the condition are broadcast together:\n",
      "    \n",
      "    >>> x, y = np.ogrid[:3, :4]\n",
      "    >>> np.where(x < y, x, 10 + y)  # both x and 10+y are broadcast\n",
      "    array([[10,  0,  0,  0],\n",
      "           [10, 11,  1,  1],\n",
      "           [10, 11, 12,  2]])\n",
      "    \n",
      "    >>> a = np.array([[0, 1, 2],\n",
      "    ...               [0, 2, 4],\n",
      "    ...               [0, 3, 6]])\n",
      "    >>> np.where(a < 4, a, -1)  # -1 is broadcast\n",
      "    array([[ 0,  1,  2],\n",
      "           [ 0,  2, -1],\n",
      "           [ 0,  3, -1]])\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(np.where)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:25.830762Z",
     "start_time": "2021-05-06T06:24:25.824566Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 2]),)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.where([True, False, True])  # returns (x-indices, y-indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is most often used when the boolean sequence is derived from a boolean expression involving an array.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:26.180571Z",
     "start_time": "2021-05-06T06:24:26.173923Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[False  True False False]\n",
      "(array([1]),)\n"
     ]
    }
   ],
   "source": [
    "A = np.array([1, 3, 5, 7])\n",
    "print(A == 3)\n",
    "print(np.where(A == 3))  # returns (x-indices, y-indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is of course also working in two dimensions. Lets see where `C` is greater than or equal to 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:26.552670Z",
     "start_time": "2021-05-06T06:24:26.545732Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 3. 5.]\n",
      " [0. 2. 4.]]\n"
     ]
    }
   ],
   "source": [
    "list1 = [1, 3, 5]\n",
    "list2 = [0.0, 2.0, 4.0]\n",
    "\n",
    "C = np.array([\n",
    "    list1,\n",
    "    list2,\n",
    "    ])\n",
    "\n",
    "print(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be also expressed as the corresponding matrix:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\mathbf{C} = \\begin{pmatrix} 1& 3& 5\\\\ 0& 2& 4 \\end{pmatrix}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:27.098310Z",
     "start_time": "2021-05-06T06:24:27.092815Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([0, 0, 1]), array([1, 2, 2]))\n"
     ]
    }
   ],
   "source": [
    "print(np.where(C >= 3))  # returns (x-indices, y-indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result are two arrays with three indices each. We get *two arrays* because `C` is *two-dimensional*. Each array has *three entries* because in total there are *three hits*. \n",
    "This means that the statement is `True` for the indices `C[0, 1]`, `C[0, 2]` and `C[1, 2]`. Lets check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:27.473754Z",
     "start_time": "2021-05-06T06:24:27.466788Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C[0, 1] = 3.0\n",
      "C[0, 2] = 5.0\n",
      "C[1, 2] = 4.0\n"
     ]
    }
   ],
   "source": [
    "# All greater than three?\n",
    "print(\"C[0, 1] =\", C[0, 1])\n",
    "print(\"C[0, 2] =\", C[0, 2]) \n",
    "print(\"C[1, 2] =\", C[1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since indices in Python are zero-based but they are 1-based in terms of linear algebra, we need to shift the indices by one to get the corresponding matrix elements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{C}_{12} = 3$\n",
    "\n",
    "$\\mathbf{C}_{13} = 5$\n",
    "\n",
    "$\\mathbf{C}_{23} = 4$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advanced"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optionally, we can pass two additional sequences to `np.where` from which values are returned\n",
    "instead of just the *true*-indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:28.408851Z",
     "start_time": "2021-05-06T06:24:28.402862Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['yes' 'no' 'yes']\n"
     ]
    }
   ],
   "source": [
    "print(np.where([True, False, True], [\"yes\", \"yes\", \"yes\"], [\"no\", \"no\", \"no\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or simply:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:28.762355Z",
     "start_time": "2021-05-06T06:24:28.756712Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['yes', 'no', 'yes'], dtype='<U3')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.where([True, False, True], \"yes\", \"no\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Special array initialisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen how to create arrays from sequences and from scratch. NumPy provides\n",
    "a handful of convenience functions to construct arrays filled with all kinds of values.\n",
    "There is for example the function `np.linspace`, which takes in principle three arguments\n",
    "`(start, stop, num)` and returns a serial array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:29.334705Z",
     "start_time": "2021-05-06T06:24:29.328841Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]\n"
     ]
    }
   ],
   "source": [
    "print(np.linspace(0, 10, 11))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have also the possibilty to create arrays full of zeros, ones, or another value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:29.700464Z",
     "start_time": "2021-05-06T06:24:29.694353Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "['x' 'x' 'x' 'x' 'x' 'x' 'x' 'x' 'x' 'x']\n"
     ]
    }
   ],
   "source": [
    "print(np.zeros(10))\n",
    "print(np.ones(10))\n",
    "print(np.full(10, \"x\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Meshgrids "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following: imagine you want to compute some function $f$ that depends on spherical coordinates $r$, $\\theta$ and  $\\phi$. If you want a resolution of 5 unit points in each dimension we might naively run the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:30.288966Z",
     "start_time": "2021-05-06T06:24:30.278389Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r =  [0.   0.25 0.5  0.75 1.  ]\n",
      "theta =  [0.         0.78539816 1.57079633 2.35619449 3.14159265]\n",
      "phi =  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n"
     ]
    }
   ],
   "source": [
    "res = 5  # set resolution\n",
    "\n",
    "r = np.linspace(0, 1, res)\n",
    "theta = np.linspace(0, np.pi, res)\n",
    "phi = np.linspace(-np.pi, np.pi, res)\n",
    "\n",
    "print('r = ', r)\n",
    "print('theta = ', theta)\n",
    "print('phi = ', phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how linspace produces evenly spaced numbers in defined ranges. \n",
    "\n",
    "However, when we want to calculate $f$ with the arrays `r`, `theta` and `phi`, the problem is that we need all possible combinations of the three input coordinates to fully describe a sphere. The solution is called `np.meshgrid` and works with two or more input arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:30.639025Z",
     "start_time": "2021-05-06T06:24:30.624074Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R = \n",
      "[[[0.   0.   0.   0.   0.  ]\n",
      "  [0.25 0.25 0.25 0.25 0.25]\n",
      "  [0.5  0.5  0.5  0.5  0.5 ]\n",
      "  [0.75 0.75 0.75 0.75 0.75]\n",
      "  [1.   1.   1.   1.   1.  ]]\n",
      "\n",
      " [[0.   0.   0.   0.   0.  ]\n",
      "  [0.25 0.25 0.25 0.25 0.25]\n",
      "  [0.5  0.5  0.5  0.5  0.5 ]\n",
      "  [0.75 0.75 0.75 0.75 0.75]\n",
      "  [1.   1.   1.   1.   1.  ]]\n",
      "\n",
      " [[0.   0.   0.   0.   0.  ]\n",
      "  [0.25 0.25 0.25 0.25 0.25]\n",
      "  [0.5  0.5  0.5  0.5  0.5 ]\n",
      "  [0.75 0.75 0.75 0.75 0.75]\n",
      "  [1.   1.   1.   1.   1.  ]]\n",
      "\n",
      " [[0.   0.   0.   0.   0.  ]\n",
      "  [0.25 0.25 0.25 0.25 0.25]\n",
      "  [0.5  0.5  0.5  0.5  0.5 ]\n",
      "  [0.75 0.75 0.75 0.75 0.75]\n",
      "  [1.   1.   1.   1.   1.  ]]\n",
      "\n",
      " [[0.   0.   0.   0.   0.  ]\n",
      "  [0.25 0.25 0.25 0.25 0.25]\n",
      "  [0.5  0.5  0.5  0.5  0.5 ]\n",
      "  [0.75 0.75 0.75 0.75 0.75]\n",
      "  [1.   1.   1.   1.   1.  ]]]\n",
      "\n",
      "\n",
      "T = \n",
      "[[[0.         0.         0.         0.         0.        ]\n",
      "  [0.         0.         0.         0.         0.        ]\n",
      "  [0.         0.         0.         0.         0.        ]\n",
      "  [0.         0.         0.         0.         0.        ]\n",
      "  [0.         0.         0.         0.         0.        ]]\n",
      "\n",
      " [[0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]\n",
      "  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]\n",
      "  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]\n",
      "  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]\n",
      "  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]]\n",
      "\n",
      " [[1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]\n",
      "  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]\n",
      "  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]\n",
      "  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]\n",
      "  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]]\n",
      "\n",
      " [[2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]\n",
      "  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]\n",
      "  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]\n",
      "  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]\n",
      "  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]]\n",
      "\n",
      " [[3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]\n",
      "  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]\n",
      "  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]\n",
      "  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]\n",
      "  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]]]\n",
      "\n",
      "\n",
      "P = \n",
      "[[[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]\n",
      "\n",
      " [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]\n",
      "\n",
      " [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]\n",
      "\n",
      " [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]\n",
      "\n",
      " [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]\n",
      "  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]]\n"
     ]
    }
   ],
   "source": [
    "R, T, P = np.meshgrid(r, theta, phi)\n",
    "print(\"R = \")\n",
    "print(R)\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"T = \")\n",
    "print(T)\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"P = \")\n",
    "print(P)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of these have a shape of `(5, 5, 5)`, that is $5^3$ elements in total. We can use them as an input for our function $f$ to evaluate it for the whole sphere (all possible combinations of the input coordinates). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:31.027393Z",
     "start_time": "2021-05-06T06:24:31.023406Z"
    }
   },
   "outputs": [],
   "source": [
    "def f(r, theta, phi):\n",
    "    \"\"\"Transform spherical to cartesian coordinates\"\"\" \n",
    "    \n",
    "    return (r * np.sin(theta) * np.cos(phi),\n",
    "            r * np.sin(theta) * np.sin(phi),\n",
    "            r * np.cos(theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:31.223552Z",
     "start_time": "2021-05-06T06:24:31.218099Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 5, 5, 5)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(f(R, T, P))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Insider tip: Outer operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already spoke about the elementwise array-arithmetics before. A really useful feature of NumPy are outer operations, instead of performing the operation elementwise. Before we will look at this, let's consider the functions `np.inner` and `np.outer` on their own and other kinds of products first, which we can connect to linear algebra again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:31.847445Z",
     "start_time": "2021-05-06T06:24:31.842111Z"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3])\n",
    "b = np.array([0, 4, 0])\n",
    "M = np.array([[3, 2, 1],\n",
    "              [1, 2, 3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inner product of two 1D-arrays (vectors) is the scalar product (tensor dot product)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:32.413670Z",
     "start_time": "2021-05-06T06:24:32.408783Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n"
     ]
    }
   ],
   "source": [
    "print(np.inner(a, b))\n",
    "# Ordinary inner product (sum over elementwise products)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:32.702305Z",
     "start_time": "2021-05-06T06:24:32.695670Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n",
      "8\n",
      "8\n"
     ]
    }
   ],
   "source": [
    "print(np.dot(a, b))\n",
    "print(a @ b)\n",
    "print(np.tensordot(a, b, axes=1))  # tensor dot product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{a} \\cdot \\mathbf{b}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use `np.dot` also for a general matrix (vector) products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:33.558824Z",
     "start_time": "2021-05-06T06:24:33.552655Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 14])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(M, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{M} \\cdot \\mathbf{a}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:34.722856Z",
     "start_time": "2021-05-06T06:24:34.716679Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 14])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(a, M.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For vectors in $\\mathbb{R}^3$ you can compute the cross product, too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:35.384045Z",
     "start_time": "2021-05-06T06:24:35.377610Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-12,   0,   4])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cross(a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{a} \\times \\mathbf{b}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The outer product (tensor product) of two 1D arrays (vectors) is a matrix of all possible elementwise products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:36.360702Z",
     "start_time": "2021-05-06T06:24:36.354095Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  4  0]\n",
      " [ 0  8  0]\n",
      " [ 0 12  0]]\n",
      "\n",
      "[[ 0  4  0]\n",
      " [ 0  8  0]\n",
      " [ 0 12  0]]\n"
     ]
    }
   ],
   "source": [
    "print(np.outer(a, b))\n",
    "\n",
    "print()\n",
    "print(np.tensordot(a, b, axes=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{a} \\otimes \\mathbf{b}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, NumPy also allows us to use the outer/inner mechanism not only for products but for any NumPy operation. Using a `.outer` for example on the `np.substract` operation will create an array containing the result of the operation performed on all possible pairs of input elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:37.543570Z",
     "start_time": "2021-05-06T06:24:37.530489Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [1 3 5 7]\n",
      "B =  [0.  2.  4.  8.5]\n",
      "\n",
      "\n",
      "np.subtract(A, B) = \n",
      "  [ 1.   1.   1.  -1.5]\n",
      "\n",
      "\n",
      "np.subtract.outer(A, B) = \n",
      "     [[ 1.  -1.  -3.  -7.5]\n",
      " [ 3.   1.  -1.  -5.5]\n",
      " [ 5.   3.   1.  -3.5]\n",
      " [ 7.   5.   3.  -1.5]]\n"
     ]
    }
   ],
   "source": [
    "A = np.asarray([1, 3, 5, 7])\n",
    "B = np.asarray([0, 2.0, 4.0, 8.5])\n",
    "\n",
    "print(\"A = \", A)\n",
    "print(\"B = \", B)\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"np.subtract(A, B) = \")\n",
    "print(\" \", np.subtract(A,B))\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"np.subtract.outer(A, B) = \")\n",
    "print(\"    \", np.subtract.outer(A,B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fix dimensionality related issues with np.newaxis "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we combine multidimensional arrays with each other, we need to take care that\n",
    "their shape can be matched, i.e. *broadcasting* must be possible. Let's have a look at these two arrays A and C:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:38.734990Z",
     "start_time": "2021-05-06T06:24:38.726454Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]] (2, 3)\n",
      "\n",
      "\n",
      "[1 2 3] (3,)\n"
     ]
    }
   ],
   "source": [
    "C = np.array([\n",
    "    [1, 2, 3],\n",
    "    [4, 5, 6], \n",
    "    ])\n",
    "\n",
    "A = np.array([1, 2, 3])\n",
    "\n",
    "print(C, np.shape(C))\n",
    "print(\"\\n\")\n",
    "print(A, np.shape(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can substract A from C without problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:39.628979Z",
     "start_time": "2021-05-06T06:24:39.622957Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0],\n",
       "       [3, 3, 3]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C - A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the cell below will give an error. That is because we want to perform an operation (a subtraction) with C.T and A, having shapes `(3, 2)` and `(3,)`, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:40.076421Z",
     "start_time": "2021-05-06T06:24:40.066052Z"
    }
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (3,2) (3,) ",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-45-8a6275f83609>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# C.T = C transposed\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mC\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: operands could not be broadcast together with shapes (3,2) (3,) "
     ]
    }
   ],
   "source": [
    "# C.T = C transposed\n",
    "C.T - A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:40.227409Z",
     "start_time": "2021-05-06T06:24:40.221750Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 4]\n",
      " [2 5]\n",
      " [3 6]] (3, 2)\n"
     ]
    }
   ],
   "source": [
    "print(C.T, np.shape(C.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are somehow trying to combine apples and oranges here, but NumPy provides a fix. `np.newaxis` can be used to extend an array in a dimension. In this case we would like A to be of shape `(3, 1)`, a one dimensional array in two dimensions so to say. Lets try:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:40.987863Z",
     "start_time": "2021-05-06T06:24:40.977212Z"
    }
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (3,2) (1,3) ",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-47-eb3ebb4fbf5a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mC\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnewaxis\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: operands could not be broadcast together with shapes (3,2) (1,3) "
     ]
    }
   ],
   "source": [
    "C.T - A[np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:41.244279Z",
     "start_time": "2021-05-06T06:24:41.238414Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]] (1, 3)\n"
     ]
    }
   ],
   "source": [
    "print(A[np.newaxis], np.shape(A[np.newaxis]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This failed, too because the new axis (or new dimension) is in the wrong place! We want to conserve the first dimension and append the new one thereafter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:42.668753Z",
     "start_time": "2021-05-06T06:24:42.662320Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 3],\n",
       "       [0, 3],\n",
       "       [0, 3]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.T - A[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-06T06:24:42.894534Z",
     "start_time": "2021-05-06T06:24:42.888647Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1]\n",
      " [2]\n",
      " [3]] (3, 1)\n"
     ]
    }
   ],
   "source": [
    "print(A[:, np.newaxis], np.shape(A[:, np.newaxis]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, our operation was succesfull."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:** Since `C` is a 2D array, we can think of it as a matrix. Using `C.T` is the transpose of this matrix. In higher dimensions array.T reverses the order of all dimensions. \n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}