smd_whirlwind_tour.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SMD Python Hands On"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%javascript\n",
    "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width: 95% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basics\n",
    "\n",
    "Not really the topic of today, these are good starting materials:\n",
    "\n",
    "* [PeP et al. Toolbox Workshop (German)](https://toolbox.pep-dortmund.org) \n",
    "* [Scientific Python Notebooks](https://github.com/maxnoe/scientific_python_notebooks)\n",
    "* [A Byte Of Python](https://python.swaroopch.com/)\n",
    "\n",
    "Especially for the Data Mining part, but a generally a very good resource:\n",
    "* [The Python Data Science Handbook](https://jakevdp.github.io/PythonDataScienceHandbook/)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Style Guide and Linter\n",
    "\n",
    "* For python code, there is a generally accepted style guide: pep8. Follow it!\n",
    "\n",
    "* The `ruff` tool checks for both the syntactical errors and style guide violations\n",
    "    * `mamba install ruff`\n",
    "    * `flake8 test.py`\n",
    "    \n",
    "* Autoformatter\n",
    "    * `ruff format` reformats files so they follow the styleguide → not have to think about it anymore\n",
    "    \n",
    "* All good editors support linter, autocompletion, etc., recommendation would be VS Code (or VS Codium):\n",
    "  * VS Code: https://code.visualstudio.com/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## conda environments\n",
    "\n",
    "* Different projects might need different versions of python and libraries\n",
    "* Best practice: isolated environemnt with fixed versions\n",
    "* Using conda: conda environments\n",
    "* Conda environments can be specified in text files listing needed packages and their versions \n",
    "\n",
    "\n",
    "* conda can be slow when installing or updating packages. Mamba is a much better and faster alternative:\n",
    "    * For a fresh installation (recommended!), download and install `miniforge3` here:  \n",
    "      https://github.com/conda-forge/miniforge#miniforge3\n",
    "    * For an existing conda installation, you can try:\n",
    "      ```\n",
    "      conda install -n base -c conda-forge mamba\n",
    "      ``` \n",
    "\n",
    "Create a new environment using the definition file:\n",
    "\n",
    "```\n",
    "$ mamba env create -f environment.yml\n",
    "```\n",
    "\n",
    "To use the environment:\n",
    "```\n",
    "$ conda activate smd\n",
    "```\n",
    "\n",
    "Your python / ipython / jupyter should now come from this environment:\n",
    "```\n",
    "$ which python\n",
    "/home/maxnoe/.local/conda/envs/smd/bin/python\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## numpy\n",
    "\n",
    "Short reminder of the most important numpy features for SMD\n",
    "\n",
    "* Python is a dynamic, interpreted language. \n",
    "    * Easy to use\n",
    "    * Powerful\n",
    "    * Slow (large overhead, especially for simple numerical operations)\n",
    "    \n",
    "* Numpy is a library providing an efficient array implementation and access to fast code through compiled C++/C/Cython/Fortran code  \n",
    "   ⇒ Application to large arrays as a whole or element-wise (vectorization)  \n",
    "   ⇒ Rule of thumb: no for-loops over numpy arrays\n",
    "   \n",
    "* Many features for data analysis, random numbers, linear algebra, ...\n",
    "\n",
    "* \"Naive\" python code being slow is one of the main criticisms of the language, especially for its use in science.\n",
    "\n",
    "* Despite that, it is now the most commonly used language for data analysis.\n",
    "\n",
    "* Most \"hot\" code is not actually python, but compiled code like numpy. Python is \"the glue\".\n",
    "\n",
    "* There are several ways of making python code faster (e.g. numba, cython, ...) or replace it completely (Julia, C++, ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Types\n",
    "\n",
    "Numpy supports many different data types, most important are probably these three\n",
    "\n",
    "* bool: True / False\n",
    "* int64: 64-Bit signed integer\n",
    "* float64: 64-Bit floating point number\n",
    "\n",
    "More on data types and how they work in the next lecture, \"Numerical Foundations\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# numpy arrays from python lists\n",
    "floats = np.array([1.0, 3.14, 1e3])\n",
    "ints = np.array([1, 2, 3])\n",
    "bools = np.array([True, False, True])\n",
    "\n",
    "print(floats.dtype, ints.dtype, bools.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic properties of numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "array2d = np.array([\n",
    "    [1, 2, 3],\n",
    "    [4, 5, 6]\n",
    "])\n",
    "\n",
    "def array_info(a):\n",
    "    print(f'{len(a)=}, {a.size=}, {a.ndim=}, {a.shape=}, {a.dtype=!s}')\n",
    "\n",
    "array_info(floats)\n",
    "array_info(array2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexing & Masks\n",
    "\n",
    "Numpy arrays can be indexed using (collections of) integers, slices or boolean masks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1.0, -3.5, 42, -5])\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[0], a[-1], a[1:-1], a[::2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a > 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[a > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parentheses are needed here\n",
    "# | = or\n",
    "# & = and\n",
    "# ~ = not\n",
    "\n",
    "a[(a > -10) & (a < 10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The axis keyword\n",
    "\n",
    "Important for aggregating operations (e.g. `np.sum`, `np.mean`, `np.prod`, `np.std`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.arange(12).reshape(4, 3)\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(X, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting\n",
    "\n",
    "Docs: https://numpy.org/doc/stable/user/basics.broadcasting.html\n",
    "\n",
    "Using \"broadcasting\", numpy can apply element-wise calculations to arrays of different sizes, if the shapes are *compatible*\n",
    "\n",
    "* The last, overlapping dimensions must be *compatible*\n",
    "\n",
    "* Compatible means\n",
    "\n",
    "    * both dimensions are the same\n",
    "    * or one dimension is 1\n",
    "\n",
    "\n",
    "Examples:\n",
    "* shape `(3, 2, 2)` is compatible with shape `(2, 2)`, the last dimensions are the same \n",
    "* shape `(3, 2, 2)` is compatible with shape `(1, 2)`, since the last dimensions are either the same, or one of them is 1\n",
    "* shape `(3, 2, 2)` is not compatible with shape `(3, 2)`, as last dimensions are note the same (the first ones are though) \n",
    "* shape `(3, 2, 2)` is compatible with `(3, 2, 1)` since all dimensions are equal or 1\n",
    "\n",
    "A a new dimension of size `1` can easily be inserted into any array using `np.newaxis`, making shapes compatible, see below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(12).reshape(4, 3)\n",
    "b = 5\n",
    "c = np.arange(3)\n",
    "d = np.arange(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape, c.shape, d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a - b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a - c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this error is expected\n",
    "a - d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can add a new dimension of size 1 at the and to make the shapes compatible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.shape, d.shape, d[..., np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a - d[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance example: Loops vs. Numpy vs. Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Task: Find closest point to a reference point in a list of points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point = [0, 1]\n",
    "\n",
    "points = [\n",
    "    (0, 0),\n",
    "    (0.5, -0.5),\n",
    "    (1, -1),\n",
    "    (0, 2),\n",
    "    (0, 1.1),\n",
    "    (-2, 3),\n",
    "    (5, 1),\n",
    "    (10, 4),\n",
    "    (-4, 2),\n",
    "    (-3, 0),\n",
    "] * 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pure python using a for-loop over a list of tuples and the math module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "\n",
    "def distance(p1, p2):\n",
    "    return math.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)\n",
    "\n",
    "\n",
    "def find_closest(points, point):\n",
    "    min_distance = math.inf\n",
    "    min_idx = None\n",
    "    \n",
    "    for i, other in enumerate(points):\n",
    "        d = distance(point, other) \n",
    "        \n",
    "        if d < min_distance:\n",
    "            min_distance = d\n",
    "            min_idx = i\n",
    "    \n",
    "    return min_idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = find_closest(points, point)\n",
    "print(idx, points[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -o \n",
    "find_closest(points, point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeit_python = _"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using numpy arrays and numpy methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = np.array(points)\n",
    "point = np.array(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_closest_numpy(points, point):\n",
    "    distances = np.linalg.norm(points - point, axis=1)\n",
    "    idx = np.argmin(distances)\n",
    "    return idx\n",
    "\n",
    "\n",
    "idx = find_closest_numpy(points, point)\n",
    "print(idx, points[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -o\n",
    "find_closest_numpy(points, point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeit_numpy = _ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now using numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import njit\n",
    "\n",
    "# a decorator is a function that changes a function\n",
    "# we apply the njit decorator to our python functions:\n",
    "distance = njit(distance)\n",
    "find_closest = njit(find_closest)\n",
    "\n",
    "\n",
    "\n",
    "# normally used like this:\n",
    "@njit\n",
    "def foo():\n",
    "    return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = find_closest(points, point)\n",
    "\n",
    "print(idx, points[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -o\n",
    "find_closest(points, point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeit_numba = _"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = timeit_python.average / timeit_numpy.average\n",
    "print(f'Numpy is {factor:.1f} times faster than python')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = timeit_numpy.average / timeit_numba.average\n",
    "print(f'Numba is {factor:.1f} times faster than numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = timeit_python.average / timeit_numba.average\n",
    "print(f'Numba is {factor:.1f} times faster than python')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What kind of sorcery is this?\n",
    "\n",
    "* Python is a very *dynamic*, interpreted language. That makes it very powerful, but also has a large overhead, especially for simple numerical calculations\n",
    "\n",
    "* Numpy is a library providing generic methods in compiled, optimized Cython, C, Fortran, ... which has much less overhead. But since every operation works on the whole array, intermediate arrays are created for the intermediate results\n",
    "\n",
    "* Numba takes the python code and \"just-in-time (Jit)\" compiles it to machine code, completely eliminating the need for intermediate arrays and thus memory allocations in this case.\n",
    "\n",
    "\n",
    "The speed-ups can be immense and numba compiled python code can be as-fast or even faster than hand-optimized C / C++ code and has trivial interoperability with python and numpy.\n",
    "\n",
    "\n",
    "For SMD, numpy will almost always be \"good enough\", numba won't be needed to make exercises run in less than a couple of seconds to minutes.\n",
    "\n",
    "But: numba is very well suited in general and especially for algorithms which cannot be vectorized, e.g. where the results depend on the previous iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pseudo Random Numbers\n",
    "\n",
    "Docs: https://numpy.org/doc/stable/reference/random/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import default_rng\n",
    "\n",
    "# create a new random generator with a fixed seed\n",
    "# this avoids \"evil\" global state when using `np.random.seed` or `np.random.<function>`\n",
    "rng = default_rng(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1-D Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uniform = rng.uniform(-5, 5, 1000)\n",
    "normal = rng.normal(0, 1, 1000)\n",
    "poisson = rng.poisson(3, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### N-D Normal Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = [2, 1]\n",
    "cov = [[2, 1],\n",
    "       [1, 4]]\n",
    "\n",
    "normal_2d = rng.multivariate_normal(mean, cov, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Important for reproducibility: setting the seed.\n",
    "\n",
    "Also important for parallel calculations or for resuming simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng.normal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = default_rng(0)\n",
    "rng.normal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matplotlib\n",
    "\n",
    "Short intro into the two most important types of plots for SMD:\n",
    "\n",
    "* Histograms\n",
    "* Scatter plots\n",
    "\n",
    "### Histograms\n",
    "\n",
    "Histograms count the number of samples in specific intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# for interactive plots in the notebook (using the ipympl package).\n",
    "# in ipython just use %matplotlib\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# relatively new feature of matplotlib, use this instead of `fig.tight_layout()`\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "ax.hist(normal, bins=20, range=[-5, 5])\n",
    "\n",
    "None # hides matplotlib objects in output, just to not mess up the notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When comparing different samples, it's important to use the same \"binning\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = 20\n",
    "limits = [-5, 5]\n",
    "\n",
    "combined = np.concatenate([normal, uniform])\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "# define common options in a single location\n",
    "hist_options = dict(\n",
    "    bins=bins,\n",
    "    range=limits,\n",
    "    histtype='step',\n",
    "    lw=2,\n",
    ")\n",
    "\n",
    "ax.hist(normal, label='Normal', **hist_options)\n",
    "ax.hist(uniform, label='Uniform', **hist_options)\n",
    "ax.hist(combined, label='Combined', linestyle=':',color='k', **hist_options)\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "None  # just to not mess up the notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color: crimson; font-weight: bold; font-size: 2rem\">\n",
    "    For discrete values (integers), always use integral-width bins centered around the values\n",
    "</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.hist(poisson, bins=15)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "ax.hist(\n",
    "    poisson,\n",
    "    bins=np.arange(15) - 0.5,  # bins can be either number of bins or bin edges\n",
    "    edgecolor='w',\n",
    "    lw=2\n",
    ")\n",
    "\n",
    "\n",
    "np.arange(15) - 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert bin edges to bin centers and bin widths\n",
    "bins = np.arange(8) - 0.5\n",
    "\n",
    "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n",
    "bin_widths = np.diff(bins)\n",
    "\n",
    "print(bins)\n",
    "print(bin_centers)\n",
    "print(bin_widths)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "ax.scatter(normal_2d[:, 0], normal_2d[:, 1])\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normal_2d = rng.multivariate_normal(mean, cov, 10000)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "# smaller dots, no border and some transparency for this many points\n",
    "ax.scatter(\n",
    "    normal_2d[:, 0],\n",
    "    normal_2d[:, 1],\n",
    "    s=5,\n",
    "    alpha=0.2,\n",
    "    linewidth=0,\n",
    ")\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more advanced example:\n",
    "\n",
    "* Points can be colored using a third array\n",
    "* Using ListedColorMap, you can get a discrete colormap for each class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "\n",
    "iris = load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "n_classes = len(iris.target_names)\n",
    "\n",
    "# new ColorMap with `n_classes` discrete colors \n",
    "# using the standard color rotation C0, C1, ...\n",
    "cmap = ListedColormap([f'C{i}' for i in range(n_classes)], name='iris')\n",
    "cmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# automatically adjust spacing\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "# scatter plot with colors per class\n",
    "scat = ax.scatter(\n",
    "    iris.data[:, 0],       # x-values, first column\n",
    "    iris.data[:, 1],       # y-values, second column \n",
    "    c=iris.target,         # data to determine color\n",
    "    cmap=cmap,             # colormap, converts data in c to an actual color\n",
    "    vmin=-0.5,             # minimum value for the color axis\n",
    "    vmax=n_classes - 0.5,  # maximum value for the color axis\n",
    ")\n",
    "\n",
    "\n",
    "# follow SI conventions (divide by unit)\n",
    "ax.set_xlabel(iris.feature_names[0].replace('(cm)', '/ cm'))\n",
    "ax.set_ylabel(iris.feature_names[1].replace('(cm)', '/ cm'))\n",
    "\n",
    "# colorbar with ticklabels\n",
    "bar = fig.colorbar(scat, ticks=[0, 1, 2], ax=ax)\n",
    "bar.set_ticklabels(iris.target_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D-Histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "hist, xedges, yedges, plot = ax.hist2d(\n",
    "    normal_2d[:, 0],\n",
    "    normal_2d[:, 1],\n",
    "    bins=[21, 21],\n",
    "    range=[[-3, 7], [-7, 9]],\n",
    "    cmap='inferno', # has more contrast than the default viridis\n",
    ")\n",
    "fig.colorbar(plot, ax=ax)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes helpful: Logarithmic scale for the color map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# two different distributions, one significantly smaller\n",
    "normal_2d = rng.multivariate_normal(mean, cov, 100000)\n",
    "normal_2d_2 = rng.multivariate_normal([-1.5, 4.5], [[0.5, 0], [0, 0.5]], 500)\n",
    "\n",
    "normal_2d_both = np.concatenate([normal_2d, normal_2d_2], axis=0)\n",
    "normal_2d_both.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "*_, plot = ax.hist2d(\n",
    "    normal_2d_both[:, 0],\n",
    "    normal_2d_both[:, 1],\n",
    "    bins=[50, 50],\n",
    "    range=[[-3, 7], [-7, 9]],\n",
    "    cmap='inferno'\n",
    ")\n",
    "fig.colorbar(plot, ax=ax)\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "*_, plot = ax.hist2d(\n",
    "    normal_2d_both[:, 0],\n",
    "    normal_2d_both[:, 1],\n",
    "    bins=[50, 50],\n",
    "    range=[[-3, 7], [-7, 9]],\n",
    "    norm=LogNorm(),\n",
    "    cmap='inferno'\n",
    ")\n",
    "fig.colorbar(plot, ax=ax)\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also often important: fixed aspect ratio of the x and y axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "hist, xbins, ybins, plot = ax.hist2d(\n",
    "    normal_2d_both[:, 0],\n",
    "    normal_2d_both[:, 1],\n",
    "    bins=[50, 50],\n",
    "    range=[[-3, 7], [-7, 9]],\n",
    "    norm=LogNorm(),\n",
    "    cmap='inferno',\n",
    ")\n",
    "\n",
    "fig.colorbar(plot, ax=ax)\n",
    "\n",
    "# force x and y to have same aspect ratio, also takes numbers\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy\n",
    "\n",
    "For SMD-A, `scipy.stats` will be most important\n",
    "\n",
    "### scipy.stats\n",
    "\n",
    "Many statistical distributions with many properties\n",
    "\n",
    "Docs: https://docs.scipy.org/doc/scipy/reference/stats.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "mean = 5\n",
    "std = 2\n",
    "\n",
    "normal_distribution = norm(mean, std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# draw random samples using our generator\n",
    "samples = normal_distribution.rvs(size=1000, random_state=rng)\n",
    "\n",
    "x = np.linspace(mean - 5 * std, mean + 5 * std, 250)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "\n",
    "# plot histogram\n",
    "ax.hist(\n",
    "    samples,\n",
    "    bins=100,\n",
    "    range=[x.min(), x.max()],\n",
    "    density=True,  # normalize histgram to an area of 1\n",
    "    label='Normalized Histogram',\n",
    ")\n",
    "\n",
    "# plot pdf and cdf\n",
    "ax.plot(x, normal_distribution.pdf(x), label='Probability Density Function', lw=2)\n",
    "ax.plot(x, normal_distribution.cdf(x), label='Cumulative Distribution Function', lw=2)\n",
    "\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rng.normal(5, 2, 100)\n",
    "\n",
    "# maximum likelihood fit \n",
    "mean_fit, std_fit = norm.fit(x)\n",
    "\n",
    "mean_fit, mean, std_fit, std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pandas\n",
    "\n",
    "Library for tabular data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a pandas.DataFrame from numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = pd.DataFrame({\n",
    "    'x': rng.normal(2, 0.5, 1000),\n",
    "    'y': rng.uniform(0.5, 1, 1000),\n",
    "    'N': rng.poisson(50, 1000),\n",
    "    't': rng.exponential(5, 1000),\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First look at the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input/Output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal.to_csv('data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = pd.read_csv('data.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# HDF5 is a fast, binary data format better suited for large datasets\n",
    "signal.to_hdf('data.hdf5', key='signal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "background = pd.DataFrame({\n",
    "    'x': rng.uniform(-4, 4, 10000),\n",
    "    'y': rng.uniform(-4, 4, 10000),\n",
    "    'N': rng.poisson(30, 10000),\n",
    "    't': rng.exponential(10, 10000),\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# you can store more than one dataset in the same file\n",
    "background.to_hdf('data.hdf5', key='background')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "background.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_hdf('data.hdf5', key='signal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the first / last values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The titanic dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('titanic.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many valid values in each column?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.count().sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop columns with too many missing values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# axis=1: drop columns\n",
    "# inplace=True: directly change df\n",
    "df.drop(['cabin', 'boat', 'body', 'home.dest'], axis=\"columns\", inplace=True) \n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many men/women on the titanic?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['sex'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Very powerful operation: GroupBy → Aggregate\n",
    "\n",
    "The dataset is split into multiple groups, aggregated values calculated per group.\n",
    "\n",
    "\n",
    "Example: mean survival rate by sex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.groupby('sex')['survived'].agg('mean')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`DataFrame` also supports mask indexing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['child'] = df['age'] < 9\n",
    "\n",
    "df[df['child']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.groupby('child').survived.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classes\n",
    "\n",
    "The overaching task this year is a simplified detector simulation and analysis.\n",
    "\n",
    "Much of the structure will be predefined, with exercises \"filling in blanks\".\n",
    "\n",
    "\n",
    "To understand the structure, some understanding of python classes and tests is helpful\n",
    "\n",
    "Docs: https://docs.python.org/3/reference/datamodel.html\n",
    "\n",
    "This is a **very** short primer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Detector:\n",
    "    \n",
    "    # special method to initialize a new instance of a class\n",
    "    def __init__(self, width, height):\n",
    "        self.width = width\n",
    "        self.height = height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detector = Detector(width=10, height=5)\n",
    "\n",
    "print(detector.width, detector.height)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One possible exercise could be to implement a method checking if a particle hits the detector.\n",
    "\n",
    "Let the coordinate system be defined as follows:\n",
    "* Lower left corner of the detector is at (0, 0)\n",
    "* width along x\n",
    "* height along y\n",
    "\n",
    "\n",
    "The given structure could then be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Detector:\n",
    "    def __init__(self, width, height):\n",
    "        self.width = width\n",
    "        self.height = height\n",
    "        \n",
    "    def is_inside(self, x, y):\n",
    "        # this is to make the function work without error,\n",
    "        # but it does not give the correct answer yet\n",
    "        inside = False\n",
    "        ###### Implement your solution between here #######\n",
    "        \n",
    "        ###### and here #######\n",
    "        return inside"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will provide tests for most of these, where you can check if your implementation\n",
    "is correct. For this, we use pytest:\n",
    "\n",
    "* Docs: https://docs.pytest.org/en/stable/\n",
    "* Lecture on testing with pytest from the Escape School 2021: https://www.youtube.com/watch?v=pGg97d8TQuY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_is_inside():\n",
    "    # a test function to be used with `pytest`\n",
    "    d = Detector(5, 5)\n",
    "\n",
    "    assert d.is_inside(2, 2), 'This point should be inside'\n",
    "    assert not d.is_inside(10, 10), 'This point should not be inside'\n",
    "    print(\"test passed\")\n",
    "    \n",
    "test_is_inside()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Python, we can implement \"special\" operations on our class using *magic* methods,\n",
    "which are defined in the data model and for example make operator overloading possible.\n",
    "\n",
    "These methods always start and end with `__` (double-underscore, \"dunder\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Detector(10, 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Detector:\n",
    "    def __init__(self, width, height):\n",
    "        self.width = width\n",
    "        self.height = height\n",
    "\n",
    "    def __repr__(self):\n",
    "        '''This methods enables conversion to a representation str, e.g. for printing'''\n",
    "        return f'{self.__class__.__name__}(width={self.width}, height={self.height})'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = Detector(10, 5)\n",
    "print(d)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}