numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "# <img width=400 src=\"http://www.numpy.org/_static/numpy_logo.png\" alt=\"Numpy\"/>\n",
    "\n",
    "\n",
    "## Why do we need numpy?\n",
    "\n",
    "* You may have heard \"Python is slow\", this is true when it concerns looping over many small python objects\n",
    "* Python is dynamically typed and everything is an object, even an `int`. There are no primitive types.\n",
    "* Numpy's main feature is the `ndarray` class, a fixed length, homogeniously typed array class.\n",
    "* Numpy implements a lot of functionality in fast c, cython and fortran code to work on these arrays\n",
    "* python with vectorized operations using numpy can be blazingly fast\n",
    "\n",
    "See: [Python is not C](https://www.ibm.com/developerworks/community/blogs/jfp/entry/Python_Is_Not_C?lang=en)\n",
    "\n",
    "But the most important reason:\n",
    "\n",
    "* More beautiful code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More beautiful code through vectorisation\n",
    "\n",
    "pure python with list comprehension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voltages = [10.1, 15.1, 9.5]\n",
    "currents = [1.2, 2.4, 5.2]\n",
    "\n",
    "resistances = [U / I for U, I in zip(voltages, currents)]\n",
    "resistances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = np.array([10.1, 15.1, 9.5])\n",
    "I = np.array([1.2, 2.4, 5.2])\n",
    "\n",
    "R = U / I\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the point with the smallest distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def euclidean_distance(p1, p2):\n",
    "    return math.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)\n",
    "\n",
    "point = (1, 2)\n",
    "points = [(3, 2), (4, 2), (3, 0)]\n",
    "\n",
    "min_distance = float('inf')\n",
    "for other in points:\n",
    "    distance = euclidean_distance(point, other)\n",
    "    if distance < min_distance:\n",
    "        closest = other\n",
    "        min_distance = distance \n",
    "\n",
    "print(min_distance, closest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point = np.array([1, 2])\n",
    "points = np.array([(3, 2), (4, 2), (3, 0)])\n",
    "\n",
    "distance = np.linalg.norm(point - points, axis=1)\n",
    "idx = np.argmin(distance)\n",
    "\n",
    "print(distance[idx], points[idx])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small example timings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "\n",
    "def var(data):\n",
    "    '''\n",
    "    knuth's algorithm for one-pass calculation of the variance\n",
    "    Avoids rounding errors of large numbers when doing the naive\n",
    "    approach of `sum(v**2 for v in data) - sum(v for v in data)**2`\n",
    "    '''\n",
    "    \n",
    "    n = 0\n",
    "    mean = 0.0\n",
    "    m2 = 0.0\n",
    "    \n",
    "    if len(data) < 2:\n",
    "        return float('nan')\n",
    "\n",
    "    for value in data:\n",
    "        n += 1\n",
    "        delta = value - mean\n",
    "        mean += delta / n\n",
    "        delta2 = value - mean\n",
    "        m2 += delta * delta2\n",
    "\n",
    "    return m2 / n "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(range(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = list(range(1000))\n",
    "a = np.array(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit var(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit np.var(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic math: vectorized\n",
    "\n",
    "Operations on numpy arrays work vectorized, element-by-element\n",
    "\n",
    "**Lose your loops**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a numpy array from a python a python list\n",
    "a = np.array([1.0, 3.5, 7.1, 4, 6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2 * a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a**a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.cos(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Attention: You need the `cos` from numpy!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "math.cos(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most normal python functions with basic operators like `*`, `+`, `**` simply work because\n",
    "of operator overloading:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def poly(x):\n",
    "    return x + 2 * x**2 - x**3\n",
    "\n",
    "poly(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly(np.e), poly(np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Useful properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.ndim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arbitrary dimension arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# two-dimensional array\n",
    "y = np.array([[1, 2, 3],\n",
    "              [4, 5, 6],\n",
    "              [7, 8, 9]])\n",
    "\n",
    "y + y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## since python 3.5 @ is matrix product\n",
    "y @ y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Broadcasting, changing array dimensions to fit the larger one\n",
    "\n",
    "y + np.array([1, 2, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trick for column broadcasting\n",
    "\n",
    "y = np.ones((3, 5))\n",
    "\n",
    "(y.T + np.arange(3)).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create new axis to enable broadcasting for unmatched shapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(3)\n",
    "b = np.arange(4, 7)\n",
    "\n",
    "a[:, np.newaxis] * b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reduction operations\n",
    "\n",
    "Numpy has many operations, which reduce dimensionality of arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.normal(0, 1, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.prod(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standard Deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample Standard Deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(x, ddof=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the numpy functions are also methods of the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.mean(), x.std(), x.max(), x.min()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Difference between neighbor elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.arange(10)**2\n",
    "diff_z = np.diff(z)\n",
    "\n",
    "print(z)\n",
    "print(diff_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reductions on multi-dimensional arrays\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "array2d = np.arange(20).reshape(4, 5)\n",
    "\n",
    "array2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(array2d, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.var(array2d, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "Write a function that calculates the analytical linear regression for a set of\n",
    "x and y values.\n",
    "\n",
    "Reminder:\n",
    "\n",
    "$$ f(x) = a \\cdot x + b$$\n",
    "\n",
    "with \n",
    "\n",
    "$$\n",
    "\\hat{a} = \\frac{\\mathrm{Cov}(x, y)}{\\mathrm{Var}(x)} \\\\\n",
    "\\hat{b} = \\bar{y} - \\hat{a} \\cdot \\bar{x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_regression(x, y):\n",
    "    a = np.nan\n",
    "    b = np.nan\n",
    "    return a, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy_linear.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 1, 50)\n",
    "y = 5 * np.random.normal(x, 0.1) + 2  # see section on random numbers later\n",
    "\n",
    "a, b = linear_regression(x, y)\n",
    "a, b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helpers for creating arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ones((5, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.full(5, np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.empty(10)  # attention, uninitialised memory, be carefull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linspace(-2, 1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# like range() for arrays:\n",
    "np.arange(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arange(2, 10, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.logspace(-4, 5, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.logspace(1, 4, 4, base=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numpy Indexing\n",
    "\n",
    "* Element access\n",
    "* Slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![index1d](images/Indexing1D.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0, 10)\n",
    "\n",
    "# like lists:\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# like lists:\n",
    "x[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all elements with indices ≥1 and <4:\n",
    "x[1:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# negative indices count from the end\n",
    "x[-1], x[-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combination:\n",
    "x[3:-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step size\n",
    "x[::2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trick for reversal: negative step\n",
    "x[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![index2d](images/Indexing2D.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([x, x + 10, x + 20, x + 30])\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only one index ⇒ one-dimensional array\n",
    "y[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# other axis: (: alone means the whole axis)\n",
    "y[:, 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inspecting the number of elements per axis:\n",
    "y[:, 1:3].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Changing array content"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[:, 3] = 0\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using slices on both sides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[:,0] = x[3:7]\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transposing inverts the order of the dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.T.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Masks\n",
    "\n",
    "* A boolean array can be used to select only the element where it contains `True`.\n",
    "* Very powerfull tool to select certain elements that fullfill a certain condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.linspace(0, 2, 11)\n",
    "b = np.random.normal(0, 1, 11)\n",
    "\n",
    "print(b >= 0)\n",
    "print(a[b >= 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[[0, 2]] = np.nan\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[np.isnan(a)] = -1\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random numbers\n",
    "\n",
    "* numpy has a larger number of distributions builtin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.uniform(-1, 1, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.normal(0, 5, (2, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.poisson(5, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating pi through monte-carlo simulation\n",
    "\n",
    "* We draw random numbers in a square with length of the sides of 2\n",
    "* We count the points which are inside the circle of radius 1\n",
    "\n",
    "The area of the square is\n",
    "\n",
    "$$\n",
    "A_\\mathrm{square} = a^2 = 4\n",
    "$$\n",
    "\n",
    "The area of the circle is\n",
    "$$\n",
    "A_\\mathrm{circle} = \\pi r^2 = \\pi\n",
    "$$\n",
    "\n",
    "With \n",
    "$$\n",
    "\\frac{n_\\mathrm{circle}}{n_\\mathrm{square}} = \\frac{A_\\mathrm{circle}}{A_\\mathrm{square}}\n",
    "$$\n",
    "We can calculate pi:\n",
    "\n",
    "$$\n",
    "\\pi = 4 \\frac{n_\\mathrm{circle}}{n_\\mathrm{square}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_square = 100000\n",
    "\n",
    "x = np.random.uniform(-1, 1, n_square)\n",
    "y = np.random.uniform(-1, 1, n_square)\n",
    "\n",
    "radius = np.sqrt(x**2 + y**2)\n",
    "\n",
    "n_circle = np.sum(radius <= 1.0)\n",
    "\n",
    "print(4 * n_circle / n_square)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "1. Draw 10000 gaussian random numbers with mean of $\\mu = 2$ and standard deviation of $\\sigma = 3$\n",
    "2. Calculate the mean and the standard deviation of the sample\n",
    "3. What percentage of the numbers are outside of $[\\mu - \\sigma, \\mu + \\sigma]$?\n",
    "4. How many of the numbers are $> 0$?\n",
    "5. Calculate the mean and the standard deviation of all numbers ${} > 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy_gaussian.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Monte-Carlo uncertainty propagation\n",
    "\n",
    "* The hubble constant as measured by PLANCK is\n",
    "$$\n",
    "H_0 = (67.74 \\pm 0.47)\\,\\frac{\\mathrm{km}}{\\mathrm{s}\\cdot\\mathrm{Mpc}}\n",
    "$$\n",
    "\n",
    "* Estimate mean and the uncertainty of the velocity of a galaxy which is measured to be $(500 \\pm 100)\\,\\mathrm{Mpc}$ away\n",
    "using monte carlo methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/numpy_hubble.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple io functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.arange(100)\n",
    "x = np.random.normal(0, 1e5, 100)\n",
    "y = np.random.normal(0, 1, 100)\n",
    "n = np.random.poisson(20, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx.shape, x.shape, y.shape, n.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\n",
    "    'data.txt',\n",
    "    np.column_stack([idx, x, y, n]),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head data.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load back the data, unpack=True is needed to read the data columnwise and not row-wise\n",
    "idx, x, y, n = np.genfromtxt('data.txt', unpack=True)\n",
    "\n",
    "idx.dtype, x.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problems\n",
    "\n",
    "* Everything is a float\n",
    "* Way larger file than necessary because of too much digits for floats\n",
    "* No column names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numpy recarrays\n",
    "\n",
    "* Numpy recarrays can store columns of different types\n",
    "* Rows are addressed by integer index\n",
    "* Columns are addressed by strings\n",
    "\n",
    "Solution for our io problem → Column names, different types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for more options on formatting see\n",
    "# https://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html\n",
    "\n",
    "data = np.savetxt(\n",
    "    'data.csv',\n",
    "    np.column_stack([idx, x, y, n]),\n",
    "    delimiter=',', # true csv file\n",
    "    header=','.join(['idx', 'x', 'y', 'n']),\n",
    "    fmt=['%d', '%.4g', '%.4g', '%d'],  # One formatter for each column\n",
    "    comments='',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head data.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.genfromtxt(\n",
    "    'data.csv',\n",
    "    names=True, # load column names from first row\n",
    "    dtype=None, # Automagically determince best data type for each column\n",
    "    delimiter=',',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['n']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear algebra\n",
    "\n",
    "Numpy offers a lot of linear algebra functionality, mostly wrapping LAPACK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# symmetric matrix, use eigh\n",
    "# If not symmetric, use eig\n",
    "mat = np.array([\n",
    "    [4, 2, 0],\n",
    "    [2, 1, -3],\n",
    "    [0, -3, 4]\n",
    "])\n",
    "\n",
    "eig_vals, eig_vecs = np.linalg.eigh(mat)\n",
    "\n",
    "eig_vals, eig_vecs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.inv(mat)"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}