Lesson 6 - Classes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 6 - Classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Authors:**\n",
    "* Yilber Fabian Bautista \n",
    "* Niklas Becker\n",
    "* Sean Tulin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the end of this lecture you will be able to:\n",
    "* Understand the concept of  **namespaces**, and differentiate the  4 different **namespaces** in python\n",
    "* Understand the meaning of  **Scope**  of a code \n",
    "* Use **classes** to efficiently organize a code\n",
    "* Differentiate between **class variables and instances variables**\n",
    "* Use **inheritance** to define **classes** from other existing **classes**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Namespaces\n",
    "\n",
    "Formally, a **namespace** is a mapping from names to objects. When you type the name of a function or variable, the python interpreter has to map that string to the corresponding object. Depending on the **namespace**, the same name can lead to different objects. For example, there is the **namespace** of built-in functions, such as `abs()` and `max()`. There is also a **namespace** for included libraries, and one for each function invocation.\n",
    "Some examples of the same name being resolved differently depending on the **namespace**:\n",
    "For instance,  although the `abs()` and `np.abs()` functions have the same effect on numbers, they belong to two different **namespaces**.  If we run the lines\n",
    "```py\n",
    "import numpy as np\n",
    "\n",
    "print(abs(-1))       # This is the python built-in absolute value function\n",
    "print(np.abs(-1))    # This is the abs function from the numpy library \n",
    "print(np.abs == abs)\n",
    "\n",
    "```\n",
    "we will get the  output\n",
    "```py\n",
    "1\n",
    "1\n",
    "False\n",
    "```\n",
    "That is, the functions do effectively the same, but correspond  different implementations in python, as seen from the last line in the output.  There is no relation between names in different **namespaces**! You can have two functions with the same name but in different **namespaces** do completely different things.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "There are 4 types of **namespaces**:\n",
    "- Built-In\n",
    "- Global\n",
    "- Enclosing\n",
    "- Local\n",
    "\n",
    "When the execution encounters a **name**, the **namespaces** are searched from **bottom to top** in this list until a name can be resolved.\n",
    "\n",
    "#### Built-in\n",
    "The first type of namespace was already shown in the example above. There are built-in functions that are available at all times. To see them in a list, run the following line\n",
    "```py\n",
    "dir(__builtins__)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Global\n",
    "The global **namespace** is created at the start of any program and contains all names defined at the level of the main program. This is where the user can come in and define **variables** and **functions**. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "moritz\n"
     ]
    }
   ],
   "source": [
    "def max(a, b):\n",
    "    return \"moritz\"\n",
    "\n",
    "print(max(5, 10)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the defined function has the same name as the **built-in**  `max` function, our definition is in the **global namespace** and as already mentioned, it  is found before than the  the  Built-in `max` function in the **Built-in  namespace** \n",
    "\n",
    "The interpreter also creates **namespaces** for all modules that are imported, such as **numpy** in the example above. \n",
    "When importing modules, the names can  be made available at the global **namespace** of the program using the syntax\n",
    "```py\n",
    "from  library import function \n",
    "```\n",
    "For instance, if we import the `max` function from **numpy** as\n",
    "\n",
    "```py\n",
    "from numpy import max\n",
    "```\n",
    "and run \n",
    "```max([5, 10])\n",
    "```\n",
    "we will get as output `10`, whereas `max(5, 10)`, defined above, will fail. The reason for this is that although the two functions are  defined in the **global namestace**,  the interpreter searches the more recent definition of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Local and Enclosing\n",
    "Every function has its own **namespace** associated with it. The **local namespace** is created when the function is executed and 'forgotten' afterwards. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "local bar a: 10\n",
      "enclosing b: 6\n",
      "local foo a: 5\n",
      "local foo b: 6\n"
     ]
    }
   ],
   "source": [
    "def foo():\n",
    "    a = 5 #creates a new variable in the local namespace of foo\n",
    "    b = 6 #creates a new variable in the local namespace of foo\n",
    "    \n",
    "    def bar():\n",
    "        a = 10   # creates a new variable in the local namespace of bar\n",
    "        print('local bar a:',a) # finds the a from the local namespace and not the enclosing\n",
    "        print('enclosing b:',b) # find the b from the enclosing namespace\n",
    "    \n",
    "    bar()\n",
    "    print('local foo a:',a)\n",
    "    print('local foo b:',b)\n",
    "\n",
    "foo()\n",
    "#print(a) would fail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `bar` is defined inside of `foo`, so `foo` is the **enclosing namespace** for `bar`. A name lookup will first search the local **namespace**, then the enclosing **namespace**, then the **global**, and finally the **built-in** one, as mentioned above.\n",
    "\n",
    "### Scope\n",
    "The **scope** of a name refers to the region of a code where it has meaning, where it can be associated with the object. Once a variable is out of scope, it is forgotten and no longer accessible.\n",
    "```py\n",
    "def foo():\n",
    "    c = 5 # c is defined\n",
    "    def bar():\n",
    "        print(c)  \n",
    "    bar()  \n",
    "    print(c)  # until here 'c' has meaning\n",
    "\n",
    "foo()\n",
    "print(c)   # outside of the scope of c\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above discussion show the importance of  to understanding  the concepts of **namespace** and **scope**, in order to have a full control, and understand the functionality  of your code!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classes\n",
    "### What is a class?\n",
    "A class is a 'blueprint' for creating objects that binds together data (variables) and manipulation of such data (via *methods*). Classes are useful ways of organizing the code, in addition to providing objects with code that can be reused. \n",
    "### Class Definition\n",
    "A class is defined using the following syntax:\n",
    "\n",
    "```py\n",
    "class MyClass:\n",
    "    classVariable = \n",
    "    def classMethod(self):\n",
    "        \n",
    "```\n",
    "A class can have variables and methods (functions) associated to them. These are collectively called **attributes** of the class. The class definition must be executed before they can be used. The class defines its own **namespace**. It needs to be accessed  with the syntax such as\n",
    "\n",
    "```py\n",
    "MyClass.classVariable\n",
    "MyClass.classMethod()\n",
    "```\n",
    "\n",
    "Let us define an specific example of a class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyClass():\n",
    "    classVariable = \"This is a class variable\"\n",
    "    \n",
    "    def classMethod(self):\n",
    "        return \"Hello Class World\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is a class variable\n",
      "Hello Class World\n"
     ]
    }
   ],
   "source": [
    "print(MyClass().classVariable)\n",
    "print(MyClass().classMethod())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class Instances\n",
    "In order to use a class **attributes**, one has to specify the entries of the class, and it is cumbersome to do it all the time we want to use the class. For that reason, one can  create  **instances** of a class. An **instance**  of a class is created  with the following function notation:\n",
    "```py\n",
    "my_instance = MyClass()\n",
    "```\n",
    "This object now has the methods and variables associated with that class\n",
    "```py\n",
    "print( my_instance.classMethod() )\n",
    "print( MyClass.classMethod(my_instance) )\n",
    "```\n",
    "Classes can access their own methods and variables using the **self** keyword. **self** represents the instance of the class. By using the **self** keyword we can access the attributes and methods of the class in python. It binds the attributes with the given arguments. These **attributes** can also be accessed from the outside"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Instance Variables\n",
    "Every instance of a class can also have their own variables - **instance variables** - associated to them. They are created similar to local variables and only associated with that instance. Consider the following example of a class\n",
    "\n",
    "```py\n",
    "class MyClass:\n",
    "    classVariable = \"This is a class variable\"\n",
    "    def startCounter(self):\n",
    "        self.counter = 0\n",
    "        b = 10\n",
    "    \n",
    "    def increaseCounter(self):\n",
    "        self.counter += 1\n",
    "\n",
    "classInstance = MyClass()    # creates an instance of MyClass\n",
    "\n",
    "```\n",
    "If we want to access  `classVariable` from `classInstance`, we simply do `classInstance.classVariable`, which will have as output the string ` \"This is a class variable\"`. We might wonder whether the variables `counter` and  `b` defined inside the **method** `startCounter()` can be accessed in the same way. If we did `classInstance.counter` we will get the error: \n",
    "```py\n",
    "AttributeError: 'MyClass' object has no attribute 'counter'\n",
    "```\n",
    "similar if we typed  `classInstance.b`. To access those  variables we first need to initialize the `startCounter()` method. \n",
    "\n",
    "```py\n",
    "classInstance.startCounter() # This function creates the counter variable of the instance\n",
    "print(classInstance.counter) # now it can be accessed from the outside\n",
    "```\n",
    "where now the output will be `0`. We might wonder whether the same procedure is true to access the value for the `b` variable. If we did \n",
    "```py\n",
    "classInstance.startCounter() # This function creates the counter variable of the instance\n",
    "print(classInstance.b) \n",
    "```\n",
    "the  output  will be the error  message\n",
    "\n",
    "```py\n",
    "AttributeError: 'MyClass' object has no attribute 'b'\n",
    "```\n",
    "The reason for this error is that unlike `counter`, `b` is not an **instance variable**. We learned then that instance variables are defined with the *self* keyword. Here `b`  corresponds to a local variable inside the `startCounter` method. \n",
    "\n",
    "Now, if we initialize the second method `increaseCounter`\n",
    "\n",
    "```py\n",
    "classInstance.increaseCounter()\n",
    "print(classInstance.counter)\n",
    "```\n",
    "we will see that  `counter` is still an instance variable but has increased its value by one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, these instance variables can differ between instances. We can see this in the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "0\n",
      "1\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "class MyClass:    \n",
    "    def startCounter(self):\n",
    "        self.counter = 0\n",
    "    \n",
    "    def increaseCounter(self):\n",
    "        self.counter += 1\n",
    "\n",
    "        \n",
    "# Create a class instance\n",
    "classInstance_1 = MyClass() \n",
    "\n",
    "# Start its counter\n",
    "classInstance_1.startCounter() \n",
    "print(classInstance_1.counter)\n",
    "\n",
    "# Increase the counter\n",
    "classInstance_1.increaseCounter() \n",
    "print(classInstance_1.counter)\n",
    "\n",
    "# Creates a second instance, unrelated to the first\n",
    "classInstance_2 = MyClass()\n",
    "\n",
    "# and starts its own counter\n",
    "classInstance_2.startCounter() \n",
    "print(classInstance_2.counter) \n",
    "\n",
    "# Increase its counter\n",
    "classInstance_2.increaseCounter()\n",
    "print(classInstance_2.counter)\n",
    "\n",
    "# Check that the  counter of the original instance was NOT affected\n",
    "print(classInstance_1.counter)     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialization\n",
    "The most important instance variables are usually assigned on initialization. For this purpose, there is a special function for classes, the *\\__init\\__*(self,...) function. The *\\__init\\__*(self,...) function is called when the instance is first created. The *\\__init\\__* function is called the **constructor** and can take a list of parameters that need to be passed on during initialization.\n",
    "\n",
    "As an example, let us imagine we want to create a **class** that has the profile information of a given person. We start by assigning the name of that person to the profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hey Michael\n",
      "Hey Janet\n"
     ]
    }
   ],
   "source": [
    "class Person:\n",
    "    def __init__(self, name, age):\n",
    "        self.name = name \n",
    "        self.age = age\n",
    "    \n",
    "    def call(self):\n",
    "        print(\"Hey \" + self.name)\n",
    "\n",
    "bill = Person(\"Michael\", 53)\n",
    "janet = Person(\"Janet\", 35)\n",
    "\n",
    "bill.call()\n",
    "janet.call()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class Variable vs Instance Variable\n",
    "Care must be taken to differentiate between class variables and instance variables. Class variables are copied on each initialization  of a class, while instance variables are created from scratch. This can have unintended effects for mutable variables, because their copies are *shallow*, i.e. they are copied by reference. To see this, let us continue with the class of our previous example, and add a `hobbies list` to our person profile as a **class variable**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Yoga', 'Marbles']\n",
      "['Yoga', 'Marbles']\n"
     ]
    }
   ],
   "source": [
    "class Person:\n",
    "    hobbies = []   # here, hobbies is defined as a class variable\n",
    "    \n",
    "    def __init__(self, name):\n",
    "        self.name = name\n",
    "        \n",
    "    def addHobby(self, hobby):\n",
    "        self.hobbies.append(hobby)\n",
    "\n",
    "\n",
    "michael = Person(\"Michael\")\n",
    "michael.addHobby(\"Yoga\")\n",
    "\n",
    "janet = Person(\"Janet\")\n",
    "janet.addHobby(\"Marbles\")\n",
    "\n",
    "print(michael.hobbies)\n",
    "print(janet.hobbies)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if  the `hobbies list` was assigned as an **instance variable** instead "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Yoga']\n",
      "['Marbles']\n"
     ]
    }
   ],
   "source": [
    "class Person:\n",
    "    \n",
    "    def __init__(self, name):\n",
    "        self.name = name\n",
    "        self.hobbies = []  # declare hobbies as an instance variable instead\n",
    "\n",
    "    def addHobby(self, hobby):\n",
    "        self.hobbies.append(hobby)\n",
    "\n",
    "\n",
    "michael = Person(\"Michael\")\n",
    "michael.addHobby(\"Yoga\")\n",
    "\n",
    "janet = Person(\"Janet\")\n",
    "janet.addHobby(\"Marbles\")\n",
    "\n",
    "print(michael.hobbies)\n",
    "print(janet.hobbies)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which generates two different outputs. Generally, it is **advisable** to use class variables only for **values that stay constant between all instances** and mostly use instance variables instead.\n",
    "\n",
    "\n",
    "\n",
    "#### Printing\n",
    "To simplify output, the classes can define a function that returns a string to represent that object. This is the *\\__str\\__* function. For instance\n",
    "\n",
    "```py\n",
    "class Person:\n",
    "    def __init__(self, name, age):\n",
    "        self.name = name      # creates an instance variable\n",
    "        self.age = age\n",
    "    \n",
    "    def __str__(self):\n",
    "        return self.name + \", \" + str(self.age)\n",
    "\n",
    "chihiro = Person(\"Chihiro\", 10)\n",
    "\n",
    "print(chihiro)\n",
    "str(chihiro)\n",
    "```\n",
    "\n",
    "### Why use classes?\n",
    "- Enables encapsulation: Bind together data and code that manipulates the data in the same place\n",
    "- Brings structure into the code, and allows **real-world mapping**\n",
    "- Polymorphism -> See below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1\n",
    "\n",
    "Make a **class** modeling a car with the attributes `model`, `cost`, `fuel efficiency` (km/liter), `condition`, and `mileage`. The `condition` is supposed to model the status of the car in terms of functionality, which of course diminishes over time. Implement a sensible `__init__` function and additional functions:\n",
    " - `Drive`: Takes distance and fuel price as an input. The distance is added to the mileage. It also prints an estimate of the fuel cost and adds it to the total cost. The condition of the car deteriorates at a rate of (1%/1000km). If this drops below 0% the car stops working.\n",
    " - `Repair`: The car is taken to the garage. The condition is set back to 100% at a cost of 100€/1%. The cost is printed and added to the total cost.\n",
    " - `LifetimeCost`: Outputs the mileage and total cost of the car so far\n",
    " - Honk: ;)\n",
    " \n",
    "* If you were to buy a new green VMW Imaginaris (cost `20k€`) with a fuel efficiency of `0.05L/km` what would be the approximated cost after driving for `10^5 km`?\n",
    " \n",
    "* And for the newer version VMW Imaginaris^2 for 25k€ with 0.04L/km?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the followng skeleton:\n",
    "\n",
    "class Car:\n",
    "    \n",
    "    def __init__(self, model, cost, fuel_efficiency, condition=1., mileage=0.):\n",
    "        pass \n",
    "    \n",
    "    def drive(self, distance, fuel_price):\n",
    "        pass\n",
    "        \n",
    "    def repair(self):\n",
    "        pass\n",
    "\n",
    "    def lifetimeCost(self):\n",
    "        pass\n",
    "\n",
    "ToyotaIm = Car(\"VMW Imaginaris\", 2e4, 0.05)\n",
    "ToyotaIm.drive(1e5, 2)\n",
    "ToyotaIm.repair()\n",
    "ToyotaIm.lifetimeCost()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inheritance\n",
    "In object-oriented programming, inheritance describes the process of basing a class on another class' implementation, thus **inheriting** its properties. \n",
    "* The **child class** inherits all variables and methods from the parent class and can expand on them. Let us see this explicit in the following example\n",
    "\n",
    "```py\n",
    "class A:\n",
    "    text = \"This is set in A\"\n",
    "    \n",
    "    def f(self):\n",
    "        print(self.text)\n",
    "    \n",
    "\n",
    "class B(A):   # This tells python that B inherits from A\n",
    "    text = \"This is set in B\"\n",
    "    \n",
    "    def g(self):\n",
    "        print(\"A new function\")\n",
    "\n",
    "a = A()\n",
    "a.f()     # calls f in A, which accesses the text variable in A\n",
    "b = B()\n",
    "b.f()     # calls f defined in A, but accesses the text variable from B\n",
    "b.g()     # calls the newly defined function, which would not work for a\n",
    "```\n",
    "with output \n",
    "\n",
    "```\n",
    "This is set in A\n",
    "This is set in B\n",
    "A new function\n",
    "\n",
    "```\n",
    "\n",
    "Thus we see that even though `f()` is not defined in `class B`, the latter **inherits** the definition from `class A`. Furthermore, the `g()` in `class B`, expands  the methods in `class A`\n",
    "\n",
    "* The child class can also *override* functions from their parent class:\n",
    "\n",
    "```py\n",
    "class A:\n",
    "    def f(self):\n",
    "        print(\"Parent function\")\n",
    "    \n",
    "\n",
    "class B(A):\n",
    "    def f(self):\n",
    "        print(\"Child function\")\n",
    "\n",
    "a = A()\n",
    "b = B()\n",
    "\n",
    "a.f()\n",
    "b.f()  \n",
    "\n",
    "A.f(b)        # This calls the function in the parent class for b\n",
    "\n",
    "```\n",
    "\n",
    "* The overwritten function can still be accessed from inside the child class using the `super()` function\n",
    "\n",
    "```py\n",
    "class A:\n",
    "    def f(self):\n",
    "        print(\"Parent function\")\n",
    "    \n",
    "\n",
    "class B(A):\n",
    "    def f(self):\n",
    "        super().f()    # first, A.f is called \n",
    "        print(\"Child function\")  # then new stuff is done\n",
    "\n",
    "\n",
    "b = B()\n",
    "b.f()\n",
    "```\n",
    "This is especially useful when initializing an object \n",
    "```py\n",
    "class A:\n",
    "    def __init__(self, text):\n",
    "        self.text = text\n",
    "    #...\n",
    "class B(A):\n",
    "    def __init__(self, text, moreText):\n",
    "        super().__init__( text)\n",
    "        self.moreText = moreText\n",
    "\n",
    "b = B(\"abc\", \"def\")\n",
    "print(b.text, b.moreText)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it yourself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why use inheritance?\n",
    "- Code Reusability\n",
    "- Enables Subtyping (IsA relation): Dog > Mammal > Animal\n",
    "- Logical Hierarchy\n",
    "- Use Polymorphism: Calling code is agnostic to the specific implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Dark Matter Halos\n",
    "In this example we will look at different analytic Dark Matter halo profiles. Since they all share the same idea of describing the density and mass profile of a halo, it makes sense to use classes and inheritance features. To this end, we define the abstract basis class DMHalo, which doesn't describe a specific halo, but tells us how what properties a halo class should have.\n",
    "\n",
    "The 3 profiles we will look at are\n",
    "1. a constant density profile\n",
    "\\begin{equation}\n",
    "\\rho_{static}(r) = \\rho_0\n",
    "\\end{equation}\n",
    "2. an NFW profile\n",
    "\\begin{equation}\n",
    "\\rho_{NFW}(r) = \\frac{\\rho_0}{\\frac{r}{r_s} (1 + \\frac{r}{r_s})^2} \n",
    "\\end{equation}\n",
    "3. a NFW profile with a central density spike\n",
    "\\begin{equation}\n",
    "\\rho_{SpikedNFW} = \\begin{cases}\n",
    "               \\rho_{spike} \\left(\\frac{r_{spike}}{r}\\right)^{\\alpha}, \\quad r <= r_{spike} \\\\\n",
    "               \\rho_{NFW}(r), \\quad  r > r_{spike}\n",
    "            \\end{cases}\n",
    "\\end{equation}\n",
    "with the parameters chosen such that $\\rho_{spike} = \\rho_{NFW}(r_{spike})$\n",
    "\n",
    "The constant and NFW profiles can inherit from this base class, and the SpikedNFW profile can inherit from the NFW and modify its density in the appropriate region. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# abstract base class for the DM halo\n",
    "class DMHalo:\n",
    "    def __init__(self, rho_0):\n",
    "        self.rho_0 = rho_0\n",
    "    \n",
    "    def density(self, r):\n",
    "        pass\n",
    "    \n",
    "    def __str__(self):\n",
    "        return \"DMHalo\"\n",
    "    \n",
    "    def mass(self, r):\n",
    "        pass\n",
    "    \n",
    "# Constant DM halo class, inherits from the abstract DMHalo\n",
    "class ConstHalo(DMHalo):\n",
    "    def __init__(self, rho_0):\n",
    "        super().__init__(rho_0)\n",
    "        \n",
    "    def density(self, r):\n",
    "        return self.rho_0 * np.ones(np.shape(r))\n",
    "    \n",
    "    def __str__(self):\n",
    "        return \"ConstHalo\"\n",
    "\n",
    "# NFW halo\n",
    "class NFW(DMHalo):\n",
    "    def __init__(self, rho_0, r_s):\n",
    "        super().__init__(rho_0)\n",
    "        self.r_s = r_s\n",
    "    \n",
    "    def density(self, r):\n",
    "        return self.rho_0 / (r/self.r_s) / (1. + r/self.r_s)**2\n",
    "    \n",
    "    def __str__(self):\n",
    "        return \"NFW\"\n",
    "    \n",
    "# Spiked NFW halo class    \n",
    "class SpikedNFW(NFW):\n",
    "    def __init__(self, rho_s, r_s, r_spike, alpha):\n",
    "        super().__init__(rho_s, r_s)\n",
    "        self.alpha = alpha\n",
    "        self.r_spike = r_spike\n",
    "        self.rho_spike = rho_s * r_s/r_spike / (1.+r_spike/r_s)**2\n",
    "    \n",
    "    def density(self, r):\n",
    "        return np.where( r < self.r_spike, \\\n",
    "                           self.rho_spike * (self.r_spike/r)**self.alpha, \\\n",
    "                        super().density(r))\n",
    "\n",
    "    def __str__(self):\n",
    "        return \"SpikedNFW\"\n",
    "\n",
    "# This is a function that plots the density for any object that has a density attribute\n",
    "def plotHalo(halo):\n",
    "    r = np.geomspace(1, 1e5)\n",
    "    l, = plt.loglog(r, halo.density(r), label=str(halo))\n",
    "    #plt.loglog(r, halo.mass(r), linestyle='--', color=l.get_c())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAELCAYAAADKjLEqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd1xW5f/H8dfFEFAEB4oDFVQcuAVFnGDlXqg58mdallnqN800zYblqqzMhpmmWWaaqThwl6AWONPceyRuzYUbuH5/HBwZKCL3fe775vN8PM4jz+G+z/2+xPhwnXOd61Jaa4QQQoiHcTI7gBBCCPsgBUMIIUSGSMEQQgiRIVIwhBBCZIgUDCGEEBkiBUMIIUSGSMEQQgiRIS5mB7AkHx8f7e/vn6n3XrlyhVy5cmVtIBsnbc4epM2O73Hbu2nTprNa6wL3H7ebgqGUKg58CZwF9mqtP3jYe/z9/dm4cWOmPi82Npbw8PBMvddeSZuzB2mz43vc9iqljqR13NRLUkqpKUqp00qp7fcdb6KU2qOU2q+UGpx6uAywSGv9PBBk9bBCCJHNmX0PYyrQ5N4DSiln4CugKUZh6KyUCgI2A52UUiuBGCvnFEKIbE+ZPZeUUsofiNZaV0zdDwOGaa0bp+4PSX3pLWC91nq1Umq21rp9OufrCfQE8PX1DZ45c2amciUmJuLp6Zmp99oraXP2IG12fI/b3oiIiE1a65D7j9viPYyiwNF79hOAUGACMEwp9QxwOL03a60nAhMBQkJCdGav42W3a54gbc4ubLnNt27dIiEhgevXr2fpeb29vXF3d8/Sc9qyjLbX3d0dPz8/XF1dM3ReWywYKo1jWmu9HUizVyGEcAwJCQnkzp0bf39/lErrR0HmXL58mdy5c2fZ+WxdRtqrtebcuXMkJCQQEBCQofOafQ8jLQlAsXv2/YDjJmURQljR9evXyZ8/f5YWC5E2pRT58+d/pN6cLRaMDUCgUipAKZUD6AQssGaAvef3cujGIWt+pBAilRQL63nUv2uzh9XOAOKBskqpBKVUD611EtAHWAbsAmZprXdYM9f4LeP59OSnvPPHO/xz/R9rfrQQQtgsUwuG1rqz1rqw1tpVa+2ntZ6cenyx1rqM1rqU1nqktXONrDuSJ7yeYOGBhbSIasGM3TNITkm2dgwhhElOnjxJp06dKFWqFEFBQTRr1oy9e/dmybnnzZvHzp077+x3796d2bNn/+s1GRnhZMaoL1u8JGW6XK65aJO3DbNbzSYoXxCj1o2i06JObDm9xexoQggL01oTGRlJeHg4Bw4cYOfOnYwaNYpTp05lyfnvLxj2xBZHSdmMUnlKManRJJYdWcaYDWPouqQrrUq1on9wf3w8fMyOJ4RDe2/hDnYev5Ql50pOTsbZ2ZmgIl6827LCA18bExODq6srvXr1unOsatWqaK0ZOHAgS5YsQSnFW2+9RceOHYmNjWXYsGH4+Piwfft2goOD+fHHH1FKMXjwYBYsWICLiwuNGjWibdu2LFiwgFWrVjFixAjmzJnzwCyJiYm0bt2a8+fPc+vWLUaMGEHr1q3/9RqtNYMGDfpXrmbNmmX+L+sBpGA8hFKKJv5NqF+0Pt9s/YYfdv5AzN8x9K7Wm45lO+LiJH+FQjiS2z/07zd37ly2bNnCX3/9xdmzZ6lRowb169cHYPPmzezYsYMiRYpQp04d/vjjD4KCgoiKimL37t0opbhw4QJ58uShVatWtGjRgvbt7z4lMHDgQEaMGPGfz3R3dycqKgovLy/Onj1LrVq1aNWq1b9uVqeVq3r16hYZRiw/7TIop2tO+gf3p3Xp1oxeN5oP1n/A3H1zGRo6lOq+1c2OJ4TDeVhP4FFkxXMYv//+O507d8bZ2RlfX18aNGjAhg0b8PLyombNmvj5+QFGb+Tw4cPUqlULd3d3XnjhBZo3b06LFi3SPfeYMWP+VUBu35/QWvPmm2+yevVqnJycOHbsGKdOnaJQoUIPzPXnn38SGBj4WO1Ni9zDeEQlvUsy8amJfNLgEy7dvES3pd14c82bnL121uxoQogsUKFCBTZt2vSf4w+aRsnNze3On52dnUlKSsLFxYX169fTrl075s2bR5MmTdJ9f3qmT5/OmTNn2LRpE1u2bMHX1/c/z01Yc3onKRiZoJSikX8j5reezwuVXmDJ4SW0jGrJtJ3TSEpJMjueEOIxNGzYkBs3bjBp0qQ7xzZs2EDevHn5+eefSU5O5syZM6xevZqaNWume57ExEQuXrxIs2bN+Oyzz9iyxRg0kzt3bi5fvpyhLBcvXqRgwYK4uroSExPDkSP/nXW8fv36/8mV1iW1rCAF4zHkdM3Jq9VfJapVFFUKVOGjDR/RIboDG09mbg0OIYT5lFJERUWxYsUKSpUqRYUKFRg2bBjPPPMMlStXpkqVKjRs2JCPPvroX5eG7nf58mVatGhB5cqVadCgAWPHjgWgU6dOjBkzhmrVqnHgwIEHZunSpQsbN24kJCSE6dOnU65cuf+8JjIy8j+5fH19H+8vIR2mz1ZrSSEhIdpaCyhprVn590o+3PAhJ66coHnJ5gwIHkCBnP9ZtMpm2fKkdJYibbYtu3btonz58ll+XplLKn1p/Z0rpdKcrVZ6GFlEKcUTJZ5gfpv5vFjpRZYfXk7LeS35fsf33Eq5ZXY8IYR4bFIwspiHiwf/q/4/olpHUa1gNT7e+DEdFnZgw8kNZkcTQojHIgXDQkp4lWD8E+MZFzGOa0nXeH7Z8wxaPYjTV0+bHU0IITJFCoYFKaVoWLwh81rPo1eVXvx25DdaRrVk6vapcplKCGF3pGBYgbuLO72r9mZe63nUKFSDTzZ9QvsF7Vl/Yr3Z0YQQIsOkYFhRMa9ifPnEl3zZ8EtuJt+kx/IeDFw1kFNXsmZSMyGEsCQpGCZoUKwBUa2jeKXKK8QcjaHlvJZM2T6FW8lymUoIsymlGDBgwJ39jz/+mGHDhgEwbNgwihYtStWqValatSqDBw9m/vz5tGnT5s7rR48eTenSpe/sL1y4kFatWlktvyVJwTCJu4s7L1d9mXmt5xFaOJSxm8bSbmE71p5Ya3Y0IbI1Nzc35s6dy9mzaU/3079/f7Zs2cKWLVv44IMPqF27NvHx8Xe+Hh8fj5eXF6dPGwNc4uLiqFOnjlWyW5pMPmgyv9x+fNHwC1YnrOaD9R/w4vIXaVSiEQNrDKRQrvSfIhXC4S0ZDCe3ZcmpPJKTwNkFClWCph888LUuLi707NmTsWPHMnLkw9dvK1CgAN7e3uzfv5/SpUtz7Ngx2rVrR1xcHG3atCEuLi7NmWjtkfQwbER9v/pEtY6iT9U+rE5YTat5rfh227dymUoIE/Tu3Zvp06dz8eLF/3xt7Nixdy5JLVu2DIDatWsTFxfHnj17CAwMpFatWsTFxZGUlMTWrVupUaOGtZtgEdLDsCFuzm68VOUlWpRqwUfrP2Lcn+OYv38+Q2oOoXbR2mbHE8K6HtITeBTXHnFqEC8vL5599lk+//xzPDw8/vW1/v378/rrr//rWJ06dYiLiyM5OZmwsDBq1qzJ+++/z+bNmylbtizu7u5Z0g6zSQ/DBhX1LMq4huP4+smvSdEpvPTrS/SP6c+JxBNmRxMi2+jXrx+TJ0/mypUrD33t7R5GXFwcYWFh5M6dm+vXrxMbG+sw9y9ACoZNq1u0LlGto+hbrS+/H/udVvNaMXHrRG4m3zQ7mhAOL1++fHTo0IHJkyc/9LVBQUEcP36cNWvWUK1aNcBYSGnChAnUru04VwfspmAopZyUUiOVUl8opbqZncdacjjnoGflnsxvM5+6RevyxeYviJwfyZqENWZHE8LhDRgwIN3RUvdSShEaGoqPjw+urq4AhIWFcfDgQYcqGKbew1BKTQFaAKe11hXvOd4EGAc4A99qrT8AWgNFgX+ABBPimqqIZxHGRowl7lgco9eP5pXfXiGiWARv1HyDop5FzY4nhMNITEy882dfX1+uXr16Z//28xhpWbRo0b/2u3fvTvfu3bM6nqnM7mFMBf61bqFSyhn4CmgKBAGdlVJBQFkgXmv9GvCylXPajNpFazOn1Rxerf4qa0+spfW81nz919fcSL5hdjQhhIMzfQElpZQ/EH27h6GUCgOGaa0bp+4PSX3pUeCm1nqWUupnrXXHdM7XE+gJ4OvrGzxz5sxM5UpMTLyzELutOp90nqjzUWy+uhkfFx/a5W1HxZwVH/7GdNhDm7OatNm2eHt7/+sp6aySnJyMs7Nzlp/XVj1Ke/fv3/+f4cMRERFpLqBki8Nqi2IUh9sSgFCMS1RfKKXqAavTe7PWeiIwEYwV9zK7spgtr0p2r0giiT8ez6h1o/jmzDeE+4UzqOYgiuUu9sjnspc2ZyVps23ZtWuXRVbGkxX30ufu7n7nRv3DmH1JKi0qjWNaa31Va91Da91Xa/2V1VPZsLAiYcxtNZf+wf1Zd3Idbea1YfyW8VxPum52NCGEA7HFgpEA3PvrsR9w3KQsdsPV2ZXnKz7PgjYLaFi8IV//9TVt5rch9mis2dGEEA7CFgvGBiBQKRWglMoBdAIWmJzJbhTKVYgxDcYwudFk3J3d6buyL71/683RS0cf/mYhhHgAUwuGUmoGEA+UVUolKKV6aK2TgD7AMmAXMEtrvcPMnPaoZuGa/NLqF14PeZ2NJzfSZn4bvtz8JdeSrpkdTQibN3LkSCpUqEDlypWpWrUq69atS/e1EyZM4IcffgAgPDycjRs3Zuozu3fvzuzZs++cJyTk7j3njRs33rnvFBsbi7e39535rJ588kkuXLhA/vz5uT2Iad26dSilSEgwnkC4ePEi+fLlIyUlJVPZbjP1prfWunM6xxcDi60cx+G4OrnSrUI3mgY05ZONn/DN1m+IPhjNoBqDiCgWgVJp3S4SInuLj48nOjqaP//8Ezc3N86ePcvNm+nPrtCrVy+L5Dh9+jRLliyhadOm//lavXr1iI6O/texQoUKsWvXLoKCgli/fj3VqlUjLi6ODh06sHbtWkJDQ3Fyerw+gi2OkhJZrGDOgnxY/0Pal2nPqHWjeDXmVeoWrcuQmkMo7lXc7HhCpOnD9R+y+5/dWXKu28NMy+Urxxs133jga0+cOIGPjw9ubm4A+Pj4AODv70/Hjh2JiYkB4KeffqJ06dIMGzYMT0/Pf01ImJKSwnPPPUexYsV47733GDx4MLGxsdy4cYPevXvz0ksvobWmb9++rFy5koCAAO5/xGHgwIGMGDEizYKRltsTIAYFBbFu3Tr69+9/p2DExcVlyRPntngPQ1hIjUI1mNVyFoNqDGLz6c20md+Gz//8XC5TCXGPRo0acfToUcqUKcMrr7zCqlWr7nzNy8uL9evX06dPH/r165fm+5OSkujSpQtlypRhxIgRTJ48GW9vbzZs2MCGDRuYNGkShw4dIioqij179rBt2zYmTZpEXFzcv84TFhaGm5vbnQJ1rzVr1ty5JHV7zY7bEyACHD58mKeffvrO5bGsWsRJehjZjKuTK12DutLEvwmfbvqUSdsm3blM5ayzz4NNwvY9rCfwKB7luQRPT082bdrEmjVriImJoWPHjnzwgTHVeufOne/8t3///mm+/6WXXqJDhw4MHToUgOXLl7N169Y79ycuXrzIvn37WL16NZ07d8bZ2ZkiRYrQsGHD/5zrrbfeYsSIEXz44Yf/Op7WJak6derwwQcfcOjQIYoXL467uztaaxITE9m0aRM1a9bMUPsfRHoY2VSBnAUYXW80U5tMxTOHJ/1j+zP+9HgOXzxsdjQhTOfs7Ex4eDjvvfceX375JXPmzAH4132/9O4B1q5dm5iYGK5fN56D0lrzxRdf3FnW9dChQzRq1OiB57itYcOGXL9+nbVrH750c2BgIOfPn2fhwoV3ikNwcDDfffcdAQEBWfJ0vxSMbC7YN5hZLWYxuOZgDt84TOSCSMb9OY6rt64+/M1COKA9e/awb9++O/tbtmyhRIkSAPz88893/hsWFpbm+3v06EGzZs14+umnSUpKonHjxnz99dfcumWsnrl3716uXLlC/fr1mTlzJsnJyZw4cSLNS08AQ4cO5aOPPspQ9rCwMMaNG3enYISFhfHZZ59l2Yy5cklK4OLkQpfyXch9PDfr3Nbx7bZviT4YzcCQgTxV4ikZTSWylcTERPr27cuFCxdwcXGhdOnSTJw4kejoaG7cuEFoaCgpKSnMmDEj3XO89tprXLx4ka5duzJ9+nQOHz5M9erV0VpToEAB5s2bR2RkJCtXrqRSpUqUKVOGBg0apHmuZs2aUaBAgQxlr1OnDosXL74z1UeWT7GutXbYLTg4WGdWTExMpt9rr263+c9Tf+r2C9rrilMr6heWvaAPXDhgbjALys7fZ1u0c+dOi5z30qVLj32OEiVK6DNnzmRBGst7lPam9XcObNRp/EyVS1LiP6oVrMaM5jMYUnMIO87uoN2CdozdNFYuUwmRzUnBEGlycXLhmfLPsDByIS1KtmDK9im0mteKpYeX/me8uBDZweHDh+88k5FdScEQD5TfIz/D6wxnWtNp5HPPx8BVA3lxxYscvHDQ7GjCQckvJNbzqH/XUjBEhlQtWJUZzWfwVuhb7Dy3k3YL2vHJxk+4cuuK2dGEA3F3d+fcuXNSNKxAa825c+dwd3fP8HtklJTIMGcnZzqW68hT/k8x7s9xTN0xlcUHFzMgZABNA5rKaCrx2Pz8/EhISODMmTNZet7r168/0g9Ge5fR9rq7u+Pn55fh80rBEI8sn3s+3qv9Hu0C2zFy3UjeWPMGs/fN5s2ab1I6b9YvrymyD1dXVwICArL8vLGxsRleVc4RWKq9cklKZFrlApX5qdlPvF3rbfae30v7he0Zs2EMiTcTzY4mhLAAKRjisTg7OdOhbAcWtllIZGAk03ZOo+W8lkQfjJbr0EI4GCkYIkvkdc/Lu2Hv8lPznyiUsxBD1gzhuWXPsff8XrOjCSGyiBQMkaUq+lRkevPpDAsbxoELB+iwsAMfrv+Qyzcvmx1NCPGYpGCILOeknGhXph3RkdG0C2zH9F3TaRnVkgUHFshlKiHsmBQMYTHebt68HfY2M1rMoKhnUYb+PpTuS7uz5589ZkcTQmSCFAxhcRXyV2Bas2m8X/t9Dl08RIfoDoxeN5pLNy+ZHU0I8QikYAircFJORAZGsjByIR3KdGDmnpm0jGrJ/P3zSdEpZscTQmSAXRUMpVQupdQmpVQLs7OIzPF282ZoraHMbD6TYrmL8dYfb9FtSTd2/7Pb7GhCiIcwtWAopaYopU4rpbbfd7yJUmqPUmq/UmrwPV96A5hl3ZTCEsrnL88PTX9geJ3h/H35bzpGd2Tk2pFcvHHR7GhCiHSY3cOYCjS594BSyhn4CmgKBAGdlVJBSqkngZ3AKYunOrsPj6snQEb0WJSTcqJN6TYsjFxIp7KdmLV3Fq3mtSJqX5RcphLCBplaMLTWq4F/7jtcE9ivtT6otb4JzARaAxFALeAZ4EWllOWyrx5D6PpeMLYiRL0MW2bAxWMW+7jsziuHF0NCh/Bzi58p4VWCd+LeoeuSruw8t9PsaEKIeyizx8UrpfyBaK11xdT99kATrfULqftdgVCtdZ/U/e7AWa11dDrn6wn0BPD19Q2eOXPmI2fyuHqcnCfX43ttD3nPb8M1yXjo7KpHES7kqcT5vJW4kKcSt3LkeeRz27LExEQ8PT1NzaC1Zv2V9cw/P5/ElETqeNahRZ4W5HLOZZHPs4U2W5u02fE9bnsjIiI2aa1D7j9ui7PVpjVH9p2qprWe+qA3a60nAhMBQkJCdHh4eKZCxMYWoVJ4OKSkwOkdcGg1OQ+tJufhPyhyYpnxooJBEFDf2ErUAQ/7LiCxsbFk9u8rK0UQwcs3X2b8lvHM2D2D7be20696PyIDI3HK4o6lrbTZmqTNjs9S7bXFgpEAFLtn3w84blIWcHKCQpWMLaw3JCfBib/g0Co4tBo2fQ/rJoBygsJVjOLhXx+K1wK37PMbTVbzyuHF4JqDiSwdyah1oxgWP4zZe2cztNZQKvpUNDueENmSLRaMDUCgUioAOAZ0wrhvYRucXcAv2NjqvQZJN+DYJji4Cg6vgfjx8Mc4cHKBoiF3eyB+NcA1+yzgklXK5ivL1CZTiT4YzaebPuWZRc/Qrkw7Xq32Knnc7btHJ4S9MbVgKKVmAOGAj1IqAXhXaz1ZKdUHWAY4A1O01jtMjPlgLm5QoraxMQRuXoWja43ex6HVsOZjWP0RuLhDsVAIqAcBDaBINXB2NTu9XVBK0bJUSyKKRTD+r/H8tOsnVhxZwf+q/Y92ge1wdnI2O6IQ2YKpBUNr3Tmd44uBxVaOkzVy5IRSDY0N4PpFOBJ/t4CsHAGMgByeRpG53QPxrWRc/hLp8szhyaAag4gsHcnIdSMZvnY4c/bN4a3Qt6hUoJLZ8YRweLZ4ScqxuHtD2SbGBnDlnHHp6nYB2bfcOO6RF/zrGr2PgPrgUwZkjew0BeYN5LvG37H40GI+2fgJXRZ3oW1gW16t/ip53fOaHU8IhyUFw9py5YcKbYwN4NJxOHS7gKyCXQuN456+d3sfAQ0gbwnzMtsgpRTNSzangV8DJvw1gem7pt+5TNW+THu5TCWEBUjBMJtXEajS0di0hvOH7/Y+Dq6Cbb8Yr8tT4m7xCKgHuQuZGttWeObw5PUarxMZaIymGrFuBHP2zWForaFUKVDF7HhCOBQpGLZEKcgXYGzB3YwCcmbPPb2PBbB5mvFan7JQssHdZ0By5jM3u8lK5SnFt42+ZdnhZYzZOIb/W/x/tCndhn7V+5HfI7/Z8YRwCFIwbJlSULCcsYX2hJRkOLn1bg9k84+wfiKgoHDl1N5HAygRBjks82S0LVNK0SSgCfX86vHNX98wbec0fvv7N/pW60uHMh3kMpUQj0kKhj1xcjaG4xapBnVehaSbxjMgh1YZ90HWTYC4z8HJFfxC7t5A9wsxhv9mE7lcc/FayGu0Kd2GUetHMWrdKObum8vQ0KFULVjV7HhC2C0pGPbMJYfRmygRBuGD//0MyMFVxvMfqz4AFw/jNbfvgRSuYhQfB1cyT0kmPTWJ5UeWM2bDGLou6Urzks15oeILlM5b2ux4QtgdKRiO5P5nQK5dgCN/3C0gvw4zjrt7g3/qA4QlGzj0EF6lFI39G1OvaD0mbp3IT7t/YtHBRTTwa8DzFZ83O54QdkUKhiPzyAPlmhsbwOVTxjMgB2ONy1i7Uyf89SwEJRtQ6IYvXCwN3n6mRbaUnK456Rfcj24VujFz90x+2v0T3ZZ2I8AtgJS/UwgvFp7lExsK4WikYGQnuX2hUntjA/jn0N1JFA+spNyVM7Dnc8hXKnUEVuo9EAcagZXXPS8vV32ZbhW6MW//PL7Z9A2vxrxKCa8S1Ctaj1qFaxFSKIRcrtlv0IAQDyMFIzu7M4S3O2jNhkXfUyP/VaOIbP0FNk4BlDFTb8lwo4gUr21c+rJzOV1z8kz5Z/A96cuNEjeYu38uv+z9hR93/YiLcqGiT0VCC4dSq3AtqhSogqvM+yWEFAyRSimuePpDWDiEvQLJt+D4ZuPex8HYuyOwnHOAX83UAhKeOomi/f4zclbONCvZjGYlm3Ej+QZbTm9h3Yl1rD2xlknbJvHN1m/wcPGgum91wgqHEVo4lDJ5y8jlK5Et2e//6cKynF2hWE1jazDQGIH1d/zd+x8xIyFmBLh5GXNglQw3Nju+ge7m7EZo4VBCC4fyP/7HpZuX2HByA+tOrGPdiXV8vPFjAPK65aVm4ZrUKlyLsCJhFPUsanJyIaxDCobImBw5ofQTxgapkyiuNgrIwVWwJ3Vy4dxF7haPkg3segoTrxxePFH8CZ4obrT51JVTrD+5nrUn1rL2+FqWHTZWXiyWuxi1CteiVuFahBYOxdvN28zYQliMFAyRObnyQ4VIYwNjDqyDsca2dyn89ZNxvEB5KBVhFJASdex6FULfXL60LNWSlqVaorXm0MVDxJ+IZ+3xtSw+tJhf9v6CQhGUP+hO76NawWrkcM5hdnQhsoQUDJE18vobN8+DuxvroJ/cereAbJwCa8cbqxD61TQKSKmGxv0PO32AUClFyTwlKZmnJF3Kd+FWyi12nN1B/PF41p5Yy/c7vmfy9sm4O7sT7Bt8p4CUyVsGZaeX7ISQgiGynpMTFKlqbHX7wa1r8Pfa1AISAzGjjHsg7t7G0N3bBSSvv9nJM83VyZWqBatStWBVXq76MlduXWHDyQ2sPbGW+OPxfLLpE9gE+d3zE1YkzNgKh1EgZwGzowuRYVIwhOW5eqQWhQjgPbhy9m7xOBBjzMILkK/k3SfV/euBu5eZqR9LLtdchBcLJ7xYOAAnr5wk/ng88SfiiTseR/RB46HJwLyB1C5cm9pFalPdtzruLrLuu7BdUjCE9eXyufsAodZwdh8cWGlsW36CDd/evXxVOrWAFK5m10vYFspViMjASCIDI0nRKez5Zw9xx+OIPx7PT7t/4vud3+Pm7EawbzC1ixgFpHSe0nL5StgUKRjCXEpBgTLGVqsXJN2Ao+vuFpCVI4wtZ/7U3scTxn9z+5qd/F+01gAZ+gHvpJwon7885fOXp0elHly9dZVNpzYRdzyOuONxd4bvFsxZkDpF6lCnaB1qFa4lo6+E6aRgCNvi4nZ3adonh0HiGePS1f7f4MBvd1cgLFQZSj9pbMVqGs+NmEgpdadoPKqcrjmp51ePen71ADiReIL4E/H8fux3fv37V6L2R+GknKjsU5k6RetQt2hdgvIHycODwuqkYAjb5lkAKncwtpQUOLUN9v9qFJC4z+H3T42HB0uGQ2Ajo4B4FTYlalZdPirsWZi2gW1pG9iWpJQktp/dzu/HfuePY38wfst4vtryFfnc81G7SG3qFa1H7SK1yeOeJ0s+W4gHsZuCoZRqAzQHCgJfaa2XmxxJWJuTk7GWR+EqUG8AXL9k3DzfvwL2/Xr35rlvJQh8Cso0Br8adjt0F8DFyeXO6Ks+1frwz/V/iDsed6eARB+Mxkk5UcmnEnWL1qWeXz3K5ysvvQ9hEaYWDKXUFKAFcFprXfGe402AcYAz8K3W+gOt9TxgnlIqL/AxIMzLgg4AAB8HSURBVAUju3P3gqBWxqY1nN4J+1YY2+3eh0deo9cR2Nh4St3OZ97N556PFiVb0KJkC5JTktlxbgdrjq3h94Tf+WrLV3y15SsKeBSgnl896vvVJ6xwGDld7X+ySGEbzO5hTAW+BH64fUAp5Qx8BTwFJAAblFILtNY7U1/yVurXhbhLKfCtYGx1+8H1i8ZN873LYd9y496HcjJGXpVtamw+ZcxO/VicnZypXKAylQtUpnfV3py7do7fj/3O6oTVLD+8nLn75uLq5EqNQjWo71efBn4N8MvteGudCOtRmb1Rl2UBlPIHom/3MJRSYcAwrXXj1P0hqS/9IHVbobX+9QHn6wn0BPD19Q2eOXNmpnIlJibi6Wm/01hkhsO2WaeQ+/J+8p/bQP5zG8mdeBCAa+6FOOFdjUuFanPROwjtZPbvT1knWSdz4MYBdlzbwfar2zmddBqAIq5FKOtSlure1Smeo3i2uXTlsP+20/G47Y2IiNiktQ65/7gtFoz2QBOt9Qup+12BUGAv0A3YAGzRWk942LlDQkL0xo0bM5UrNjaW8PDwTL3XXmWbNl9MMOa72rOUlAOxOOlbxlPngY2N1QlLP2nXc16l5cilI8QejSX2aCx/nvqTFFLw8fChgV8DGhZvSGjhUNyc3cyOaTHZ5t92qsdtr1IqzYJhi79SpTXURGutPwc+t3YY4YC8/aDGC1DjBf74dQn1iibB7sVGEdk2C5zdjKfSyzWHMk2NkVp2roRXCbpV6Ea3Ct1Y9NsidIAm9mgsSw8vZc6+OeR0MYb2PlH8CeoWrUvuHLnNjixskC0WjASg2D37fsBxk7IIB5fs4gHlw6F8S0hOgqNrYfciY73zvUsBBSVqQ/lWxmu87X/ti1zOuQgvGU6Lki24mXyTdSfW8dvfvxFzNIZlh5fh4uRCaOFQnir+FA2LNySve16zIwsbYYsFYwMQqJQKAI4BnYBnzI0ksgVnF2MxKP+60HgUnNoOu6KN4bpL3zA2vxoQ1NooIHlLmJ34seVwznHnocG3U95m69mt/HbkN377+zeGxQ9j+NrhhBQKoVGJRjQs3hAfDx+zIwsTmT2sdgYQDvgopRKAd7XWk5VSfYBlGMNqp2itd5gYU2RHKnUt80KVIGKIMd/VzvnGtvwtYytSDSq2M9YE8bb/0UfOTs5UK1iNagWrMSBkAHvO72H54eWsOLKC4WuHM2LtCIJ9g2nk34hGJRqR3yO/2ZGFlZlaMLTWndM5vhhYbOU4QqTPJxDqv25s/xwyeh07ou4Wj2KhUKEtVGhj16sM3qaUoly+cpTLV46+1fqy/8J+VhxZwfLDyxm1bhQfrP+A0EKhNA1oSsPiDWWeq2zCFi9JCWHb8gVAnVeN7dwBo3Bsn5t62WqwcUmrckfjgUJ3+/9BqpQiMG8ggXkDeaXqK+w7v48lh5aw5NAS3ol7h/fXvk/donVpFtCM8GLheLh4mB1ZWEiGCoZSyglw0lonWTiPEPYlf6m7PY8ze2D7HOMhwQV9YNEA4wHByh2NoboujrFU6+3i0bdaX3ac28GSQ0tYengpsUdjyeWaiyeLP0nLUi2pUahGtnnOI7t4aMFIvZ/wLnBTKXUW+EJr/a3FkwlhbwqUhYg3IXwIHNsEW382eh475xlTlFRoC9W6QJHqxj0SO6eUoqJPRSr6VGRAyAA2ntzIwoMLWXFkBfMPzMc3py/NSzanZcmWlM5b2uy4IgtkpIcxAKiktT6plCoMjFJK+Wmth1k2mhB2SinwCzG2xqOMVQW3zoQt02HjZChQ3igclTuCZ0Gz02YJJ+VEzcI1qVm4Jm+Gvsmqo6tYeHAh3+/4ninbp1DJpxJtSrehaUBTecbDjmWkYCQCpwG01ieUUj2ALcAwC+YSwjE4u0KZRsZ2/aLR49j8o3Gj/NdhxpTs1f7PeMrc2TFuKXq4eNAkoAlNAppw7to5Fh1cRNT+KIavHc6YDWN4qsRTRAZGEuIbIisK2pmM/Av9GvhFKfWG1no/UBy4atlYQjggd28Iec7YTu82ehx/zYQ9iyF3YajWFao/C3mKPfxcdiK/R36erfAsXYO6suPcDubum8uSQ0tYeHAhxXIXo32Z9rQp3YZ87vY9i3B28dA7Ulrr8cB04Ful1D/AfmCPUupppVSgpQMK4ZAKloNGw+G1ndDpJ/CtCKvHwLjKML0D7FkCKclmp8wyt+93vBP2Dis7rGRU3VEUzFmQsZvG8uQvT/LG6jfYdGpTplctFNaRoT6w1nouMFcp5QIEAdWA2sDLQEPLxRPCwTm7GnNWlWsO54/Anz/A5mkwoxN4+Rm9keDukMtxnrD2cPGgZamWtCzVkgMXDvDL3l9YsH8Biw8tpnSe0jxd5mlalWqFZw7HmgDSETzSmDetdZLWeqvW+nutdX+ttRQLIbJK3hLwxNvQfwd0mAY+pWHlcPg0COb1hhN/mZ0wy5XKU4rBNQfz69O/8n7t93F3dmf0+tE88csTjF43miOXjpgdUdzDMe6yCeFInF3vriR4ejesnwh/zYAtP0LxMAh9Ccq1dJib5AA5XXMSGRhJZGAkO87uYPqu6czaO4sZu2dQz68eXcp3IaxwmNwkN5k8VSOELStYDlp8Cq/tgkYj4fIJ+KU7fF4V4sfDjctmJ8xyFXwqMKreKFa0X0GvKr3YfnY7L614icj5kczZO4ebyTfNjphtScEQwh545IHafaDvn9BpBngXg2VDYGwF+PU9uHzS7IRZzsfDh1eqvsKK9isYUWcEOZxzMCx+GI3nNGbytslcunnJ7IjZjukr7llSZlfce2/hDuJ2/k2ePHkskMp2XbhwQdpsR0rf3E3LK7OpeT2OZJxZ4xHBQs+nOe7y4GG59tpmjeaK027OOS/jivNOnLQ7eZPrky/pCVx58Jod9trmzPJKucSklxtn+v32tOKeECID9ucox9gcb+GbdJzmV+YScXU54ddWsM69LlGenTniWtLsiFlKofBMKY9nSnmuJx3lrMsyzjn/yjnn38iTXAuf5Kbk0I7x5LytkoKRhndbViA29xnCw8PMjmJVxjrA0mb71A6unIW14wlbN5Gws2ugbDOo9zr4Bf/rlY7R5jCgAwmXE/hh5w/M3TeXgynxNA9ozouVXyTAO+Bfr3aMNmdcbGysRc4r9zCEcBS5fOCJd6D/NogYCkfi4NuGMC0SjsSbnc4i/HL78Wbomyxtt5Rng57l179/pfW81gxaPYj95/ebHc/hSMEQwtF45IUGg6D/dnhyGJzYCt81MQpHwiaz01mEj4cPA0IGsLTdUp6v+Dyrjq6i7YK2DIgdwIELB8yO5zCkYAjhqNxyQ93+0G8bPDUcjm+BbxtScdtIOLnN7HQWkc89H/2C+7Gs3TJerPwifxz/g7YL2jLt7DSOXj5qdjy7JwVDCEeXIyfU+R/02woRb5Hnwg6YUNd4nuPMHrPTWUQe9zz0rdaXpW2X0i2oG5uvbqZVVCuGxw/n1JVTZsezW3LTW4jswi03NBjI2pvlqeu0GdZ+DTvnQ9UuxsJPXkXMTpjl8rjn4bWQ1yh1oRTbPLcxZ+8c5h+YT6eynXih0gvkcc8+Q22zgvQwhMhmklw9jTmr+m2F0F7GFOufV4ff3jfW7HBA3i7evFXrLRZELqCxf2Om7ZpGs7nN+G77d9xIvmF2PLthNwVDKZVLKfW9UmqSUqqL2XmEsHu5fKDJaOizwZgtd80n8Hk1WPcNJDnm9BvFchdjZN2RzG45myoFq/Dppk9pFdWKRQcXkaJTzI5n80wtGEqpKUqp00qp7fcdb6KU2qOU2q+UGpx6uC0wW2v9ItDK6mGFcFT5AqD9ZOgZCwWDYMkg+KqmcbnKQWeCCMwbyNdPfs2kRpPwcvNi8JrBPLPoGTac3GB2NJtmdg9jKtDk3gNKKWfgK6ApxtobnZVSQYAfcHuYg+OsLCOErShSDbothC6zwdUDZj0L37c0huU6qFqFa/Fzi58ZWXckZ6+d5fllz/Pqylc5eklGVKXF9LmklFL+QLTWumLqfhgwTGvdOHV/SOpLE4DzWutopdRMrXWndM7XE+gJ4OvrGzxz5sxM5UpMTMTTM3st4CJtzh4y0maVkkzhE8sJODQdl6REThRuxKGAZ7iVwz5vEmekzTdTbhJ7OZZlF5eRolOI8IqgkXcj3J3crZQy6zzuv+uIiIg055JCa23qBvgD2+/Zbw98e89+V+BLIBfwHcYa410ycu7g4GCdWTExMZl+r72SNmcPj9Tmq/9ovWSw1u/l03qUn9Z/fK71rRsWy2Ypj9LmU1dO6SGrh+iKUyvqiJ8j9Pz983VySrLlwlnA4/67BjbqNH6mmn1JKi1prZCitdZXtNbPaa1f1lpPt3oqIbIjj7zGjfGX46FYKCx/C76uDQdizE5mMQVzFmRUvVH82OxHfHP6MvT3oXRd3JVtZxzzYcdHYYsFIwG4d35mP+C4SVmEEAAFysD/zYZnfoGUJJjWBmZ1g4vHzE5mMVUKVGF68+kMrzOcY4nH6LK4C+/Fv8fFG4459DgjbLFgbAAClVIBSqkcQCdggcmZhBAAZRrBK2uNyQ33LoUva8Af4xx2GK6TcqJN6TZER0bzf0H/R9S+KFpGtSRqX1S2HIZr9rDaGUA8UFYplaCU6qG1TgL6AMuAXcAsrfUOM3MKIe7h6m5Mbth7HZRsACveMaYaObTa7GQW45nDk0E1BvFzi5/x9/bnnbh36LakG3v+ccypVdJjasHQWnfWWhfWWrtqrf201pNTjy/WWpfRWpfSWo80M6MQIh15/aHzDOj8MyTfMIbgRvWCK+fMTmYxZfOVZWqTqQyvM5wjl47QMbojH234iKu3rpodzSps8ZKUEMKelG1iXKaq9zps+wW+DIHN0x32ob/bl6kWRi6kXWA7ftz5I63nt2bV0VVmR7M4KRhCiMfn6mHMT9Xrd/AJhPmvGD2Os467iJG3mzdvh73ND01/wNPVkz4r+zAgdgBnrp4xO5rFSMEQQmSdguXhuaXQ4jPjCfGva8Oqjxz2pjhA1YJVmdViFn2r9SX2aCyt57Xml72/OORNcSkYQois5eQEIc8ZkxqWbwExI2FiAzjmmKv9Abg6u9Kzck/mtJpD+fzleT/+fZ5b+hwHLx40O1qWkoIhhLCM3L7QfopxU/zaBfj2SVg2FG467g1if29/vm30LcPrDGf/hf08veBpJm+bTFJKktnRsoQUDCGEZZVtAr3XQvVuEP+lcZnq0BqzU1mMUoo2pdswv8186vvV57M/P6PL4i4OMQRXCoYQwvLcvaHlZ9At2tj/vgUsfBWuXzI3lwX5ePgwNmIsnzT4hJNXTtIpuhNfbfmKW8m3zI6WaVIwhBDWE1APXo6D2n3hzx9gfBgcWGl2Kotq5N+I+a3n0ySgCRP+mkCH6A7sOGufzyJLwRBCWFeOnNBoBDy/3BiOOy0SFvaDG5fNTmYxedzzMLreaL564isu3bxEl8Vd+GLzF3bX25CCIYQwR7Ea0GuN0dvYNBXG14aDsWansqj6fvWJah1F85LNmbh1Ip0WdWL3P7vNjpVhUjCEEOZx9UjtbSwDZ1f4oTVEvwY3Es1OZjFeObwYWXckn0d8zj/X/6FzdGe+/utrbqXYfm9DCoYQwnzFQ42nxGv1ho1TYEId+Hut2aksKqJ4BFGtomjk34jxW8bTZVEX9p3fZ3asB5KCIYSwDTlyQpNR8NxiYx6qKU1gxbuQdMPsZBaTxz0PH9b/kE/DP+XU1VN0jO7I9zu+t9mnxKVgCCFsS4na8PIfUP1Z+OMzmBgBJx17tbunSjzF3FZzqVu0Lh9v/Jgey3pwLNH2FqeSgiGEsD1uuaHV5/DMLLhyxigaaz6FlGSzk1lMfo/8jIsYx/u132fXP7tot6Ad8/bPQ9vQrL9SMIQQtqtMY2Pq9HLN4Lf34LtmcP6w2aksRilFZGAkc1rNoVy+crz9x9v0i+nHuWu2scaIFAwhhG3LlR+e/h7aToLTu+DrurDlJ4ddbwOgqGdRpjSewushr7Pm2BraLmjL6gTzVzSUgiGEsH1KQeUO8PLvULgKzHsZfukGV/8xO5nFOCknulXoxswWM8nvkZ/ev/VmxNoRXEu6Zl4m0z5ZCCEeVZ7i0G0BPPke7F5sTC2y/zezU1lUmbxlmNF8Bs8GPcvPe36mY3RHdp7baUoWKRhCCPvi5Ax1+8GLvxmTGv7YFpYMhlvXzU5mMW7ObgysMZBJjSZx5dYVuizqwrfbviXZyoMApGAIIexT4Srw0iqo+RKs+xomNYRT5vzmbS21Ctdibqu5NCzekHF/jqPH8h6cSDxhtc+3q4KhlGqjlJqklJqvlGpkdh4hhMlcPaDZR9BlNlw5DZMiYN1Eh74h7u3mzccNPmZEnRHsOreLdgvbsezwMqt8ttUKhlJqilLqtFJq+33Hmyil9iil9iulBj/oHFrreVrrF4HuQEcLxhVC2JPAp4xp0/3rwZKBMKMTJJ4xO5XFKKVoXbo1s1vOxt/Ln9dXvc47f7zD1VuWXc3Qmj2MqUCTew8opZyBr4CmQBDQWSkVpJSqpJSKvm8reM9b30p9nxBCGDwLQpdfoOlHcCDGWNlv369mp7KoYl7F+L7p97xY6UXm7Z9n8bU2lDWfIlRK+QPRWuuKqfthwDCtdePU/SEAWuvR6bxfAR8AK7TWaf5LUEr1BHoC+Pr6Bs+cOTNTWRMTE/H09MzUe+2VtDl7yA5tzpV4mKCdn5Dr6t8c9WvNtoKR5PTKa3Ysi9p3fR8/nP2BS8mXaJSzEU19muKkMtcniIiI2KS1DvnPF7TWVtsAf2D7PfvtgW/v2e8KfPmA9/8P2ARMAHo97POCg4N1ZsXExGT6vfZK2pw9ZJs237yq9aLXtX7XS18aU1XrM/vMTmRxF65f0P1j+uuKUyvqmL9jMn0eYKNO42eq2Te9VRrH0u3yaK0/11oHa617aa0nWDCXEMLeuXpAszHQ6Sfcr5+Bb+rD5ukOf0P8kwaf0Ltgbxr4Ncjy85tdMBKAYvfs+wHHTcoihHBE5ZqzocY4KFod5r8Cc3rA9Ytmp7IYpRTlPMphXMHPWmYXjA1AoFIqQCmVA+gELDA5kxDCwdx0yw/PzoeGb8OOeTChLhzdYHYsu2PNYbUzgHigrFIqQSnVQ2udBPQBlgG7gFlaa8vd4hdCZF9OzlD/dXh+qbH/XRP4fSyk2OZiRbbIxVofpLXunM7xxcBia+UQQmRzxWrCS2tg4avw6zA4uAraTjSG5YoHMvuSlBBCWJ9HHnh6KrT4DP6Oh6/rwIGVZqeyeVIwhBDZk1IQ8hy8GAM588O0SGMN8eRbZiezWVIwhBDZm28QvLgSgrsba4h/1xQu/G12KpskBUMIIXLkhJbjoP13cGaPMYpqV7TZqWyOFAwhhLitYltjyvS8AfBzF1g8CJJumJ3KZkjBEEKIe+UrCT2WQ61XYP03MPkpOHfA7FQ2QQqGEELcz8UNmoyGTjPg/BFjWpFts81OZTopGEIIkZ5yzaDX7+BbwZhSZMH/4NY1s1OZRgqGEEI8SJ5i0H0R1O0Pf34Pk56AM3vNTmUKKRhCCPEwzq7w5DDoMgcST8LEcPgrc2vt2DMpGEIIkVGBTxqXqIpUhaiXYF5vuGnZZVFtiRQMIYR4FF5F4NkFUH8gbJkOkxrC6d1mp7IKKRhCCPGonF2g4VvQdS5cPQuTImDLDLNTWZwUDCGEyKxSDY2Zb4tUh3m9YL5jX6KSgiGEEI/Dq7CxOFP9QcYSsJMaGtOLOCApGEII8bicXaDhUOMS1ZUzDjuKSgqGEEJklVINU0dRVTdGUc3v41AP+knBEEKIrHT7ElW912HzNONBv7P7zE6VJaRgCCFEVnN2gSfehv+750E/B5iLSgqGEEJYSuknjVFUhSoZc1FF94db181OlWlSMIQQwpK8i0K3hVDnVdg4xa6nS7ergqGUyqWU2qSUamF2FiGEyDBnV3jqfej8s7H868Rw2Dnf7FSPzCoFQyk1RSl1Wim1/b7jTZRSe5RS+5VSgzNwqjeAWZZJKYQQFla2CfRaAz5lYNazsGQwJN00O1WGWauHMRVocu8BpZQz8BXQFAgCOiulgpRSlZRS0fdtBZVSTwI7gVNWyiyEEFkvT3F4bomxot+6r+G7Jkavww4orbV1PkgpfyBaa10xdT8MGKa1bpy6PwRAaz06nfePBHJhFJdrQKTWOiWN1/UEegL4+voGz5yZuYdnEhMT8fT0zNR77ZW0OXuQNtsOnzNxlNv9BVo5sbtcP8751MiS8z5ueyMiIjZprUP+8wWttVU2wB/Yfs9+e+Dbe/a7Al9m4DzdgRYZ+czg4GCdWTExMZl+r72SNmcP0mYbc+6A1l/X1fpdL62Xv6110q3HPuXjthfYqNP4mWrmTW+VxrGHdne01lO11tEWyCOEENaXryT0WAHBz8Ef4+D7lnDpuNmp0mRmwUgAit2z7wfY5t+SEEJYkqs7tPwM2k6CE3/BhHpwIMbsVP9hZsHYAAQqpQKUUjmATsACE/MIIYS5KneAnjGQywemRULsh5CSbHaqO6w1rHYGEA+UVUolKKV6aK2TgD7AMmAXMEtrvcMaeYQQwmYVKAsvroTKHSF2FExvD1fOmp0KABdrfIjWunM6xxcDi62RQQgh7EaOXBA5AUqEweJBxiWqp7+D4rVMjWVXT3oLIUS2oRQEd4cXfgUXN/iuGcR9AVZ6FCItUjCEEMKWFa4ML62Csk1h+Vvw8//BtQumRJGCIYQQts7dGzr+CI1Hwd6lMLEBHN9i9RhSMIQQwh4oBWG9oftiSL4FkxsZs99a8RKVFAwhhLAnxUONNTb86xrra0S9BDevWOWjpWAIIYS9yZUfusyGiKGwdRZMaghn9lj8Y6VgCCGEPXJyggaDoGuU8ZzGxAjY+otlP9KiZxdCCGFZpSKMNTYKV4a5L0D0a6iUWxb5KCkYQghh77yKGMvA1u4LGydTbfNgOH8kyz9GCoYQQjgCZ1doNAI6/YTbjXNw43KWf4RVpgYRQghhJeWasy7UlfqFKmb5qaWHIYQQDibFOYdFzisFQwghRIZIwRBCCJEhUjCEEEJkiBQMIYQQGSIFQwghRIZIwRBCCJEhUjCEEEJkiNImLvdnaUqpM8ARwBu4eM+X7t1P788+QFasvH7/Zz/Oa9P7+oPal9ax7NDmjH7P7aXNGTlmr23O6Pc4rWPS5rTb/LjtLaG1LvCfo1prh9+AientP+DPGy3x2Y/z2vS+/qD2Zdc2P8L33C7anJFj9trmjH6Ppc0Zb3NWtff+Lbtcklr4gP30/mypz36c16b39Qe1L61j2aHNGf2eZxVLtzkjx+y1zRn9Hqd1TNps+Tbf4dCXpB6HUmqj1jrE7BzWJG3OHqTNjs9S7c0uPYzMmGh2ABNIm7MHabPjs0h7pYchhBAiQ6SHIYQQIkOkYAghhMgQKRhCCCEyRAqGEEKIDJGCkUFKqVxKqe+VUpOUUl3MzmMNSqmSSqnJSqnZZmexFqVUm9Tv8XylVCOz81iaUqq8UmqCUmq2Uupls/NYS+r/z5uUUi3MzmINSqlwpdSa1O91eGbPk60LhlJqilLqtFJq+33Hmyil9iil9iulBqcebgvM1lq/CLSyetgs8iht1lof1Fr3MCdp1nnENs9L/R53BzqaEPexPWJ7d2mtewEdALt9TuER/18GeAOYZd2UWesR26yBRMAdSMj0h1ri8XF72YD6QHVg+z3HnIEDQEkgB/AXEAQMAaqmvuYns7Nbo833fH222blNaPMnQHWzs1ujvRi/AMUBz5id3RptBp4EOmH8UtDC7OxWarNT6td9gemZ/cxs3cPQWq8G/rnvcE1gvzZ+u74JzARaY1Rlv9TX2O3f2yO22SE8SpuV4UNgidb6T2tnzQqP+j3WWi/QWtcG7PZS6yO2OQKoBTwDvKiUssv/nx+lzVrrlNSvnwfcMvuZLpl9owMrChy9Zz8BCAU+B75USjXHwvO1mCDNNiul8gMjgWpKqSFa69GmpLOM9L7PfTF+A/VWSpXWWk8wI5wFpPc9Dse43OoGLDYhlyWl2WatdR8ApVR34Ow9P0wdQXrf57ZAYyAP8GVmTy4F479UGse01voK8Jy1w1hJem0+B/SydhgrSa/Nn2P8cuBo0mtvLBBr3ShWk2ab7/xB66nWi2I16X2f5wJzH/fkdtkVs7AEoNg9+37AcZOyWIu02fHbnN3aC9JmyOI2S8H4rw1AoFIqQCmVA+Pm2AKTM1matNnx25zd2gvS5ixvc7YuGEqpGUA8UFYplaCU6qG1TgL6AMuAXcAsrfUOM3NmJWmz47c5u7UXpM3WarPMViuEECJDsnUPQwghRMZJwRBCCJEhUjCEEEJkiBQMIYQQGSIFQwghRIZIwRBCCJEhUjCEEEJkiBQMIUyUOjuu/H8o7IL8QxXCypRS/kqpXUqp8cCf/HvuHyFsljzpLYSVKaX8gYNAba31WnPTCJFx0sMQwhxHpFgIeyMFQwhzXDE7gBCPSgqGEEKIDJGCIYQQIkPkprcQQogMkR6GEEKIDJGCIYQQIkOkYAghhMgQKRhCCCEyRAqGEEKIDJGCIYQQIkOkYAghhMiQ/wcbJnjNB1GmeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotHalo(ConstHalo(5.))\n",
    "plotHalo(NFW(10, 1e3))\n",
    "plotHalo(SpikedNFW(10., 1e4, 1e2, 7./3.))\n",
    "plt.xlabel(\"r\"); plt.ylabel(r\"$\\rho$\")\n",
    "plt.legend(); plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2\n",
    "\n",
    "Implement a mass function inside the classes and plot it alongside the density profile. Recall the mass can be computed as the volume integral of the DM density.\n",
    "\n",
    "There are two possible approaches to this:\n",
    "\n",
    "1. Implement the analytical mass function for each class individually\n",
    "2. Implement a numerical integration scheme in the base class and let inheritance do the rest  (tip: use  solve_ivp)\n",
    "\n",
    "* What are the pros and cons of each approach?\n",
    "\n",
    "*  Implement both and compare the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3\n",
    "\n",
    "In this exercise we will automatize  Exercise 3 in Lecture 5 for all data present in the `Rotation curves` directory, using classes. As usual, we  will divide the exercise in several steps:\n",
    "\n",
    "\n",
    "1. Load the pandas, and pyplot libraries; they will be used in your code \n",
    "2. Define a new class `rotation_curve`, and initialize it. The initialization function should take the input  `Galaxy_ID`. \n",
    "\n",
    "3. In the `__init__` function, define the instance variable `df_circular`, which creates a  **DataFrame** for the given `Galaxy_ID`. (This corresponds to step 2 in Exercise 3 Tutorial 5.)  Create a test class instance of `rotation_curve`, for a given galaxy ID. Check that for instance, for `Galaxy_ID = 'IC2574'` the following lines\n",
    "```py\n",
    "test = rotation_curve('IC2574')\n",
    "test.df_circular.head()\n",
    "```\n",
    "are produced as output \n",
    "\n",
    "![](./Figures/circ_velocity.png \"Circular velocity\")\n",
    "\n",
    "\n",
    "4.  Define the class method `update_df_circular` which takes as input the variables `key` and `data`. The method adds a new column  `data` to the **DataFrame** `df_circular`, with keyword `key`. \n",
    "\n",
    "5. Define the class method `add_M_DM`, which will add the column `DM Mass (M_sol)` to   the **DataFrame** `df_circular`. (This corresponds to step 3 in Exercise 3, Tutorial 5. Hint, use the `update_df_circular` defined above.) Check that when you run the lines\n",
    "\n",
    "```py\n",
    "test = rotation_curve('IC2574')\n",
    "test.add_M_DM()\n",
    "test.df_circular.head()\n",
    "```\n",
    "you get as output \n",
    "\n",
    "![](./Figures/rot_curve_m.png)\n",
    "\n",
    "6. Now you can test that the `rotation_curve` class works well with all of the other `Galaxy_IDs` in the  `Rotation curves` directory. (You will probably have troubles loading the data for galaxy `'F563-1'`.)\n",
    "\n",
    "7. We will now proceed to do define a function outside of our class, that plots two columns of  the **DataFrame** `df_circular`, for all of the Galaxy_IDs in the  `Rotation curves` directory. First we need to have a list with all of the `Galaxy_IDs` in the  `Rotation curves` directory. \n",
    "\n",
    "* Use the **os** library, and the `os.listdir()` function, to create such a list, and remove `'F563-1'` from the list if you got error message when loading it in step 6.\n",
    "* Define the function `plot` that takes as entries two keywords `x` and `y`. It will loop over all of the elements of `Galaxy_IDs`, creating an instance of the class `rotation_curve` for every iteration,  and plot the columns `x` vs `y` in  the  `df_circular` class function.  For instance, if you use that function to plot `'radius (kpc)'` vs `'DM Mass (M_sol)'`, one should get a plot similar to:\n",
    "![](./Figures/M_DM_vs_r.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}