day_5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Programming for Scientists - Day 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "In the fifth (and last!) day we will consider performance, i.e. \"how to make my code faster\"!\n",
    "\n",
    "* `multiprocessing`\n",
    "* `threading`\n",
    "* `numba` - \"just in time\" compilation of python code\n",
    "* `pybind11` - using external c/c++/fortran code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [0] Parallel Programming\n",
    "\n",
    "One way to make a given task execute faster is to use **parallel programming**, i.e. having more than one worker working on the task at once, reducing the overall time to completion.\n",
    "\n",
    "There are two important concepts: **processes vs threads**:\n",
    "* If you start a program on a computer, its running instance could be called a process.\n",
    "* A given process can start one, or many, threads.\n",
    "* Both processes and threads can be thought of as workers which execute a series of commands (e.g. a block of python code).\n",
    "\n",
    "The main difference is:\n",
    "* Threads (of the same process) run in a shared memory space, i.e. they can easily share variables and memory. Threads are lightweight (created very fast, low memory overhead).\n",
    "* Processes run in separate memory spaces: they cannot share variables or memory, unless they explicitly \"communicate\" (i.e. send and receive) information between themselves. Processes are heavy (slow to start, memory overhead).\n",
    "\n",
    "You can either use multiple processes, or multiple threads (or both!) to make a program parallel.\n",
    "\n",
    "By definition all threads (of a process) run together on a single computer. On the other hand, multiple processes can  all run on the same system, or they can be spread across multiple computers.\n",
    "\n",
    "This leads to the three main models for parallel computing:\n",
    "\n",
    "* if multiple threads or processes run on the same computer, they can use **shared-memory parallelism**.\n",
    "* if multiple processes are spread across different computers, this is **distributed-memory parallelism** or **message passing**, in which case communication must occur across a network.\n",
    "* if the problem is **embarrassingly parallel** it can be trivially split into a number of sub-problems, each being fully independent, and the final result can be assembled after.\n",
    "\n",
    "We will mostly look at the first case (one computer).\n",
    "\n",
    "If you want to write a program that needs distributed-memory parallelism (i.e. high-performance computing, analysis of datasets >100s GB which cannot be otherwise split), the low-level package to learn is [mpi4py](https://mpi4py.readthedocs.io/en/stable/). A high-level, very modern option is [dask](https://dask.org/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [1] Multiprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A straightforward way to speed up a computation in Python is to have multiple Python processes working together. In this case we use the built-in `multiprocessing` library. This approach is good if:\n",
    "* there is lots of computation to do, and not much data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing as mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine we have a function which is very expensive to run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to run it for a list of inputs values, and obtain the outputs, we can just do a loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in [1,2,3]:\n",
    "    print(f(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three calls to `f()` run **in serial** (one after another). We can instead run them **in parallel**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with mp.Pool() as p:\n",
    "    args = [1,2,3]\n",
    "    print(p.map(f, args))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: We use the `with` syntax, just like opening a file, to let Python automatically clean up all the internal aspects of `Pool` when we are done with it.\n",
    "\n",
    "The `Pool.map` function is a helper function which:\n",
    "* (i) starts a number of independent \"child\" Python processes\n",
    "* (ii) distributes (\"maps\") the set of arguments between these processes\n",
    "* (iii) runs the function `f` with the argument(s) which each process is responsible for\n",
    "* (iv) collects the results by sending them back to the \"parent\" process\n",
    "* (v) shuts down all the child processes\n",
    "\n",
    "Note: Initializing `mp.Pool()` will use all the available CPU cores on the current machine (1 process per core). Instead, `mp.Pool(4)` will start only 4 processes, so these will occupy only 4 cores.\n",
    "\n",
    "> If you are sharing a system (non-exclusive), you will want to avoid using all the CPU cores.\n",
    ">\n",
    "> If you are on a system with only 1 CPU core, there is no point to multiprocessing.\n",
    "\n",
    "Note: the `Pool.starmap` function is similar, except that each element of `args` is extended to also be a list, which is then unpacked and passed to `f` as a list of arguments.\n",
    "\n",
    "> For example, an iterable of `args = [(1,2), (3, 4)]` results in `[func(1,2), func(3,4)]`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a convenient helper which abstracts away much of the complexity. For instance, what if we want to run `f()` on 10 different arguments, but only have 4 cores (and so 4 processes)? This is automatically handled:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with mp.Pool(4) as p:\n",
    "    args = np.arange(10)\n",
    "    print(p.map(f, args))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One important note to keep in mind: data passed between the processes (the arguments, and the return) is pickled (\"serialized\", using the same `pickle` library we saw earlier to save/load an arbitrary Python object). This is **only ok for small data sizes**, but extremely slow/problematic for large data (e.g. > GB)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is it actually faster? Let's make a more expensive function, then benchmark it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "args = [1,2,3]\n",
    "\n",
    "def f(x):\n",
    "    j = 0\n",
    "    for i in range(10000000):\n",
    "        j += x*x\n",
    "    return j"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Profiling timing (benchmarking)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "for i in args:\n",
    "    print(f(i))\n",
    "    \n",
    "print(f'Took {time.time() - start_time:.2f} seconds.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can avoid having to do the timing ourselves with a useful helper (in IPython or a notebook):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "for i in args:\n",
    "    print(f(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare that serial result, to the parallel result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with mp.Pool() as p:\n",
    "    print(p.map(f, args))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ideally, the second case finished in a shorter amount of time. If the \"scaling\" was perfect, then running `N` processes on `N` cores, the total time to finish the computation should be `N` times lower. In practice, it will never be quite as good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: The `%%time` times the whole cell (and must be the first line!). If you have a single-line command, you can do e.g. `%time func()`, with just one `%`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Write a loop to benchmark a function like `f()` above, testing multiprocessing pools with different numbers of processes (e.g. `1, 2, 4, 8`). Plot the time to solution versus number of processes. Plot the \"speedup\" (relative to `N=1`) versus number of processes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2] Threading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Python standard library provides `threading`, which provides a clean interface to work with threads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import threading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the same function as above, as an example of something which takes a bit of time to run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "args = [1,2,3]\n",
    "\n",
    "def f_save(x, results, results_index):\n",
    "    j = 0\n",
    "    for i in range(10000000):\n",
    "        j += x*x\n",
    "        \n",
    "    results[results_index] = j"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of returning the answer, we pass our work function a `results` array, as well as a `results_index` indicating where in that array it should store its result.\n",
    "\n",
    "Our plan is that: each thread will compute one answer as a number, and all threads will store their results in the same `results` array (in different places)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make threads, start them manually, and collect the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "threads = []\n",
    "\n",
    "# allocate an array to store the results\n",
    "results = np.zeros(len(args), dtype='int64')\n",
    "\n",
    "# create a list of Thread objects, each set to run f(arg) when it starts, one thread per element of args\n",
    "for i, arg in enumerate(args):\n",
    "    t = threading.Thread(target=f_save, args=(arg,results,i))\n",
    "    threads.append(t)\n",
    "    \n",
    "    # spawn this thread, tell it to start, and \"detach\" it so it can run independently (likely on a different CPU core automatically)\n",
    "    t.start()\n",
    "    \n",
    "# now the threads are all working\n",
    "\n",
    "# wait for each to finish (re\"join\" the master thread, i.e. this one)\n",
    "for thread in threads:\n",
    "    thread.join()\n",
    "    \n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subclassing Thread to store a return value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative to preparing a separate `results` array (or list) ahead of time, and then filling a specific index with the result from each thread, we can instead create a custom sub-class of `Thread` (using inheritance), which uses an attribute (class variable) to store the return of `f()`, which we can then access after all threads have completed their work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyThread(threading.Thread):\n",
    "    def __init__(self, target, *args, **kwargs):\n",
    "        # call the __init__() method of the parent class\n",
    "        super(MyThread, self).__init__()\n",
    "        \n",
    "        # save the target function, any arguments\n",
    "        self.target = target\n",
    "        self.args = args\n",
    "        self.kwargs = kwargs\n",
    "        \n",
    "    # when Thread.start() is called it in turn calls run(), which we override, saving the return of target() into an attribute variable\n",
    "    def run(self):\n",
    "        self.result = self.target(*self.args, **self.kwargs)\n",
    "        \n",
    "    # when join() is called, call the join() method of our parent Thread class, then return our result\n",
    "    def join(self):\n",
    "        super(MyThread, self).join()\n",
    "        return self.result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `MyThread` instead of `threading.Thread`, with its added functionality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threads = []\n",
    "\n",
    "# create a list of MyThread objects, each set to run f(arg) when it starts, one thread per element of args\n",
    "for arg in args:\n",
    "    t = MyThread(f, arg)\n",
    "    threads.append(t)\n",
    "    \n",
    "    # spawn this thread, tell it to start, and \"detach\" it so it can run independently (likely on a different CPU core automatically)\n",
    "    t.start()\n",
    "    \n",
    "# wait for each to finish (re\"join\" the master thread, i.e. this one), and capture the return of .join() into a list\n",
    "results = [t.join() for t in threads]\n",
    "\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Race conditions\n",
    "\n",
    "Because **threads share memory**, it is very easy to get into trouble.\n",
    "\n",
    "You should **never write the same piece of memory from more than one thread**.\n",
    "\n",
    "You should also **never read a piece of memory from one thread which has been written to by another**.\n",
    "\n",
    "* When multiple threads are running, they are effectively racing against each other.\n",
    "* You don't know which will reach a particular code block first, so you don't know if such a \"common\" memory space has been written by another thread already, or not.\n",
    "* As a result, the value of that memory space (variable) can be unexpected.\n",
    "* This leads to **race conditions** (multiple threads interfere with each other) and **deadlock** (due to such interference, a thread may end up waiting forever and never finishing).\n",
    "* These types of problems do not necessarily happen each time you run the same program, but can happen only occasionally: they are not easily reproducible, and can be **non-deterministic**. Very hard to fix. Better to avoid at all costs.\n",
    "\n",
    "Let's make a very simple example to understand the problem. We will create a race condition which happens every time, to be clear:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "class VisitorCounter():\n",
    "    def __init__(self):\n",
    "        self.count = 0\n",
    "\n",
    "    def update(self, thread_num):\n",
    "        print('Thread %d, starting update.' % thread_num)\n",
    "\n",
    "        local_value = self.count\n",
    "        local_value += 1\n",
    "        time.sleep(2.0)\n",
    "        self.count = local_value\n",
    "\n",
    "        print('Thread %d, finished update.' % thread_num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `update()` function emulates what might happen in a database, for example. `value` could be a visitor counter on a website, and each time someone visits, we want to increase the counter by one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threads = []\n",
    "counter = VisitorCounter()\n",
    "\n",
    "# create and start two threads\n",
    "for i in range(2):\n",
    "    t = threading.Thread(target=counter.update, args=(i,))\n",
    "    threads.append(t)\n",
    "    \n",
    "    t.start()\n",
    "    \n",
    "# wait for both to finish\n",
    "for thread in threads:\n",
    "    thread.join()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Exercise\n",
    "\n",
    "What is the final value of `counter.count`? What happened? Copy the code block from above and change it to run the two threads in serial, one after another. What happens to the final counter value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synchronization between threads (Lock)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the example above may seem silly and easy to avoid, such problems actually happen and can be very hard to identify.\n",
    "\n",
    "There are a number of ways to avoid or solve race conditions. In the [threading documentation](https://docs.python.org/3/library/threading.html) you will find \"synchronization primitives\" called `Lock`, `RLock` (reentrant lock), `Semaphore`, and so on. These are all parallel programming concepts on how multiple threads can coordinate with each other. As are: `Mutex`, `Events`, `Signals`, `Condition variables`.\n",
    "\n",
    "We will only look quickly at the simplest, which is a `Lock` (very similar to what a `Mutex`, or \"mutual exclusion\", is in other languages).\n",
    "\n",
    "To solve the race condition above, we need to allow only one thread at a time into the read-modify-write \"section\" of the code. A `Lock` lets us do this:\n",
    "* a lock is an object which is like a hall pass in school.\n",
    "* only one thread at a time can have the lock.\n",
    "* any other thread which wants the lock must wait until its current owner gives it up (i.e. back to the teacher).\n",
    "* the basic methods are `.acquire()` and `.release()`.\n",
    "\n",
    "> Note: if one thread gets the lock, but never gives it back, your program will never finish: a deadlock.\n",
    "\n",
    "Let's update our `VisitorCounter` class from above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class VisitorCounterWithLock():\n",
    "    def __init__(self):\n",
    "        self.count = 0\n",
    "        self._lock = threading.Lock()\n",
    "\n",
    "    def update(self, thread_num):\n",
    "        print('Thread %d starting update.' % thread_num)\n",
    "\n",
    "        with self._lock:\n",
    "            # before entering this block, the function must acquire the lock\n",
    "            print('Thread %d has lock.' % thread_num)\n",
    "            local_value = self.count\n",
    "            local_value += 1\n",
    "            time.sleep(2.0)\n",
    "            self.count = local_value\n",
    "\n",
    "        print('Thread %d finished update, releasing lock.' % thread_num)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threads = []\n",
    "counter = VisitorCounterWithLock()\n",
    "\n",
    "# create and start two threads\n",
    "for i in range(2):\n",
    "    t = threading.Thread(target=counter.update, args=(i,))\n",
    "    threads.append(t)\n",
    "    \n",
    "    t.start()\n",
    "    \n",
    "# wait for both to finish\n",
    "for thread in threads:\n",
    "    thread.join()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Python's \"Global Interpreter Lock\": the bane of multithreaded programming.\n",
    "\n",
    "Before finishing with threading, let's test our performance gain by profiling our slow function `f()` which we have sped up."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Call `f(arg)` for each `arg` in `args` (just use a for loop, like you did in the exercise above). This will run the computation in serial. Time it: what is the total time? How does this compare to our first result (where we used `%%time` above)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have achieved absolutely no speedup!\n",
    "\n",
    "The problem is that, due to the [Global Interpreter Lock](https://docs.python.org/3/glossary.html#term-global-interpreter-lock), **only one thread can execute Python code at once**.\n",
    "\n",
    "> This is quite a bummer. (Many people would like to change this in Python, but not yet).\n",
    ">\n",
    "> So why would you ever use multiple threads?\n",
    ">\n",
    "> If a thread started in Python \"lowers\" out of the Python interpreter, e.g. starts to run pure C or C++ code, then it can release the GIL and circumvent this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3] Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But, writing C/C++ code is hard. We rarely want to do this. The best solution is the [numba](https://numba.pydata.org/) library.\n",
    "\n",
    "> Numba translates Python functions to optimized machine code \"just-in-time\" (JIT), i.e. at the moment you run the code (using the \"LLVM\" compiler).\n",
    ">\n",
    "> Numba-compiled numerical algorithms in Python can approach the speeds of C or FORTRAN (or be faster).\n",
    ">\n",
    "> You don't need to change your code, replace the Python interpreter, run a compiler, or even have a C/C++ compiler installed. Just apply the Numba decorator to a Python function, and Numba does the rest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import njit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's look at a simple example, a method of calculating the value of $\\pi$ using a Monte Carlo approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day5_mc_pi.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea is to use the equation for the area of a circle: $A = 2 \\pi r^2$ relative to the area of a square $A = r^2$. The ratio between the two is $2 \\pi$. If we imagine throwing random darts at a square board, the number which fall within the circle, relative to the total, give us a \"Monte Carlo\" estimate of this number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def monte_carlo_pi(nsamples):\n",
    "    acc = 0\n",
    "    for i in range(nsamples):\n",
    "        x = random.random()\n",
    "        y = random.random()\n",
    "        if (x ** 2 + y ** 2) < 1.0:\n",
    "            acc += 1\n",
    "    return 4.0 * acc / nsamples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: you should **never write a for-loop based code in Python, when you can use numpy instead!** We show this only for demonstration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time monte_carlo_pi(10000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Re-write the `monte_carlo_pi()` function above using numpy, and no for loops. Compare the speed. (This has nothing to do with numba, it is just a check that we remember Day 2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the numba 'jit' decorator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@njit\n",
    "def monte_carlo_pi_numba(nsamples):\n",
    "    acc = 0\n",
    "    for i in range(nsamples):\n",
    "        x = random.random()\n",
    "        y = random.random()\n",
    "        if (x ** 2 + y ** 2) < 1.0:\n",
    "            acc += 1\n",
    "    return 4.0 * acc / nsamples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time monte_carlo_pi_numba(10000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first time you run the function, numba compiles it behind the scenes. This can add some time, so for benchmarking, always run it a few times.\n",
    "\n",
    "> In practice, the time for the function to run will be much, much greater than the compilation time.\n",
    "\n",
    "Although for this simple algorithm we could use numpy, many more complex algorithms and analyses cannot easily be re-written without loops. **Use numba to write algorithms with C-like loops.**\n",
    "\n",
    "> The catch: **numba JIT'ed code must be kept simple**: only numpy arrays, (most) functions, and simple operations. No file I/O! No fancy python: avoid dictionaries, classes, etc. Keep it simple..\n",
    ">\n",
    "> See [supported python features](https://numba.readthedocs.io/en/stable/reference/pysupported.html) and [supported numpy features](https://numba.readthedocs.io/en/stable/reference/numpysupported.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Exercise\n",
    "\n",
    "If `nsamples=1e9` (one billion), what is the peak memory usage of `monte_carlo_pi_numba` (don't code, just think about it)? Peak memory means the largest amount of memory used, at any point in the program. What about `monte_carlo_numpy`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Profiling memory usage\n",
    "\n",
    "We can actually use a \"memory profiler\" to measure how much memory is used during a given operation.\n",
    "\n",
    "Remember on Day 2 we installed it with `pip install --user memory_profiler` (otherwise: need to run this command now, on the terminal).\n",
    "\n",
    "Then, we can activate this functionality in the notebook by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can profile the memory usage of all the code in a cell by using `%%memit`` (similar to ``%%time``):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%memit\n",
    "pass # total memory currently used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%memit\n",
    "x = np.zeros(int(1e8), dtype='float32') # 100 million numbers, 4 bytes each = 4e8 bytes = 400 MB\n",
    "x += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Re-consider the previous exercise. Actually profile the memory usage of the numpy vs. numba functions (use `%memit`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automatic parallelism with numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba can also automatically run code in parallel on multiple CPU cores. It can also run code automatically on [GPUs](https://numba.readthedocs.io/en/stable/cuda/overview.html) (outside of our scope). Let's take a quick look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import njit, prange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: decorating with `@njit` is the same as with `@jit(nopython=True)`.\n",
    "\n",
    "You always want to use this \"no-python\" mode, which means that the code must compile successfully and run outside of python, outside an error is thrown.\n",
    "\n",
    "If you remove this requirement, then numba will simply run code it cannot compile as normal, but then you will silently have slow code. Better to:\n",
    "* (i) extract the compute intensive code from your existing Python into a new function,\n",
    "* (ii) make it accept simple arguments and return simple outputs,\n",
    "* (iii) make it use only simple python and numpy functionality, an\n",
    "* (iv) make sure numba can compile it.\n",
    "\n",
    "> `prange` is like `range`, but \"parallel\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@njit(parallel=True)\n",
    "def prange_test(A):\n",
    "    s = 0\n",
    "    # without \"parallel=True\" in the jit-decorator, the prange statement is equivalent to range\n",
    "    for i in prange(A.shape[0]):\n",
    "        s += A[i]\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 100, 100000000)\n",
    "%time print(prange_test(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"Wall time\" is the actual time (on a clock) which has gone by.\n",
    "The \"CPU time\" is the amount of work that the computer has done.\n",
    "\n",
    "For example, if you have 2 cores working for 1 second each, at the same time, the CPU time is 2 seconds.\n",
    "\n",
    "In this case, the ratio between the two gives us the speed-up, i.e. about the number of cores used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Write a function `stddev(x)` which computes the [standard deviation](https://en.wikipedia.org/wiki/Standard_deviation) $\\sigma$ of an array of input numbers `x`. \n",
    "\n",
    "> Hint: First compute the mean $\\mu = \\frac{1}{N} \\sum_{i=1}^{N} x_i$. Then compute $\\sigma^2 = \\frac{1}{N} \\sum_{i=1}^{N} (x_i - \\mu)^2$. Then, take the square root to get $\\sigma$.\n",
    "\n",
    "Do this the slow way, i.e. use loops! Benchmark it using a large array for `x`. Then, accelerate it with numba, and compute the speed-up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] Calling C or C++ (or fortran) code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We said \"writing C/C++ code is hard, we rarely want to do this\" - which is true. But:\n",
    "\n",
    "* sometimes you have no choice (complex algorithm, need advanced C/C++/fortran features)\n",
    "* sometimes you have existing code (from colleagues)\n",
    "* sometimes you want to use an existing library (which only exists in C/C++/fortran)\n",
    "\n",
    "The easiest approach is to use the **pybind11** library. This is a lightweight header-only library that exposes C++ types in Python and vice versa. Its main use is to create **Python bindings (i.e. wrappers) of existing C/C++ code**.\n",
    "\n",
    "It works like this. First, we would start with a file `example.cpp` with a function we would like to use:\n",
    "\n",
    "```c++\n",
    "    int add(int i, int j) {\n",
    "        return i + j;\n",
    "    }\n",
    "```\n",
    "\n",
    "We would like to be able to call this function from python. The most basic syntax is to change the source code file by adding the following lines:\n",
    "\n",
    "```c++\n",
    "    #include <pybind11/pybind11.h>\n",
    "\n",
    "    int add(int i, int j) {\n",
    "        return i + j;\n",
    "    }\n",
    "\n",
    "    PYBIND11_MODULE(example, m) {\n",
    "        m.doc() = \"pybind11 example\"; // optional module docstring\n",
    "\n",
    "        m.def(\"add\", &add, \"A function that adds two numbers\");\n",
    "    }\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "1. Right-click in the file explorer on the right, select \"New File\". Rename it to `example.cpp`. Open it, and paste the code above.\n",
    "2. Open a terminal (File -> New -> Terminal), and compile this source file with the command:\n",
    "\n",
    "`c++ -O3 -Wall -shared -std=c++11 -fPIC $(python3 -m pybind11 --includes) example.cpp -o example.so\n",
    "`\n",
    "\n",
    "3. A file `example.so` should have been made - this is a \"shared library\", i.e. a binary file with functions which are meant to be called from other programs.\n",
    "4. Try to run `import example` in Python, then get its documentation with `?example`, and finally, try to use the `add` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numpy arrays with pybind11\n",
    "\n",
    "You can also pass numpy arrays into your C++ function. (Also possible to create new numpy arrays in the C++ code, but a bit more complicated). Let's update our source code to:\n",
    "\n",
    "```c++\n",
    "    #include <pybind11/pybind11.h>\n",
    "    #include <pybind11/numpy.h>\n",
    "\n",
    "    namespace py = pybind11;\n",
    "\n",
    "    int add(int i, int j) {\n",
    "        return i + j;\n",
    "    }\n",
    "\n",
    "    float multiply_vecs_then_sum(py::array_t<float> a, py::array_t<float> b) {\n",
    "        float sum = 0.0;\n",
    "        auto buf_a = a.unchecked<1>(); // request direct memory pointer with no checks\n",
    "        auto buf_b = b.unchecked<1>(); \n",
    "\n",
    "        for(auto i = 0; i < a.size(); i++)\n",
    "            sum += buf_a[i] * buf_b[i];\n",
    "        return sum;\n",
    "    }\n",
    "\n",
    "    PYBIND11_MODULE(example, m) {\n",
    "        m.doc() = \"pybind11 example\"; // optional module docstring\n",
    "\n",
    "        m.def(\"add\", &add, \"A function that adds two numbers\");\n",
    "        m.def(\"my_f\", &multiply_vecs_then_sum, \"A pretty strange function!\");\n",
    "    }\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Update the `example.cpp` source file with the code above, re-compile the `.so` library, and re-import (you may need to Kernel -> Restart Kernel). Create two numpy arrays and test our new function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One glance at the [docs for using numpy with pybind11](https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html) shows how things get complicated quickly.\n",
    "\n",
    "But! It's a very powerful technique. You can pass data from python to C/C++ (without making any copies, reading it directly from memory), use advanced algorithms and techniques in C/C++, and return the result. From the point of view of the Python user, it's always just as simple as `import library` and `library.compute()`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid #bbb; margin: 30px 0\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The end!\n",
    "\n",
    "All done? You're free! I hope you enjoyed the course. \n",
    "\n",
    "> Don't forget to email me your homework solutions by the due date!\n",
    "\n",
    "  -- Dylan Nelson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid #bbb; margin: 30px 0\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 5 Practice Problem - Parallelized $\\pi$ Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will approximate $\\pi \\simeq 3.1415$ using one of the many series expressions:\n",
    "\n",
    "$$ \\pi = \\frac{4}{N} \\sum_{i=1}^N \\frac{1}{1 + \\left(\\frac{i-0.5}{N}\\right)^2} $$\n",
    "\n",
    "The accuracy increases as $N$ gets larger.\n",
    "\n",
    "## Task A\n",
    "\n",
    "Write a function `series_pi(N)` which computes the above approximation to $\\pi$. Run it for values of $N$ spanning several orders of magnitude, and compare the result to `np.pi`, as a function of $N$, with a y-log plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task B\n",
    "\n",
    "Use the `multiprocessing` module to create a pool of python workers. Each worker should use `series_pi` to compute a **subset** of the needed terms for the requested $N$, then the results should be combined. Use `map`. Modify the function as needed. Check the performance benefit versus the serial version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task C\n",
    "\n",
    "Switch to `threading` to solve the same problem. Make the target function for each thread numba-accelerated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task D\n",
    "\n",
    "What if you wanted to calculate $\\pi$ to a certain accuracy (with respect to `np.pi`)? Use a lock to protect a shared variable which stores the current estimated value, and is periodically updated as calculations progress. Once the accuracy reaches a goal (tolerance), the threads should stop working."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid #bbb; margin: 30px 0\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 5 Challenge Problem - Bitcoin Historical Data with Pandas\n",
    "\n",
    "A very popular and potentially powerful library is `pandas`. This is particularly useful for loading/analyzing/saving tabular data, e.g. large CSV files, especially data with dates, or anything which looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "display.Image(\"images/day5_pandas.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you're interested, read and work through the [10 minutes to pandas tutorial](https://pandas.pydata.org/pandas-docs/stable/user_guide/10min.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task A\n",
    "\n",
    "Use `pandas` to load [data/day5_btc_eur_daily.txt](data/day5_btc_eur_daily.txt). This file contains historical data on the price of bitcoin, specifically the BTC/EUR exchange rate, on a daily cadence, since the start of 2020.\n",
    "\n",
    "You can use `pd.read_csv()`, followed by `.describe()` to get an overview of the structure of the loaded data (and make sure pandas interpreted it correctly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task B\n",
    "\n",
    "Use pandas to resample the daily historical data to a weekly cadence. Plot the resulting time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task C\n",
    "\n",
    "Compute the running average ([SMA](https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average)), with a time window of one month. Overplot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task D\n",
    "\n",
    "\"Technical indicators\" are statistics computed on the price of e.g. a stock, or other financial instrument. One example is the [RSI](https://en.wikipedia.org/wiki/Relative_strength_index). Compute it and include it in your plot as a sub-plot below the main price panel. What is your interpretation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}