08b_neighbor_search_methods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# ASTR 21100 \n",
    "\n",
    "# \"Computational Techniques in Astrophysics\"\n",
    "\n",
    "## Instructor: Andrey Kravtsov\n",
    "\n",
    "### office: ERC 415; email: kravtsov@uchicago.edu\n",
    "### office hours: Tue, 10:30am-noon (unless noted otherwise)\n",
    "\n",
    "#### Teaching Assistants: \n",
    "#### Dimitrios Tanoglidis (dtanoglidis@uchicaago.edu)\n",
    "#### Georgios Zakharegkas (gzakharegkas@uchicago.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Efficient ways to carry out near neighbor searches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In many applications one needs to find neighbors of a point in some parameter space (or actual regular space, we are used to). The problem were this need arises is estimation of densities, which require defining of a scale $R$ corresponding to say a sphere of such radius and then finding the number of points around a local point that fall within such radius. \n",
    "\n",
    "If we have $N$ points, density estimate at each point location via brute-force neighbor search would scale as $\\mathcal{O}(N^2)$ and thus is extremely slow for large $N$, as we've already seen in direct $N$-body integration. \n",
    "\n",
    "However, efficient ways of finding nearest neighbors exist and you should be aware of at least the most basic of these, so that when you are faced with this problem, you know that the $\\mathcal{O}(N^2)$ search can be easily avoided. These efficient methods employ specialized data structures that encode information about spatial locations and structure of the data point distribution.\n",
    "\n",
    "Perhaps the simplest data structure of this kind is <a href=\"https://en.wikipedia.org/wiki/Linked_list\">linked list</a> - a linear collection of elements, called nodes, which consist of some data and a pointer pointing to the next node in the list. The last node's pointer is either zero in the standard linked list, or points to the first (aka head) node in a <i> circular linked list</i>. Linked list data structure was developed in 1955-1956 by researchers at the RAND corporation and has since then was the main workhorse data structure for lists in computer science and programming. \n",
    "\n",
    "Here is an example of a linked list class with functions to construct and manipulate it adopted from <a href=\"https://www.codefellows.org/blog/implementing-a-singly-linked-list-in-python/\">codefellows</a>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Node(object):\n",
    "\n",
    "    def __init__(self, data=None, next_node=None):\n",
    "        self.data = data\n",
    "        self.next_node = next_node\n",
    "\n",
    "    def get_data(self):\n",
    "        return self.data\n",
    "\n",
    "    def get_next(self):\n",
    "        return self.next_node\n",
    "\n",
    "    def set_next(self, new_next):\n",
    "        self.next_node = new_next\n",
    "\n",
    "\n",
    "class LinkedList(object):\n",
    "    \"\"\"\n",
    "    class for manipulating the nodes in the linked list\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, head=None):\n",
    "        self.head = head\n",
    "\n",
    "    def insert(self, data):\n",
    "        new_node = Node(data)\n",
    "        new_node.set_next(self.head)\n",
    "        self.head = new_node\n",
    "         \n",
    "    def size(self):\n",
    "        current = self.head\n",
    "        count = 0\n",
    "        while current:\n",
    "            count += 1\n",
    "            current = current.get_next()\n",
    "        return count\n",
    "\n",
    "    def search(self, data):\n",
    "        current = self.head\n",
    "        found = False\n",
    "        while current and found is False:\n",
    "            if current.get_data() == data:\n",
    "                found = True\n",
    "            else:\n",
    "                current = current.get_next()\n",
    "        if current is None:\n",
    "            print(\"%s is not in the list\"%data)\n",
    "            return None\n",
    "        return current\n",
    "\n",
    "    def delete(self, data):\n",
    "        current = self.head\n",
    "        previous = None\n",
    "        found = False\n",
    "        while current and found is False:\n",
    "            if current.get_data() == data:\n",
    "                found = True\n",
    "            else:\n",
    "                previous = current\n",
    "                current = current.get_next()\n",
    "        if current is None:\n",
    "            print(\"%s is not in the list\"%data)\n",
    "            return None\n",
    "        if previous is None:\n",
    "            self.head = current.get_next()\n",
    "        else:\n",
    "            previous.set_next(current.get_next())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some examples of how a linked list can be constructed and manipulated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jisheng\n",
      "Georgios\n",
      "Rebecca\n",
      "Dimitrios\n",
      "Nora\n",
      "Yuxi\n",
      "-------\n",
      "Andrey\n",
      "Jisheng\n",
      "Georgios\n",
      "Rebecca\n",
      "Dimitrios\n",
      "Nora\n",
      "Yuxi\n",
      "-------\n",
      "Jisheng\n",
      "Georgios\n",
      "Rebecca\n",
      "Dimitrios\n",
      "Nora\n",
      "Yuxi\n",
      "-------\n",
      "data for Rebecca: Rebecca\n",
      "Albert is not in the list\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "classlist = ['Yuxi', 'Nora', 'Dimitrios', 'Rebecca', 'Georgios', 'Jisheng']\n",
    "\n",
    "class_list = LinkedList()\n",
    "\n",
    "def print_class_list(class_list):\n",
    "    next_node = class_list.head\n",
    "    while next_node != None:\n",
    "        print(\"%s\"%next_node.get_data()) \n",
    "        next_node = next_node.get_next()\n",
    "    print('-------')\n",
    "    \n",
    "for name in classlist:\n",
    "    class_list.insert(name)\n",
    "\n",
    "print_class_list(class_list)\n",
    "class_list.insert('Andrey')\n",
    "print_class_list(class_list)\n",
    "class_list.delete('Andrey')\n",
    "print_class_list(class_list)\n",
    "\n",
    "print(\"data for Rebecca: %s\"%(class_list.search('Rebecca').get_data()))\n",
    "class_list.search('Albert')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The definition of a linked list via class is good for general linked lists. However, it also obfuscates the simplicity of linked lists, which can be constructed very simply and cheaply for specific purposes.\n",
    "\n",
    "For example, if we just need to construct linked lists for set of $N$ points, we just need an array of size $N$. Let's say we have a collection of $N$ points with $x_i$ values from 0 to 10 ($i\\in[1,N]$), but the values are random and are not correlated with the point index $i$, and we want to easily find all points with $x\\in[0,1)$, $x\\in [1,2)$, etc. We can achieve this by constructing 10 linked lists using single array of size $N$, as shown in the code below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14 4.868376465100422\n",
      "15 4.193175134882519\n",
      "27 4.502007728072269\n",
      "32 4.423568439749088\n",
      "42 4.61058298299315\n",
      "45 4.114409979987523\n",
      "47 4.417233952129211\n",
      "69 4.361769211074628\n"
     ]
    }
   ],
   "source": [
    "N  = 100 # number of points\n",
    "x  = np.random.uniform(high = 10., size=N)\n",
    "ll = np.zeros(np.shape(x), dtype=int)\n",
    "lh = np.zeros(10, dtype=int)\n",
    "\n",
    "# construct linked lists\n",
    "for i, xd in enumerate(x):\n",
    "    ich = int(xd) # list number to which x belongs\n",
    "    lcurr = lh[ich] # head point in the list\n",
    "    if lcurr: # if the list is not empty\n",
    "        lnext = ll[lcurr] # get next point in the list \n",
    "        while lnext: # while next point exists \n",
    "            lcurr = lnext      # walk down the list\n",
    "            lnext = ll[lcurr]\n",
    "        ll[lcurr] = i # when lnext=0, insert i as the last entry in the list\n",
    "    else: # if list is still empty\n",
    "        lh[ich] = i; ll[i] = 0 # insert the current point as its head\n",
    "        \n",
    "# now we can easily look at points with, say, 4 < x < 5 \n",
    "# -> i.e., we find points within)\n",
    "lcurr = lh[4]; lnext = ll[lcurr]\n",
    "while lnext:\n",
    "    print(lnext, x[lnext])\n",
    "    lnext = ll[lnext]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simple 1-d example illustrates how linked lists help with neighbor searches. If I want to find all points within $4.5\\pm 1$, I will need to examine also chunks 3 and 5 and computed distances for the points in those chunks to 4.5 and extract points which fall within specified distance. \n",
    "\n",
    "Analogously to splitting into 1-d chunks in the code above, if we want to find neighbors in 3D space within a certain distance $R$ from any point among $N$ points. If we split the volume of interest into $N_c$ cubes of size $2R$. Consider two arrays: <tt>lh</tt> of size $N_c$ and <tt>ll</tt> of size $N$. If we go through all the points \n",
    "and assign each entry in <tt>lh</tt> corresponding to a given cube the pointer to the first particle that's located in the cube, then assign pointer to the next particle located in the same value to the <tt>ll</tt> entry for that particle or 0 if there are no particles in the cube. This constructs a linked list of particles located in a given cube. \n",
    "\n",
    "If we do this for all cubes, we will have the ability to find all particles located in any cube, by going to its <tt>lh</tt> array entry to get the head particle in its linked list, and then following the sequence of pointers in <tt>ll</tt> until we hit 0 to get all the particles in the cube. \n",
    "\n",
    "Then to find neighbors around a given point we just need to find all the nearby cubes that could contain particles within $R$ and then scan all linked lists of these cubes and examine distances of the particles in the linked lists to the current particle selecting all those with distance $\\leq R$. \n",
    "\n",
    "Thus, having constructed linked lists of particles in all cubes in $\\mathcal{O}(N)$ operations, we can then use the linked lists to search for neighbors for each particle in the set in $\\mathcal{O}(N_l)$ operations, where $N_l$ is the number of particles in the linked lists of neighboring cells. This number will depend on the distribution of points. For uniform distribution, $N_l\\approx N/N_c$. In such case, we would have reduced the complexity of neighbor searching from $\\mathcal{O}(N^2)$ to $\\mathcal{O}(N[1+1/N_c)]$, i.e. to $\\approx \\mathcal{O}(N)$ for large $N_c$! \n",
    "Not bad for a bit of coding! This is an illustration how a clever data structure, like linked list can simplify algorithmic complexity of a problem. \n",
    "\n",
    "Construction of the linked lists is simple, but their use for neighbor searchers requires specifying tesselation (i.e., partitioning) of space into cells. If cells are equal volume cubes and are covering the entire spatial region of interest uniformly, the number of cells will scale with the number of spatial dimensions $D$, while most cells will be empty for non-uniform particle distributions. This is inefficient in terms of memory and memory requirement will become prohibitive for large $D$. \n",
    "\n",
    "For such cases, tree data structures are used to find neibors efficiently. Trees can be used to support tesselation of space that adhers to the point distribution and thus eliminates or elmost eliminates empty cells.  \n",
    "\n",
    "The trees are data structures allowing hierarchical organization of data. They consist of a number of parent and child nodes and links between them that support the tree data structure and make it traversable and searchable. There is a large variety of tree data structures and the way that they are implemented. \n",
    "The trees that are commonly encountered are binary, quad, oct, etc. tree. Such trees start with a root node and root node is split according to a specified rule, into 2 (binary tree), 4 (quad tree), 8 (oct tree), etc. child nodes. Binary, quad, and oct trees are natural choices for tesselation of 1-d, 2-d, and 3-d space, where all nodes are equal-size segments, squares, or cubes. However, the number of cells with each subsequent tree level grows exponentially as $2^l$, $4^l$, $8^l$ for the binary, quad, and oct tree, respectively. If we would continue this strategy for tesselating space into larger number of dimensions, each level would require creation of $2^{dl}$ cells. So for $d=10$, to go to level 2, requires $2^{20}\\approx 10^6$ cells. This is one of the manifestations of what's known as <i>the curse of dimensionality</i>. \n",
    "\n",
    "Thus, in practice for work in $D>3$ spaces, binary trees are preferred to minimize the rate of growth of the number of cells with increasing level of the tree. Even in this case, the curse is minimized, not eliminated. \n",
    "\n",
    "Here we will consider specific instances of <i>binary trees</i> - the k-d tree and ball tree. \n",
    "\n",
    "#### <a href=\"https://en.wikipedia.org/wiki/K-d_tree\">$k$-d tree:</a> \n",
    "\n",
    "short for \"$k$-dimensional tree\" is a useful binary tree data structure for nearest neighbor searches introduced in 1975 by Jon Louis Bentley along with a heuristic algorithm for its construction, an american computer scientist who later became also well known for his \"Programming pearls\" books. \n",
    "\n",
    "There are different ways to construct a $k$-d tree. In the simplest and most illustrative method, one starts with a root node representing the entire domain of interest and encloses all data points in the $k$-dimensional space. Then 2 child nodes are defined as sub-volumes of the root that are split by the hyperplane in the first dimension, with the position of the hyperplane corresponding to the median of the point distribution, or coordinate close to it. \n",
    "Each of the two child nodes is then split into further two child nodes along the next dimension, using the points corresponding to median coordinate of points in each child node along that dimension etc. The process continues until there are no more points to be used for splitting nodes or if a given specified level of refinement is reached. \n",
    "\n",
    "The tree construction can be achieved in $\\mathcal{O}(N\\log N)$ operations.\n",
    "This is because finding a split along first dimensions requires scanning through all points, i.e. $\\mathcal{O}(N)$ operations. Then on the next level, this will have to be done twice for $N/2$ particles, so again $\\mathcal{O}(N)$ and so on, so that for the tree level $l$, the number of operations is $\\mathcal{O}(lN)$. \n",
    "If we want to construct a tree with nodes that contain a single or a few particles the level of the required tree will be $l=\\mathcal{O}(\\log_2 N)$ (i.e. to for 2 points clearly $l=1$ always, for 4 points, $l=2$, etc). Thus, the total number of operations required for the binary tree construction is $\\mathcal{O}(N\\log_2 N)$. This is quite efficient, especially if we note that tree needs to be constructed only once and it can then be used to find neighbors for all points. \n",
    "\n",
    "Code example below is from the AstroML library and shows an example implementation of $k$-d tree class and shows how a $k$-d tree is constructed for a 2-d distribution of points. The implementation here uses recursion to make the tree construction easy to code up and transparent. Perhaps an even simpler implementation is shown on the <a href=\"https://en.wikipedia.org/wiki/K-d_tree\">$k$-d tree wikipedia page.</a>\n",
    "\n",
    "In astrophysics entire codes, such as N-body code <a href=\"http://www-hpcc.astro.washington.edu/faculty/trq/brandon/pkdgrav.html\"><tt>PKDGRAV</tt></a>, are built based on the k-d tree concept. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time to construct example KDTree = 0.000997304916381836\n",
      "time to construct scipy.spatial.cKDTree = 0.001994609832763672\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAFgCAYAAADTt046AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaIElEQVR4nO3dfZAkdX3H8c939+6Ip4QsdwgcC7OeJ6ceGrw9w1kaDRExJljogQbqyjxYgFaBRIlJNFqHdYVlnrSwEtQoFqZSHIXmDkVQAhhKDeESdjbKg0+Bc0c2CMLdkpxC2If55Y/pPpq5edzpnunv9PtVtQU7j907n/vMr3/T020hBAGAVyODXgAA6AUlBsA1SgyAa5QYANcoMQCuUWIAXKPEALhGiQFwjRJzxMweMrNXtLnNF8zsin4t07AwsxkzO2PQy4HuUWJOmNmYpOMlfb/Hx+n7P9boOZ8ys58nfv6un8uA4bVi0AuAjr1M0gMhhP/L8knMbEUIYTGDh35zCOH2DB4XBcdIzI+XS7pPksxstZntMrM9ZvbrZjZtZgfN7HpJv9TsAczsHyWdJOmr0WjoT6PLZ8zsz8zsHkm/MLMVZrbOzHab2WNm9mMzuzTxOE2v64aZvdDMDpjZ5sTjPm5mvxH9/gEzezBat++Z2Vvr7j9jZn9iZveY2S/M7PNmdqyZfT26z+3RCDa+7Qejx5kzs2vMrOHfKq31Q5+EEPhx8CPps5Iul/QCSdPR/x8hqSLpfZJWSjpX0oKkK1o8zoykMxpc9h1JJ0p6jmpvbmVJOyStkrRe0j5Jb2x1XTfPmbjuQtU2kVdL+mdJf5O47m2S1kXP+buSfiHp+LrH3SvpWEknSPpZ9Ld5RfS3+RdJlydue1+0jkdLujP5d4qXcTnrx8+A/20MegH46fCFku6S9BlJP5Z0dnTZayU9LMkSt/u3ZZbYOxO/nybpJ3W3+aCka1pd1+Y5fy7picTPhYnrb5R0r6R7JB3R4nG+E6974nG3J37fLenTid/fI+nLidu+O3Hdb0t6sP7vspz142ewP8yJOWBmJukU1UYFnwghfCW6ap2k/w7Rv7RIJXG/7ZL+Pvr12yGEN7V4mocS/1+StM7MnkhcNirp222ua+Utofmc2OdUK7KLQghPJ5b/9yRdJmkiuuh5ktbW3ffRxP8/1eD35yV+T65jRbW/X73lrh8GhBLz4QXRf8+Q9A0z+0YIYUrSTyWdYGaWKLKTJD0oSSGEayVdW/dYzQ4gl7z8IUk/DiG8qP5GZvaqZtcth5k9T9KVkj4v6SNmtjuEcMDMSqqV2+sl3RVCWDKz70iyHp7uxMT/n6TaKLZe03VHPjGx78PLJd0TQrhX0kWSbjCz41XbxFyUdGk0Gb9N0q+1eaxHVRvRtfIfkv43mux/jpmNmtkpZvbKNtctxycllUMIF0i6WbVNZkl6rmrF+pgkmdkfqjYa7cXFZjZuZkdL+nNJ1ze4Tdrrh4xRYj68TLX5IoUQvqzaJP+XVXv9tkn6A0lzqk1+72nzWB+T9GEze8LM3t/oBiGEJUlvlnSqanNwj0u6WtJRra5r87zxJ6Lxzw1mdrak35L07ug2l0nabGbbQwjfk/Rx1Yr60ehvcGeb52hnl6RbVZuo3yfpsJ2Ce1g/DIg9ezoFGE5mNiPpghbzcnCqcCOxfuyxzld/iolsDUbhSmzQzOx4M7vRzB42s2BmE4NeJgwHM/sdM/vXaKrgETP7nJkdOejlyhol1n9VSbdIOmfQC1IkIYSJAmxKHqXaPN86SS+RNC7prwe6RH1Q6BIzs5HEV1v2m9kXo0+uZGa3mNkldbf/bvQJoMzsxWZ2W/S1mR+a2ds7ec4QwqMhhE9Jujv1FUJuDChbu0IIt4QQngwhzKm2i8qr0163vCl0iUm6VNJbJL1OtXevOUlXRdftknR+fEMze6lqO0LebGbPlXRbdJvnR7f7lJlt6t+iI+fykK3XSrp/uSvgRdFL7F2SPhRCmI32FP+IpHPNbIWkGySdGu10KUnbJe2JbneWpJkQwjUhhMUQwrRqX3k5t/+rgJwaaLbM7A2Sfl+174AOtaKXWEm1HUefiL5m8n1JS5KODSEcVG3ny/Oi256nZ/Z+L0k6Lb5fdN/tko7r7+IjxwaWLTPbqtpI7twQwo/SWZ38KvrXjh5S7YvPzXaivE7S5Wb2LdWO7nBH4n7fDCG8oQ/LCJ8Gki2rHfn3xui5v7Gcx/Cm6COxz0j6aDysN7Njor3IY19T7Z1xp6TrQwjV6PKbJJ1sZu8ws5XRzyvN7CWdPKnVjmN1RPTrEdbkuFZwre/ZMrNTVPvk+z0hhK+mujY5VvQS+6Rq71q3mtlB1Y5NdVp8ZTRHsUe1L17vSlx+UNKZqm0GPCzpEUl/qWeKqZ2nVDs0jST9IPodw2UQ2fpjScdI+nzi611DP7HP144AuFb0kRgA5ygxAK5RYgBco8QAuNbVfmJr164NExMTGS0KhkG5XH48hHBMt/cjW2inWba6KrGJiQlNTU2lt1QYOmZWaX+rw5EttNMsW2xOAnCNEgPgGiUGwDVKDIBrPR/FYmJiQpXKsuZy4VipVNLMzEymz1HUbI2OjmppaWnQizEw3War5xKrVCri+5fFY9bLibg7U9RsmVkh1zvWbbbYnATgGiUGwDVKDIBrlBgA1ygxAK5RYgBco8QAuEaJAXAt0xIrV+Z01R0PqFyZy/JpUEBkC7HMTp5brsxp+9V7Nb9Y1aoVI7r2gq2aLI1l9XRIUbkyp7379mvr+jW5fM3Ill9ZZCuzEtu7b7/mF6uqBmlhsaq9+/YTNAc8FATZ8imrbGW2Obl1/RqtWjGiUZNWrhjR1vVrsnoqpKhRQeQN2fIpq2xlNhKbLI1px1mb9PX7fqo3nXI875ROxAWxsFjNbUGQLZ+yylamc2I7b7pf84tV3T1zQBuPO5KwOTBZGtO1F2zN/ZwY2fInq2wxJ4bDTJbGcv1akS2/ssgWc2Jwh2whKdM5sbxvlsAnsoWkzEpMyv9mCfwiW4jxtSMArlFiAFyjxNAQ301EVtLOVqZzYvDJw1eP4FMW2WIkhsN4+OoRfMoiW5QYDsN+WMhKFtlicxKHYT8sZCWLbFFiOKT+WE+UF9KSZbYoMahcmdPu6Vn9U3lWi0tM5iM9/cgWJVZw8adFTy9UFaLL+FI10tCvbDGxX3Dxp0VxyExM5iMd/coWI7GCSx6obnTE9LYtJ2rb5nFGYehZv7JFiRUcn0QiK/3KFiUGPolEZvqRLebEALhGiQFwjRID4BolBsA1SgyAa5QYANcoMQCuUWIAXKPEALhGiQFwjRID4BolBsA1SgyAawMpMU7MiqyQreLp+6F4ODErskK2iqnvIzFOzIqskK1i6nuJcWJWZIVsFVPfNydbHbK2/tx0QDfIVjEN5PDUjQ5Zy3wG0kC2iic3u1jUz2fsnp7lUyakIpmt+YWqrrz9R+RqiOTmRCH1p3fibNRIS5yt+YWqqpLufOBx3T1zgFwNidyMxOL5jMvO3Ki3bTlRi0t8yoR0xNl69YvWasREroZMbkpMqoXt4tM3aNvmcT5lQqomS2N67xknk6shlJvNySRO6IoskKvhlMsSkzihK7JBroZPrjYnAaBblBgA1ygxAK5RYgBco8QAuEaJAXCNEgPgGiUGwDVKDIBrlBgA1ygxAK5RYgBco8QAuEaJAXCNEgPgGiUGwDVKDIBrlJhT5cocp7RDJrxlK7eHp0ZznAwWWfGYLUZiDtWfaJhTjyEtHrNFiTkUnwyWU48hbR6zxeakQ8lTj42tXnXo3TLvw37kX5ytPdOzCoNemA5RYk7FheVt/gI+7J6e1fxiVXumZ3OfKzYnHfM4f4H885YrSswxj/MXyD9vuWJz0rHk3NjW9WtyPeSHH95yRYk5N1kay33I4I+nXLE5CcA1SgyAa0NVYt6+8wU/yFZ+uZwTK1fmDpt09PidL+QP2fLHXYmVK3M6/7N3aWEpaOWo6bqLXqXJ0ljDfVsIGrpBtnxytzm5e3pW80tBQdL8UtDu6VlJ/vZtQf6QLZ/cjcSsye/e9m1B/pAtn9yV2LbN4/pSeVYLi1WtXDGibZvHD13nad8W5A/Z8sldiU2WxnTdhbwrIn1kyyd3JSbxrojskC1/3E3sA0ASJQbANUoMgGuUGADXKDEArlFiAFyjxHKAIyQgK0XIlsv9xIYJR0hAVoqSLUZiA+btzDLwoyjZosQGjCMkICtFyVYhNicbHeguLzhCgm9ka/CGvsQGNS/QTbj5vp5PZCsfhr7Eejkq53LfZYsyoVp0ZCsfhr7E4nmB+BhRnc4L9BIWDmdcDGQrH4a+xJY7L9BLWJqFO8/zJ+ge2cqHoS8xqfG8QLsXvT4sY6tX6ao7Huh4HqI+3N28+xY5kN6QrcErRInV6/RF37Z5XCZp07qjtPOm+7sOycWnbzh0eafvvsx5+Nbu9YvzseOsTZp7cl5jq1cdlq1Wj022DlfIEmv3ote/2JJ6DsnY6lUaMZMUWm4GMOfhW6vXr1E+Ot0hlWw1V8gSazchW/9iB6mjCdxmISlX5rTzpvu1VA0aHTHtOGtT082A5U4WIx9avX6N8tHp6022mitkibWbkK1/sc/ZPK5zNo+3nUtoFpI4gEFSCEFzT84/6/JkMC8+fUMhdlAcVq2y1SgfnX44QLaaK2SJSa13AqwPlqSOXvhmgWwWwGaXF2EHxWHW7PVrNilPtnpjIYSOb7xly5YwNTX17AcwUzeP4c1yJ0Prw9ksrF4/LWr2uptZOYSwpdvHI1vPZKvdepOtZyvsSKxTy5kMbRbOZu/OngKG9JCtdHAUizaWcySAohwCBb0hW+lgJNZGJxOv9cP2Yf0UCOkiW+mgxDrQaljebHg/jJ8CIX1kq3eUWI+azWsUZT4C2SFbnWFOrEPNTrhQlKNnIjtkqzeMxDrQajcLhvfoRaNsxchWZyixDrT7KJzhPZar3aeNZKs9Nic7wLAeWSFbvWMk1gGG9cgK2eodJdYhhvXICtnqDZuTAFyjxAC41vPmZKlUkpmlsSxwpFQq9eU5ipit0dHRQq53rNts9VxiMzMzvT4E0BDZQifYnATgGiUGwDVKDIBrlBgA17o6xr6ZPSapkt3iYAiUQgjHdHsnsoUONMxWVyUGAHnD5iQA1ygxAK5RYgBco8QAuFa4EjOzGTM7I+Pn+IKZXZHlcyB/yNZgFK7EBs3MTjeze83sCTPbb2Y3mNkJg14uDBczu8bMgpltGPSyZI0S67/vSXpjCOFXJK2T9F+SPj3YRcIwMbPXSHrhoJejXwpdYmY2YmYfMLMHo1HRF83s6Oi6W8zskrrbf9fMtkX//2Izu83MDpjZD83s7Z08Zwjh0RDCw4mLliQN/btl0QwiW9F9V0j6W0mXtLvtsCh0iUm6VNJbJL1OtVHRnKSrout2STo/vqGZvVRSSdLNZvZcSbdFt3l+dLtPmdmmTp7UzE4ysyckPSXp/ZL+KpW1QZ4MJFuS3ifpWyGEe9JYCQ+KXmLvkvShEMJsCOFpSR+RdG70bnaDpFPNLD5C23ZJe6LbnSVpJoRwTQhhMYQwLWm3pHM7edIQwk+izcm1kj4s6QeprhXyoO/ZMrMTo+fdkf7q5FfRS6wk6YZokv0JSd9XbfPu2BDCQUk3Szovuu15kq5N3O+0+H7RfbdLOq6bJw8hHJD0D5K+EoUbw2MQ2bpS0s4Qwv+kuSJ5V/R/OA9JemcI4c4m118n6XIz+5ak50i6I3G/b4YQ3pDCMqxQbbPhlyUdSOHxkA+DyNbrJb3GzJLTE3eZ2R+FEHYt4/FcKPpI7DOSPhoP683sGDM7O3H911R7Z9wp6foQQjW6/CZJJ5vZO8xsZfTzSjN7SbsnNLNtZrYxmvg9RtInJP1nNCrD8Oh7tiSdLOlXJZ0a/UjSm1XbfB1aRS+xT0q6UdKtZnZQ0l5Jp8VXRnMUeySdodpEa3z5QUlnqrYZ8LCkRyT9paQjOnjOEyTdIumgpHslVSW9NYV1Qb70PVshhJ+FEB6Jf6KLHw8hPJXOKuUTh+IB4FrRR2IAnKPEALhGiQFwjRID4FpX+4mtXbs2TExMZLQoGAblcvnx5ZwohGyhnWbZ6qrEJiYmNDU1ld5SYeiY2bLOWES20E6zbLE5CcA1SgyAa5QYANd6/gL4xMSEKpVinrh5dHRUS0tLg16MgSiVSpqZmcn0ObLIlofXzMMyxrJY1m6z1XOJVSoVFfWrS2ZW6HXPWhbZ8vCaeVjGWBbL2m222JwE4BolBsA1SgyAa5QYANcoMQCuUWIAXKPEALhGiQFwjRID4FqmJVauzOmqOx5QuTKX5dOggMgWYpmdPLdcmdP2q/dqfrGqVStGdO0FWzVZGsvq6ZCicmVOe/ft19b1a3L5mpEtv7LIVmYltnfffs0vVlUN0sJiVXv37SdoDngoCLLlU1bZymxzcuv6NVq1YkSjJq1cMaKt69dk9VRIUaOCyBuy5VNW2cpsJDZZGtOOszbp6/f9VG865XjeKZ2IC2JhsZrbgiBbPmWVrUznxHbedL/mF6u6e+aANh53JGFzYLI0pmsv2Jr7OTGy5U9W2WJODIeZLI3l+rUiW35lkS3mxOAO2UJSpnNied8sgU9kC0mZlZiU/80S+EW2EONrRwBco8QAuEaJoSG+m4ispJ2tTOfE4JOHrx7BpyyyxUgMh/Hw1SP4lEW2KDEchv2wkJUsssXmJA5JHiaF/bCQpiyzRYlBUuO5iotP3zDoxcIQyDpbbE5C5cqcrrz9R8yDIXX9yBYjsYKL3yWfXqgqSBphHgwp6Ve2KLGCiz8tCqoNy1+9Ya3ee8bJzIOhZ/3KFpuTBZf8tGjVyhEKDKnpV7YYiRUcR4RAVvqVLUoMHBECmelHtticBOAaJQbANUoMgGuUGADXKDEArlFiAFyjxAC4RokBcI0SA+AaJQbANUoMgGuUGADXBlJinNMQWSFbxdP3o1hwTkNkhWwVU99HYpzTEFkhW8XU9xLjnIbICtkqpr5vTnIkUWSFbBXTQI7s2uxoj8kTbBJALEe7bGH45Obw1PWTsjvO2qS5J+cpNPQsmS1J+tAN92rb5nFyNSRyU2LJSdn5hap2fOU+VUPgUyb0LJktSdr17z/R7ulZcjUkcrOza3JSdmTEVA2BT5mQijhbFv0eRK6GSW5KLJ6UvezMjdp59il8yoTUxNk6/7STJIlcDZncbE5Kz56U3XjckUzyIzVxtj4m6bIzN5KrIZKrEkviXIjIysWnbxj0IiBFudmcBIDloMQAuEaJAXCNEgPgGiUGwDVKDIBrlBgA1ygxAK5RYgBco8QAuEaJAXCNEgPgGiUGwDVKDIBrlBgA1ygxAK5RYgBco8QAuEaJOVWuzOmqOx5QuTI36EXBkPGWrdweYx/N1Z9omPMnIi0es8VIzKHkyWA5fyLS5DFblJhDyRMNc/5EpMljtticdCg+Gezeffs1tnrVoXfLvA/7kX9xtvZMzyoMemE6RIk5FReWt/kL+LB7elbzi1XtmZ7Nfa7YnHTM4/wF8s9brigxxzzOXyD/vOWKzUnHknNjW9evyfWQH354yxUl5txkaSz3IYM/nnLF5iQA1ygxAK65LLFm3+3y9p0v5A/Z8sfdnFi5MqfzP3uXFpaCVo6arrvoVZosjbn8zhfyhWz55G4ktnt6VvNLQUHS/FLQ7ulZSf72bUH+kC2f3JWYNfnd274tyB+y5ZO7zcltm8f1pfKsFharWrliRNs2j0vyt28L8ods+eSuxCZLY7ruwsaB8rRvC/KHbPnkrsQkAoXskC1/3M2JAUASJQbANUoMgGuUGADXKDEArvX86eTo6KjM6ncTBHqXVbbI63DpucSWlpYUgpdTCqSLfwzZGuZsmVnTdSNX3WFzMgc4QgKyUoRsudzZdZhwhARkpSjZYiQ2YBwhAVkpSrYosQHjCAnISlGyVYjNyXJlLrdHIOAICb6RrcEb+hIb1LxAN+HmS8c+ka18GPoSazQv0OmLutx32aJMqBZdL9mSpKvueIBspWDo58SWOy8Qh+Xjt/5Q26/e29VH1EWZUC26XrIliWylZOhHYsudF+j0XbaTHRMv+Qvpkq6XPN9GR0cHvQgD10u2JLXMVqc7vOYhW2nvnNtttoa+xKTG8wLtNhXjd9n4UMVjq1c1HP432+u6/vE73QzwsrnAXuU1y82WpEMjuEbZavVNheTjb5k4WlMzBwqdrUKUWL1OX8xtm8dlkjatO0o7b7q/oxc/GbCLT99w6PI907N6eqGqoNYju17nWTBY7bIV50OSLjtzo8ZWrzosW60em2wdrpAl1u7FrA+ipI5e/GYBLlfm9KWphxS/t46OPjN/Uv+uXT8CHNZ9e4ZVq2wl8yHVRmSdznG1KseiZ6uQJdbuxawPVpA6evGbBXjvvv1arNZiZpLOnRxvuYlZhH17hlWrbCXzEf/eabG0KseiZ6uQJdbuxawP1jmbx3XO5vG2L36zQDZ6PKl5MIuwb8+wapWtZA7i3zstllZlV/RsWTeHOtmyZUuYmpp69gO0OKSIZ8mhuKSGIWu07s0mdRtdHr9bxgHM60RrI81edzMrhxC2dPt4RclWnINLfvNFCiE0zEWz9W5226mZA4XOFiXWRqu5iFbr3umOsnn+2korlFhv4vJplK12613/6WSvGcybbrNVyM3JbiznE51uPsoexuE9OpNGtlopSraGfo/9Xi1nr2z2qkYn0sgWGIm11cnEa1E+yka60sgWKLGOtBqWF+mjbKSv12xtuaLPC5xDlFiPivRRNvqLbHWG8WiHmp1woShHz0R2yFZvGIl1oNWnjWw6oheNshUjW52hxDrQ7qNwhvdYrnafZJOt9tic7ADDemSFbPWOkVgHGNYjK2Srd5RYhxjWIytkqzdsTgJwjRID4FrPm5OlUqmwx1sfHR0t7LqXSqW+PMew/n1bZaebXA1jBrvNVs8lNjMz0+tDAA2RLXSCzUkArlFiAFyjxAC41tXhqc3sMUmV7BYHQ6AUQjim2zuRLXSgYba6KjEAyBs2JwG4RokBcI0SA+AaJQbANUoMgGuUGADXKDEArlFiAFyjxAC49v+jM3dxdiSe+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.spatial import cKDTree\n",
    "from time import time\n",
    "\n",
    "# We'll create a KDTree class which will recursively subdivide the\n",
    "# space into rectangular regions.  Note that this is just an example\n",
    "# and shouldn't be used for real computation; instead use the optimized\n",
    "# code in scipy.spatial.cKDTree or sklearn.neighbors.BallTree\n",
    "class KDTree:\n",
    "    \"\"\"Simple KD tree class\"\"\"\n",
    "\n",
    "    # class initialization function\n",
    "    def __init__(self, data, mins, maxs):\n",
    "        self.data = np.asarray(data)\n",
    "\n",
    "        # data should be two-dimensional\n",
    "        assert self.data.shape[1] == 2\n",
    "\n",
    "        if mins is None:\n",
    "            mins = data.min(0)\n",
    "        if maxs is None:\n",
    "            maxs = data.max(0)\n",
    "\n",
    "        self.mins = np.asarray(mins)\n",
    "        self.maxs = np.asarray(maxs)\n",
    "        self.sizes = self.maxs - self.mins\n",
    "\n",
    "        self.child1 = None\n",
    "        self.child2 = None\n",
    "\n",
    "        if len(data) > 1:\n",
    "            # sort on the dimension with the largest spread\n",
    "            largest_dim = np.argmax(self.sizes)\n",
    "            i_sort = np.argsort(self.data[:, largest_dim])\n",
    "            self.data[:] = self.data[i_sort, :]\n",
    "\n",
    "            # find split point\n",
    "            N = self.data.shape[0]\n",
    "            N2 = int(N/2)\n",
    "            N21 = int(N/2 - 1)\n",
    "            split_point = 0.5 * (self.data[N2, largest_dim]\n",
    "                                 + self.data[N21, largest_dim])\n",
    "\n",
    "            # create subnodes\n",
    "            mins1 = self.mins.copy()\n",
    "            mins1[largest_dim] = split_point\n",
    "            maxs2 = self.maxs.copy()\n",
    "            maxs2[largest_dim] = split_point\n",
    "\n",
    "            # Recursively build a KD-tree on each sub-node\n",
    "            self.child1 = KDTree(self.data[N2:], mins1, self.maxs)\n",
    "            self.child2 = KDTree(self.data[:N2], self.mins, maxs2)\n",
    "\n",
    "    def draw_rectangle(self, ax, depth=None):\n",
    "        \"\"\"Recursively plot a visualization of the KD tree region\"\"\"\n",
    "        if depth == 0:\n",
    "            rect = plt.Rectangle(self.mins, *self.sizes, ec='k', fc='none')\n",
    "            ax.add_patch(rect)\n",
    "\n",
    "        if self.child1 is not None:\n",
    "            if depth is None:\n",
    "                self.child1.draw_rectangle(ax)\n",
    "                self.child2.draw_rectangle(ax)\n",
    "            elif depth > 0:\n",
    "                self.child1.draw_rectangle(ax, depth - 1)\n",
    "                self.child2.draw_rectangle(ax, depth - 1)\n",
    "\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Create a set of structured random points in two dimensions\n",
    "np.random.seed(0)\n",
    "\n",
    "X = np.random.random((30, 2)) * 2 - 1\n",
    "X[:, 1] *= 0.1\n",
    "X[:, 1] += X[:, 0] ** 2\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Use our KD Tree class to recursively divide the space\n",
    "t1 = time()\n",
    "KDT = KDTree(X, [-1.1, -0.1], [1.1, 1.1])\n",
    "t2 = time()\n",
    "print(\"time to construct example KDTree =\",t2-t1)\n",
    "\n",
    "t1 = time()\n",
    "ckdt = cKDTree(X)\n",
    "t2 = time()\n",
    "print(\"time to construct scipy.spatial.cKDTree =\",t2-t1)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot four different levels of the KD tree\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "fig.subplots_adjust(wspace=0.1, hspace=0.15,\n",
    "                    left=0.1, right=0.9,\n",
    "                    bottom=0.05, top=0.9)\n",
    "\n",
    "for level in [1,2,3,4]:\n",
    "    lev = level\n",
    "    if lev > 4 : lev = 4\n",
    "    ax = fig.add_subplot(2, 2, lev, xticks=[], yticks=[])\n",
    "    ax.scatter(X[:, 0], X[:, 1], s=9)\n",
    "    KDT.draw_rectangle(ax, depth=level - 1)\n",
    "\n",
    "    ax.set_xlim(-1.2, 1.2)\n",
    "    ax.set_ylim(-0.15, 1.15)\n",
    "    ax.set_title('level %i' % level)\n",
    "\n",
    "# suptitle() adds a title to the entire figure\n",
    "fig.suptitle('$k$d-tree Example')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the above 2-d example, $k$-d tree produces cells that tend to be elongated along direction that undergoes splitting later. If $d=10$, for example, by the time we will split along the 10th dimension for the first time, $2^{10}\\approx 10^3$ nodes will already have been created. This means that splitting along later dimensions will not be efficient.\n",
    "\n",
    "This can be rectified by choosing the dimension on which to split, not in the order of dimensions, but by some measure such as the dimension along which the extent of points is largest. This, however, brings up a problem of dimensionality if our space is not metric (i.e., different dimensions have different units): how to we compare extents along different dimensions meaningfully? \n",
    "\n",
    "The problem of measuring distances in non-metric spaces, however, is a general problem of neighbor searches because distance needs to be defined to define what a neighbor is. Solution of this problem depends on specific case. Here I would just note that definions of distance that side-step the fact that space is non-metric do exist (e.g., <a href=\"https://papers.nips.cc/paper/1581-classification-in-non-metric-spaces.pdf\">Weinshall et al.</a>). Also note that in the $k$-d tree one can compute local density by taking the number of particles in a tree cell and dividing by cell volume. This density estimate requires no distance estimates, but does estimate local density (see <a href=\"http://adsabs.harvard.edu/abs/2005MNRAS.356..872A\">Ascasibar & Binney 2005</a> for details).\n",
    "\n",
    "Another type of tree - the <i>ball tree</i> - rectifies the assymetry of the cells in the $k$-d tree, which may increase efficiency of neighbor searches for many point distributions. \n",
    "\n",
    "\n",
    "Neighbor searches can be done using $k$-d tree efficiently, because we can walk the tree in $O(N\\log_2 N)$ steps. Starting with a root node we can examine the child nodes and select only those nodes whose extent falls within the range of distance within which we are searching for neighbors (you can see detailed algorithm <a href=\"https://en.wikipedia.org/wiki/K-d_tree#Nearest_neighbour_search\">here</a>). \n",
    "\n",
    "#### Another type of tree: <a href=\"https://en.wikipedia.org/wiki/Ball_tree\">ball tree</a>\n",
    "\n",
    "A ball-tree is a binary tree, in which every node defines a d-dimensional hypersphere, or \"ball\", containing a subset of the points for which the tree is constructed. Each internal node of the tree partitions the data points into two disjoint sets which are associated with different balls. \n",
    "While the balls themselves may intersect, each point is assigned to one or the other ball in the partition according to its distance from the ball's center. Each leaf node in the tree defines a ball and enumerates all data points inside that ball.\n",
    "\n",
    "As with k-d trees there are many different possible algorithms for construction of ball tree. One of the simplest ways is similar to the k-d tree construction: the tree is constructed recursively starting with a sphere that encloses all points. The points are then split into two subsets along the successive dimensions, with partition coordinate defined by the median value of all points along that dimension or a position close to the median. Two balls are then defined with the centers given by the mean position of points in the subsets and the radius that corresponds to the distance from the centroid to the most distant point in the subset.  \n",
    "\n",
    "As for the $k$-d tree, ball tree construction for $N$ points requires $\\mathcal{O}(N\\log_2 N)$ operations. Below is an example code constructing ball tree recursively from the AstroML package in 2-d for the same point distribution as above, which illustrates its properties and differences from the $k$-d tree nicely. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAFgCAYAAADTt046AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd1hUx9fHvwMsipUqRRTEgl1pogkqYoti7NEYo8ZujBqjxt6iJsYUTYyxx/KzxwImVozYkijSxBo1Cog0EWlK3d3z/rEsLxpWdpe7e3dxPs+zj7I795xz2cP3zsydOZcRETgcDsdYMRE7AA6Hw6kIXMQ4HI5Rw0WMw+EYNVzEOByOUcNFjMPhGDVcxDgcjlHDRawSwhhzZYwRY8ys+OfzjLFxYsdlTLz6O+QYLlzEDBTGWBxjLI8x9pwxlsEYO84Yq6cjX8QYa6QL2yr87WCMFRafm/IVoy//nMoFFzHD5l0iqgHAEUAqgJ/ECEJHvZFviKhGqVcbHfjgvAFwETMCiCgfwCEAzZXvMcYCGWPRjLFsxlgCY2ypNrYZYxeL/xtT3CMayhjzZ4w9ZozNYYylANhe3LYPY+waYyyTMfY3Y6x1KTtOjLHDjLE0xlgsY2yalvEMZYw9ZIzVKv65F2MshTFmV/zzj8Xnm80Yi2SMdSx17FLG2EHG2G7GWA5j7AZjrAljbB5j7EnxcT1KtT/PGFvJGLvKGMtijB1ljFmriKs2Y+wXxlgyYyyRMbaCMWaqzTlyhIWLmBHAGKsGYCiAK6XefgFgJABLAIEAPmaM9dfUNhF1Kv5vm+Ie0YHinx0AWANwATCBMeYJYBuAiQBsAGwC8BtjrApjzATA7wBiANQF0BXAdMZYTy3iOQDgMoC1jDEbAL8AGEdEacVNwgG0LY5tL4CDjLGqpUy8C2AXACsA0QBOQ5HndQEsK467NCMBjAHgBEAKYK2K0HYWf94IgAeAHgD4PKMhQET8ZYAvAHEAngPIhOKPJwlAq9e0/wHAmuL/uwIgAGbFP5+HQghUHUsAGpX62R9AIYCqpd7bAGD5K8fdBdAZgC+AR698Ng/AdhX+dgDILz435Wtnqc8tATwCcAPApnJ+TxlQCDAALAVwptRn7xb/Dk2Lf65ZfK6WpX4vX5dq37z4vE1L/w4B2AMoAGBRqu0wAOfEzhP+IvA7L4ZNfyL6o3jY0g/ABcZYcyJKYYz5AvgaQEsA5gCqADhYnkHG2C0oelcA0IuILqlomkaKYawSFwCjGGNTS71nDkUPRgbAiTGWWeozUwCqbAPAd0S0sKwPiCiTMXYQwAwAg16JfyYUPSAnKESmFgDbUk1SS/0/D8BTIpKV+hkAakAhnACQUKp9PADJK/YAxblLACQzxpTvmbxyLEck+HDSCCAiGREdgUIs/Irf3gvgNwD1iKg2gI0AmAoTpW21oP+fTH+dyLxa3iQBwJdEZFnqVY2I9hV/FvvKZzWJqLem5woAjLG2UAzx9qHU8K54/msOgCEArIjIEkAW1Djv11D6jm99AEUAnr7SJgGKnphtqfOrRUQtKuCXIxBcxIwApqAfFPM8d4rfrgngGRHlM8baAfigAi5SAbiV02YLgEmMMd/ieKoX31yoCeAqgOziGwEWjDFTxlhLxpiPpoEUz2/tBjAfwGgAdRljk4s/rgnF0DoNgBljbDEUPbGK8CFjrHnxvOMyAIdK9dwAAESUDCAEwPeMsVqMMRPGWEPGWOcK+uYIABcxw+Z3xthzANkAvgQwiohuFX82GcAyxlgOgMUAfq2An6UAdhbfdRxSVgMiigAwHsA6KOah/gXwUfFnMijmn9oCiIWiJ7MVQO3X+Jz9yjoxZe9nJYDHRLSBiAoAfAhgBWOsMRST9CcB3INi6JePig/pdkExR5cCoCoAVXdVR0IxfL4NxfkfgmLpC0dkGBEvish5M2GMnQewm4i2ih0LR3t4T4zD4Rg1XMQ4HI5R88aJGFPsSeymYx87GGMrdOmDU3GIyF/IoSTPLXF440RMbBhjjoyx3xhjSUyx8dpV7Jg4lYPiu8V/Ft+gSWGMbSm+e1yp4SKmf+QATuGVRZwcjgDUBrACioXAzQA4A/hW1Ij0wBstYsXrfeYyxh4wxtIZY78qNwAzxk4xxqa80j6GMTaw+P9NGWNnGGPPGGN3VS1NeBUiSiWi9VDsAeRUUkTKrb1EdIqIcokoA4q1fW8LfW6GxhstYlCsCeoPxf4/JyjW//xc/NleKPbHAQAYY82h2H5ynDFWHcCZ4jZ1itutZ4zxFdwcJYaQW50A3Cq3lZHzpovYRAALiOhx8cLKpQAGM0X9rCAAbRljyn2GwwEcKW7XB0AcEW0nIikRRQE4DGCw/k+BY6CImluMse4ARkGxELpS86aLmAuAoOKJ0EwotvTIANgTUQ6A4wDeL277PoA9pY7zVR5XfOxwKMrXcDiAiLnFGGsPRU9uMBHdE+Z0DJc3vYpFAoAxRPSXis/3AVjCFIUDLQCcK3XcBSLqrocYOcaJKLnFGPOAojDAGCI6q40NY+NN74ltBPClslvPGLMr3mit5AQUV8ZlAA4Qkbz4/WMAmjDGRjDGJMUvH8ZYM3WcFm9yrlL8YxX2clE/TuVA77nFGGsJxZ3vqUT0u6BnY8C86SL2IxRXrZDijdRXoCjwBwAonqM4AqAbFN1z5fs5UFT2fB+KYoUpAFbh/4WpPPKgKNYHAP/g/+tccSoPYuTWTAB2AH4ptbG+0k/s8w3gHA7HqHnTe2IcDsfI4SLG4XCMGi5iHA7HqOEixuFwjBqN1onZ2tqSq6urjkLhVAYiIyOfEpGdpsfx3OKUh6rc0kjEXF1dERERIVxUnEoHYyxem+N4bnHKQ1Vu8eEkh8MxariIcTgco4aLGIfDMWq4iHE4HKOGixiHwzFquIhxOByjhosYh8MxariIcTgco+ZNr+yqNkSE7OxsPH/+HFKpFFKpFIwxSCQSSCQSWFtbw9zcXOwwOZw3Di5ipZBKpbh16xbCw8Nx9+5dJCUlITExEYmJiUhKSgJjDLVq1YKZmRnMzMxARJBKpSgsLERGRgYsLS3h5OSEunXrwsnJCfXq1UObNm3g4+MDJycnsU+Pw6mUvNEilpaWhjNnzuDq1asIDw9HTEwM6tWrBx8fHzRv3hytW7cuEaS6deuiZk3VD1OWyWRIS0srEb6kpCTEx8dj48aNGDNmDKpWrQofHx/4+Pjg7bffhp+fH8zM3uhfP4cjCG/cX1FcXByCg4MRFBSEa9euISAgAO3bt8eKFSvg5eWFWrVqaWXX1NQUDg4OcHBwgKen50ufERFiY2MRHh6OiIgIzJo1C3FxcejTpw/69++PHj16oFq1akKcHofzxvFGiNjTp0/xyy+/4MCBA0hISEDfvn3x+eefo2vXrrCwsNC5f8YY3Nzc4ObmhqFDhwIAEhISEBwcjLVr12LUqFHo2rUrRo8ejd69e8PU1FTnMXE4lYVKfXcyJiYGY8eORePGjXH37l2sWbMGycnJ+OWXX9CnTx+9CJgq6tWrh6lTpyI0NBQPHz5Enz59sGzZMjRp0gQ//PADsrKyRIuNwzEmKp2IyeVyBAUFwd/fH4GBgWjYsCHu3buHbdu2oXPnzgY5D2VjY4MxY8bg6tWr2LNnD65evYoGDRpgypQpuH//vtjhcTgGTaUSsQsXLsDX1xdffvklJk+ejNjYWMyfPx92dhrX6BMFxhjat2+PvXv34ubNm7CyskKHDh0wdepUpKWliR0eh2OQVAoRu3v3Lvr3749Ro0ZhxowZuHr1KoYMGQKJRCJ2aFrj5OSE5cuX459//gFjDM2aNcOqVauQl8cfUcnhlMaoRSwrKwtTpkyBn58f/Pz88M8//2DYsGEwMTHq03oJW1tbrF27Fn///TfCwsLQtGlTHDhwAPx5oRyOAqP9aw8NDUXr1q1RVFSEO3fuYNasWahatarYYemMJk2a4MiRI9i9ezeWLVuGoUOHIj09XeywOBzRMToRy8vLw2effYYRI0Zg48aN2LRpE2xtbcUOS2907NgRERERqFu3Llq3bo1Tp06JHRKHIypGJWLR0dHw9vZGYmIirl+/jl69eokdkihYWFhgzZo12LVrFyZOnIjJkyfjxYsXYofF4YiC0YjYrl270KNHD8yfPx8HDhyAjY2N2CGJTkBAAGJiYpCdnY0OHTrg0aNHYofE4egdgxcxuVyORYsWYfHixTh//jyGDx8OxpjYYRkMlpaW2LVrFz766CN06NAB4eHhYofE4egVw1v5WYq8vDyMHj0a8fHxCAsLQ506dcQOySBhjGHGjBlo2LAhevfujY0bN2LQoEFih8Xh6AWDFbEnT56gX79+cHFxQWhoqKhbhIyFfv36oV69eujbty/u37+POXPm8F4rp9JjkMPJ1NRU+Pv7o0uXLti7dy8XMA3w9PREWFgY9u7diwULFvD1ZJxKj8H1xFJTU9GlSxcMHToUS5YsETsco6Ru3boIDQ1FQEAAGGNYsWIF75FVEJlMhidPniA9PR1SqRRyuRxmZmawsLCAo6MjatSoIXaIryU/Px85OTmQSqUgIkgkElSvXr1SlIAyKBHLyMhA9+7d8d5773EBqyC2trY4e/YsAgICULVqVSxatEjskIyC/Px8xMTEIDw8HOHh4bhz5w6SkpLw5MkTWFlZwcbGBubm5mCMQSaT4cWLF0hOToaZmRmcnJxQv359eHp6ol27dvDx8YGzs7PeLiAymQy3b98uqVv38OHDkgKdOTk5qF27dkkBBKlUipycHFStWrWk6Kerqyu8vLzg4+ODVq1aGU25dYMRsRcvXiAwMBDdunXD0qVLxQ6nUmBnZ4eQkBB07NgR1tbW+OSTT8QOySB59OhRSaHMq1evokmTJvDx8YGfnx8mTZoEZ2dnODg4qNyLS0TIyspCUlISYmNjERkZiW3btuHjjz+GmZkZ3n33XQwYMABdunQRXBiUsR89ehRXr16Fk5MTfHx84O3tjT59+pQIlI2NzX+24xERMjIySoTuwYMHuHr1Kn7++Wc8fPgQbdu2Rd++fTFgwAA0adJE0LgFhYjUfnl5eZEukMlk1KdPHxo9ejTJ5XKd+HiTefjwITk7O1NQUJDOfQGIIA1yinScW6rIzs6mn376iTw9PcnGxoY++ugjOnr0KL148UIwH3K5nO7evUvffPMNdejQgSwtLWnEiBF05cqVCtlNT0+nb7755qXYg4ODKSMjQ6DIiXJycuj06dP08ccfk6OjIzVv3pyWLFlCiYmJgvnQFFW5ZRCJtnTpUurYsSMVFhbqxD6H6OrVq2RnZ0d37tzRqR9DF7H79+/T9OnTydramgYPHkx//PEHFRUV6cV3UlISfffdd9SgQQPy9fWlPXv2UEFBgdrHX79+ncaPH0+WlpY0atQoCg0N1UvsMpmMrly5QpMnTyYrKysaNmwYXb58We8dDoMVsd9//52cnZ0pOTlZcNucl9m6dSs1bdqUsrOzdebDUEUsMTGRxo4dS7a2tjR37lyKj4/Xqb/XIZVK6ejRoxQQEEANGjSg/fv3v1YQbt68SYGBgeTk5ETLly+n1NRUPUb7MpmZmbRmzRpyc3MjPz8/unz5st58G6SI3bt3j+zs7Ojvv/8W1C5HNRMmTKCBAwfq7CpqaCKWk5NDS5YsIWtra5o9e7agQy4hCA0NJU9PT/L19aVLly699FlycjJNmDCB7OzsaM2aNZSfny9SlP9FKpXS9u3bqW7dujRkyBB68OCBzn0anIjl5ORQy5YtacOGDYLZ5JRPfn4+tWvXjlauXKkT+4YkYufPnydXV1caNmwYxcbGCm5fKGQyGe3atYvq1atHEydOpOzsbNq4cSPZ2NjQrFmz6NmzZ2KHqJIXL17Q8uXLydramlasWKHT4a3BidiQUePJp1s/Co9NF8wmRz0SEhLIzs6Odh87T+tC71NEnHB/JGKLWETcM1pz6iaNmDCFHB0d6dixY4LY1QeZmZk0ZMgQqlatGrm7u9Pt27fFDkltHj16RF27dqUOHTrQ/fuKnNJXbomyYn/bkVM4fPgwUlu8jw9/CUNkfIYYYbyxODs7Y/LsJRg7bhy+O3Ubw7deqRTfQWR8Bt77NghzRvRB0IUo/O/YBQQGBoodltpcvHgRFy5cQO/evZGZmYmgoCBFT8MIqFevHkJCQjB06FB4t/NFnynL8H3IXb3klt5FrKCgALM//QTWXSeAVa2JIqkcVx7yCqX6xtazB1jVGsi8GlxpvoO9v53Box2zUL1NT9j2n4e7WcaxS4GI8P3332PSpEkICgrCwYMHERUVhaCgIIwaNQoFBQVih6gWJiYm+PTTTzHl+914+tevSA/dhsIiqc5zS+8iNnX2IuRXqwML97cBAKamJmjvxmuD6ZsODW3hGDgV2WGHQdkpRv8d7Nu3D1uXTIZNr09Ry+tdEBisqhn+ivOioiJMmjQJ//vf/3D58mV06NABgOJBMRcuXMCLFy/QvXt3PH36VORI1WdAQHu4jlmDwuR7SAteiTYOui0br1cRu3nzJvbt3Arr7pPAGAMDMNjLGV4uVvoMgwPAy8UKBz8fgH6jJsMmeic861uKHZLWrFu3DnPnzsUn3+5EtYbeABSJnZFbKG5g5ZCXl4e+ffsiMTERf/75J+rXr//S59WqVcPBgwfRoUMH+Pn5ISUlRaRINcPLxQr7p3bH11sPoEsrV3w2cgCePXumM396FbFZs2bh4xlzUc2qDkwZUEVigkGezvoMgVMKLxcr/PrTckjzchAUFCR2OFqxfv16fPfdd7hw4QIGdXsLVSQmMGWAucSwe/j5+fno168frKysEBwcjJo1a5bZzsTEBKtWrcLw4cPRpUsXpKam6jlS7fByscKnPZrj6K974O/vj+7duyMjQ0dzY2XN9qt6VeQO0sWLF8nV1ZUKCgp0cueCox0Rcc9o0lebya1JU5JKpRW2Bz3endy1axc5OzvTw4cPS94zhtwqLCykd999l4YOHarRkoTFixdT27ZtDW6tmyqU30V4bDpNnz6d2rdvTzk5OVrbU5VbehExuVxOnTt3pu3bt2t9AhzhiYh7Ru4LT5DrnN+pat2mtOLHLRW2qS8RCw0NJXt7e7p161aFY9YncrmcPvroI+rdu7fG2+zkcjlNnTqVAgIC9LZVSluUudVg7jFyX3iCwmPTacyYMdSjRw+tY1eVW3oZToaGhiI5ORkffvihPtxxVBAZn4Gfz/2LvWGP8PO5f3Ek6jEKpXLFJHjHD7H2268glUrFDrNcEhIS8MEHH2D37t1o3ry52OFoxLp16xAdHY2DBw9q/IR6xhjWrFkDU1NTzJ8/X0cRaoeq3JITUCSVIyz2GTZt2gS5XI4FCxYI6lvnpXiICAsXLsTSpUtLahlx9E9kfAaGb72CgiI5CIAJA8xMGMxMTSCTyVGzoQdqxDph9+7d+Oijj8QOVyX5+fkYNGgQPvvsM3Tr1k3scDTi0qVLWLFiBS5fvqx1MUJTU1Ps27cPPj4+8PHxwXvvvSdwlJpTXm5JzBTzk2ZmZiWxe3t7Cxa7zlXl7NmzyM7OxtChQ3XtivMarjxML+51KZATIJMThrarh7qWFmjvZoPcnisxevRojBw58j+1pwyFGTNmwMXFBZ9//rnYoWhEUlIS3n//fezcuRNubm4VsmVjY4PDhw+jR48eaN26Ndzd3QWKUjvUyS3lCgRbW1scPnwYPXv2RIsWLQTpSetcxDZu3IipU6ca7B/Fm0J7NxuYm5mgsEgOORRXS4mZ4u6wMsGovh9q1qyJs2fPonv37uIGXAbnzp3D77//jps3bxpVuW0iwqRJkzBmzBi88847gtj08PDAokWLMH78eJw/f17Uvy91cqs0np6e+PLLLzFq1Chcvny54iO0sibKVL00nXxNTk4mS0tLysrK0moijyMsyrtFe67Eq7yDt2HDBho8eLDWPqCjif3c3Fxq1KgRHT16VOvYxOLAgQPUvHlzwatQSKVSat++vUEUUVAnt0ojl8upS5cu9P3336vtQ1Vu6VTEVq5cSePGjdPoGI64ZGVlkaWlJaWkpGh1vK5EbM6cOTRkyBCtYhKT9PR0cnBw0Fm5qZs3b5KtrS09fvxYJ/Z1yf3798nGxkbtMj6qcktnfVC5XI4tW7ZgwoQJunLB0QG1atXCoEGDsGPHDrFDKeHBgwfYunUr1q5dK3YoGjNv3jwMHjy4ZDuR0LRo0QITJ07EvHnzdGJflzRq1AizZ8/GlClTKmaoLGVT9dKkJ3bmzBlq27Ytr5lvhISFhZGbmxvJZDKNj4UOemIjR46kJUuWaByL2Ny7d49sbGx0Xg8sMzOTbG1tdV56XBcUFBSQq6srXbx4sdy2qnJLZz2xAwcOYOTIkUY1ActR4OPjA4lEgqioKLFDwT///IMTJ07gs88+EzsUjfniiy/w6aefwspKt3uDa9eujZkzZxrlU8LMzc2xePFiLFq0SDG/pQU6ETEiwvHjx9GnTx+NjlMumKsMta2MGcYY+vTpg2PHjokdCr744gvMnDkTtWvXrpAdfefW7du3ERISgk8//VQv/qZOnYrz58/j+vXrevEnJCNGjEBSUhJCQ0O1Ol4nIhYdHY0aNWqgcePGah+jXDCnr0JqnNcTGBiI48ePixpDYmIiTp8+XeHnZYqRW19//TVmzJiBWrVq6dwXAFSvXh3Tpk3DTz/9pBd/QmJmZoYlS5ZgxYoVWh2vExE7fvy4xhU1lQvmlNsUKkORPmPGz88P//77r6jlX7Zt24ahQ4eqrPCgLvrOrWfPnuG3337D+PHjdernVcaMGYNDhw4hOztbr36F4L333sPt27dx7949jY81GBFTLpgzLV4o197N5qUhAB9q6heJRIJu3brh5MmToviXyWTYunUrJk6cWGFb+s6tXbt2ITAwEDY2+i0F5ODggK5du2Lfvn169SsE5ubm+Oijj7BlyxaNj2WaTKZ5e3tTRETEa9s8ffoUDRs2xJMnT1ClShWNgomMz8CVh+kldaCGb72CQqkcZqYmABGkcoK5mQn2jGvPCynqgR07duD48eM4ePCg2scwxiKJyFtTX6/m1smTJ7FkyRJcvXpVU1Nloq/cIiK0bNkSGzZsQKdOnQSJXRNCQkIwb948REZG6t13Rfn333/x1ltvISEhoUztUJVbgvfEIiIi4O3trbGAAYpCap90aQQvF6v/DAGKZMSHmnrGz88PYWFhovg+cuQIPvjgA8Hs6Su3oqKiUFhYiI4dOwoUuWZ069YNCQkJiI+PF8V/RWjUqBFatGihce9fcBGLioqCh4dHhe28OgSQmLKXhgMc3ePm5oasrCy913cnIpw4cUJnTyrSZW4dO3YMffv2FW1pkYmJCXr16iX6TRlt6devn8axC74BPDo6GgMGDKiwHS8XK+wZ1/6lIYDy/3woqR9MTEzQtm1bREdH63VD+LVr11C9enWN7m5rgi5z6/jx4/j666+FCFNrAgMDsXPnTkyePFnUOLQhMDAQ33zzDYhI7QuB4CIWFRWFZcuWCWLLy8XqpaTi4qV/PDw89C5iuuyFKdFFbqWkpODevXvw8/OrsK2K0KNHD4wbNw55eXmwsLAQNRZNady4MWrUqIHo6Gh4enqqdYygw8msrCykpqaiSZMmQprliIhSxPTJ33//jc6dO+vVpxCcPXsWAQEBMDcX91FxlpaWaNKkCa5duyZqHNrSq1cvhISEqN1eUBG7fv06WrVqBVNTUyHNckREDBETal5V30RGRqJdu3ZihwFA8b0ZwrYxbWjXrp1GsQsqYo8ePUKDBg2ENMkRmQYNGiAhIUHrfW2akpycjMLCwv88g9EYiI6ONhjx9fT01PvFRyg0vXAKKmLJyclwcHAQ0iRHZGrUqAEAeP78uV78KYXA2AoHEJFBiZgx98Tc3d2RnJyMrKwstdoLKmIpKSlwdHQU0iRHZBhjcHR0RHJysl78/fvvv2jatKlefAlJfHw8qlevjjp16ogdCgCgadOmuH//vthhaIWpqSlatmyp9mZ23hPjlIuDg4Pe9lAmJycb5YUwISEBrq6uYodRQu3atSGVSvHixQuxQ9EKV1dXPH78WK22vCdmRIi1f9TR0VFvImasOZSSkmJQF3BlD1qd780Q9yVrcuEUdJ1YSkoK7O3thTTJKUZZTqZQKtf7/lEHBwe9DSeNtTdviD1I5ffWsGFDlW3EzKvXockUhqA9sby8PFSvXl1Ik5xiSu/3KyiS40iUel1tIahWrRry8vL04is9PR22trZ68SUkhtYTAwA7O7tyt4yVzqtCqRw//HHPIHpkmlw4BRUxqVTKn/KtI9q72SgqLgAgAAcjEvSWbGZmZpDJZHrxVVRUJPpiUW3IycmpcPVZoZFIJJBKpa9to9xHasIUD7398/5TgyhKWqtWLeTk5KjVVnAR4wtddYOXixUGezlDufBAJie9VfMwMzNDUVGRXnzJZDKjfNCyIV7ATUxMyhUx5T7StxvZgkFxgTSESjGaXDgFzRZDWttjiJOVFWWQpzOqSEz0Xs2DiPQmLKampnrr9WlLWbmlr8XAmiCTydQSVi8XK0zv1kSU3FIFYwxyuVyttoJeOtTpvuqCvWGPcPJmMnq1dMQHvvUNdrKyorxafUFf5ySVSlG1alW9+JJIJCgsLNSLL3VQN7f0OeRWl6KiIkgkErXaipVbqtCkZyuoiOlz2KFkb9gjzA+6AQC4dF8xiZmRW/ifmupifylC8Wr1BX1QVFSkt2kCOzs7pKWl6cVXeWiSW1WrVkVubq6Y4f6HJ0+ewM7OTu32YuSWKvLy8tS+cAo6RrC0tERGhn6HbydvJv/n57JqqnO0JyMjA9bW1nrxpc81aeWhSW7Z29sjNTVVjDBVYqxr7gAgNTVV7bu9gvbE9Lk9RUmvlo4lV0nlz4bWNTZ29Ll2S59r0spDk9xydHQ0qA3XRGSQyz7URZN1d4KKmBgJ+IGvotpB6XkLwLC6xsaOPhdyOjo64tatW3rxVR6a5JY+t2apQ0ZGBqpUqWJ0RRGVpKSkoFmzZmq1FbwnJsYX+YFv/ZIE4wiPPq/oTZo0waFDh/TiSx3UzS0nJye19/rpg9u3b8Pd3V3sMLQmMTERTl3qBMEAACAASURBVE5OarUVdE7MkIYCxoAxLAORyWRIS0vT23YyDw8PXLt2Te3b64ZCo0aN8PjxY4PZcH307F8wd2ho0LmlCiLC9evX0bJlS7XaCypiYsyJGSvKW/Xfh9w1iBXSqkhLS0Pt2rXVvlVfUezs7FCzZk3ExsbqxZ9QmJubo1mzZmqXj9ElkfEZ2HzkD/xTaGPQuaWKxMTEkg3s6iCoiDVu3Bh37twR0qTg6KP3o46PV599KPYKaVXcuXNH78MST09PoyvoFxmfgSoODXH07F8696NObuWnPIDEvqFB55YqNC2MKeicWNOmTZGYmIjs7GzUqlVLSNOCoM0i2NJPjlbnRoG6PpS36oukcoNeBiJGtVI/Pz+cO3cO7733nl79aovyO39aZIPbh87gvRFj1c4VXeSWe22CNCsVFnUaGHRuqULTZywI2hMzMzNDy5YtERMTI6RZwSir9/O6K5s2Qz5VPaxX/Shv1c/o4Y4949oDgEHOj0VFRan96Cyh6N27N44fP26QW3nKQvmdm9drheex0bj8QLEsQ6zcSr51BZ06d8HM3i0NOrdUERISAn9/f7XbC75jVfmAArEe4/46Xu39WFUzf+2V7UjUYxQUyUEACosUZUqmd2ui8qoZGZ+BxMw8mJmaQCb7/x7W3rBHWHz0JuREL/lRvgx5m1R0dDRmzpypV5/NmzcHYwy3bt1Se3JXTJR5Bbv6YIzBrugJIuNtVX6nkfEZ+OGPe/9fAkfg3Dp+/DiGv9cf47o0MujcKov09HTcvHlTo0f2CS5iHh4euHz5stBmBeHVhYplXdlKJ9rBiAQo+wJyAH/9+xThcc/KTITSyWJmwvB+u/oY6OkMAFh89CakcoWlwlJXUHXiEJPc3Fw8fPgQLVq00KtfxhgCAwNx/PhxoxCx0nl1Ibkv4q79iafm9mV+p8o8UV4cGYTNrb/up+LkqdNoNWhqyXDVEHNLFadOnUKXLl002qsreGkCT09PRERECG1WMLxcrPBJl0bwcrF67fakKw/TS5JDibJwXFkTpaWTRSYnOFlawMvFClcepkNWyo4JYyU9QOVQwqqauUFuk4qJiUHTpk1Fqe81ePBg7N6922iGlMq8Gv3+QPz+++8qc0uZJwTFH5+LTTWFkAmUW3HRf6KwpiN+ico06NxSxbFjxzR++rvgPbE2bdrg8ePHBlmu91Vetz2p9NCTMUBavGxJToBVtf/+UZdub2pqgsTMPETGZ6C9mw2qSExQWCSHiQnDsn4t/7OJOCO30CC3SYWEhKBLly6i+Pb390dBQQGuXLmCDh06iBKDNgQEBGDUqFGoXZRe5nf66pTGOy0csPHiQwDC5NZXn36DGm3eMfjcKouMjAycPHkSa9eu1eg4psmVztvbm9TpZQ0ZMgQ9e/bE2LFjNQrG0FB2x5My87A37FHJ1XNmT3d80qXRS22UV7gjUY9xMCIBUvn/z1EAeCmJlMMDZSIb6jxFu3btsHLlSnTt2lXtYxhjkUTkramvsnLru+++w61bt7B9+3ZNzYnKrFmzYGZmhq+//rrMz0vnzJWH6fju9F1BcsuO5aB1m7aoM3EbZExi0LlVFuvWrcOff/6J/fv3l/m5qtzSSSnKPn36IDg42OhFrPTE++GoxyiUykuGg0DZt7ydLC0gldNLcxDK4Wtpu4Z+dUxNTcW9e/dEvUEzatQoNG7cGKtXr4aVleH9jlQxfvx4dO7cGcuWLStzKP7q3ssqEhNBcmvRotUYOeJDjPq4s0HnVlkQETZv3owffvhB42N1Uq6zV69eOHv2LAoKCnRhXu94uVhhcZ8WMGEMMjlh2bFbKidN1S0DVHpuzhA5efIkunXrJmq9ezs7OwwePBirV68WLQZtcHd3R7NmzRAUFFRuW6FyKyMjAxs3bsTUqVMNPrfK4vLly8jLy9Nq+kInImZnZ4dmzZrh0qVLujAvChm5hZATvVSDvKykenX9lzElUmmOHz+u8QSrLli0aBHWr19f7lN7DI2ZM2dixYoVau0BFSK3Vq9ejb59+6Jx48Y6PCvdsWzZMsyYMUO7EvdEpPbLy8uL1GX58uX0ySefqN3e0ImIe0buC0+Q29xj5L7wBEXEPSt5f13o/ZKfKwMvXrwgS0tLSk5O1vhYABGkQU6RGrk1efJk+vzzzzWORUzkcjn5+vrS/v37y21b0dxKS0sjGxsbio2NFSJ0vXPp0iVydXWlgoKC17ZTlVs6E7G4uDiytramFy9eaHpOBktlFKyy2LFjB/Xu3VurY3UhYo8fPyZra2uKi4vTKiaxCAkJIXd3dyoqKiq3bUVy69NPP6WPP/5YmxBFRy6XU+fOnWnbtm3lttW7iBER9erVi3bu3KnRMRzxeeuttyg4OFirY3UhYkREX375Jb3zzjskl8u1iksM5HI5derUiTZt2qQzH2FhYWRvb09paWk686FLTpw4QY0bN1ZL6EURseDgYHr77bc1OoYjLjdu3CAnJye1kqosdCVihYWF1Lp1a9q9e7dWcYnFtWvXyM7OjpKSkgS3XVBQQK1ataI9e/YIblsf5OTkkIuLC50+fVqt9qpyS6cPEwwMDERsbKzBlBvmlM+WLVswduxYg3sQrEQiwdatWzFz5kyDeRqSOrRp0wbjx4/HlClTBLf9zTffoF69ehg2bJjgtvXBggUL0LlzZ/To0aNihspSNlUvTXtiREQLFy6kadOmaXwcR//k5uaSjY1NheaeoKOemJJ58+ZR9+7dSSqVah2jvsnLyyN3d3c6fPiwYDYvXLhA9vb2FB8fL5hNfXL58mVycHCgp0+fqn2MqtzSuYglJCSQtbU1paamanwsR7+sXr2a+vfvXyEbuhaxoqIi6tq1K82bN69CceqbP//8k+zt7QW5g/j48WNycnJSexhmaKSnp1ODBg3o0KFDGh0nmogREU2ZMoVmzJih1bEc/ZCTk0P29vYUExNTITu6FjEioidPnpCLi4ugPRt9sGbNGvLw8KDc3FytbeTn51P79u3pq6++EjAy/SGVSqlnz540c+ZMjY8VVcSSkpLIysqKHj9+rNXxHN2zcuVKGjp0aIXt6EPEiIjCw8PJzs6Orly5UuGY9YVcLqcPPviARo4cqdVdVplMRsOHD6cBAwYY1V3a0ixYsID8/f21unEkqogREc2aNYsmT56s9fEc3ZGZmUm2trZ0586dCtvSl4gRER07dozs7e3p+vXrFY5bX7x48YJat25NX3/9tUbHyeVy+uSTT6hTp04V6smJya5du6h+/fpaTy2JLmJPnjwha2tro11VXJlZsmQJjRw5UhBb+hQxIqL9+/eTo6Mj3bhxo+LB64nHjx+Tq6srbdy4Ua32crmcpk2bRj4+PpSZmanj6HTD0aNHyd7enm7duqW1DdFFjIho0aJFNHjw4ArZ4AhLfHw82djY0L///iuIPX2LGBHR3r17yd7eni5dulSx4PXI/fv3ydnZmX755ZfXtisoKKDRo0eTj48PZWRk6Ck6YTl27BjZ2dnR1atXK2THIEQsNzeXGjVqpPVqcI6wyOVy6t27Ny1fvlwwm2KIGBHR6dOnyc7Ojnbt2lUhO/rk7t27VLduXVqzZk2Zc1zp6enk7+9Pffv2pefPn4sQYcXZt2+fYHOXBiFiRETnzp0jZ2dno+0WVyb27t1LLVu2LHfjrSaIJWJERDdv3iRXV1datGiR0awji42NpRYtWtDHH3/80mT33bt3qUmTJjRz5kyjOZfSyOVy+uKLL6h+/foVvuOtxGBEjIho/PjxNGnSJEFscbQjLS2N7O3tKSwsTFC7YooYEVFKSgp17tyZOnXqZDTzr5mZmdSzZ0/q0aMHZWRk0IYNG8jW1pY2b94sdmhakZeXR8OHDycfHx+tKqGowqBELCMjg5ycnOjixYuC2ONozogRI2j69OmC2xVbxIgUa5G+/fZbsrW1pW3bthnFcoSioiIaM2YMWVhYkLu7uyB3isUgJiaGWrVqRUOGDBG8go1BiRiRYnN4gwYNKD09XTCbHPXYt28fubm5UU5OjuC2DUHElFy/fp3atGlDXbt2pejoaMHtC0V+fj599913ZGtrS++//z45ODjQ7NmzKT8/X+zQ1EYqldKqVavI1taWduzYoZMLh8GJGBHRjBkzqGfPnkY55jdWrl+/Tra2tjr7ozYkESNSVL9Yv3492dvb06hRoyghIUEnfrRBLpfT/v37qUGDBvTuu+/S7du3iUixHGnAgAHUokULOnPmjMhRlk9YWBi99dZb1LlzZ50O4Q1SxIqKisjf358WLFggqF1O2WRkZFDDhg11egfP0ERMSVZWFi1YsICsrKxo9OjRovbMXrx4QZs3b6aWLVuSp6cnhYaG/qeNXC6nw4cPU6NGjeidd94xyHVwDx8+pPfff5+cnJxo69atJJPJdOrPIEWMiCg1NZXq1atHQUFBgtvm/D8ymYwCAwNp6tSpOvVjqCKm5MmTJ/TVV1+Rs7MzdezYkfbs2UNZWVk69yuXy+n69es0e/ZssrW1pXfffZfOnDlT7rCroKCAfvzxR7Kzs6Phw4dXeK2VENy6dYsmTZpE1tbW9MUXX+ht+YfBihiRojtqZ2dX0p3mCM+SJUvIz89P0OUUZWHoIqaksLCQfv31V3rnnXeoZs2a1KtXL9q8ebOgw838/Hy6dOkSzZw5k9zc3MjFxYVmzpyp1cLijIwM+vbbb8nFxYXat29Pe/fu1eucWVFREf3222/UvXt3sre3p8WLFwt651EdDFrEiBR13evXr0+PHj3SmY83lU2bNpGrq6tOqou+irGIWGmysrJo3759NGTIELKxsSFHR0fq27cvLV++nIKDgyk8PJwSExNVblrOycmhf/75h0JDQ2nbtm00efJk8vb2JgsLC2rbti0tXryYoqKiBJnslkqlFBQURAEBAVS7dm0aMmQI7du3Tye9ydzcXAoODqZRo0aRjY0N+fr60q5du0S74aAqt3TyBHBtWb16NTZv3oyLFy+iTp06OvPzJnHgwAHMmDEDFy5cQKNGjXTuT8gngIsBESEuLg7h4eEIDw/HnTt3kJSUhKSkJKSnp6N27dqQSCQwMTGBVCpFbm4uZDIZ6tatCycnJ9SrVw9eXl7w8fGBh4cHqlevrrNYU1NT8fvvvyMoKAiXLl1CixYt4OPjU/Jq2LAhJBKJWrZkMhni4+NLzjs8PBzR0dHw9vbGgAED0K9fP9SvX19n56IOqnLLoEQMAJYsWYKgoCCEhobC1tZWp74qO0FBQfj4448REhKC1q1b68WnsYvY6ygqKkJmZiakUilkMhkkEgmqVq2KWrVqafe8RAF5/vw5IiIiSgQoIiICCQkJsLGxKRFYS0vLkrLjUqkUOTk5SEpKQmJiItLS0mBvb18iwMqXIT15XVVuGVYhdQBLly5FUVERunbtirNnz3Ih0xKlgJ08eVJvAlbZkUgksLOzEzuMMqlRowb8/f3h7+9f8p5UKsWTJ09KhCo7OxtSqRREBIlEgurVq5cInIODg9q9NkPD4ESMMYYvv/wSjDH4+fnh2LFjehkGVSY2bdqEJUuW4MSJE/Dw8BA7HI5ImJmZwcnJCU5OTvD21rhzbDQYnIgB/y9k9erVg5+fHw4ePIiOHTuKHZbBI5PJMHv2bBw7dgx//vknF3/OG4FOH9lWUSZNmoT//e9/GDRoEHbv3i12OAbN8+fPMXDgQERHR+Py5ctcwDhvDAYtYgDQo0cPnDt3DosWLcKiRYsgl8vFDsngePz4MTp16gRbW1ucOnUK1tbWYofE4egNgxcxAGjRogXCwsJw8eJFBAQEID4+XuyQDIb9+/fDy8sLw4YNw9atW2Fubi52SByOXjEKEQOAOnXqIDQ0FL1794a3tzd27NgBTZaHVDaePXuGYcOG4YsvvsDx48fx+eefi36bn8MRA6MRMQAwNTXF7Nmz8ccff2D16tUYNGiQUT3SXihCQkLQpk0b1KlTB1FRUZX6zhOHUx5GJWJK2rRpg/DwcDRu3BitWrXC+vXrUVRUJHZYOufRo0cYMWIExo0bh+3bt+PHH3+EhYWF2GFxOKJilCIGAFWqVMGqVatw6tQpHDlyBK1atcJvv/1WKYeYWVlZmDdvHjw8PNCgQQPcunUL3bp1EzssDscgMFoRU9K2bVucOXMGa9aswfz58+Hv74/w8HCxwxKEwsJCrFu3Du7u7khNTcX169exbNky1KxZU+zQOByDwehFDFAsju3VqxeuXbuGESNGoH///ggICEBwcDBkMpnY4WnMkydPsGLFCjRo0ADHjh3D6dOnsW3bNtStW1fs0Dgcg6NSiJgSMzMzjBs3DrGxsRg/fjxWrVqFRo0a4bvvvkNGRobY4ZVLVFQUPvroI7i7uyM+Ph4nT57EqVOn0KZNG7FD43AMlkolYkrMzc0xbNgwXL58GQcOHEBMTAzc3NwwbNgw/Prrr8jJyRE7xBLu37+Pb775Br6+vujfvz+aNWuGf//9F1u2bOEbtzkcNaiUIlaadu3aYdeuXbh79y66dOmC7du3o27duggMDMTWrVuRmpqq13jkcjkiIyOxcOFCtGjRAp07d0ZsbCyWL1+Ohw8fYs6cObCxsdFrTByOMWNw9cT0QVZWFk6ePImgoCCcOnUKNjY2JfWT2rVrB09PT9SoUaPCfogIjx49eqnQXGRkJBwcHNC/f38MGDAA7dq1g4lJ5bmWVOZ6YhxxMZqiiPpGLpfj7t27LwnNjRs3YG9vDycnp5J6S8p/a9asCTMzM0gkEsjlckilUhQWFiItLa2kAmhiYiKSkpLw6NEjmJmZvVRkzsfHp1LXSOMixtEVRlMUUd+YmJigWbNmaNasGUaOHAlAsbTh0aNH/xGliIgIPH/+HFKpFEVFRTAxMYFEIoGZmRns7Ozg5OQEDw8PBAYGwsnJCc7OzrC3t+fbgTgcHfLGi1hZmJubo1GjRrycDYdjBFSeyRgOh/NGwkWMw+EYNVzEOByOUcNFjMPhGDVcxDgcjlHDRYzD4Rg1XMQ4HI5Ro9GKfcZYGgD+lA7O63AhIo0fk81zi6MGZeaWRiLG4XA4hgYfTnI4HKOGixiHwzFquIhxOByjhosYh8Mxat44EWOMxTHGdPq8M8bYDsbYCl364BgePLfE4Y0TMbFhjHVhjN1gjGUyxtIZY0GMMf4YI46gMMa2M8aIMVbp60lxEdM/twH0JCJLAE4A7gPYIG5InMoEY8wPQEOx49AXb7SIMcZMGGNzGWMPintFvzLGrIs/O8UYm/JK+xjG2MDi/zdljJ1hjD1jjN1ljA1RxycRpRJRUqm3ZAAq/dXyTUOM3Co+1gzATwCmlNe2svBGixiAaQD6A+gMRa8oA8DPxZ/tBTBM2ZAx1hyAC4DjjLHqAM4Ut6lT3G49Y6yFOk4ZY/UZY5kA8gDMAvCNIGfDMSREyS0AnwG4SETXhTgJY+BNF7GJABYQ0WMiKgCwFMDg4qtZEIC2jDGX4rbDARwpbtcHQBwRbSciKRFFATgMYLA6TonoUfFw0hbAQgD/CHpWHENA77nFGKtX7Hex8KdjuLzpIuYCIKh4kj0TwB0ohnf2RJQD4DiA94vbvg9gT6njfJXHFR87HICDJs6J6BmAnQCOFic3p/IgRm79AGAZEWUJeSKGzpv+h5MAYAwR/aXi830AljDGLgKwAHCu1HEXiKi7ADGYQTFsqAXgmQD2OIaBGLnVFYAfY6z09MRlxtinRLRXC3tGwZveE9sI4Etlt54xZscY61fq8xNQXBmXAThARPLi948BaMIYG8EYkxS/fBhjzcpzyBgbyBhzL574tQOwGkB0ca+MU3nQe24BaAKgDYC2xS8AeBeK4Wul5U0XsR8B/AYghDGWA+AKAF/lh8VzFEcAdINiolX5fg6AHlAMA5IApABYBaCKGj7rAjgFIAfADQByAAMEOBeOYaH33CKiJ0SUonwVv/2UiPKEOSXDhJfi4XA4Rs2b3hPjcDhGDhcxDodj1HAR43A4Rg0XMQ6HY9RotE7M1taWXF1ddRQKpzIQGRn5VJsHhfDc4pSHqtzSSMRcXV0REREhXFScSgdjTKsnFvHc4pSHqtziw0kOh2PUcBHjcDhGDRcxDodj1HAR43A4Rg0XMQ6HY9RwEeNwOEYNFzEOh2PUcBHjcDhGDRcxDodj1HAR43A4Rg0XMQ6HY9RwEeNwOEYNFzEOh2PUcBHjcDhGDRcxDodj1HAR43A4Rg0XMQ6HY9RwEeNwOEYNFzEOh2PUcBHjcDhGjUYPCuGIS1FREW7duoXIyEg8ffoU+fn5yMvLg1QqRdWqVVG1alVUq1YNDRs2hI+PD5ycnMQOmcPROVzEDBS5XI779+8jPDy85HX9+nXUr18f3t7esLe3h4WFBWrXrg0zMzPk5+cjNzcXT58+RUhICMaMGYOqVavCx8en5OXt7Q1ra2uxT43DERQuYgbGs2fPsHXrVvz8889gjKFdu3bw8fHBgAED4OnpiVq1aqllh4gQGxtbIoBffvkloqOj0bFjR0ybNg3du3eHiQmfTeBUAohI7ZeXlxdxdMP169dp/PjxZGlpSSNHjqTw8HDBfeTm5tK2bduoTZs25O7uTuvWraPs7GxBfQCIIA1yinhucdREVW7xS7GIEBF+++03BAQEoGfPnqhXrx7++ecf7Ny5E97e3oL7s7CwwOjRoxEdHY3Nmzfj3LlzcHFxwfTp0xEXFye4Pw5HH3ARE4nU1FT069cPCxYswLhx4xAXF4dFixbB3t5e574ZY+jUqRMOHTqEa9euoUqVKvD29sbGjRuhuOBxOMYDFzERCA4ORtu2bdGyZUtERkbigw8+gLm5uSix1K9fH6tWrcKlS5ewdetWBAYGIjk5WZRYOBxt4CKmR7KzszF69GjMmjULhw4dwldffSWaeL1Ks2bNcPnyZXh7e6Nt27Y4ePCg2CFxOGrBRUxPXLx4Ea1bt4ZEIsG1a9fw9ttvix3Sf5BIJFi2bBl+++03LFiwAB9++CEyMzPFDovDeS1cxPTAwYMHMWTIEKxbtw6bN29GjRo1xA7ptfj6+iI6OhrVqlWDv78/0tLSxA6Jw1EJFzEds2fPHkybNg2nT59Gnz59xA5HbapXr45NmzahT58+8Pf3R0pKitghcThlwhe76pDdu3djzpw5OHv2LJo3by52OBrDGMOKFStQpUoVdO7cGRcvXtTL3VMORxMMSsSys7MRGRlZsso8JiYG2dnZyM/PB2MMFhYWsLOzg5eXV8lWmtatWxvM5HhpTp48iVmzZuHcuXNo1qyZ2OFUiEWLFkEqlaJ37944f/48atasKXZIHE4Jog8ni4qKcODAAfj5+cHJyQkLFy5EUlIS+vfvj+DgYERHRyM2NhYPHjxAeHg4tm3bBl9fX0RGRmLMmDGwsrLCe++9h0uXLhnMGqewsDCMGjUKwcHBRi9gSpYuXQpvb28MGjQIhYWFYoejkpycHGzatAlDhw6Fm5sbqlSpAolEAnt7e/Tq1QvLly/Ho0ePxA6TIyRlLeNX9RJya8iZyHvUZ8xnVMfBkfz9/enIkSNUWFiosZ2srCz66aefqEmTJuTh4UHbt2+nvLw8weLUlKdPn5KjoyP9/vvvosWgLhFxz2hd6H2KiHumVvuioiLq27cvffrppyrbQKRtR7m5uTR37lyqbWlFbTr2oKXfr6fbt29TXl4eFRQU0OPHjyk4OJimTp1K1tbWNHjwYIqPj6+QT45qNM0tdVCVW3oXMblcTou+WUumFrWoRuse5DJ+nSAnKpPJ6OTJk9SzZ09q2LAh/fXXXxW2qQ0jR46k6dOni+JbEyLinpH7whPUYO4xcl94Qu3v4OnTp+Tg4EB///13mZ+LIWJRUVHUrFkz6h7Yn9ym7ij3nLKzs2nZsmVkZ2dHe/bs0dov52UKCgro0KFDNG3eUqrzzsfkPGmrRrlVHqpyS6/DySdPnmDAgAHY9PM62L//JWx6TYOpjSuuPEyvsG0TExO88847OHXqFL799lsMHDgQCxYs0OvQJyQkBBcvXsTy5cv15lNbrjxMR6FUDjkBRVK52t+BjY0NfvjhB4wbN84ghpV//fUXevbsiQULFqDfzG9B1W0hJ6DwNedUs2ZNLFq0CKdPn8bChQuxevVqPUdd+di0aRNcXFywfv163Pw3AbnJD5C0cwYSDixDSMQ/OvWtNxE7ceIE2rRpA0tHV9Qc+g0kdRoAAExNTdDezUZQXwMGDEBMTAxu3LiB9u3b4+7du4LaL4vnz59j4sSJ2Lhxo8GvAwOA9m42MDczgSkDJGaafQdDhgxBw4YNsXLlSh1GWD7379/HwIEDsWvXLgwfPhxW1cwhL54WlRNgVe31N3w8PDxw8eJFrF27Fvv27dNDxJWTb7/9Fj/88ANCQkJw9uxZfPPtt6j77nTUn7wNFvau2Dh7JNLTK95RUUlZ3TNVL227/Hv27CEHBwe6ePEirQu9Tw3mHiOXOcfIdc4xmnfkegU7maqRy+W0YcMGcnR0pJiYGJ35ISL67LPPaMSIETr1ITQVmbdISEggW1tbunXr1kvvQ0/DSZlMRp06daI1a9aUvLcu9D65zlHkVoM5x2hd6H21bEVHR5OdnR0lJydrFAOH6P79+2Rra0uPHz9+6f3SuTVlyhSaOHFihX2pyi2d98T27NmDzz//HH/88Qc6duz4Ug+gisQEgzyddeabMYZJkybhxx9/RI8ePXD9+nWd+Ll69Sr27t1rdMMSLxcrtHezwZWH6YiMz9DoWGdnZyxbtgzjxo2DTCbTUYSq2blzJwoKCjB16tSS99q72aCKRJFb5hL1e5dt27bFuHHjMHPmTF2FW2nZunUrRo0ahbp16770func6v3hJBw4cADPnz/XTRBlKZuql6ZXy1OnTlGdOnX+c7XWxZ2L8ti/fz85OztTXFyc4LZ9fX1p165dgtvVNdpO7iuRyWT01ltv0f/+97+S96CHt5N7KQAAIABJREFUnphcLqdWrVrRmTNnyjwnbXIrOzubLC0tKSkpSaPj3nTatGlDV65c+c/7r+aWp+/bdPr06Qr5UpVbOlvs+uDBA3z44YcIDg7+z2p1LxcreLlY6cp1mQwdOhQpKSno06cPIiIiUKVKFUHsRkVFITk5GcOGDRPEnj6IjM/AlYfpSMrM+8/kvibfi4mJCebMmYNVq1ZhxIgROoz4ZcLDw5Gbm4uAgID/fKZtbtWsWRNDhgzB9u3bMX/+fCHCfCPIzc1F7dq1S35WlVtFplWRm5urkxh0MpwkIkyYMAFz5swxqGoN06ZNQ4MGDfD1118LZnPLli0YN24cTE1NBbOpKyLjMzA/6AaGbbmC70Pu4mBEAsxMtZvcV9K7d2/ExcXh5s2bOoi4bEJDQ9G3b1/BnxHQr18/nDt3TlCblZ26devi3r17r80tM1OG7JRH/xlyCoVOemLbt29HVlYWpk+frgvzWsMYw/r16+Hh4YHBgwejRYsWFbL3/PlzHDhwADdu3BAoQt0RGZ+B4VuvoKBIDuW+BpmcMLRdPdS1tEB7NxutejBmZmYYO3YsNm/ejLVr1wobtAqio6Px7rvvCm7Xw8MDUVFRICIwxgS3XxkZMWIEvvtpA560l6jMrarp9/HtIROdlFwHdNATS0lJwdy5c/HLL7/AzMygtmYCUExIL1++XJAJ6QMHDqBjx446u8IIiXJdmDLJGBS9r0GezvikS6MKDe/Hjh2LPXv26Gy48CoPHjxAkyZNBLfr6OiIwsJCZGVlCW67sjJ06FDciLmGZzculJlbH3raYd2X8zFjxgydXRgEF7H58+dj7NixaNOmjdCmBWPChAkwMzPDrl27KmRn8+bNmDhxokBR6ZbSd4XNTRk+8K2PPePaCzI36eLiAl9fXxw6dEiASMunsLBQsDnNV6lSpQqKiop0YrsyUr16dfy84wCend2MzHPbwJ6n4QPf+tg5ygsPws7g7bffRseOHTF+/HjdBVHWbL+qV3l3kNLT06l27dqUlpZWobsQ+uD48ePk4+Oj9fHR0dFUr149kkqlAkalW3R5VzgoKIjeeustvdyd9PX1pT///FPYEyDF3VYLCwvKyckR3HZl59jf1yngvTFUq7Yl2draUrVq1ahTp0505MgRksvlgvhQlVuC9sR2796N3r17w9bWVkizOqFnz55ISUlBdHS0VscfPXoU77//vlFM6CvxcrGq8NBRFYGBgXqb3G/WrBlu3boluN3Y2FjY2NgYxY4LQyOwQyuc/fUXPE17glu3biElJQUXLlzAgAEDdD6/KJiIERE2b96MCRMmCGVSp5iammLcuHHYsmWLVsdHR0fDx8dH4KiMF4lEgrZt2+rFl4eHByIiIgS3GxkZCQ8PD8HtvklIJBLUqVNHrzXnBBOxy5cvo6ioCJ07dxbKpM4ZM2YM9u/fjxcvXmh8bHR0NE/4V9DX76N3794IDg5GQUGBoHb3799vVCXEOQoEE7GQkBAMGjTIqG5NOzs7o3nz5ggLC9PouPT0dGRmZsLNzU1HkRknnp6eevHTqFEjtG7dGocPHxbMZmJiIs6fP29Ui5Y5CgQTsaioKHh5eVXIRmR8Bn4+96/G+/gqgpeXF6KiojQ6Jjo6Gm3bthV8saWxo8+e6dy5czF//ny19+OVl1szZ87EhAkTeOltI0Swv8KKDq+UizG/D7mL4Vuv6E3IPDw8NJ7c50PJstFnKe5u3bqhc+fOmDt3brlty8utoKAgREREYPHixboKl6NDBBGxtLQ05OTkoEGDBlrbKKtInz56ZtqIWFRUFBexMtD34uY1a9bgjz/+KLd6yOty65fDJzFx4kTs3r0b1apVe+m4oqIi5ObmQi6X6/I0OBVEkKyLiYlB27ZtKzQfplyMWSSVQ2JmAqtq5hi+9QoKpXKYm5lgcZ8WyMgt1Hp7jCqaN2+OuLg45OXlwcLCQq1j7t27hxkzZggWA0c7rK2tERISgoCAACQkJGDlypWoWrXqf9qVlVsfbLmMZ9f+QMa5XzDrq3WIzLNFwh+XcevP07h69SoiIiKQnp4OiUSCoqIiNGzYED4+PnjrrbfwwQcfwMpKvwUMOKoRpCeWnp5e4ecRerlYYc+49pjRwx17xrVHRm5hydWzUCrH4qM3dTLUlEgkqF27NjIzM9U+Jjc3l68lMhDq16+PsLAwJCYmwsvLC0ePHv3PdrJXcyvmWhQeH1iGrKtBsB+6Av+7HI9ZY97DsAGBeJjyDOPHj0dUVBSkUiny8/NRUFCAI0eOoFu3bvj777/h5uaGSZMm4eHDhyKdNac0gvTE8vPzy7wCasqrZVSUV0/GGGRyAkG7kjHlUbVqVeTn56vdXqjz5QiDjY0NDhw4gKCgIKxcuRLTp09Hnz594OPjAzc3N5iYmODJkydICg/HlNBQxCU8RvWmPVCt+TQ8O78deQk3UdtvOGo19YNv75YY0KXRS/bNzMzQsmVLtGzZEqNHj0ZKSgo2bNgAX19frFixAhMmTDCqu/KVDUFETCaTCX6nTnn1vPIwHVbVzLHs2K2S4YDQNflNTU0hlUrVbp+Xl8dFzMBgjGHgwIEYOHAgIiMjcf78eZw4ceL/2rvzuCirNQ7gvwFZHWUZhkUQ0FiVZBvWG35QU+G6IIhLGSKmQlpKZuXVTMHUxGtYblcM0hTMSk0DE68KXTJRwOuGqCiLu4KhCMM+z/2DywTC6AzMwEye7+czfzTOe94zw9Mz73vmnOfg9u3baGpqAo/Hg7u7O5YvX47XX38dB3/NQ+RboRB4+eF24GY0qWlJHVumpqaIiYnBlClTEBYWhvT0dKSkpLCY6CFySWLa2tpyn3gItL0yszftg+yiR3IfEwOar6ykHQ9jlJ+7u/tzp/tcvXoV74UFY/3naxAeHi4u5CdrbA0aNAinTp3C9OnTMWHCBBw6dEgpdqMnIpSWliInJwc5OTkoLS1FXV0dNDU10a9fPwgEAnh4eMDW1vYvMU1ILu9AR0dHcfWz/0+R6/6qq6tlSmI6Ojoy3X4yyqOiogIBAQH47LPPEB4eDqBrsaWpqYndu3dDW1sbUVFR8u6uTGpqapCYmAgXFxf4+Phg9+7d0NPTQ3BwMMLDwxEaGgpzc3McOnQIo0ePhq2tLeLj42UaD1ZKHa0Kl/SQVGkgPz+fbGxs5LJSvbvdvXuXeDyeTCvtHRwc6PLlywrslepCD+0ALq1Zs2bRO++8I/d2q6qqaODAgZSamir3tl9EJBLRjh07iM/n09ixYyk9PZ2amppeeMypU6fojTfeIAMDA4qLi1P6iiySYksut5P29va4d+8eKisr0bdvX3k02W1aJq7KMjCrq6ur8CtPRv4yMjKQnp6ukGobvXv3RkJCAiIiIpCfn99tM//Ly8sRGRmJwsJCHD16VOpF+BwOB97e3vD29kZxcTFmzJiBn3/+GTt37uzSfM+eIJfbSXV1dTg5OeH8+fPyaK5bdWb2vZ2dHQoKChTUI0ZRVqxYgbi4OIV90Y4YMQI+Pj5ITExUSPvPKikpgZeXFwYMGICcnJxOVxEZMGAAMjIyMH78ePj4+Mi8DK+nyW1Uz83NTeXePNA8+17WhcsttdiZtmT5hbe7XblyBdeuXcPEiRMVep558+YhISEBzXc/inPv3j0MHz4c0dHR+Oc//9nlSrdqampYtGgRtm7disDAQIXUa1MUuSUxLy8vZGZmyqu5btHQ0ICTJ0/C09NTpuPc3Nw6XUzxr0yZr06TkpIwY8YMaGhoKPQ8fn5+EIlEyM7OVtg5GhsbMX78eLz99tttNg+Wh+DgYHzxxRcYM2aMyuw1ILckFhwcjMzMTNy/f19eTSpcWloabG1tZS6p4+rqinPnznX7mrqeqPIhC2VO7FlZWQgMDFT4eTgcDgICAvDbb78p7Bzr16+HgYGB3PbHfDaupk2bhlGjRkm1uF4ZyC2J9e3bF6GhodixY4e8mlS4zlai5fF40NPTQ3FxsQJ61bGeqvIhC2VNYk1NTbh48WK3Vp5V1HDD9evXsW7dOmzbtk0uqwQkxVVcXBwOHTqErKysLp9D0eQ6023OnDnYvn27Sqz6Ly0txZkzZxAaGtqp47t7XKx1JYa6BhH2n73dbeeWlrKOE964cQPGxsbd9su5s7OzwvYiXbVqFaKjo+X2C2LruKpvFGHDseaNcPX19REXF6cS5YnkmsQEAgH69u2L9PR0eTarEP/6178wbdq0Ts/Ud3NzQ05Ojpx7JZn3QB56qTf/uQjAD7m3lOpqrKGhAefOnevpbnSosrKyW6tOGBgYoLKyUu7tPn78GAcOHJDrPhYtFT7UOICIgN8Ky8VXZJMmTcLly5dx7do1uZ1PEeSaxDgcDlasWIEFCxYo9Yz2K1euYPv27V0qpxMUFITvvvuuyxvwSsvdygCh7hZouYFoEhGyix51y7mlkZaWBicnp57uhlLgcDgKuRtJTk7G6NGjYWxsLLc2W9Yo/83GCBygTZEFTU1NzJgxAwkJCXI7nyLIfeFUUFAQnJ2dERsbK++m5UIkEmH27NlYvnw5rKysOt2Oi4sLTE1Nu/Wqc6KbBbQ0mjfAVcRC+K5Q5p2uuFyuQq6MJHny5IlCJrump6d3evjjedytDBD9ul2HsTV58mTlv7PqaBq/pIe0S0Pu3btHfD6fzp07J991BxJ0tCmspI1it27dSt7e3nJZYrF9+3YKCgrqcjuyUOQGuJ1VUlJChoaGJBQKlXLZUX19Peno6FBVVZXMx8oSWy2Sk5MpNDS00/2VxNzcnG7cuCH3dlt09L5qa2tJR0eHhEKhws4rLUmxpZB6wqampvj8888xc+ZM/PbbbwqtEJFXWoE3tmeLy/Tsme0NAG2qwibP8oa7lQGKioqwbNkyZGZmymXT26lTp+LDDz/EnTt3YG5u3uX2pPFszTVlkJiY2KXxRUXT0NDAoEGDcP78efj6+kp9nCyx1Zoiypc/fPgQVVVVCl0S1FFsaWlpwc7ODpcuXVLafVYVVocjIiIC9vb2mDZtmkLHjfafvY36RhEIzb+u7D97u8Oa6mVlZQgICEBMTAwGDx7crh0iwsmTJ7FmzRqEhobC398fw4cPR0REBLZs2dJhFU8ul4spU6YgKSlJYe9P2TU2NiIxMRGzZ8/u6a48l6+vL44dOybTMdLGVmtEhOPHj8PHx0eOvQdu3ryJV155pUeKL9rY2KCkpKTbzysthSUxDoeDb775BpWVlZgzZ47Cpl08u7iD8OcvLi3394N46ggMDMTkyZMxd+7ctq8nwp49e+Di4oJZs2ahvLwcEydOxPLly7FkyRJ4eXkhNzcXXl5eGD9+fLu5UJGRkfj66687TNR37tzB4cOH8f333yM1NRU3b95U+HKU7nb48GFYWVnh1Vdf7emuPFdERAQSExNl+kKVJraeHZfMy8vDkydP5L6JdG1tbZeXFnWWrJWPu11H95iSHp0Zt3j69Cn5+/tTWFgYNTQ0dPGuuL3ckj/IdkkaWX+cSrZL0sT38y3398fOFpKLiwstXLiwXbmdsrIymjBhAg0ZMoTS09OfW46nurqatmzZQnw+n1auXNlmTM3Ly4t2795NRETl5eW0evVqsra2JiMjIxo5ciSFhobSqFGjyNjYmCwsLCgmJoYePHgg98+iuzU1NZGvry99++234ueghGNiLQQCAR06dEjq178otjoaE5s5cyatXr1abn1ukZ2dTZ6ennJvVxqTJ0+mPXv29Mi5W5MUW90SaNXV1RQYGEi+vr50/fr1Lr6V9iQF1bFjx8jS0pKWLVvWLkHdvXuX7O3taeHChVRbWyv1uW7dukV+fn701ltviRPZ6dOnycTEhBITE8nY2JhmzJhBeXl57c4pEonowoULNHv2bOLz+bRr1y6Z6pgpm61bt5KPj0+bhK7MSeynn34iW1tbqqmpkfoYWX5IaYmD8vJyIiKqqamhBw8eyGVQ/OrVq2Rtbd3ldjpj2LBhlJ6e3iPnbq1HkxhR87d2fHw8GRkZ0bZt2xT6P69QKKQFCxaQhYVFhx9+TU0Nubq6UkxMTKfb9/f3p0WLFhFRc3Ly9PSkPn360KlTp6RqIy8vjxwdHSk6OlolE9mtW7fIyMiILl261OZ5ZU5iREShoaH0j3/8Q+7t1tXVkZOTE+3evZv+/e9/07hx40hLS4t4PB5paWnRyJEj6eeff35hsUJJGhsbicvlUkVFhZx7/nwikYgMDAzo/v373XrejvR4EmuRn59Prq6uNGbMGLp48WKX22tNJBJReno6OTg40NSpU+nRo0cdvm7ZsmU0YcKELiWPR48ekZmZGZ08eZLWr19PgwYNov79+9ORI0ekbiPjfDH1t3OiyPcXd7ofPUEkEtG4ceNoxYoV7f5N2ZPYvXv3yMzMjA4cOCC3NkUiEc2cOZPGjRtH4eHhZG9vTwkJCVRdXU1EzdMUdu7cSc7OzhQcHCzTlWBrvr6+lJGR0e75/Px8ioqKIj6fTxoaGmRsbEyTwt6mT3Yc6fJUnOLiYjIzM+tSG/KiNEmMqPlba/Xq1WRmZkbDhw+ngwcPdmne1tOnT2nLli3k6OhIQ4YMoe+//17ia4uKisjIyIju3LnT6fO12Lt3L9nZ2RGPx6PS0lJKT08na2trevr06QuPzS35g+w/OUz9391F6r31aXdqRpf701327t1LgwYN6vA2XNmTGBFRTk4O8fl8uZSSbmpqonfffZc8PT0pKiqKRowYIXE+Wm1tLQUHB9O0adM69QX6wQcf0NKlS8X/LRKJaNWqVWRsbEzLly+nkpISqqmpoYNZ/yXea2+Qmq4eGQ2P6FIi+/rrrykkJKTTx8uTUiWxFnV1dZScnExeXl40YMAA+uijj+iHH36gkpKS5/6RGxsb6eLFi5SUlERz5swhQ0NDCgkJoczMzBcGR8Tc92n4pJlymSwqEolIV1eX3n//ffFz06dPpzdmRr1wHGXTiUIasDiVrD5OJf6Y92nAINcu96c7lJeXk6mpKf3+++8d/rsqJDGi5oFyExMT+uSTT6i+vr5Tbdy+fZtGjhxJfn5+tOvQCeLqG1LmhZLnHiMUCsna2prOnDkj8/kuXrxIfBMz+vLfBZRb8gdt3LiRHB0d6e7du21e1xJb5nN3kIahBU2J7tywCVHzj1Y9sW9AR5QyibV2+vRpio2NpXHjxpGpqSnx+XwaNWoUTZ48mcLDwyksLIxCQ0PptddeIy6XS7a2tvTmm29SfHw8FRcXS3WOnOJHpM41JPNZW8n+k8NSJbLnDexevXqVuFwuRUREiJ879t/rpM41JJPQ5c89R8uV2MDFqWS35GcyNbegs2fPSvU+ekpjYyMFBQXR/PnzJb5GVZIYUfOPO4GBgeTm5kZHjx6V+uqoqqqKNm/eTMbGxhQbG0vZ1x+SgfdE0vMKkSqu1qxZQxEREfTxynVkNsCOtHV0ycrKilasWCFxCISoOWa0LRzJOGQp2X60j/rq6VNhYWGHr2uJrQHvJJCevoH41lYW586dIwsLC6XZQERSbClkxn5neHp6iiusEhHu3LmD8+fPo7KyEjU1NVBTU4OOjg6MjIzg5ubWqaoEv2RfAkiEXrz+bSYqStpzsKXWUn2jCL3UOJgk6I8QNwvx61JSUhAUFNSmmsWVCoJxyFI8+DEW6jrLkF1kC3crg3Z7G7beHNh7IA/bqt7Ah59vxtq4dUo3Ix9o/pvMnTsX1dXVWLduXU93Ry7MzMyQlpaGXbt24YMPPkBjYyPefvtt+Pj4wNXVtc0KhLKyMuTk5OD48ePYuXMn/Pz8cPjwYbi7u2NzxnVUXTsNo/EftdmhXtJ+lqGhoXAcNAia5oPR12cGzPvZwr8fkPv7fvj6+iIrKwt8Pr9df7OLHqGvIAgVWckQVVfA2kkAGxubF8SWLxZc+QGRMRsRPXeO1LFFRPj4448xf/58uaxuUSSlSWKtcTgcWFhYwMLCQq7t9q66DW3TV8QTFQ10NZ+7hGT/2duoa/j/jO0mQsrpm9h39rb4dbm5uQgLC8O+fftQW1uL/Ac1uPO4Btz+jqAx7+PB/lUwCf8bUk5r4tODlyAiEp8H+DN5AsAvZX1Rnn0E077O7nApS0+LiYlBbm4uMjMzlWKDWHnhcDiYPn06wsLC8Ouvv+K7775DSkoKCgoKoK+vD01NTQiFQtTX10MgEMDX1xc5OTltlv94D+RBVFcFDa6+eAJs6y/AZ2MrPmEnGhsbYTE5BsRRRyOA4w8BLccweBnxsWDBAqSkpABo3kty69atSN67DxX1AMfqb1Dvw0NlXiqil32ClNM3XxhbN7iv4tyhNORoDpE6tpKTk3H//n1ER0fL+ROXP6VMYopi2Kse/s62CBxlD++BvA6XkLT8gfNKK/BD7q02s7Zblp+0vK6goACurq7gcrnIyr+J9366Ib5qm/lGCLR8TDB3Wgi0x38KNcP+wP+P33f2tnhJi2YvNUx0swDHsD/qy2+164cyWLlyJVJSUpCVldVtW5F1Nw6HA39/f/j7+wNoniFfXl6OhoYG6OjowNjYWOJu2e5WBjA11MN0gTEm+HvC3coAmzOudxhbuSV/ICEhARwNbRBHXVz+piW2nMaEY1XYcJSVlcHQ0BDjxo1DHanj4SuBqK+pxuPf92KwiwAXT5zHk8cV+PTgJTSKmqNUUmyRFhdN9UKpY6u0tBSLFi1CamqqwvckkIeXKolxOBwY9tbAvGE24uc0e6mJF/i2XkKSXfRIHBytiQgw0G2+EqmpqUHv3r1BRDh767E4aBtFhJt/CBEdMgU3yqqRvHEpjCfFQNN4ANQ4HHCANgFOALR0dEFN9UpVYoeIsGzZMhw4cAC//vorTExMerpL3UZbW1umO4HA0SOhe/cs3K1GA/hzedKzsZV5sRSNVRXQecUDagAG8nvjelk1gObYMjc2wpAhQ5Cfn4+nT5/iyZMnmL5mN+KPX4c6AdrWLshPmAMbRydsilsFwzfXoZdhc/EBSbGFmidQ1+ZKFVsPHz5EQEAAFi9eDIFAIOvH1iMUtnZSGZmbm7dZyNoydrBwlH27y+zWa+R6qXHExQjVAFQI6wE0b5h6PO8KnlYL0Y9v1GGFzNFBk2AyKhIPv1+GuqIcxAY5IcTNos36u4luFtgQ4oA+vXsrza1kdXU1oqKikJqaiszMTJiZmfV0l5RaS2n2hoYGAJJjy+sVI6CpEX2HjISmhhq8BvLaxZZQKISGhga+O/gLzF39weujI44tjhYXvfq/ivua5jA0McWD75ag4f519FLjSIwty4r/YvLE4BfGVnFxMfz9/TF58mSVuI1s8VIlsY52KXK3MsC8YTbt/ritgzA2yElcME5T489vM/MBdli4/htwDC2x6shVfDp2cLsKmRXCehz64kO8+9kmcH5Pwq9Jq2DP02gf4H/chKebs1IksNOnT8PV1RVCoRCZmZkdDjIzbQkEAjg5OeGdd94Rx9ezsUVEOJj0Jbh9+mDsIB6SZ3kjpFWhS00NNZirPcHNmzehxn8F6YWVyDhXiNjU/Dax1VRdAU0rFzx5VIa5776HyoOxCGjMxhSBebvkyXlUjJLCAmxdGiUxtogISUlJ8PT0RFRUFGJiYrrrY5OLlyqJ8Xg8WFlZSb2dVksQvullieRZ3pjqaYkQtz9vMfpaOuDp1WxoWQwWJ6yOKmS6WxlgQ/Q0XLl8CfX19XBxcUHdnYI2AZ6VldXj9ZoaGhrw6aefYvz48Vi9ejV27doFfX39Hu2TquBwOEhOTsa1a9cwduxYnDx5snkO0//l5eVh0qRJOHHiBLZs3oTM5C9hplEjTjpTPS0xbrAR1iz9AHPnzsXZO1XQcRiKqksZqH54SxxbDaV5aHx8H/r2nlgcuxY/JSch9af9uJz7G3x8fLBv3z44m/fBvGE24NY+xIQJE7Bx48YOK2AQEY4ePYrRo0fjq6++QkZGBubPn9+dH5t8dDTvQtKjJ+byyNuGDRvozTfflPm4lrk3AxaniucC7TuRQwCHzCM2tpkf9KJFw/v37ycTExNasmQJ1dXVUX19PZmamtLly5e79N66oqCggNzd3SkwMLDd5ElZQIXmiSmCUCikL7/8kmxsbMjOzo6GDh1Kjo6OZGlpSXFxcVRZWUlERJ9//jkZGRnR4sWLae3WHcQfMZN66ZlQX6dhdPpGmTje+H9fQGraXBo2eiyNGDGCeEbGtPCrveLYSkhIIAMDA4qMjKTVq1eTQCAgY2NjcnZ2Ji6XS+vWrWsz/+3hw4eUlpZGK1asIAcHBxoyZAht376d6urqeuTzkoWk2OIQtR+8lkQgEFBubq7CEmp3qKiogIODg3iOj7Q2Z1zH+qNXISJAnQMsHGWPJ6e+x/r4DfALDseypUtkuhV88OABZs+ejZKSEri7u6O0tBQnTpzozFvqklu3bmHz5s1ITEzEypUrERkZ2aXCexwOJ4+IZB4R/ivEVmsikQgXLlxARUUF9PT0MGTIEPTq1fZ3tGvXriEhIQHHTl9AUZU6eju9Dh1ze3ww2gHzhv05/8tBH3hQcAY6OjoICAhoV0H3/v372LZtG3788Uc8fvwYWlpa6N+/P9TV1XHhwgU8efIE2traqK+vh5aWFgQCATw8PBAQEIChQ4f2SKHFzpAYWx1lNkmPv8q35bfffkvOzs4ylUhpPQva/pPDtPfoSeLxeJSVlUU8Ho8KCgpk7odIJKJt27aRhoYG8fl8io2N7ZZqASKRiP7zn/9QaGgoGRoa0oIFC6Re9fAieMmvxDrj2diS5/4JIpGIhEIhVVRUUHV1tUpWTGkhKbZeykATiUQ0depUGj9+vEyX0S23iYeyzpGVlZW4GGBiYiINHDiQSktLZerH3bt3ycHBgb766itxnTF9fX2aPn065ebmytSWNGpqaigpKYlcXFzIzs6ONm3aJL69kReWxDpHGTeAUTYsiT2jrq6OJkyYQK+99hoVFRVJfdwvv/wpoZbvAAADX0lEQVRC/fr1o40bN7Z5Pj4+nvr160dpaWlStXPs2DHq378/ffbZZ22eLy8vp7Vr15KlpSVZW1vTpEmTKC4ujjIyMmRKOCKRiIqKimjv3r20aNEiGjp0KPXp04f+/ve/05EjRzpd1+pFWBJjFEVSbL10Y2KtNTU1IT4+HmvXrsXcuXMRFRUlcT5UXl4evvjiC2RlZeGbb77BiBEj2r0mMzMTM2bMgJOTE9577z34+/u3+VWovr4eWVlZ2LRpE3JycpCYmIjRo0d3eD6RSITCwkLk5OSIHxcuXIClpSU8PDxgYmICHR0daGtrQ0NDA7W1taipqYFQKMSVK1eQm5sLTU1NeHh4iB8CgQCGhoby+fAkYGNijKJIiq2XOom1uHbtGjZs2IA9e/bAzs4OHh4eMDMzQ1NTEwoLC3HmzBnU1tZi3rx5iIyMhJ6ensS2hEIhUlJSsG3bNuTn58Pe3h56enp4+vQprly5Ant7e8yePRthYWHgcrky9bOhoQH5+fnIy8tDeXm5OHG1LI3R1taGjo4ObGxs4OHhgX79+nX1o5EZS2KMorAkJoXq6mqcPXsWOTk5KC8vh5qaGqytreHp6YnBgwfLvJpfKBQiPz8fVVVV0NXVxeDBg2VOXKqGJTFGUSTF1ku1dvJFevfuDT8/P/j5+cmlPV1d3R6fwMowf3Uv1Yx9hmH+elgSYxhGpbEkxjCMSmNJjGEYlcaSGMMwKo0lMYZhVBpLYgzDqDSWxBiGUWksiTEMo9JYEmMYRqWxJMYwjEpjSYxhGJXGkhjDMCqNJTGGYVQaS2IMw6g0lsQYhlFpLIkxDKPSWBJjGEalsSTGMIxKY0mMYRiVxpIYwzAqjSUxhmFUGktiDMOoNJbEGIZRaSyJMQyj0lgSYxhGpbEkxjCMSmNJjGEYlcaSGMMwKo0lMYZhVBpLYgzDqDSWxBiGUWksiTEMo9JYEmMYRqWxJMYwjEpjSYxhGJXGkhjDMCqNJTGGYVQaS2IMw6g0lsQYhlFpLIkxDKPSWBJjGEalcYhI+hdzOGUAShXXHeYvwIqI+LIexGKLkUKHsSVTEmMYhlE27HaSYRiVxpIYwzAqjSUxhmFUGktiDMOoNJbEGIZRaSyJMQyj0lgSYxhGpbEkxjCMSmNJjGEYlfY/n8woB7D4sLEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.patches import Circle\n",
    "\n",
    "# We'll create a BallTree class which will recursively subdivide the\n",
    "# space into circular regions.  Note that this is just an example and\n",
    "# shouldn't be used for real computation; instead use the optimized\n",
    "# code in scipy.spatial.cKDTree or sklearn.neighbors.BallTree\n",
    "class BallTree:\n",
    "    \"\"\"Simple Ball tree class\"\"\"\n",
    "\n",
    "    # class initialization function\n",
    "    def __init__(self, data):\n",
    "        self.data = np.asarray(data)\n",
    "\n",
    "        # data should be two-dimensional\n",
    "        assert self.data.shape[1] == 2\n",
    "\n",
    "        self.loc = data.mean(0)\n",
    "        self.radius = np.sqrt(np.max(np.sum((self.data - self.loc) ** 2, 1)))\n",
    "\n",
    "        self.child1 = None\n",
    "        self.child2 = None\n",
    "\n",
    "        if len(self.data) > 1:\n",
    "            # sort on the dimension with the largest spread\n",
    "            largest_dim = np.argmax(self.data.max(0) - self.data.min(0))\n",
    "            i_sort = np.argsort(self.data[:, largest_dim])\n",
    "            self.data[:] = self.data[i_sort, :]\n",
    "\n",
    "            # find split point\n",
    "            N = self.data.shape[0]\n",
    "            N2 = int(N/2)\n",
    "            N21 = int(N/2 - 1)\n",
    "            split_point = 0.5 * (self.data[N2, largest_dim]\n",
    "                                 + self.data[N21, largest_dim])\n",
    "\n",
    "            # recursively create subnodes\n",
    "            self.child1 = BallTree(self.data[N2:])\n",
    "            self.child2 = BallTree(self.data[:N2])\n",
    "\n",
    "    def draw_circle(self, ax, depth=None):\n",
    "        \"\"\"Recursively plot a visualization of the Ball tree region\"\"\"\n",
    "        if depth is None or depth == 0:\n",
    "            circ = Circle(self.loc, self.radius, ec='k', fc='none')\n",
    "            ax.add_patch(circ)\n",
    "\n",
    "        if self.child1 is not None:\n",
    "            if depth is None:\n",
    "                self.child1.draw_circle(ax)\n",
    "                self.child2.draw_circle(ax)\n",
    "            elif depth > 0:\n",
    "                self.child1.draw_circle(ax, depth - 1)\n",
    "                self.child2.draw_circle(ax, depth - 1)\n",
    "\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Create a set of structured random points in two dimensions\n",
    "np.random.seed(0)\n",
    "X = np.random.random((30, 2)) * 2 - 1\n",
    "X[:, 1] *= 0.1\n",
    "X[:, 1] += X[:, 0] ** 2\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Use our Ball Tree class to recursively divide the space\n",
    "BT = BallTree(X)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot four different levels of the Ball tree\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "fig.subplots_adjust(wspace=0.1, hspace=0.15,\n",
    "                    left=0.1, right=0.9,\n",
    "                    bottom=0.05, top=0.9)\n",
    "\n",
    "for level in range(1, 5):\n",
    "    ax = fig.add_subplot(2, 2, level, xticks=[], yticks=[])\n",
    "    ax.scatter(X[:, 0], X[:, 1], s=9)\n",
    "    BT.draw_circle(ax, depth=level - 1)\n",
    "\n",
    "    ax.set_xlim(-1.35, 1.35)\n",
    "    ax.set_ylim(-1.0, 1.7)\n",
    "    ax.set_title('level %i' % level)\n",
    "\n",
    "# suptitle() adds a title to the entire figure\n",
    "fig.suptitle('Ball-tree Example')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computational efficiency of linked lists and trees\n",
    "\n",
    "In practical applications involving large number of points and dimensions, one has to be concerned with computational efficiency of the particular data structure implementation. Small changes in the implementation may result in significantly different performance. \n",
    "\n",
    "Check out exploration of efficiency of different list implementations for range queries (another fancy term for neighbor searches that one often encounters) by Phil Mansfield in his blog <a href=\"https://floatingpointastronomy.wordpress.com/2016/10/06/high-performance-lists-i-lists-and-you/\">here</a>, <a href=\"https://floatingpointastronomy.wordpress.com/2016/10/10/high-performance-lists-ii-naive-linked-lists/\">here</a>, and <a href=\"https://floatingpointastronomy.wordpress.com/2016/10/20/high-performance-lists-iii-array-lists/\">here</a>.\n",
    "\n",
    "In this <a href=\"https://jakevdp.github.io/blog/2013/04/29/benchmarking-nearest-neighbor-searches-in-python/\">blog post</a> Jake van der Plas compares performance of different implementations of $k$-d and ball trees in scipy and scikit-learn. The differences in some regimes are substantial and one needs to consult these benchmarks to decide which implementation is best for the job, if the job is very computationally intensive. \n",
    "\n",
    "Here we will concern ourselves with the primary application of tree structures for neighbor searches and density estimates. The first example is estimate of density of MCMC samples of the posterior distribution to estimate the Bayes factor (see notes). The second example concerns searches for pairs of halos in a halo catalog separated by a specified separation, $R$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example of how k-D tree can be used to estimate density\n",
    "\n",
    "In this particular example the density estimate in multi-d parameter space from the MCMC samples of the posterior distribution is used to estimate Bayes factors of a model for model comparison in a Bayesian analysis. How the Bayes factors are used in Bayesian model selection is well-illustrated and contrasted with the frequentist approach in this <a href=\"https://jakevdp.github.io/blog/2015/08/07/frequentism-and-bayesianism-5-model-selection/\">blog post</a> by Jake van der Plas (look towards the end). There he has simply computed evidence integrals needed for estimate of the Bayes factor via standard integration techniques because problems he considered were 2- and 3-d. However, in multiple dimensions such integration is not feasible, while Monte Carlo integration converges too slowly. \n",
    "\n",
    "Ivezic et al. in their book \"Statistics, Data Mining, and Machine Learning in Astronomy\" introduced a neat trick (S 5.8.4), which allows one to estimate the evidence integral by using a density estimate of MCMC samples near the peak of the posterior distribution. This is because computing the evidence for a model $M$ involves integrating over unnormalized posterior distribution $\\pi(\\vec{\\theta}\\vert \\vec{d}, I)$:\n",
    "$$ E(M) = \\int \\pi(\\vec{\\theta}\\vert \\vec{d}, I) d^k\\vec{\\theta}.$$\n",
    "By the nature of the MCMC algorithm, the density of MCMC samples in a converged chain at a given location should be proportional to $\\pi(\\vec{\\theta}\\vert \\vec{d}, I)$ at that location:\n",
    "$$\\rho(\\vec{\\theta})=C N\\pi(\\vec{\\theta}), $$\n",
    "where $C$ is an unknown constant of proportionality and $N$ is the total number of MCMC samples in the chain. If we integrate this equation over the entire parameter space we get: \n",
    "$$\\int \\rho(\\vec{\\theta})d^k\\vec{\\theta}=N= CN \\int \\pi(\\vec{\\theta})d^k\\vec{\\theta}= CN E(M).$$\n",
    "Thus, $E(M)=1/C$ and we can estimate at any point $\\vec{\\theta}$ in the parameter space the evidence as: \n",
    "$$E(M)\\approx N\\frac{\\pi(\\vec{\\theta})}{\\rho(\\vec{\\theta})}.$$\n",
    "\n",
    "Of course, it makes sense to do this at a point where density of the points is highest - i.e., near the peak of the posterior. \n",
    "This is the best we can do, because in that region the estimate of density will be the most accurate. Note, however, that measurement of density in the non-metric space rears its head again. If models with very different parameters and parameter units are compared, this estimate may not make much sense. For dimensionless parameters (such as say $\\Omega_{\\rm m}$ or any logged variable), however, this will work well. \n",
    "\n",
    "The routine from AstroML package below shows how scikit-learn's BallTree routine is used to estimate density of the MCMC samples near the peak of the posterior. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clever estimate of the Bayes factor from the AstroML library\n",
    "# http://www.astroml.org/book_figures/chapter5/fig_model_comparison_mcmc.html#book-fig-chapter5-fig-model-comparison-mcmc\n",
    "from sklearn.neighbors import BallTree\n",
    "from scipy.special import gamma\n",
    "\n",
    "def estimate_bayes_factor(traces, logp, r=0.05, return_list=False):\n",
    "    \"\"\"Estimate the bayes factor using the local density of points\"\"\"\n",
    "    N, D = traces.shape\n",
    "    \n",
    "    # compute volume of a D-dimensional sphere of radius r\n",
    "    Vr = np.pi**(0.5*D) / gamma(0.5*D + 1) * (r**D)\n",
    "\n",
    "    # use neighbor count within r as a density estimator\n",
    "    bt = BallTree(traces)\n",
    "    count = bt.query_radius(traces, r=r, count_only=True)\n",
    "    BF = logp + np.log(N) + np.log(Vr) - np.log(count)\n",
    "\n",
    "    if return_list:\n",
    "        return BF\n",
    "    else:\n",
    "        p25, p50, p75 = np.percentile(BF, [25, 50, 75])\n",
    "        return p50, 0.7413 * (p75 - p25)\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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}