day_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Programming for Scientists - Day 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "In the fourth day we will return to some more advanced fundamentals in Python:\n",
    "\n",
    "* classes\n",
    "* object-oriented programming\n",
    "* decorators\n",
    "* error handling\n",
    "\n",
    "We will then look at one other important data analysis library:\n",
    "\n",
    "* [scipy](https://scipy.github.io/devdocs/tutorial/general.html) (for many mathematical and numerical algorithms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [1] Classes and OO (object-oriented) programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classes are the epitomy of \"object-oriented (OO) programming\": they bundle data (variables) and functionality (methods) together. \n",
    "\n",
    "A class provides the blueprint for a new object. You can then create new \"instances\", as many as needed, of a class. Each class instance has its own attributes (data), and share common functions/methods.\n",
    "\n",
    "Python classes have similar functionality and structure to other languages, such as C++.\n",
    "\n",
    "> Note: It is **typically not necessary to use classes** for scientific data analysis. Loading/creating, analyzing, and plotting data is likely more clear and also understandable to others without fancy Python techniques: **keep it simple!**. However, more complex programs and large code bases will start to benefit from OO techniques, and you will also encounter them in Python code you are reading/needing to fix/update/use.\n",
    "\n",
    "The simplest form of a class definition looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyClass:\n",
    "    \"\"\"A simple example class\"\"\"\n",
    "    i = 12345\n",
    "\n",
    "    def f(self):\n",
    "        return 'hello world'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: a general Python naming convention (PEP8) is that normal functions `should_be_lowercase_with_underscores()` while classes `ShouldBeCapsLikeThis`.\n",
    "\n",
    "We can then create a single instance of this class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = MyClass()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.f()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: We called `.f()` without an argument, even though it was expecting an argument named `self`. What actually happens here?\n",
    ">\n",
    "> The special aspect of class methods is that the instance object (`x`) is passed as the first argument, always. It is common practice to call this argument `self` (although the name is arbitrary).\n",
    "\n",
    "When we created the instance and named it `x`, we instantiated an empty object. It is common to perform some startup tasks in this \"init\", as so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyClass:\n",
    "    \"\"\"A simple example class\"\"\"\n",
    "    def __init__(self):\n",
    "        print('__init__() called')\n",
    "        self.data = []\n",
    "\n",
    "    def f(self):\n",
    "        return 'hello world'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = MyClass()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `__init__()` method may also accept arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyClass:\n",
    "    \"\"\"A simple example class\"\"\"\n",
    "    def __init__(self, value1, value2):\n",
    "        print('__init__() called')\n",
    "        self.data = [value1, value2]\n",
    "\n",
    "    def f(self):\n",
    "        return 'hello world'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = MyClass(1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classes have two types of attributes:\n",
    "1. Data (variables)\n",
    "2. Methods (functions)\n",
    "\n",
    "Data variables in classes do not need to be declared ahead of time, nor do function methods. Everything in **Python is dynamic**! Just like normal local variables, they \"spring into existence\" when first assigned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.more_data = []\n",
    "x.counter = 3\n",
    "\n",
    "for i in range(x.counter):\n",
    "    x.more_data.append(i)\n",
    "\n",
    "print(x.more_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note also an important feature of Python, which applies to class methods as well as all normal functions: you can give them another name by which they can be called."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_function = x.f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_function()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Write a new class `Lottery` with a state variable `winning_number`, which upon instantiation is set to a random integer between `0` and `5`. Create a class method `.draw_ticket()` which selects a random integer between `0` and `5`, and if it equals the previously set winning number, prints \"Bingo!\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instance variables vs. class variables\n",
    "\n",
    "Instance variables are for data unique to each instance and class variables are for attributes and methods shared by all instances of the class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Dog:\n",
    "    kind = 'canine'         # class variable shared by all instances\n",
    "\n",
    "    def __init__(self, name):\n",
    "        self.name = name    # instance variable unique to each instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1 = Dog('Max')\n",
    "d2 = Dog('Lassie')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1.kind, d2.kind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1.name, d2.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to default function arguments, **be careful when using mutable objects as (shared!) class variables**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Dog:\n",
    "    tricks = []             # class variable which will be common to all instances: just don't do it!\n",
    "\n",
    "    def __init__(self, name):\n",
    "        self.name = name    # instance variable unique to each instance\n",
    "        \n",
    "    def add_trick(self, trick):\n",
    "        self.tricks.append(trick)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1 = Dog('Max')\n",
    "d2 = Dog('Lassie')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1.add_trick('play dead')\n",
    "d2.add_trick('roll over')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1.tricks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Methods can call other methods using `self`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyCollection:\n",
    "    def __init__(self):\n",
    "        self.items = []\n",
    "\n",
    "    def add_item(self, x):\n",
    "        self.items.append(x)\n",
    "\n",
    "    def add_twice(self, x):\n",
    "        self.add_item(x)\n",
    "        self.add_item(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = MyCollection()\n",
    "x.add_twice('hi')\n",
    "print(x.items)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a new class `Vehicle` which accepts two initialization variables: `fuel_kg` and `mileage_km_per_kg`. Write method `.drive(km)` which (i) decreases the available fuel given the requested distance and the mileage, (ii) prints the total distance driven so far, and the remaining fuel. If no more fuel is available, driving further should not be possible. Write a method class `.out_of_fuel()` returning a bool to indicate this state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inheritance\n",
    "\n",
    "A class can [inherit](https://en.wikipedia.org/wiki/Inheritance_(object-oriented_programming)) from a \"base class\" (or \"parent class\"), meaning that the new class is based on the base class, i.e. has similar structure and functionality, but likely with additions. The resulting class is called a \"derived class\".\n",
    "\n",
    "Inheritance is a foundation of object-oriented programming: it allows code to be structured into a hierarchy of models (i.e. classes).\n",
    "\n",
    "The syntax is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Animal:\n",
    "    def __init__(self, name):\n",
    "        self.name = name\n",
    "\n",
    "class Dog(Animal):\n",
    "    kind = 'canine'\n",
    "\n",
    "class Cat(Animal):\n",
    "    kind = 'feline'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Dog('fluffy')\n",
    "y = Cat('whiskers')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derived classes may override methods of their base classes (by using the same name)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Animal:\n",
    "    \"\"\" Parent class. \"\"\"\n",
    "    def __init__(self, name):\n",
    "        self.name = name\n",
    "        \n",
    "    def speak(self):\n",
    "        print('Animals make a noise.')\n",
    "\n",
    "class Dog(Animal):\n",
    "    \"\"\" First derived class. \"\"\"\n",
    "    def speak(self):\n",
    "        # derived class method overrides parent\n",
    "        print('Dogs bark!') \n",
    "\n",
    "class Cat(Animal):\n",
    "    \"\"\" Second derived class. \"\"\"\n",
    "    # derived class does not have a speak(), so parent method will be used\n",
    "    pass\n",
    "    \n",
    "class Bird(Animal):\n",
    "    \"\"\" Third derived class. \"\"\"\n",
    "    def speak(self):\n",
    "        # we want to keep parent behavior, so call the parent class method\n",
    "        Animal.speak(self)\n",
    "        # but we also want to add some additional functionality\n",
    "        print('In this case, a bird-like noise.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = Dog('fido'), Cat('tiger'), Bird('iago')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.speak()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.speak()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z.speak()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few more notes about classes:\n",
    "\n",
    "1. Multiple inheritance is allowed: `class MyDerivedClass(BaseClass1, BaseClass2, BaseClass3)`.\n",
    "\n",
    "2. \"Private\" variables, which can be accessed only by class methods, and not from outside the class, don't exist in Python. But convention is that a name prefixed with an underscore, such as `_my_var` or `_my_func()` should be treated as non-public."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The concept of looping over \"containers\" of objects with `for` loops is pervasive throughout Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for element in [1, 2, 3]:\n",
    "    print(element)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Behind the scenes, the for loop is calling the [iter()](https://docs.python.org/3/library/functions.html#iter) function on the container object. This returns an \"iterator\" object, with a method called `__next__()`, which accesses elements in the container one at a time. When there are no more elements, a \"StopIteration\" exception is raised, telling the for loop to stop.\n",
    "\n",
    "It is easy to add this behavior to a class we write."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Reverse:\n",
    "    \"\"\"Iterator for looping over a sequence backwards.\"\"\"\n",
    "    def __init__(self, data):\n",
    "        self.data = data\n",
    "        self.index = len(data)\n",
    "\n",
    "    def __iter__(self):\n",
    "        return self\n",
    "\n",
    "    def __next__(self):\n",
    "        if self.index == 0:\n",
    "            raise StopIteration\n",
    "            \n",
    "        self.index = self.index - 1\n",
    "        return self.data[self.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rev = Reverse('science!')\n",
    "\n",
    "for c in rev:\n",
    "    print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Copy the `Reverse` class from above, rename it `Forward`, and modify it so that it iterates \"normally\" i.e. in the forward direction through its data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2] Decorators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In python, functions are objects, just like everything else. We saw above how we can create a second name for a function (`my_function = x.f`). Functions can also be passed around and used as arguments, just like anything else."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def square_op(x):\n",
    "    return x**2\n",
    "\n",
    "def cube_op(x):\n",
    "    return x**3\n",
    "\n",
    "def do_operation(data, op_func):\n",
    "    return op_func(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "do_operation(3, cube_op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also:\n",
    "* define functions inside other functions (often a nice way to re-use some plotting/analysis snippet multiple times - don't go overboard).\n",
    "* return functions from a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def children_greetings(child_num):\n",
    "    \n",
    "    # define an \"inner function\"\n",
    "    def first_child():\n",
    "        return \"Hi, I am Sophia.\"\n",
    "    \n",
    "    # define another inner function\n",
    "    def second_child():\n",
    "        return \"My name is Liam.\"\n",
    "    \n",
    "    # actual statments of this function\n",
    "    if child_num == 1:\n",
    "        return first_child\n",
    "    else:\n",
    "        return second_child"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first = children_greetings(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **decorator** uses these concepts, and provides a simple syntax with the **@ symbol** for calling \"higher-order functions\". A higher-order function takes another function, and extends its behavior, without explicitly modifying it.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_decorator(func):\n",
    "    def wrapper():\n",
    "        print(\"Before func() is called.\")\n",
    "        func()\n",
    "        print(\"After func() is called.\")\n",
    "        \n",
    "    return wrapper\n",
    "\n",
    "@my_decorator\n",
    "def print_hi():\n",
    "    print(\"Hello!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_hi()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The more general case should accept arguments (any number), and also make sure a return from the wrapped function is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_decorator(func):\n",
    "    def wrapper(*args, **kwargs):\n",
    "        # do something before\n",
    "        value = func(*args, **kwargs)\n",
    "        # do something after\n",
    "        return value\n",
    "        \n",
    "    return wrapper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A real world example: let's create a decorator which times the execution of any function it is applied to, and prints the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "def timer(func):\n",
    "    def wrapper(*args, **kwargs):\n",
    "        start_time = time.time()\n",
    "        value = func(*args, **kwargs)\n",
    "        elapsed_time = time.time() - start_time\n",
    "        print(f'Executing {func.__name__} took [{elapsed_time:.2f}] sec.')\n",
    "        return value\n",
    "        \n",
    "    return wrapper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@timer\n",
    "def make_calculation():\n",
    "    for i in range(100000):\n",
    "        j = i**3 + np.sqrt(i) + i\n",
    "    return j\n",
    "        \n",
    "@timer\n",
    "def another_calculation():\n",
    "    for i in range(100000):\n",
    "        j = i**2 + np.exp(i/1000) + i*5\n",
    "    return j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "make_calculation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "another_calculation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other common uses for decorators:\n",
    "* Debug code (e.g. print out the function called and its arguments, then call it)\n",
    "* Slow down code (e.g. sleep one second inbetween calls to an online API)\n",
    "* Registry (e.g. just save the name of the function to a dict, return original function without wrapping it)\n",
    "* Verify state (e.g. `@is_logged_in` before running a function which requires user authentication)\n",
    "* Validate output (e.g. check that a returned numpy array has no NaN values, or that a web API return is valid JSON)\n",
    "\n",
    "More advanced aspects of decorators:\n",
    "* In classes, can also decorate methods in the same way, and/or decorate the entire class!\n",
    "* Decorators can be nested (e.g. more than one applied).\n",
    "* Decorators can accept arguments which modify their behavior.\n",
    "* Decorators can have state (e.g. count the number of times a function has been called).\n",
    "* You can use a class, instead of a function, as a decorator.\n",
    "* You can create a function cache (e.g. save the results of an expensive calculation) - see also [lru_cache](https://docs.python.org/3/library/functools.html#functools.lru_cache).\n",
    "* You could add important metadata as an attribute of the wrapped function, e.g. physical units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Keeping track of units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_unit(unit):\n",
    "    \"\"\"Register a unit on a function\"\"\"\n",
    "    def decorator_set_unit(func):\n",
    "        func.unit = unit\n",
    "        return func\n",
    "    return decorator_set_unit\n",
    "\n",
    "@set_unit(\"cm^3\")\n",
    "def volume(radius, height):\n",
    "    return np.pi * radius**2 * height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "volume(3,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "volume.unit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several python packages which are focused on units. In particular, they all interface with numpy, allowing numpy arrays to also carry information about their units. They can also perform automatic (or explicitly requested) unit conversion, via symbolic algebra. Two are: [pint](https://pint.readthedocs.io/en/stable/) and [unyt](https://unyt.readthedocs.io/en/stable/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pint\n",
    "ureg = pint.UnitRegistry()\n",
    "vol = volume(3, 5) * ureg(volume.unit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vol.to(\"cubic inches\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] Error handling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've mostly ignored error messages until now.\n",
    "\n",
    "There are (at least) two kinds of errors: syntax errors and exceptions.\n",
    "\n",
    "## Syntax errors\n",
    "\n",
    "Syntax errors, also known as parsing errors, are very common while you are learning Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "while True print('Hello world')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parser repeats the offending line and displays a little ‘arrow’ pointing at the earliest point in the line where the error was detected. The error is caused by (or at least detected at) the token preceding the arrow: in the example, the error is detected at the function `print()`, since a colon (`':'`) is missing before it. File name and line number are printed so you know where to look in case the input came from a script."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exceptions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even if a statement or expression is syntactically correct, it may cause an error when it is executed.\n",
    "\n",
    "Errors detected during execution are called \"exceptions\". They are sometimes fatal (cause the program to stop), but do not need to be. By \"handling\" them, you can control the behavior of what happens.\n",
    "\n",
    "Most exceptions are not handled, however, and result in error messages as shown here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "10 * (1/0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "4 + spam*3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'2' + 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line of the error message indicates what happened. Exceptions come in different types, and the type is printed as part of the message: the types in the example are `ZeroDivisionError`, `NameError` and `TypeError`. The string printed as the exception type is the name of the [built-in exception](https://docs.python.org/3/library/exceptions.html#bltin-exceptions) that occurred."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling exceptions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can handle certain exceptions, if we want to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    x = 5\n",
    "    z = 0\n",
    "    y = x / z # error occurs, y is not assigned\n",
    "except ZeroDivisionError:\n",
    "    print('Not good! Division by zero in the try clause.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `try`,`except` syntax is similar to error handling in other languages. The `try` statement works by:\n",
    "\n",
    "1. First, the try clause (the statement(s) between the try and except keywords) is executed.\n",
    "\n",
    "2. If no exception occurs, the except clause is skipped and execution of the try statement is finished.\n",
    "\n",
    "3. If an exception occurs during execution of the try clause, the rest of the clause is skipped. Then, if its type matches the exception named after the except keyword, the except clause is executed, and then execution continues after the try/except block.\n",
    "\n",
    "4. If an exception occurs which does not match the exception named in the except clause, it is passed on to \"outer\" try statements; if no handler is found, it is an unhandled exception and execution stops with an error message.\n",
    "\n",
    "You can handle multiple exceptions, listed in a tuple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    x = 5\n",
    "    z = 0\n",
    "    y = x / z\n",
    "except (ZeroDivisionError,TypeError, NameError):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raising exceptions\n",
    "\n",
    "You can intentionally generate an exception in code by using `raise`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute(x):\n",
    "    if x == 0:\n",
    "        raise Exception(\"This is a failure!\")\n",
    "    return x*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compute(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compute(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also define new exceptions, by making a class which derives from the `Exception` base class.\n",
    "\n",
    "## finally: clean-up actions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `try` statement has an optional clause `finally`, which is intended to define clean-up actions that absolutely must be executed under all circumstances. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    y = 2 / 0\n",
    "finally:\n",
    "    print('Clean-up e.g. close open files, write a checkpoint file.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a `finally` clause is present, the `finally` clause will execute as the last task before the `try` statement completes. The `finally` clause runs whether or not the `try` statement produces an exception."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Debugging\n",
    "\n",
    "There are two classic ways to debug a problem in code:\n",
    "\n",
    "1. Put lots of `print()` commands everywhere, to try to identify when and where variables have unexpected values which lead to a problem.\n",
    "2. Pause the program at the moment an error occurs (a `breakpoint` or `tracepoint`), in order to inspect its state (e.g. variables) interactively.\n",
    "\n",
    "The first is fairly self-explanatory - always a good idea when you are developing new code. You can always remove (or comment out) these later.\n",
    "\n",
    "For the second, if an error has just occured, you can start up a debugger at that point, by entering (i.e. immediately in the next cell):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%debug"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also request to **automatically** start the debugger whenever an error occurs. Entering `%pyb` in a cell toggles this option between 'on' and 'off':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are using the ipython interpreter, instead of the Jupyter notebook, you can add the following contents to your `~/.ipython/profile_default/ipython_config.py` file:\n",
    "\n",
    "    c = get_config()\n",
    "\n",
    "    # lines of code to run at IPython startup.\n",
    "    c.InteractiveShellApp.exec_lines = ['%autoreload 2']\n",
    "\n",
    "    # A list of dotted module names of IPython extensions to load.\n",
    "    c.InteractiveShellApp.extensions = ['autoreload']\n",
    "\n",
    "    # Automatically call the pdb debugger after every exception.\n",
    "    c.TerminalInteractiveShell.pdb = True\n",
    "    \n",
    "The last turns on the automatic debugger. The first two automatically watch for changes to any files, and if any module, file, function, etc is changed (i.e. in a text editor), then it is immediately reloaded, so that you always run this 'current' version."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Run the following code, and debug the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n",
    "\n",
    "def subtract_mean(x1, x2):\n",
    "    \"\"\" Return the first input array after subtracting it by the number of elements in the second. \"\"\"\n",
    "    x1_new = x1 - x2.shape(0)\n",
    "    return x1_new\n",
    "\n",
    "data1 = np.array([1,2,3,4,5])\n",
    "data2 = np.array([2,1,4])\n",
    "\n",
    "subtract_mean(data1, data2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [5] Scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to numpy and matplotlib, the [Scipy](http://www.scipy.org) library is the third most used in scientific data analysis.\n",
    "\n",
    "Scipy is basically a very large library of algorithms (numerical and mathematical) that you can use. The [documentation](https://docs.scipy.org/doc/scipy/tutorial/index.html) has a complete list of functions, in the following categories:\n",
    "\n",
    "* Special functions (scipy.special)\n",
    "* Integration (scipy.integrate)\n",
    "* Optimization (scipy.optimize)\n",
    "* Interpolation (scipy.interpolate)\n",
    "* Fourier Transforms (scipy.fft)\n",
    "* Signal Processing (scipy.signal)\n",
    "* Linear Algebra (scipy.linalg)\n",
    "* Sparse eigenvalue problems with ARPACK\n",
    "* Compressed Sparse Graph Routines (scipy.sparse.csgraph)\n",
    "* Spatial data structures and algorithms (scipy.spatial)\n",
    "* Statistics (scipy.stats)\n",
    "* Multidimensional image processing (scipy.ndimage)\n",
    "* File IO (scipy.io)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We will just look at one or two examples.\n",
    "\n",
    "## Fitting a model (function) to data\n",
    "\n",
    "One method of analyzing data is to fit a theoretical model. A model is typically a function of one or more parameters, and \"fitting\" is the process of finding the optimal parameters such that the resulting model is a good description of the data -- this is the \"best fit model\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are a number of related routines in scipy for fitting: we will use the simplest, [curve_fit](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the typical simplest example: fitting a straight line to a set of `x,y` points in two dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We first create a \"fake\" dataset, representing a linear relation with some added random noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "x = np.random.uniform(0, 10, N)\n",
    "y = 4.0 * x + 1.0\n",
    "\n",
    "rng = np.random.default_rng()\n",
    "y_with_noise = y + rng.normal(loc=0.0, scale=2.0, size=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first visualize the data, and the \"ground truth\" model (which in a real experiment would be unknown):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x, y_with_noise, 'o', label='data')\n",
    "ax.plot(x, y, '-', label='truth')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If this was real data, we would start by defining a model which we believe represents the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line(x, a, b):\n",
    "    return a * x + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The arguments to the function should be `x` (i.e. the independent variable), followed by one or more parameters. The scipy `curve_fit()` function then uses the [non-linear least-squares](https://en.wikipedia.org/wiki/Non-linear_least_squares) method to find the best-fit parameters.\n",
    "\n",
    "> Note: This is overkill for this problem! You can use the simpler \"linear least-squares\" for fitting a line, e.g. with [np.polyfit()](https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html). However, the scipy `curve_fit()` is a useful optimization method, since it can work with *any* function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_opt, params_cov = curve_fit(line, x, y_with_noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The `curve_fit()` function has two returns. The first is an array of the optimal values for the parameters, such that the sum of the squared residuals `f(x,*params_opt)-y` is minimized. The second is a 2D array giving the estimated covariance: the diagonal elements give the variance (errors) of the parameter estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: the covariance (error estimate) doesn't mean much if we don't give input errors on the y data points.\n",
    "\n",
    "Let's visualize the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "y_fit = line(x, params_opt[0], params_opt[1])\n",
    "\n",
    "ax.plot(x, y_with_noise, 'o', label='data')\n",
    "ax.plot(x, y, '-', label='truth')\n",
    "ax.plot(x, y_fit, '--', label='best fit')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's now try and fit the data assuming each point has an error (standard deviation) of $2.0 \\pm 0.5$ (uniform)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_err = 2.0 + (rng.uniform(size=N)-0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.errorbar(x, y_with_noise, yerr=y_err, fmt='none');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now redo the `curve_fit()`, including the uncertainties on `y` as inputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "params_opt, params_cov = curve_fit(line, x, y_with_noise, sigma=y_err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now `params_cov` will contain a reasonable estimate of the variance and covariance of the parameters. The simplest information to extract are the standard deviations for each parameter, which are given by the square root of the diagonal elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_stddev = np.sqrt(np.diag(params_cov))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for param, stddev in zip(params_opt, params_stddev):\n",
    "    print(f'Best fit {param = :0.3f} with {stddev = :.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more complex models (i.e. high dimensionality of the parameter space, hard to find minima), more sophisticated techniques may be required for fitting, but `curve_fit()` is a good starting point, and will work well for simple cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: carefully consult the documentation of `curve_fit()`, and any scipy function to use, to make sure arguments and options are being set correctly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Given the following random data points `x,y`, with associated errors `y_err`, fit both linear ($y = mx+b$) and quadratic ($y=ax^2 + bx + c$) functions using `curve_fit()`, and overplot the data and the two best-fit models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(424242)\n",
    "\n",
    "N = 80\n",
    "x = rng.uniform(0.0, 20.0, N)\n",
    "y = np.polyval([1, 2, -3], x) + rng.normal(0.0, 15.0, N)\n",
    "y_err = rng.uniform(5, 10, N)\n",
    "\n",
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Load the same [data/day2_munich_temps.txt](data/day2_munich_temps.txt) data as before, which gives the temperature in Munich every day for several years. \n",
    "\n",
    "Fit the function $f(t) = a~\\cos{(2\\pi t + b)} + c$ with three parameters to the data, where $t$ is the time in years. What are the best-fit values of the parameters? Make a plot of the data and the best-fit model in the range 2000 to 2004. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolating data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common task is data interpolation: given a set of data points e.g. $y=f(x)$ at a given set of values $x$, we want to estimate new values $y'$ at corresponding points $x'$.\n",
    "\n",
    "For one-dimensional data, the [interp1d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d) function is perfect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, start by creating some synthetic data we can use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([0.0, 1.0, 3.0, 8.0])\n",
    "y = np.array([6.5, 4.0, 3.7, 3.1])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x,y,'s-');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use different kinds of interpolation, the simplest being [linear](https://en.wikipedia.org/wiki/Linear_interpolation), by creating an interpolating function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = interp1d(x, y, kind='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can then evaluate this function `f()` for any single value of $x$, or any numpy array $x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x,y,'s-',label='original data')\n",
    "\n",
    "xx = np.linspace(x.min(), x.max(), 30)\n",
    "ax.plot(xx, f(xx), 'o', markersize=5.0, label='interpolated values')\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, requesting interpolation outside the original range of the $x$ data values (extrapolation) raises an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can allow extrapolation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = interp1d(x, y, kind='linear', fill_value='extrapolate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of linear interpolation, we can also use higher-order methods, such as `quadratic` (second order) or `cubic` (third order).\n",
    "\n",
    "> Note: Both of these are spline-based, meaning that polynomials of the given order are constructed between each pair of input points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = interp1d(x, y, kind='cubic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x,y,'s-',label='original data')\n",
    "\n",
    "xx = np.linspace(x.min(), x.max(), 30)\n",
    "ax.plot(xx, f(xx), 'o', markersize=5.0, label='interpolated values')\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Always **be careful when interpolating**, and check against overfitting, oscillations, and other poorly behaved results.\n",
    "\n",
    "> You can also interpolate in higher dimensions: if the input points are regularly spaced (gridded), then in 2D with [interp2d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html#scipy.interpolate.interp2d) or in arbitrary dimensions with [interpn](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interpn.html#scipy.interpolate.interpn).\n",
    ">\n",
    "> For \"unstructured\" data, e.g. irregular point clouds, then [griddata](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata) uses significantly more complex methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Try quadratic interpolation on the same data. For this dataset, which of the three methods would you prefer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration\n",
    "\n",
    "Here we seek the answer for a discrete integral $\\int_a^b f(x) dx$ of some function $f(x)$ over some range $(a,b)$. By \"numerical\" we mean that the answer is a number, as opposed to symbolic integration (in a program like Mathematica), which can compute the answer in terms of variables.\n",
    "\n",
    "The main workhorse here is [scipy.integrate](https://docs.scipy.org/doc/scipy/tutorial/integrate.html), which implements many different standard approaches.\n",
    "\n",
    "In general, you can either integrate:\n",
    "\n",
    "* a python function `def f(x)` which returns values $y$ for any $x$.\n",
    "* a discrete \"function\" which is defined by two numpy arrays $x$ and $y$.\n",
    "\n",
    "In the first case, we can use the `quad(f,a,b)` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x**2 + 3*x + 1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integral, error = integrate.quad(f, 1.0, 2.0) # from a=1 to b=2\n",
    "\n",
    "print(integral)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the second case, we frequently only have measured data of a \"function\" sampled at some given points $x$.\n",
    "\n",
    "* If the sampling is even, you can use `romb` (Romberg integration).\n",
    "* If the sampling is uneven (arbitrary), you can use the `trapezoid` or `simpson` methods.\n",
    "\n",
    "Let's create samples of the same function as above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(1.0, 2.0, 200)\n",
    "y = x**2 + 3*x + 1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate.trapz(y,x) # renamed trapezoid() in more recent scipy versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate.simps(y,x) # renamed simpson() in more recent scipy versions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Reduce the number of samples: when does the numerically computed integral start to diverge from the true value? Write a loop to identify the number of samples required to obtain 1\\% accuracy. Plot the error versus the number of samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving differential equations\n",
    "\n",
    "We can compute the solution to [ordinary differential equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation) (ODEs) using scipy.\n",
    "\n",
    "> Remember that an ODE is a differential equation of functions of one variable $x$, and derivatives of those functions.\n",
    ">\n",
    "> The much more complex [partial differential equations](https://en.wikipedia.org/wiki/Partial_differential_equation) provide the foundation for our modern scientific understanding of many physical phenomena: sound, heat, diffusion, electrostatics, electrodynamics, thermodynamics, fluid dynamics, elasticity, general relativity, and quantum mechanics (Schrodinger equation, Pauli equation, etc). Such \"solvers\" become expensive and complex, especially in two or more dimensions.\n",
    "\n",
    "As an example, consider the equation for the motion of a [damped harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator#Damped_harmonic_oscillator). While in a simple harmonic oscillator the only force acting on the mass is the restoring force, in a damped harmonic oscillator there is in addition a frictional force which is always in a direction to oppose the motion, which we model as being proportional to the velocity v of the object. Balancing the forces (i.e. Newton's second law) gives: $\\ddot{x} + 2 \\gamma \\dot{x} + k x = 0$.\n",
    "\n",
    "To solve this equation we will use the [odeint()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `odeint()` solver accepts systems of first order ODEs. In our case, the second-order differential equation needs to be transformed into a system of two first-order equations:\n",
    "\n",
    "1. $\\dot{x} = v$\n",
    "\n",
    "2. $\\dot{v} = -k x - 2 \\gamma v$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def diffeq(y,t,k,gamma):\n",
    "    \"\"\" Specify first-order ODE system.\n",
    "    y is a 2-vector of [x,v] - must be the first argument.\n",
    "    t is time - must be the second argument.\n",
    "    k and gamma are parameters - must be the last arguments. \"\"\"\n",
    "    # allocate for result\n",
    "    dydt = np.zeros(2, dtype='float32')\n",
    "\n",
    "    # dx/dt = v\n",
    "    dydt[0] = y[1]\n",
    "    dydt[1] = -k*y[0] - 2*gamma*y[1]\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose constants\n",
    "k = 0.1\n",
    "gamma = 0.25\n",
    "\n",
    "# choose initial condition\n",
    "y0 = np.array([2.0,-1.3], dtype='float32')\n",
    "\n",
    "# choose requested times for solution\n",
    "times = np.linspace(0.0,30.0,60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = odeint(diffeq, y0, times, args=(k,gamma))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the solution for $x(t)$ and $v(t)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(times,solution[:,0],'s-',label='x')\n",
    "ax.plot(times,solution[:,1],'o:',label='v')\n",
    "ax.set_xlabel('t')\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a corresponding analytical solution available in this case, so we can check our numerical integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def damped_oscillator(time,y0,k,gamma):\n",
    "    \"\"\" Solution for the damped harmonic oscillator at a given time, for parameters (k,gamma), where y is a 2-vector of [x0,v0]. \"\"\"\n",
    "    omega = np.sqrt(k-gamma**2)\n",
    "    \n",
    "    A = (y0[1]+gamma*y0[0])/omega\n",
    "    B = y0[0]\n",
    "    \n",
    "    y = (A*np.sin(omega*time) + B*np.cos(omega*time)) * np.exp(-gamma*time)\n",
    "    \n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(times,solution[:,0],'-',label='numeric x')\n",
    "ax.plot(times[::2],damped_oscillator(times[::2],y0,k,gamma),'o',label='analytic x')\n",
    "ax.set_xlabel('t')\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "A critically damped harmonic oscillator has $\\gamma$ = $2 \\sqrt{k}$. Re-solve the new ODE and plot the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Change the equation system to represent a driven (or \"forced\") harmonic oscillator. Use the forcing term $f(x)= A_f \\sin (\\omega_f t)$, choosing the two (constant) parameters $A_f$ and $\\omega_f$ as you like (just add this term to the right hand side of the original equation). Solve and plot the initial (early time) behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other libraries\n",
    "\n",
    "In addition to numpy, matplotlib, and scipy, there are many other Python libraries for scientific data analysis which may be useful:\n",
    "\n",
    "* [pandas](https://pandas.pydata.org/) - with `DataFrame`, to load/manipulate/search data which is structured like tables.\n",
    "* [astropy](https://www.astropy.org/) - domain specific, for astronomy and astrophysics.\n",
    "* [emcee](https://emcee.readthedocs.io/en/stable/) - MCMC sampler (and many others).\n",
    "* [scikit-learn](https://scikit-learn.org/stable/) - \"classical\" machine learning (e.g. regression, classification, clustering with random forest, SVM, k-means, PCA, etc).\n",
    "* [pytorch](https://pytorch.org/) - modern machine learning (e.g. convolutional/deep/graph/... neural networks).\n",
    "* [tensorflow](https://www.tensorflow.org/) - the other modern machine learning library (all ML applications and research use one of these two).\n",
    "* many others - domain specific, for your field.\n",
    "\n",
    "You should always explore existing libraries and modules which may accelerate your science. On the other hand, using large complex packages as \"black boxes\" makes it very difficult to understand what is going on, e.g. at the level to convincingly explain it in a paper, talk, or to change/extend the method to your particular problem. In these cases, implementing algorithms and concepts \"from scratch\" may give you more flexibility in the long term."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid #bbb; margin: 30px 0\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 4 Practice Problem - Radioactive Decay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Radioactive decay](https://en.wikipedia.org/wiki/Radioactive_decay) is the emission of particles from an unstable atom, which causes it to lose energy. Common \"products\" of radioactivity are \"alpha particles\" (Helium nuclei) and \"beta particles\" (electrons).\n",
    "\n",
    "In this exercise we will analyze a dataset giving: the detection rate of alpha particles from a radioactive sample, which we will call \"element X\". \n",
    "\n",
    "The dataset is available in the [data/day4_radioactive_decay.txt](data/day4_radioactive_decay.txt) file:\n",
    "* First column: time, in days since the start of the experiment.\n",
    "* Second column: detection *rate*, in events per second.\n",
    "* Third column: uncertainty in the detection rate, in events per second.\n",
    "\n",
    "## Task A\n",
    "\n",
    "Load the data, and plot the rate of detection versus time, including error bars. Make sure the plot has all the required information to make it scientifically accurate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task B\n",
    "\n",
    "The rate of decay of a radioactive element can be modeled by the function $R(t) = R_0 \\exp{(-t/\\tau)}$ where $R_0$ is the rate of decay at time $t=0$, and $\\tau$ is the *mean lifetime* of the element. \n",
    "\n",
    "Define this function (three arguments: time, plus two parameters), and fit it to the data using `curve_fit()`.\n",
    "\n",
    "Print out the best-fit parameters and uncertainties. Overplot the resulting best-fit function on the data. Is the fit a good representation of the data? If not, why not?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task C\n",
    "\n",
    "Unfortunately, experimental data is never perfect. In this case, **background** radioactive decay (from sources other than those in the actual experiment) is an additional source of noise.\n",
    "\n",
    "We can include such a background, assuming it to be constant in time, by modifying our function $R(t) = R_0 \\exp{(-t/\\tau)} + R_{\\rm bg}$.\n",
    "\n",
    "Fit this new function to the data, and overplot it on the data. Draw a horizontal line to represent the background level $R_{\\rm bg}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Task D\n",
    "\n",
    "The \"goodness of fit\" can be quantified with a parameter called reduced $\\chi^2$, which is essentially the sum of the squared errors (model minus data, divided by uncertainty), normalized by the number of degrees of freedom.\n",
    "\n",
    "Specifically, we can write $\\chi_{\\rm red}^2 = \\frac{1}{N - p}\\sum_{i=1}^{N}~\\left(\\frac{y_i - m_i}{\\sigma_i}\\right)^2$ where:\n",
    "\n",
    "* $N$ is the number of datapoints\n",
    "* $p$ is the number of parameters for the model\n",
    "* $y_i$ are the data values\n",
    "* $m_i$ are the model values at the same positions\n",
    "* $\\sigma_i$ are the uncertainties on the data $y_i$.\n",
    "\n",
    "We typically consider a fit to be \"good\" if $\\chi_{\\rm red}^2 \\approx 1$.\n",
    "\n",
    "Compute $\\chi_{\\rm red}^2$ for both models above: how do the quantitative values compare to your visual impression from each plot?"
   ]
  },
  {
   "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 4 Challenge Problem - Solving N-Body Gravity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A [N-body simulation](https://en.wikipedia.org/wiki/N-body_simulation) is a numerical calculation of the force of gravity acting on $N$ particles, each with a mass $m$.\n",
    "\n",
    "The only trouble is: **every particle exerts a gravitational force on every other particle**. As a result the number of force interactions is $N \\times (N-1) \\simeq N^2$.\n",
    "\n",
    "Solving this problem \"brute force\" means computing all of these individual forces, which is why gravity is a classic example of a $\\mathcal{O}(N^2)$ computational problem.\n",
    "\n",
    "Yesterday we have seen the gravitational force law $F = \\frac{G~M_1~M_2}{r^2}$ which gives the gravitational force between two masses $M_1$ and $M_2$ separated by a distance $r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rewrite this in terms of the acceleration experienced by particle `i` due to the gravity of all the other particles:\n",
    "\n",
    "$$ \\vec{a}_i = G \\sum_{i \\ne j}^{N} m_j \\frac{\\vec{r}_j - \\vec{r}_i}{|\\vec{r}_j - \\vec{r}_i|^3} $$\n",
    "\n",
    "where $G$ is the gravitational constant. Consider the three-dimensional problem: $\\vec{a}$ and $\\vec{r}_i$ are 3-vectors, i.e. each has three components for the $x,y,z$ directions.\n",
    "\n",
    "> Note: The expression $|\\vec{r}_j - \\vec{r}_i|^3$ is the cube of the 2-norm (or [Euclidean norm](https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm)).\n",
    "\n",
    "> Note: Instead of the actual value $G = 6.6738 \\times 10^{-11}~\\rm{m}^3~\\rm{kg}^{-1}~\\rm{s}^{-2}$, instead we will take $G = 1$, i.e. we will work in a simpler unit system where G is unity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task A\n",
    "\n",
    "Write a `Particle` class to represent a single particle. It should be initialized with the mass, position, and velocity of the particle.\n",
    "\n",
    "Write a method `compute_accel` which accepts as an argument another instance of a Particle class, and computes the acceleration due to the interaction with that particle. The result should be stored as an attribute. (Hint: initialize such an attribute to zero, and add to it each time this method is called)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task B\n",
    "\n",
    "Write a `ParticleSystem` class to represent our collection of particles. It should be initialized with an integer `N` for the requested number of particles, and a float `mass` for the (constant) mass of each particle.\n",
    "\n",
    "Upon init, it should generate these particles, using a Gaussian \"standard normal\" distributions (zero mean, variance of one) for both positions and velocities.\n",
    "\n",
    "The particles (instances of the `Particle` class) should be stored as an attribute of the `ParticleSystem`.\n",
    "\n",
    "> Note: I suggest to 'shift to the center of mass frame': after creating the particles, subtract the mean velocity of the system (component-wise) from each particle velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task C\n",
    "\n",
    "Write a new method for your `ParticleSystem` called `compute_accelerations()` which should compute the total acceleration on each particle. (Hint: use a double for loop)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task D\n",
    "\n",
    "Now that we have the acceleration on each particle, we need to advance the system by some amount of time. This is called \"taking a timestep $\\Delta t$\". This will let us \"update\" the positions and velocities of the particles.\n",
    "\n",
    "One method of doing so is called the [leap-frog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration), which consists of three steps (\"KDK\"):\n",
    "\n",
    "1. First, each particle receives a half-step \"kick\": $\\vec{v}_i = \\vec{v}_i + ( \\Delta t / 2 ) \\vec{a}_i$.\n",
    "\n",
    "2. Then each particle undergoes a full-step \"drift\": $\\vec{r}_i = \\vec{r}_i + ( \\Delta t ) \\vec{v}_i$.\n",
    "\n",
    "3. Finally, each particle finishes with another half-step \"kick\" (as in \\#1).\n",
    "\n",
    "Write two methods `kick()` and `drift()` for the `ParticleSystem` which implement these updates.\n",
    "\n",
    "Then write a third method `take_timestep()` which accepts a single argument `dt`, the requested timestep. It should perform the following calculations in order: (half) kick, drift, compute new accelerations, (half) kick, and finally increase a `time` attribute by `dt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task E\n",
    "\n",
    "Add a final method `simulate()` to your `ParticleSystem`. It should start at $t=0$ and advance to $t=10$ with a constant $dt=0.01$. Try with a particle $\\rm{mass}=1.0$ and a small number e.g. $N=10$.\n",
    "\n",
    "Every few timesteps call a method `plot()` which visualizes the current particle positions (optional: and velocities, as arrows)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: we would like to \"re-draw\" (update) the same plot multiple times, in the notebook. One technique is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example code of updating a single plot\n",
    "from IPython import display\n",
    "\n",
    "x = []\n",
    "y = []\n",
    "\n",
    "for i in range(10):\n",
    "    if i == 0:\n",
    "        hdisplay = display.display(\"\", display_id=True)\n",
    "        fig, ax = plt.subplots()\n",
    "        line, = ax.plot(x,y)\n",
    "    else:\n",
    "        line.set_xdata(x)\n",
    "        line.set_ydata(y)\n",
    "        ax.set_title(i)\n",
    "\n",
    "        hdisplay.update(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task F\n",
    "\n",
    "Does it work? Does the behavior look reasonable?\n",
    "\n",
    "One common ingredient in numerical N-body integration is to add a \"softening\" to the force. This simply means adding a small constant value $\\epsilon$ to the distance denominator, in order to avoid numerical issues when the distance between two particles approaches zeros (such that the force approaches infinity):\n",
    "\n",
    "$$ \\vec{a}_i = G \\sum_{i \\ne j}^{N} m_j \\frac{\\vec{r}_j - \\vec{r}_i}{|\\vec{r}_j - \\vec{r}_i + \\epsilon|^3} $$\n",
    "\n",
    "Try the same integration with $\\epsilon = 0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task G (optional)\n",
    "\n",
    "An important method to assess the accuracy of a numerical integration is to verify an \"invariant\" of the physical problem. Meaning, some property of the system which be \"conserved\", and so not change during the numerical integration.\n",
    "\n",
    "In case of an ensemble of point masses acting under the force of gravity, the total energy of the system is given by:\n",
    "\n",
    "$$E_{\\rm tot} = \\sum_i \\frac{1}{2} m_i v_i^2 - \\sum_{i=1,i<j}^{i<N} \\frac{G m_i m_j}{|r_j - r_i|}$$.\n",
    "\n",
    "The first term is the kinetic energy, while the second term is the potential energy. Add a computation of the total energy into your `ParticleSystem`, and print out the value occasionally. How well is it conserved?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "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
}