SphericalParticleSelection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A typical use case is selection of solvent particles near a protein. Proteins can be treated as a group of atoms compressed in a form of sphere. Here, we consider a single protein at the centre of the simulation box  to compare PBC aware and non PBC data structures to evaluate the \"around\" selection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is to compare different shortlisted methods for this particular use case. The two types of data structures (hierarchical and grids) and naive brute frorce are used here. For linear grid data structures,  CellGrid ([here](https://github.com/MDAnalysis/cellgrid)) and Neighbour search implemented in FatsLiM([here](http://fatslim.github.io/)) are used. For hierarchical tree structures, Periodic KDtree implemented in MDAnalysis is used. Another tree structure, OCtree, which recursively divides a cuboid into its 8 children, is adopted from PointCloud library ([here](https://github.com/strawlab/python-pcl)). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A similar benchmarks for single query are done for the same datastructures, which can be accesed from [here](https://github.com/ayushsuhane/Benchmarks_Distance/blob/master/Notebooks/OctreeComparison.ipynb). A variant of similar benchmarks are done by Jonathan ([here](https://github.com/jbarnoud/Benchmarks_Distance)) and Sébastian' ([here](https://github.com/seb-buch/Benchmarks_Distance/blob/master/Notebooks/CythonNS.ipynb) as well. This is an extension to multiple particles to check the scaling behaviour of build and query time for multiple selections and with different cutoff radius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from cellgrid import capped_distance_array\n",
    "from scipy.special import gammainc\n",
    "from collections import defaultdict\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_protein(center,radius,nparticles,box=None):\n",
    "    \"\"\"\n",
    "    Generates set of particles distributed\n",
    "    inside a sphere\n",
    "    \"\"\"\n",
    "    r = radius\n",
    "    ndim = center.size\n",
    "    x = np.random.normal(size=(nparticles, ndim))\n",
    "    ssq = np.sum(x**2,axis=1)\n",
    "    fr = r*gammainc(ndim/2,ssq/2)**(1/ndim)/np.sqrt(ssq)\n",
    "    frtiled = np.tile(fr.reshape(nparticles,1),(1,ndim))\n",
    "    p = center + np.multiply(x,frtiled)\n",
    "    return p.astype('float32')\n",
    "\n",
    "def generate_solvent(center,radius,nparticles,box):\n",
    "    \"\"\"\n",
    "    Solvent particles are generated in the box uniformly\n",
    "    except inside the sphere\n",
    "    \"\"\"\n",
    "    points = np.random.uniform(low=0,high=1.0,size=(nparticles,3))*box[:3]\n",
    "    idx = np.where(np.sum((points-center)**2,axis=1)<radius**2)\n",
    "    points[idx] += np.sign(points[idx]-center)*radius\n",
    "    return(points[np.where(np.sum((points-center)**2,axis=1)>radius**2)].astype('float32'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CellGrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cg_select(protein,solvent,maxdist,box):\n",
    "    indx,dist = capped_distance_array(protein,solvent,maxdist,box[:3])\n",
    "    mask = np.where(dist<maxdist)\n",
    "    indices = (np.unique((indx[mask])[:,1]))\n",
    "    return solvent[indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Periodic KDtree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.lib.pkdtree import PeriodicKDTree\n",
    "def kdtree_select(protein,solvent,maxdist,box):\n",
    "    kdtree = PeriodicKDTree(box,bucket_size=10)\n",
    "    kdtree.set_coords(solvent)\n",
    "    kdtree.search(protein,maxdist)\n",
    "    indices = kdtree.get_indices()\n",
    "    return solvent[indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cython Neighbour Search - FATSLiM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from core_ns import FastNS\n",
    "def ns_select(protein, solvent, maxdist, box):\n",
    "    triclinic_box = np.array([[box[0], 0, 0],[0, box[1], 0],[0,0,box[2]]], dtype=np.float32)\n",
    "    searcher = FastNS(triclinic_box)\n",
    "    searcher.set_cutoff(maxdist)\n",
    "    searcher.set_coords(solvent)\n",
    "    searcher.prepare()\n",
    "    results = []\n",
    "    for coord in protein:\n",
    "        res = searcher.search(np.array([coord,]))[0]\n",
    "        if res.shape[0] != 0:\n",
    "            results.append(res[0])\n",
    "    if len(results) != 0:\n",
    "        return np.unique(results, axis = 0)\n",
    "    return\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Octree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pcl\n",
    "def oct_select(protein, solvent, maxdist, box):\n",
    "    cloud = pcl.PointCloud(solvent)\n",
    "    resolution = 0.2*box[0]\n",
    "    octree = cloud.make_octreeSearch(resolution)\n",
    "    octree.add_points_from_input_cloud()\n",
    "    result = []\n",
    "    for coord in protein: \n",
    "        [ind, sqdist] = octree.radius_search (coord, maxdist)\n",
    "        if len(ind) != 0 :\n",
    "            result.append(solvent[ind][0])\n",
    "    if len(result) != 0:\n",
    "        return np.unique(result, axis = 0) \n",
    "    \n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Brute Force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.lib.distances import distance_array\n",
    "def bf_select(protein,solvent,maxdist,box):\n",
    "    distance = distance_array(protein,solvent,box)\n",
    "    mask = np.where(distance<maxdist)[1]\n",
    "    return solvent[np.unique(mask)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarks\n",
    "\n",
    "Total time to select solvent particles around a spherical protein using different methods is shown below.\n",
    "\n",
    "For the first set, the scaling with number of solvent particles is shown. A radius of 6 units is chosen here for the purpose of demonstration. It is then followed by the variation of cutoff radius. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initialize\n",
    "box = np.array([100,100,100,90,90,90],dtype=np.float32)\n",
    "center = box[:3]/2\n",
    "radius = 10.0\n",
    "maxdist = 6.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm_notebook\n",
    "def timer(**params):\n",
    "    Nprotein = params.get('Nprotein', None)\n",
    "    Nsolvent = params.get('Nsolvent', None)\n",
    "    maxdist = params.get('maxdist', None)\n",
    "    centre = params.get('centre', None)\n",
    "    radius = params.get('radius', None)\n",
    "    protein = generate_protein(center,radius,Nprotein)\n",
    "    results = defaultdict(list)\n",
    "    \n",
    "    for Nsolvent in tqdm_notebook(np.unique(Nsolvent)):\n",
    "    \n",
    "        solvent = generate_solvent(center, radius, Nsolvent, box)\n",
    "        for cut in tqdm_notebook(np.unique(maxdist)):\n",
    "            \n",
    "            \n",
    "            res = %timeit -q -o kdtree_select(protein, solvent, cut, box)\n",
    "            results['kd'+'.'+str(cut)].append(res.average)\n",
    "            \n",
    "            res = %timeit -q -o ns_select(protein, solvent, cut, box)\n",
    "            results['ns'+'.'+str(cut)].append(res.average)\n",
    "            \n",
    "            res = %timeit -q -o cg_select(protein, solvent, cut, box)\n",
    "            results['cg'+'.'+str(cut)].append(res.average)\n",
    "    \n",
    "            res = %timeit -q -o bf_select(protein, solvent, cut, box)\n",
    "            results['bf'+'.'+str(cut)].append(res.average)\n",
    "\n",
    "            res = %timeit -q -o oct_select(protein, solvent, cut, box)\n",
    "            results['oct'+'.'+str(cut)].append(res.average)\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nprotein = 100\n",
    "Nsolvent = np.logspace(3 ,5 ,num = 6, dtype = np.int)\n",
    "center = box[:3]/2\n",
    "radius = 10.0\n",
    "maxdist = 6.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "    'Nprotein':Nprotein,\n",
    "    'Nsolvent':Nsolvent,\n",
    "    'maxdist':maxdist,\n",
    "    'centre':center,\n",
    "    'radius':radius\n",
    "         }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8f842e58285749c2a0639af98360cabb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=6), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "daffbfdb0f0345208686438c2cf99408",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ae0f96046d3e466e9b7fbfd4d6129cc5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1debd38d45f94f69854a7c5e9a9f3c21",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ba7f6afd76754234a2b89ec3f694f478",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b4b269e5c1d4a5f9cbc0e5f1bcb47c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "62373593575d41bf8984442139103a79",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "timing = timer(**params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff3f09f7978>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEaCAYAAAAG87ApAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4lFX68PHvyaTXSScFSCAQQgcDCBaqIAgIFqKrr/pbewdXXbEgukVdXVFX1GUtoOxKAipIUZAmUpQmvSVAAimkkZn0ZMp5/5hkSC+QSSbJ+XDlysxTzwyZ555TnvsIKSWKoiiK0lQObV0ARVEUpX1RgUNRFEVpFhU4FEVRlGZRgUNRFEVpFhU4FEVRlGZRgUNRFEVpFhU4FMVGhBCuQggphAi/jH3vF0KstkW5FOVKqcChdCpCiMIqP2YhREmV53c1su+NQoikFirH6SrnNQkhSqs8f0ZK+ZmUclpLnEtRWppjWxdAUVqTlNKz8rEQIhl4QEq5sQ3K0bNKOX4FPpRSLm3tcijK5VA1DkWpQgjhJoRYKITIEEKkCiHeFkI4CSH8ge+AHlVqBv5CiGuEEL8JIXRCiHQhxAIhxBV/IRNCPCKE2FjxuLLJ65GKmkq+EOJlIUS0EGK3EEIvhPhv1fMKIWYKIQ5VlOsXIUTfKy2TolRSgUNRqnsNGAgMAK4CxgDPSylzgZnAGSmlZ8VPLmAAngD8geuAacADNirbeGAQMBp4FfgAuB2IBIYDtwIIIa4GPgL+r6JcXwErWyKgKQqowKEoNd0FvCqlzJFSZgJ/Bf5ffRtLKXdLKfdIKU1SytPAp1gu7LbwppSyUEr5O3AKWCelTJFSXgQ2AEMqtnsYS9PXvopyLQJcsARCRbli6huIolQQQgigC5BSZXEKENbAPn2BfwJDATcsn6kdNipiZpXHJXU8r+y/6Q7MEkI8V2W9Mw28DkVpDlXjUJQK0pIq+gKWC2+lbkBa5SZ17PYfYD/QU0rpDbwOCFuWswnOA/OklNoqP+5Sym/buFxKB6ECh6JU9zXwakXHdxDwElA52ikTCBJCeFbZ3gvQSykLhRD9gAdbt7h1WgQ8KYSIFRaeQojpQgj3ti6Y0jGowKEo1c0DjgFHgQNYmp3+UbHuIPA9kFIxWskPmAM8IIQoBBYC8a1f5OqklDuAp4B/Azos/SF/oO4ak6I0m1ATOSmKoijNoWociqIoSrOowKEoiqI0iwociqIoSrOowKEoiqI0iwociqIoSrN0yDvHAwICZERERFsXQ1EUpV3Zt29fjpQysLHtOmTgiIiIYO/evW1dDEVRlHZFCJHS+FaqqUpRFEVpJhU4FEVRlGbpUIFDCDFNCLFIr9e3dVEURVE6rA4VOKSUq6WUD/n4+LR1URRFUTqsDhU4FEVROqVDCbCgP8zXWn4fSrDp6TrkqCpFUZRO41ACrH4KDCWW5/rzlucAA2fZ5JSqxqEoitIemU1Qkgc/zbsUNCoZSmDT6zY7tapxKIqitAWzCcryoVR/eT9l+Q0fX59qs6KrwFHpUIIlQutTwSccxs+zWTVPUVpVZ/zbbo3XbDZd2UW/sQs/Aly9wdWn4kcLvhFVnvuAizdsextKLtbe3Se8ZV9vFSpwQJu0ESpKq6jvb9tssvxtCwcQbT1Fegtr6ufZZLyyb/zlBY0UpI4Lv19k9Qt/Qz/OXuDQhN4Ej4DqrxfAyc0SLG2kQ84AGBsbK5uVcmRBf8sfV03OHjDkHnByBceKHyc3cHQBx4rfTm5V1rnWvVzjZJ8fTvVNtO1fs8kAhmLLh95QDOVVHtf7u+aykvqPUaprQiGEJYA4aCoCSeWPxvJ3azfrKoJcQ+uEA+z7AsrquKhrXMA/6sou/E296Dfnwt8SWujvWgixT0oZ29h2qsYB9bcFlhfBgf+BsQRM5Zd/fOFQd+CpLyBVXW4NSFW3qy9Q1djfQVN/mTpjLas5r1lKMJZdxgW7nnW1jlHx3Gxs/uvQVPx/O3tYfju5gZO75be7f/VluxfVf5wxL4I0V/yYqjw2g9l8metMlveurvWV68wmkIa61zW0X1PXGUvrfr2msopv/E0MAM6erXfhv1IDZ7Xq51bVOKD+GodPV5hzxPLYbLb8QVb+GEqqPC61BBdjWTOW1/W4YltjyaV9pfny3wgHpzoCSsXjjIOWD1JNTm7Qe/Lln7NSi9ewWuB4J9dZLtY1aZzAt0ftizyX8dmovIA71XFRr/bYHZzday+r9bvmMreGvxDU1JS/7Y6mM77mFtIpaxxCiGnAtKioqObtOH5e422EDg6WD7qze4uUtclMhjoCSn2BpmZAaiCY1RU0wLLthcNXWOgW/jLSUl9u6goaYHmPg/o0csGucfF3rmMbR1f7a5Jsyt92R9MZX3MrUzWOSvbW9m1rnfFbWWd8zdD5/rahc77mFtDUGocKHJ1VzfZ+sHwrm/ZBx/2AdcbXrCjN0NTA0U56fpQWN3CW5YLp0xUQlt8d/QLaGV+zotiAqnEoiqIogKpxKIqiKDaiAoeiKIrSLCpwKIqiKM2iAoeiKIrSLCpwKIqiKM2iAoeiKIrSLCpwKIqiKM2iAoeiKIrSLCpwKIqiKM2iAoeiKIrSLHafVl0I4QF8BJQDW6WU/23jIimKonRqbVLjEEJ8LoTIEkIcqbH8RiHESSFEkhDihYrFtwArpJQPAtNbvbCKoihKNW3VVLUYuLHqAiGEBlgITAb6AncKIfoC4UDlJAqmViyjoiiKUoc2CRxSym3AxRqLhwNJUsozUspyYBlwM5CKJXhAA+UVQjwkhNgrhNibnZ1ti2IriqIo2FfneBiXahZgCRhhwLfArUKIj4HV9e0spVwkpYyVUsYGBgbatqSKoiidmN13jkspi4D/a+tyKIqiKBb2VONIA7pWeR5esazJhBDThBCL9Hp9ixZMURRFucSeAsceoJcQIlII4QzcAXzfnANIKVdLKR/y8fGxSQEVRVGUthuO+zWwC4gWQqQKIe6XUhqBJ4D1wHEgQUp5tC3KpyiKotSvTfo4pJR31rN8HbDuco8rhJgGTIuKirrcQyiKoiiNsKemqiummqoURVFsr0MFDkVRFMX2VOBQFEVRmqVDBQ41HFdRFMX2OlTgUH0ciqIottehAoeiKIpieypwKIqiKM3SoQKH6uNQFEWxvQ4VOFQfh6Ioiu11qMChKIqi2J4KHIqiKO2cfvVqEseN53hMXxLHjUe/ut6pi1qE3c/HoSiKotRPv3o1Ga/MQ5aWAmBMTyfjlXkA+EybZpNzdqgah+ocVxSls8l6d4E1aFSSpaVkLXjPZufsUIFDdY4ritJZSJMJ/Zq1GDMy6lxf3/KWoJqqFEVR2hFpMpG/7gdyPv6Y8jNnwNERjMZa2zmGhNisDCpwKIqitAPWgPHRR5SfPYtL796Evf8+5rJSLsx7tVpzlXB1JWjObJuVRQUORVEUO1ZfwPC6YQLCwdLbIIQga8F7GDMycAwJIWjObJt1jAMIKaXNDt7aqswA+GBiYmJbF0dRFOWyWQLGOnI++tgaMAIef7xawGhpQoh9UsrYRrfrSIGjUmxsrNy7d29bF0NRFKXZagWM6GgCHn8Mrwm2CxiVmho4VFOVoiiKHZBG46WAkZyMS3Q0YR+83yoBo7lU4FAURWlD7SlgVFKBQ1EUpQ3UGTD+9QFe48fbbcCopAKHoihKK5JGI/lr11oCRkoKLn36tJuAUUkFDkVRlFbQEQJGJRU4FEVRbKiugBH+4b/wHDeu3QWMSh0qcFS5j6Oti6IoSicnjUb0a9aQ8/HHGFLO4RIT0+4DRqUOFTiklKuB1bGxsQ+2dVkURemc6gwYCz+0BAwhbHLOlb+n8fb6k6TrSgjVuvHcpGhmDAmzybmggwUORVGUttIWAQMsQWPut4cpMZgASNOVMPfbwwA2Cx4qcCiKolwBaTSiX72GnE9aN2BUenv9SWvQqFRiMPH2+pMqcFwpg8FAamoqpTUmPFFsx9XVlfDwcJycnNq6KIrS4qwB4+OPMZw7h0vfGMI/Wojn2LGtEjAADCYzabqSOtel17O8JXSawJGamoqXlxcRERGt9p/amUkpyc3NJTU1lcjIyLYujqK0GHsIGKUGE8v3pfLJ1tP1bhOqdbPZ+TtN4CgtLVVBoxUJIfD39yc7O7uti6IoLUIajei/X03OJ5+0WcAoLjfyv9/OsWjbGbIKyhjSTcuUAV1Y+msKJQazdTs3Jw3PTYq2WTk6TeAAVNBoZer9VjqCmgHDtW9fgj/6CM+xY1rtbzy/1MCXO5P5fEcyF4vKGdnDn/fiBjOypz9CCPqF+qhRVYqiKG1NGgyXAsb5820SMC4WlfP59rMs2ZVMQamRsdGBPDEuiqu6+1XbbsaQMJsGippU4GhFnp6eFBYWArBu3Tpmz57NTz/9xBdffMF//vMfAgMDKSoqYsCAAfz1r3+lb9++zJw5k7Nnz1JYWEh2dra1v+Cjjz5i1KhRbflyFKVDsoeAkZVfyqJtZ/jvb+coNZqY3L8Lj42Jon+YT6ucvzEqcNTDljfUbNq0iaeeeor169fTvXt3AObMmcOzzz4LQHx8POPGjePw4cN89913AGzdupV33nmHNWvW1HlMo9GIo6P671SUy1UrYPTrR/DHH+E5pvUCxvmLxfx722kS9qZiMktuHhTKY2N7EhXk1Srnb6oOdaVpqZQjtryhZtu2bTz44IOsW7eOnj171rlNXFwca9eu5X//+x9PP/10vccKDw/n7rvvZv369bz44osMHjyYJ554gpycHDw8PPj000/p3bs3mZmZPProo5w7dw4HBwc++OADrr766it6HYrSUVgCxvfkfPwJhtTUNgkYZ7IL+WjraVb+noYQcNtV4Twyuifd/T1a5fzN1aECR1NTjry2+ijH0vPrXf/7OR3lJnO1ZSUGE8+vOMTXu8/VuU/fUG9endavwfKVlZUxY8YMtm7dSp8+fRrcdujQoZw4caLBbQCCgoL4/fffARg7diyffvopPXv2ZMeOHTzxxBNs2LCBp556iueff56rr76a5ORkpk6dypEjRxo9tqJ0ZHUGjJdebNWAcTwjn4Vbklh7OAMXRwf+38juPHR9D0J8bDeUtiV0qMDRUmoGjcaWN5WTkxOjRo3is88+4/33329w26bOBR8XFweATqfj119/5dZbb7WuMxqNAGzcuJGTJ09al+fl5VFSUoKbm33/cSqKLdQKGP37E/zyS3iOHt1qAePAeR0fbk5i4/FMPJw1PDK6J/dfG0mAp0urnP9KdcrA0VjN4Jo3N9d5N2aY1o34h0de9nkdHBxISEhg/Pjx/P3vf+fFF1+sd9vff/+d2NhG54zHw8NSlZVSEhAQwIEDB2ptI6Vk9+7dODs7X3bZFaW9kwYD+lWryPnk320WMH47k8uHW5L4JTEHHzcnZk/oxX2jItC6t6/PZvvO7Wsjz02Kxs1JU21ZS91Q4+7uztq1a/nvf//LZ599Vuc233zzDRs2bODOO+9s8nF9fX0JCQmxdqabzWYOHjwIwIQJE1i4cKF127qCi6J0JPrVq0kcN57jMX1JHDuO9Jdf4fSNk8l4+RU0vr50/fcnRCxPwKsVmqWklPx8KpvbP9lJ3KJfOZ6Rz9zJfdjxwjhmT+jd7oIGdNIaR2MqO8BtNarKz8+PH3/8keuvv57AwEAAFixYwNKlSykqKqJ///5s3rzZuq6pli1bxqOPPsr8+fMpLy/n7rvvZtCgQSxcuJBHH32UL774AqPRyNixY6sFEkXpSPSrV5PxyjxkRV46Y0YG+hUrcAwPp+u/P8Hj+utbpYZhNkt+Op7Jh5uTOJymJ8THldem9yNuWFdca3wxbW9EU9vS25PY2Fi5d+/easuOHz9OTExMG5Wo81Lvu9LaEkePwZiZWWu5Y2gIvTZvtvn5TWbJmkPpLNySxKnMQrr7u/Po6J7cMjQcZ0f7buQRQuyTUjbaRq5qHIqitHvSZKJo+3byli+vM2gAGDMu2LQM5UYzK39P46OtSSTnFtMryJP37xjMTQNCcNTYd8BoLhU4FEVptwxpaei++Rbdt99ivHABjZ8fDh4emIuKam3rGBJikzKUGkwk7D3PJ1tPk64vpX+YN5/cfRUT+wbj4NAx87WpwKEoSrsiy8sp2LwF3YoVFO3YAYDHtdcSPHcuXmPHkL9+fbU+DgDh6krQnNktWo6iMiP//S2FRdvOklNYRmx3X/52ywDG9A7s8Ak+VeBQFKVdKDtzFt2KFehXrsR08SKOISEEPPYY2ltm4hR2aeCKz7RpAGQteA9jRgaOISEEzZltXX6l9MUGFu9M5oudZ9EVG7g2KoAnxg1hRKRfhw8YlVTgUBTFbplLSylYv5685csp2bsPHB3xGjsG7e2343HNNQhN3aOTfKZNa7FAUSm3sIzPtp/ly10pFJYZmRATxONjoxjSzbdFz9MeqMChKIrdKT1xAl3CcvSrV2MuKMCpezcC//QM2hkzcGzmMPUrdUFvyVT7v90plBnNTBkQwuNjougb6t2q5bAnKnC0ogsXLjB79mz27NmDVqslODiY9957j969e9e5fWUa9svNL/XAAw/wzDPP0Ldv32rLFy9ezN69e/nwww8v+7UoSkszFRaSv3YduuXLKT1yBOHsjNekSWhvuw334cNavRnoXG4xH/98mm/2pWKSkhmDw3h0TE+igjxbtRz2yO4DhxCiB/AS4COlvK3VTnwoATa9DvpU8AmH8fNg4KzLPpyUkpkzZ3LvvfeybNkyAA4ePEhmZma9geNKmEwmPv300xY/rqK0JCklJQcOoFuxgvwffkQWF+PSuzfBL72Ez7SpaLTaVi9TUlYBH205zaqD6WiE4PZYS6barn7urV4We9Vg4BBC1J9CtmITIENKWeeVTwjxOTAVyJJS9q+y/EbgfUADfCqlfLO+E0gpzwD3CyFWNFKWlnMoAVY/BYaKfFX685bncNnBY8uWLTg5OfHII49Ylw0aNAiAt99+m4SEBMrKypg5cyavvfZavccpLi7mvvvu48iRI0RHR5Oens7ChQuJjY3F09OThx9+mI0bN7Jw4UJefvll3nnnHWJjY/niiy9444030Gq1DBo0CBeX9pFMTemYjHl55K9ejW75csoSkxDu7vjcNAXt7bfjOmBAm3QyH03X89GW06w7koGro4b7RkXw4HU96OLj2uplsXeN1ThOSymHNLSBEOL3BlYvBj4EvqyyvQZYCNwApAJ7hBDfYwkib9TY/49SyqxGyth8P7wAFw7Xvz51D5jKqi8zlMCqJ2Dfkrr36TIAJtcb/zhy5AhXXXVVreUbNmwgMTGR3bt3I6Vk+vTpbNu2jeuvv77O43z00Uf4+vpy7Ngxjhw5wuDBg63rioqKGDFiBP/85z+r7ZORkcGrr77Kvn378PHxYezYsQwZ0uB/q6K0OGk2U7x7N7rlKyjYsAFpMOA6aCBd/vI63pOnoPFsm7kn9qXksXBLEptPZOHl4shjY3ryx2si8W8nmWrbQmOB49ZG1je4jZRymxAiosbi4UBSRU0CIcQy4GYp5RtYaidtr2bQaGz5FdiwYQMbNmywXsgLCwtJTEysN3Bs377dOrlT//79GThwoHWdRqOplla90m+//caYMWOsua/i4uI4depUS78URamTISsL/Xcr0X3zDYZz53Dw8UF7xx1ob7sV1+grTxx6OaSU7DqTy4ebk9h5OhetuxN/uqE394yKwMfNqU3K1J40GDiqXNw9gBIppVkI0RvoA/wgpTRUbtMMYcD5Ks9TgRH1bSyE8Af+BgwRQsytCDB1bfcQ8BBAt27dGi5BAzUDABb0tzRP1eTTFf5vbcP71qNfv36sWFG7tU1Kydy5c3n44Ycv67hVubq6oqlneKKitCZpNFK4fTu65Sso3LoVTCbchw8n8Mkn8bphAg6ubdP8I6Vk68lsPtySxL6UPAK9XHhpSgx/GNENDxe77/K1G01NoLINcBVChAEbgP+HpRnK5qSUuVLKR6SUPesLGhXbLZJSxkopY5ubVbaW8fPAqcYkR05uluWXady4cZSVlbFo0SLrskOHDuHt7c3nn39OYWEhAGlpaWRl1d86d80115CQkADAsWPHOHy4gSa3CiNGjODnn38mNzcXg8HA8uXLL/t1KEpDylPTyP7gA5LGTyD1kUcpOXgQ/z/+Hz1//IHuXy7BZ9rUVgkaK39P45o3NxP5wlqueXMz3+1L5YfDGUz913b+b/EeLuhL+cvN/fjl+bE8eH0PFTSaqanvlpBSFgsh7gc+klL+QwhxuZM6pAFdqzwPr1hmPyo7wFtwVJUQgu+++47Zs2fz1ltv4erqSkREBO+99x5arZaRIy0TRHl6erJ06VKCgoLqPM5jjz3GvffeS9++fenTpw/9+vXDx8enwXOHhIQwf/58Ro4ciVarrdYvoihXypICZDO65Sso2rkTAI/rriX45Zcs8104tW7Tz8rf05j77WFKDCYA0nQlPLP8IBKIDPDgH7cNZMbgMLvPVGvPmpRWvaID/DFgAXC/lPKoEOKwlHJAE/aNANZUjqoSQjgCp4DxWALGHuAPUsqjl/siqpxrGjAtKirqwcTExGrrOkp6b5PJhMFgwNXVldOnTzNhwgROnjxpt7P7dZT3Xamt7MwZdMsrUoDk5eEYGoL2llvR3noLTjZKKNgU9c3g6evuxN6Xb0DTQRMPtoSWTqv+NDAX+K4iaPQAtjShEF8DY4AAIUQq8KqU8jMhxBPAeiwjqT5viaABIKVcDayOjY19sCWOZ4+Ki4sZO3YsBoMBKSUfffSR3QYNpeMxl5SQv349uhUrLqUAGTcO7e234TFqVL0pQFpLuq6kzqABoCs2qKDRQpoUOKSU27D0c1Q+PwM81YT96pz7VEq5DljXxDIqVXh5eVFzkipFsbXS48fRLV+OfvUazAUFOHfvTtBzz+Jz8804BgS0admklOw+e5HFO5PZcKzuuTgAQrVu9a5TmqexGwDnSynnX+k2raVKU1VbF0VR2j1TYSH5a9ZaUoAcPYpwccFr0kRLCpBhrZ8CpKaSchOrDqSxeGcyJy4U4OPmxAPXRRLk5cI7609Z+zgA3Jw0PDepbYb+dkSN1TgeaOTucQHcAcxvsRJdgc7QVKUotiSlpOT3yhQgPyBLSnCJjib45ZctKUAaGYjRGlLzivnq1xTi95xHV2ygTxcv3rxlADcPDsPN2dJU5u/hwtvrT5KuKyFU68Zzk6KZMSSskSMrTdVY4PgP4NWEbRRFaceMeXnoV61Ct2IF5UmncXB3x2fqVLSzbse1f/82r11IKfn1zEUW7zzLTxXNUZP6deHeURF1zoMxY0iYChQ21NgNgPUnTbJDqqlKUZpOms0U//YbuuXLKfhpI9JgwG3QIEL++he8J0/GwaNtUoBUVVJu4rvf01iyM5mTmQVo3Z14eHRP7r66O2Gqz6LNdKi7Xuy9qaq5adUPHDhAeno6U6ZMAWD+/Pl4enry7LPPtmi55s+fzz/+8Q+Sk5Ot949UpnQH+Nvf/sb//vc/NBoNDg4O/Pvf/2bEiHpv9lfsjH716mqz4fn98Y/IwkJLCpDz5y0pQO68A+2tt+Ea3fKZmi/H+YuXmqP0JQZiQrz5x60DmT44FFcnlR2hrXWowNGS1p5Zy/v73+dC0QW6eHTh6aFPc1OPmy77eJeTVv3AgQPs3bvXGjhsKSAggH/+85+89dZb1Zbv2rWLNWvWsH//flxcXMjJyaG8vNzm5VFahn716mrzbxvT08n6618BcB8xgsCnn7akALGDbMlSSnaezmXxzmQ2Hs/EQQhurGiOGhbh2+bNZcol6tbJOqw9s5b5O+eTUZSBRJJRlMH8nfNZe+by8lRB/WnVr7vuOu655x5WrlxpXX7XXXexatUq5s2bR3x8PIMHDyY+Ph6wpBkZM2YMPXr04IMPPrDu8+6779K/f3/69+/Pe++9B0BycjIxMTE8+OCD9OvXj4kTJ1JSUvcY9z/+8Y/Ex8dz8eLFasszMjIICAiwpmEPCAggNDT0st8HpXVlvfOONWhU5RgURPcli/GZelObB43iciNLf01h4oJt3PXpb+xLyeOxMT3Z/uexLLxrKMM70Vze7UVT7xzvDXwMBEsp+wshBgLTpZR/tXUBm6Opd46/tfstTlw8Ue9xDmUfotxc+1u1s4MzAwMH1rEH9PHrw5+H/7neY37wwQecPXuWBQsW1Fr3888/s2DBAlauXIler2fw4MEkJiaydOnSajP1zZ8/nw0bNrBlyxYKCgqIjo7mwoULHDp0iPvuu49ff/0VKSUjRoxg6dKl+Pr6EhUVxd69exk8eDCzZs1i+vTp3H333dXOX9kEVlxcjMlk4rXXXrM2VRUWFnLttddSXFzMhAkTiIuLY/To0fW+zprUneOtT5pMFG7bhi4+wZJgsC5CEHP8WKuWq6ZzucV8uSuZ+L3nKSg10j/Mm3tHRjBtkGqOaitNvXO8qTWO/2C5c9wAIKU8hGUYrl2RUq6WUj7UWO6mxtQVNBpafqVGjx5NYmIi2dnZfP3119x66604OtbdinjTTTfh4uJCQEAAQUFBZGZmsn37dmbOnImHhweenp7ccsst/PLLLwBERkZac1NdddVVJCcn11uOp556iiVLllBQUGBd5unpyb59+1i0aBGBgYHExcWxePHiFnvtSssxZGaSvXAhSRNuIPXRxyg5egQHz7qnOXVso5QgUkp+Sczm/sV7GP3OFhbvTGZ070C+eXQkq5+4lttju6qg0Q40tY/DXUq5u0Z10WiD8rSKhmoGABNXTCSjKKPW8hCPEL648YvLOmd9adUr3XPPPSxdupRly5bxxRf1n6PqzH0ajQajseH/hprb19dUBaDVavnDH/7AwoULqy3XaDSMGTOGMWPGMGDAAJYsWcJ9993X4HmV1iFNJop27CCvsnZhMuFxzTUEz30Br7Fjyf/xx2p9HADC1ZWgObNbtZxFZUa+3Z/K4p3JnM4uwt/DmSfGRnHXiO5qhr12qKmBI0cI0ROQAEKI24DaV9YO4umhTzN/53xKTZc+bK4aV54e+vRlH3PcuHG8+OKLLFq0iIceegiwpFXX6/Vcd9113HfffQwfPpwuXbpvP52BAAAgAElEQVTQt29fwJJepOq3//pU7v/CCy8gpeS7777jq6++uqxyPvPMMwwbNswakE6ePImDgwO9evUCLB323bt3v6xjKy3HkJWF/ttv0SUsx5CejsbfH/8//hHtrNtx7nop+bTPtGkA1UZVBc2ZbV1ua8k5RSzZlcyKvakUlBkZGO7Du7MGcdPAEFwcVc2ivWpq4HgcWAT0EUKkAWeBuxvepf2qHD3VkqOqGkqrDhAcHExMTAwzZsyw7jN27FjefPNNBg8ezNy5c+s99tChQ62BB+CBBx5gyJAhDTZL1ScgIICZM2da+2IKCwt58skn0el0ODo6EhUVVW1OEaX1SLOZol270MUnULB5MxiNuF99NUHPP4fXuHGIepJd+kyb1mqBAsBslvySlMPiHWfZeiobjRBMGRDCfddEMKSrVnV0dwBN6hy3bmyZCdBBStn41+A20J7TqhcXFzNgwAD279/f6Pwa7Ym9v+/tgTE3F11l7eL8eTS+vvjMnInvrNtxjoho6+JZFZQa+GZfKl/uSuFMThEBni7cNaIbd43oRpC3ao5qD1o0rboQQgvcA0QAjpXfGKSUjWbIbU32fgNgfTZu3Mj999/PnDlzOlTQUC6flJLi334jLz6ego2bwGDAfdgwy30XE2/AwY5S6Z/JLuTLXSms2JdKYZmRwV21vBc3mCkDQtRkSR1UU5uq1gG/AocBs+2K0zlNmDCBlJSUti6GYgeMeXnov/0OXUIC5SkpaHx88PvDH9DGzcKlR4+2Lp6V2Sz5+VQ2i3cm8/OpbJw0gqkDQ7l3VASDu2rbuniKjTU1cLhKKZ+xaUkUpZOSUlK8Z4+l72LDBkvOqKuuIvTxx/CaNKnNb9CrKr/UwIq9qXy5K5nk3GKCvFyYM6E3d47oSpCXao7qLJoaOL4SQjwIrAHKKhdKKS/Wv4uiKA0x6XToV60iLz6B8jNncPDyQhsXh2/cLFwqRrHZi6SsQr7clcw3+1IpKjcxtJuWOTf0ZnJ/1RzVGTU1cJQDbwMvUTEkt+K3/dSdFaUdsMx38Tu6+Hjyf/gRWV5uyUj797/jPflGHNzsJ+Or2SzZcjKLxTuT+SUxB2eNA1MHhXDfqAgGhqvmqM6sqYHjT0CUlDLHloW5UiqtumKvTPn56Fd9jy4hnrLEJBw8PdHedivauDhco+1rZjp9iYHle8/z5a4Uzl0sJtjbhT/d0Js7R3QjwNN+ms2UttPUwJEEFNuyIC3B3kdVaTQaBgwYgJQSjUbDhx9+yKhRo5p1jMWLFzNx4sRmJRqcP38+//nPfwgMDATgxhtv5M0332zWeZXmk1JSevAgefEJltn0SktxHTDAMt/FlCk4uLu3dRGrScwsYMmuZL7dn0ZxuYnY7r48NymaG/t3wUmjmqOUS5oaOIqAA0KILVTv47Cr4bgtqeYcBi1xt62bmxsHDhwAYP369cydO5eff/652jZGo7HePFVgCRz9+/dvdobaOXPmXNY8HiaTCY1G3eHbHKbCQvTff48uPoGykycts+lNn442bhZu/fq1dfGqMZklm45nsmRXMjuScnF2dODmQZbRUf3D1NBwpW5NDRwrK346hbrmMMh4ZR5Ai92Bm5+fj6+vLwBbt27llVdewdfXlxMnTrBhwwamTp3KkSNHAHjnnXcoLCykf//+7N27l7vuugs3Nzd27drFsWPHeOaZZygsLCQgIIDFixcT0sQEdps2beLZZ5/FaDQybNgwPv74Y1xcXIiIiCAuLo6ffvqJ559/ntjYWB555BGys7PRaDQsX76cnj178vbbb5OQkEBZWRkzZ87ktdfa1YSRLa7k8BHy4peRv3adZa7uvjF0mT8f76lT0Xi23Wx6K39PqzX/9tjoIOL3nuPLXSmk5pUQ4uPKc5OiuWNYV/xVc5TSiCYFDinlElsXpDVd+PvfKTtef1r1koMHkTUmK5KlpWS89DK6hOV17uMS04cuL77Y4HlLSkoYPHgwpaWlZGRksHnzZuu6/fv3c+TIESIjI+tNFXLbbbfx4Ycf8s477xAbG4vBYODJJ59k1apVBAYGEh8fz0svvcTnn39ea98FCxawdOlSAN566y1Gjx7Nfffdx6ZNm+jduzf33HMPH3/8MbNnW5Lf+fv7s3//fgBGjBjBCy+8wMyZMyktLcVsNrNhwwYSExPZvXs3UkqmT5/Otm3buP766xt8DzoaU2ER+WvXoouPp/TYMYSbG943TcE3Ls4u5upe+Xsac789TInBBECaroQ/LT+IA2AwS4ZH+vHilBgm9g3GUTVHKU3UYOAQQiRIKWcJIQ5zaTSVlZSy7skp2rmaQaOx5U1Vtalq165d3HPPPdZaxfDhw4mMjGzW8U6ePMmRI0e44YYbAEuzUn21jZpNVQcPHiQyMtI6++C9997LwoULrYEjLi4OgIKCAtLS0pg5cyYArq6WsfobNmxgw4YNDBkyBLDktEpMTOw0gaP02DFL38Xq1ZiLi3GJjiZ43iv4TJuGxsurrYtn9fb6k9agUclklrg4a1j1yCj6hnq3UcmU9qyxGkdlOtipti5Ia2qsZpA4bjzG9PRayx1DQ+n+1ZctUoaRI0eSk5NDdnY2AB4el5oyHB0dMZsv3aBfWscMbmDpfO3Xrx+7du1qkTJVVbU89Z177ty5PPzwwy1+bntlLi4mf9068uITKD18GOHqivfkyfjGzcJ10KA2r11UZTZLdidfJE1Xdxr9knKTChrKZWuwbiqlrEyd/piUMqXqD/CY7YvXPEKIaUKIRXq9/oqOEzRnNsK1+l2wLT2HwYkTJzCZTPj7+9daFxwcTFZWFrm5uZSVlbFmzRrruqqp1qOjo8nOzrYGDoPBwNGjR5t0/ujoaJKTk0lKSgLgq6++qnNmPy8vL8LDw61T25aVlVFcXMykSZP4/PPPKSwsBCAtLY2srKxmvAPtR+nJk1x4/XUSrx9NxsuvYC4pJvill+j181ZC3/g7boMH20XQkFJyNF3P39cdZ9Sbm7lj0a/UV6pQrf3cL6K0P03tHL8BqDn70eQ6lrWplhqOa6s5DCr7OCrKypIlS+ocseTk5MS8efMYPnw4YWFh9OnTx7ruvvvu45FHHrF2jq9YsYKnnnoKvV6P0Whk9uzZ9GvCyB1XV1e++OILbr/9dmvneNX50Kv66quvePjhh5k3bx5OTk4sX76ciRMncvz4cUaOHAlYZgpcunQpQUFBl/PW2B1zSQn5P/yILj6ekoMHEc7OeN04Cd+4ONyGDrWLQFHpXG4xqw6ksepgOklZhTg6CEb3DuTFm2IoLTfy6vfHqjVXuTlpeG6Sfd07orQvDaZVF0I8iqVm0QM4XWWVF7BDSmmXc3LExsbKvXv3Vlum0nu3jfb2vpclJZEXn4B+1SrM+fk4R0aijZuFz80341gxCs4eZBeUsfZQOqsOpvP7OR0AwyP8uHlIKFP6h+DrcSl7bl2jqmYMCWuroit2rKXSqv8P+AF4A3ihyvICladKaY/quj/Ha+JECtavJy8+gZJ9+xBOTnhNnIg2bhbuw4bZTe2ioNTAhqOZrDqYzo6kHExmSUyINy9M7sO0QaGE1dP8NGNImAoUSotqMHBIKfWAHrizdYqjKLZT1/056S/MRcx7FVlSgnP37gQ99xw+M2fg6OfXxqW1KDOa+PlkNqsOpLPxeCZlRjPhvm48MroHNw8Oo3ew/YzgUjqPpvZxKEq7l7XgPWvQsDKZQEq6Lf4C9+HDEQ5tfy+D2Sz57exFVh1IY93hDPJLjfh5OBM3rCs3Dw5laDdfu6kFKZ2TChxKp1CemlrnEGsAWVaGx9VXt3KJapRBSo6m57PqQBqrD2ZwIb8Ud2cNk/p1YfrgUK6NClD5ohS7oQKH0mFJo5HCbdvI+3oZRdu317udYxNTtNhCck4R3x9MZ9WBNE5nF+HoIBgTHchLN8UwISYYN2eVJ0yxPypwKB2OITML3Tcr0CUsx3jhAo5BQQQ89hgOWi3Z//xnteaqlr4/pymyCkpZeyiDlQfSOXi+YkRUpB/3X9uDyf27VBsRpSj2SAWOVlSZVt1oNBITE8OSJUtwd3fnwoULzJ49mz179qDVagkODua9997D2dm5WrLD+ixatIh3330XAG9vb959912uvfbaerc/cOAA6enpTJkypUVfX1uSZjPFv/5K3tfLKNi8GUwmPK65huCXXsRr7FhERcZhR61Pi9+f0xQFpQbWH81k1YE0diTlYJbQN8SbuRUjotQNeUp7ogJHPU79doFdq05TeLEMTz8XRt7ck94julzRMavmqrrrrrv45JNPmDNnDjNnzuTee+9l2bJlgCWPVGZmJl27dm30mGvWrOHf//4327dvJyAggP379zNjxgx2795Nly51l/fAgQPs3bu3zsDRWFp3e2PMy0P/3Up08fGUp6Sg0Wrxu+9efGfNwrl791rb+0yb1iqBAiwjoracyOb7g2lsPJ5FudFMVz83HhsTxc2DQ+mlRkQp7VT7uUI0QUvNAHjqtwts+e8JjOWWfFGFF8vY8l9LNt0rDR6VrrvuOg4dOsSWLVtwcnKqdtf2oEGDAOrNklvVW2+9xdtvv01AQAAAQ4cOtSYs/Mtf/sKePXt4+umnKSoqwsXFhZ9++ol58+ZRUlLC9u3bmTt3LsePH+f06dOcOXOGbt26sXTpUl544QW2bt1KWVkZjz/+uDUnlT2kUrdMv3oAXfyyS9OvDh1K6BOP4zVxIg4ubZcW3GSW/HYml1UH0ll3JIOCUiP+Hs7cOawr0weHMbSbVo2IUtq9DhU4mppy5JeEU+ScL6x3feZZPSZj9TvqjeVmNn91nKPb6x6ZE9DVk+tm9W5SOY1GIz/88AM33ngjR44c4aqrrmrSfnU5evRorf1jY2NZsmQJ5eXlxMXFER8fz7Bhw8jPz8fd3Z3XX3+dvXv38uGHHwKWGQKPHTvG9u3bcXNzY9GiRfj4+LBnzx7Kysq45pprmDhxIomJiW2aSt1UWET+6u/JWxZvmSDJwwPtbbdVTL/atPfeFqSUHEmrGBF1KJ3M/DI8aoyIUinLlY6kQwWOllIzaDS2vKmq5qq67rrruP/++/nkk0+u6JgNOXnyJCEhIQwbNgyw9H/UZ/r06bi5WdrZN2zYwKFDh1ixYgUAer2exMTENkulXnriBHlfL7uUwrxvDF1efw2fm27CoZEsvrZ0NqeI7w+ks+pgGmeyi3DSCEb3DuKVqaGM76NGRCkdV6cMHI3VDJa8uIPCi2W1lnv6uTDzT0Mv+7xV+zgq9evXz3qBvhx9+/Zl3759jBs3zrps3759TUp0WFXVNOpSSv71r38xadKkattUTnfbGqnUzaWl5P/4I7qvl1mSDLq44D1lCr533oHrgAFt1tyTVVDKmoMZrDqQxsFUPULAiEg/HrzOMiJK665GRCkdn6o/12HkzT1xdK7+1jg6OzDy5p4tfq5x48ZRVlbGokWLrMsOHTrEL7/80qT9n3/+ef785z+Tm5sLWDq+Fy9ezGOPPUZ0dDQZGRns2bMHsEzKZDQaq6Vmr8ukSZP4+OOPMRgMAJw6dYqioqJWSaVedvYsmW+8SeLoMWS8MBdTfj7Bc1+g17afLSnMBw5s9aCRX2ogYe957v70N67++yZeX3MMo1ny4pQ+7PjzOJY9NJI7h3dTQUPpNDpljaMxlR3gLT2qqi5CCL777jtmz57NW2+9haurKxEREbz33nuApbkpPDzcuv2CBQu4/fbbrc+nT59OWloao0aNQgiBl5cXS5cutc4EGB8fz5NPPklJSQlubm5s3LiRsWPH8uabbzJ48GDmzp1bq0wPPPAAycnJDB06FCklgYGBrFy50map1KXBQMGmzeQtW0bxr7+CoyNeN0zAN+4O3EcMb5PaRanBxNaTWaw6kM6mE5YRUd383Hl8rGVEVFSQGhGl2I+1Z9by/v73uVB0gS4eXXh66NPc1OMmm52vwbTq7ZVKq24/GnrfDenp5C1fjm7FCkzZOTiGhuA7Kw7trbfgGBjYyiW1jIj69Uwuqw6k8cORCxSUGgnwdGbqwFBuHhzK4K5qRJRif9aeWcv8nfMpNV26sdVV48r8UfObHTxaKq26orQoaTJRtH07ecviKfz5Z5ASz9Gj0d4Rh+d11yHqmNiqJdWcm+LZib3pGeTJqgPprD6YTlZBxYio/l24eXAY1/T0VyOiFLtTbionoyiDtMI03vjtDcIv9GPEual4lvtS6JzHb93W8P7+921W61CBQ2kVxpwcdN98iy4hAUNaGpqAAPwfehDf22/HKax15opY+Xsac789bJ0NL01XwjMJB5GAk0YwJjqIGYPDGB8ThKuTGhGltJ1SYynpRemkF9b4qViWXZJt3TYq+ypGn7kDJ7Olj82r3I/RZ+7gZ5bZrHydKnBIKVVTQysym83IsjLSnnmG/J82gsGA+4gRBD37J7zGj0c4t15nspSSv607Xm0KVQAJaN2c+Pm5sfi4O7VaeZTOrdhQXC0QVH2cVpjGxdLq8+Q5OjjSxb0LYZ5hXBt2LSGeIYR5hhHqEcqWN85bg0YlJ7Mzo1Jvtln5O03gcHV1JTc3F39/fxU8bEyaTBgvXiQn4wKm/fsp3L4Dvz/ciTYuDpcePVqtHCXlJnaezmHj8Sw2n8gku6D2EGsAfYlBBQ2lRRWWF5JWmGZtTkovTK/2WFemq7a9s4MzIZ4hhHqEMrbrWEI9Qwn1DCXMM4wQjxAC3QLROFyqBZcWGkg5mkvylhzcy+q+P8u9tP77tq5Upwkc4eHhpKamkp2d3fjGymWR5eWYi4oxl5aAyYRGryfEzx+/bT/j4OraKmXIzC9l0/EsNh3PZMfpHEoNZjycNVzXK5Bfz+aiKzbU2kclGFSaQ0pJfnl+7RpDlef55fnV9nHVuFoCg2co/f37W2sMIR6W3/5u/jiI+vvSpJTkpheScjiX5MM5XDitR0pw83bG0UWDscxcax9PP9t95jpN4HByciIyMrKti9HhmIuL0a9di+7rZZQeO4ajuzs+U6fie0ccrjfcYPvzmy0TIG08nsmmE5kcSbN8YMO0bsTFdmV8TDAjevjh4qip1ccB4Oak4blJ0TYvp9J+SCnRlenq7GNIK0ojozCDQkP1lEVujm6WpiPPUAYFDrI+DvWw1Bz8XP2a3dJhMphJS8wj+ZAlWBTkWkZNBXT15KrJEUQMDCComxeJezKr5dYD2913Zj2+zY6sdGhliYnkLYtHv2oV5sJCXHr1InjeK/hMm4bGy7b3OJSUm9iRlMOmE5lsOp5FVkEZQsDQbr48NymaCTHB9A72rPVBnTHE0glfdVTVc5OircuVjqOh+xqklOSW5pJRmEFaUVqdHdAlxpJqx/N08rQ2HQ0LHnapGckzhDCPMHxcfFqkCbxIX0bKkVxSDudy7vhFjGUmHJ0cCI/x46obu9O9vz+evtVrEq1531mlTnMfh3LlzOXlFKzfQF78Mkr27kM4OeE1+UZ877gDtyFDbNp3lKEvYfOJLDYdz2JHUg5lRjOeLo5c3zuAcX2CGRsdiL9n22XFVezH2jNreXXnq5SZLvVpaYSGHj49MEojGYUZ1e55APB29rbWEiqbj6wd0J6heDvbpr9ASknO+UKSD+eQfCiHrBRLRgdPXxe6DwggYoA/4dG+OLZS3rMOcx+HEGIGcBPgDXwmpdzQxkXqdMrPnUOXkIDum28x5eXh1K0bQc89h88tM3H09bXJOc1myZF0PRsr+iuOpluaoLr6uXHn8G6MjwliRKQ/zo7qHovOymg2klaYRrI+mbP6syTnW34fyDqAmept/iZpIjk/mTFdx3B92PXWzufK5iRPZ89WK7eh3ETqiTySD+eQciiHIn05CAiO8GbE9B5EDPTHP6x2jdme2DRwCCE+B6YCWVLK/lWW3wi8D2iAT6WUb9Z3DCnlSmClEMIXeAdQgaMVSKORwq1byVsWb5mvW6PBa9xYtHfcgcfIkQiHlr9gF5cb2ZGUy6bjmWw+YWmCcqhogvrzjX0YHxNEryD7/kApLS+/PL9WcEjWJ5NSkILRbLRu5+fqR4R3RK2gUcloNvLumHdbq9jVFFwsJeVwDsmHc0k9mYfJYMbJRUO3vn5EDAygWz9/3L3bT64zW9c4FgMfAl9WLhBCaICFwA1AKrBHCPE9liDyRo39/yilrMyi93LFfkoL0a9eXWsaVffhw9EtX4Fu+XKMmZk4BgcT8MQTaG+/Dafg4BYvQ4a+xDoKaufpXGsT1OjegYyPCWJMdBB+ag7uDs9kNpFelG4NEGfzz1of55bmWrdzFI509e5KpHcko7uOJtInkkifSCK8I/Bx8QFg4oqJZBRl1DpHFw/btfnXZDZLspLzLU1Qh3PJTbV0pnsHuNLvulAiBgQQ2kuLpp3WmG3exyGEiADWVNY4hBAjgflSykkVz+cCSClrBo3K/QXwJvCTlHJjA+d5CHgIoFu3blelpKS04KvoePSrV5PxyjxkaZW2XgcHkBKkxOPaa/G9Iw7PMWOs83W3BLNZcihNz+bjmWw8nsWxDEsTVDc/d8bHBDEhJphhEX6qCaqDKjIUkaxP5oz+jLX2cFZ/lnP55yg3l1u383HxIdI7slpgiPSJJMwrDCeHhu+5acncTc1RXmLk3LGLpBzOIeVoLiUFBoSDIKSnDxEDAogY6I822N2ua8z23McRBpyv8jwVGNHA9k8CEwAfIUSUlLLOmY+klIuARWDpHG+hsnZYWQveqx40AMxmHDw9iPz2W5y7dWuxcxWXG9memMOm41lsPplFdkUT1FXdfXlhch8mxATRM1A1QXUUZmnmQtGFWk1LZ/VnySq5lIZfIzSEe4UT6R3JtWHXVgsSvq6X33dWGRxaI1usPrvYOlw2PVGH2SRxcXekWz9/Igb6062vP64eHe/mUrvvHJdSfgB80Nbl6CiklBTv3oMxve4pcM1FxS0SNNJ1JWw6cakJqtxoxsvFkeujA5kQE8SY3kH4qiaodq3YUExKfkqtpqWU/JRq3/a9nLyI9Ink6tCrLcGhoibR1asrThrbXFRv6nGTTQKF2WQm47Se5MO5pBzOIe9CMQC+XdwZNL4rEQMC6NLDG4cOnhizLQJHGtC1yvPwimVXTAgxDZgWFRXVEofrUIx5eei/W4kuIYHy5GQQwtIsVYNjxTwezWU2Sw6m6th8IouNx7M4XtEE1d3fnbtHdGdCTBDDIv1w6uAfqI5GSklmcWa12kPl4wtFF6zbOQgHQj1CifSJZHjIcGuAiPCJwN+1faf5KS0ycO5oLsmHczl3NJeyYiMOGkFYby39rg8jYoA/PoHubV3MVtUWgWMP0EsIEYklYNwB/KElDiylXA2sjo2NfbAljtfeSSkp2buXvPgECtavRxoMuA0ZQsgjDyPNZjJf/0u15irh6krQnNlNPn5RmZHtSTkVo6CyySm0NEHFdvdj7uQ+jI8Jpmegh11dNFp7wht70JTXXGostdQe8s9Wa1pKzk+udjOch5MHkd6RxAbHVut76ObdDRdNx7iPRkpJXkZxRcd2lfQeXk5EDg4kYoA/XWP8cHa1+wYbm7H1cNyvgTFAgBAiFXhVSvmZEOIJYD2WkVSfSymP2rIcnY1Jp0O3ciW6hOWUnzmDg5cX2lmz0M6ahWv0pfnWHZycao2q8pk2rcFjp+lKrB3bu85UNEG5OjImOojxfYIYEx1ot1Oo1uw0zSjKYP7O+QAdNnjU9Zrn7ZjHr+m/4unsaQ0O6YXpSC7VQCtrD0ODh1qbliJ8Igh0C7SrLwItxZreo6IJKj+nenqP7gP8Ce7ujXDoeK/9cnSoO8erNFU9mJiY2NbFaVVSSkr27ycvPp6CH9cjy8txGzQIbVwc3pNvxMHt8hL5mc2SA6k6Nh/PYuPxTE5csNzZGhngwfg+QYyLCWJYhP03QV0svcgtq26pNrSzkqvGlXHdxlkunBUfB1n5r+LzUXlRrfq82uM69q3cvupzy2Z17Ft1XY19a563oePWXHcq71S1ex2qcnN0I8I7ggifiGp9D928u+Hm2PETPxbnl5NyxDJc9vyxixjKTGicHOjax5fuAwLo3t8fLxsmCrRHTR1V1aECR6XOlHLEpNejX7WKvIQEypNO4+Dpic/0aZbaRZ8+De5bcza8yrxNRWVGfknMZtPxLLaczCKnsByNgyC2uy8TYoIZVzEKyh4VlheSpEu69JOXRKIusdb8BjV187IMCBBCILj0rbLqc4Go9m27rnU197Wuq7KvQFC5WZ3rau5b4zgN7ltj3S9pv9T5egWCA/ccaDAja3t26rcLtXI39RoeTE5qIcmHLMEiK9nSD+ehdSFigD8RAwII6+OLUyul97BHKnB04MAhpaTk9wPo4uPJ//FHZFkZrgMH4hs3C+/Jk3Fwb7yjrq5MsU4aQY9AT85mF1FuMuNd2QQVE8To3vbVBFVqLOWs/ixJOktgSMqzBIqqN365OboRpY2y/nx25LM6A0iIRwgbbuuYCQnquxmuI7/mU79dqJUtVjgInFwcKC8xWdN7RAzwp/uAAALC1VDwSvZ8H4dymUz5+ehXfY8uIYGyxEQcPDzwmTkD37g4XGNiGt1fSomu2ECaroTX1xyrNRuewSQ5nVXI/10Twbg+wcRG+LZ5E5TRbORc/jlLcKioQSTpkjhXcA6ztFwYHB0c6eHTgyFBQ5jlO8saKEI9Q6t9o/Z386/zxrCnhz7d6q+rtTw99OlO9ZrNZsmOb5KqBQ0AaZaYjZJx9/She/+AdpXewx51qMDREYfjSikpPXiQvPgE8n/4AVlaimv//nT5y+v4TJmCg4eHdVuzWZJdWEZqXglpuhLS8kpI0xVbnlcsKy43NXA2MJklL93U19YvqxazNJNemG5tYkrMswSKs/qzGMyWyZcchAPdvLoRpY3ixsgbidJG0Uvbi67eXRu9mxha98Ywe9HRX7PZLMlNLSTtVB5pp3RkJOkoKzdDynIAABKdSURBVK67T8doMBMzKrSVS9gxqaYqO2UqKEC/ejW6+ATKTp7Ewd0dz5tuwjDlZjKCIqoFhsrH6bpSyk3Vv2lp3Z0I07pZfnzdCPd1J0zrxiurjtQ5lWqY1o0dL4yz2euSUpJTklOteanyp+qwzxCPEEvNwdcSHKK0UUT6ROLq2Lk6K5XqGgoUPoFuhPXWcuZANqVFtYOHp58L9/79mtYucruimqraoZJyI6m79lGwPAHnbZvQlJeRHRLBzvH3sCF4IMklYF51Abh041WQlwthvm70D/NhUv8uhFcJEKFaNzxd6v4vLjWYbD4bnr5MX62DujJA6Mv01m38XP3ope3FLb1usTYx9dT2xMvZtpNBKe2DNEty0gpJO1l3oOg5JJDQ3r6E9dZaJzgK6+3b6jPidTYqcLSi/FKDpZZQ2ZRUUVPIzrpIxO+/cP3J7fTUp4PGmQ3hQ1jfYySF3XsR5uvGYK0bU30tQSFM606YrxshPq64Ol3eCJCWnA2v2FDMWf3ZarWIRF0iWcWX8hJ5OnkSpY3ihu43WJuYemp74u/mf1nlVzqmykCRfkpH2qk80hMvBQrvQDd6DAkkrLcvob209Q6VbYsZ8TqbDtVU1Zb3cUgpuVhUbg0Glf0Ml/obiskvrV597pufxi1puxl2Zi/OhjIKwiMpnDgNzylTCQsPINjLBUc7uj/CYDKQnJ9crQ8iSZdEakGq9R4CF40LPXx60Mu3l7UG0cu3F8HuwWrkilJLY4EirLe20UChtBw1HLeZfRz13dNQyWSWZBWUWmsL1QJDXjHputJao5Q8XRwJ07oRbq0puBHuKum6/xdc16/GdOI4ws0N7ymTLSOjBgxo1YtrfakoTGYTaYVptfohkvXJGKXlQ60RGiK8I4jyjbLWIKJ8owj3DEfj0HnHwSsNk2ZJbnohaSdVoLBHKnA0I3DUdU+Do4NgaDctGgcH0nQlZOhLMJiqv1d+Hs41Op6rPNa64+3maA0EJUePootPIH/NGszFxbhER6ONm4XPtGlovFq3Pd9gMrAyaSVv7Xmr1rzMwe7BXCy9WG34ZrhneLVO6ijfKCK8I3DWqCGNSsMaDBQBroRV9E+E9vZVgcIOqM7xZnh7/UkMbnvx6LYe4aRDGrSU/f/27jxGzvq+4/j7O3vNzsweM+O1d7w2PoKBkHCFK5AoohKpsSi4RSlOeqVJlKhRowSkRNBTqdKmh0oTtaSNLKVxEkjAKghMoFCamoSWoOIC5jSGImhsr1m8h71z7e7M/PrH8+zsMbOzO7M7e35e0sje53lmnmdWz85nfs/z+31/7+7k0NuXcMlZUS7e3Mn1FyaKobA52srGzlZCzZV/fYVUiqFHHmHo3v1kX3oJCwZp37WL6J6bCV500YytC+ccI/kRMrlM2Uc6lyYzVn7djI9J24+3GqbLuzz92X72nLuneIlpe8d2Qk1rq/Kn1K4kKN4YYiQ1ERTbL+5SUKwCqyo4ah3H0Vd4imDifizgjRew5iGCifsZAe77/B+XfU7BFUiPpb0P8Wkf0mNHjtL80BOEDz5LQ2aE1OY4v/jkh3nzg5sZbimQPvNDMge/M+MHezafLQ5um6vWxtayj/ZQe9nldz5/Z9nXGc2P8pXLv1LVvmXt8oIi5YXE0SGOvz44NSguUlCsRqsqOGotq9664d9wfmiMs8AYwcR9fPqxN8t+u598KQegZdRx9auOa58rsKMXRhvhP88zHr+kgaM9QwQC/01r74slH+ChxhDxYLx0eVNoxjCY/gg2BquuOXTf6/ct+bzMsvLMFhTbxoNiRyft8dVfKHGtWlXBUSvXOFR+RSBHvpCno6WD7nD31A/sJu/f6LEzbHj8MB0HnyeQzpLf2kPhS79M5Pqd7I6t5+P+ts2B5mXVq2itlaKQ2igopBwFB5AId89YCO57u75XsryQyXDmXx9l6N57yRw+jDU303bdTqJ79tD6gQ8sq4CYyWovRSG1cQXHQG+qODL7xNEhsimvNa6gkHEKDub+7Tt79ChD9+7n9IEDFIaHad6+nQ1/cDvtN95IYzS62Ic9b/Wal1mWl3IlxscHw1UKirZ4kK0XxovdY9vXKSjEo+Cg8rfvQjbLmUcfZeje/WSeew5raqJt506ie26m9bLLVkTrQtau6SXGkwMj/MddR/i/V/sZGykoKKQmqyo45lMd98MvFzjnH/PkenM0JvJEP9HLyR99ndMPPkjhzBmat25l/W230fGru1dk60LWpqce+N+SEuP5sQKvPf0ObTEFhdRGAwCB0w89RO+f/CkuO7WnFIEA7dddR+eePYSuuFytC1m28vkCQ++k6T+epP9Yyvv3eJLkYGkF5HG//+36VUGWlUkDAKvQ941vloYG0NjVRc/f3bEERyRSnnOO9JnRKQFx6niSwZMpCjnvS2AgYEQTITbu6OStF/sZzZQvMS5SKwUHkOst7VEFkOvrK7tcZDHkRvMM9KaKIXHKb0VkkxNjjsIdzcQ3RTjr/BjxngjrNkXo3BCiodEb11NuGlWVGJf5UnAAjYkEuRMnyi4XqTfnHMP92eLlpVN+S+J0X5rxK8mNTQFiG8Nsu2idFxA9EeI9EYKRyjMfqsS41IOCA1h/6y0l9zgsGGT9rbcs4VHJajSSydF/PMnA8SSnjqfoP5ak/0SSsexEgc32rlbiG8Ocfdn6YkC0d7USCNR2j+2cK7sVFLKgFBxAxw03AN69jlxvL42JBOtvvaW4XKRahXyBob5MsRXR74fE8MDEl5OWUCOxjWHOu7Kb+CYvIGIbwzQH9Wcpy5t6VYnMU/Fm9aSQGDiRIp/z7itYwIh2h4hvDBcDIt4TIRJtUU89WVbWZK+q+YzjEJlNbizPYG+62JNp/HJT5sxocZtQu3ez+oJreiZaEd1hGpqWz0yOIvOlFoesKZXKb4xzzpEcHKH/WLLYk6n/eIqhd9K4gvf30tDo3ayO94S9FsSmCPGNEULtmtxKVq412eIQqaRc+Y2Ddx1h6N00ofaWKZeaJo99aIsHifdE2H7xumKX146uVgLLaD54kcWk4JA1ITea57/uf6Ok/EZurMAzP34LgKZgA+t6Ipxz+YZiSyLWE6GlVX8mIpPpL0JWldFsjsGTaQZ7UwyeTDHQm2agN8WZUxmocFX2t//8KtriQd2sFpkDBYesSNnUGIO9KQZ6U8WgGOhNTanNFGgwOjeE6NrcxrlXbODFnx6fMup6XCTWogJ/IlVQcMiy5ZwjMzzmhYP/GPBbEZN7MjU2Bejs9mozRRNhYt1hookQ7V2tNEy6D9G5PqTyGyILQMEhS268F1O5FsRIeuImdVOwgVgizJb3x4l2h4glwsQSYdpiQWwOo6pVfkNkYSg4ZNEUCo7h/gwDvemJFoQfFGMjEyU3guEmookQZ1+6nmi3Fw7RRJhw5/znbVf5DZH5U3DIgsvnC5zuy5TcoB56J01+bOIyUaijmVgizHlXJ4h1h7zLTIkwrW0aCyGynCk4pGa5sTxD76T9exATl5dO92UoFCa6MLXFgkQTYTadFy1eXop2h2gJVa7sKiLL06oKDpUcqc5cRlHDpC6uJ8cvL3lhMXwqUyz7beZVdY0lwmy7qItYwmtBdG4IqWifyCqjkiNrVLkJfhqaAlzy0c20xVqLPZkGTqZIDpR2cY36PZfGWxAd61tpbGpYirciIgtEJUekRD5XIDk4QnIwy5P7j5aMos6PFTj0yNuAFyLR7hAbz+6cdIO6tIuriKw9Co5VwhUc6eFRkgMjDA9kSQ5mSQ54IeH9PEJ6eLTi6Olxv/U1bxR1rRMHicjqpuBYAZxzjGZyDPtBkBwsDYfk4AiF/NRUaGwO0BYLEom2ePM/+P9viwX5932vkD49WrKvSKyFji6NohaRmSk4loHcWH5S62AiCJIDWYb9fyePcwAIBIxwZwuRWAsbtnVw9qUtRKLBKeHQEmqccdzDh246W6OoRaQmCo46KxQc6dMjE62EgRGGB7Mk/ctHycEsmeHS+kmtbU1eN9YNITafF50SCJFokFBH87wuJWkUtYjUSsHhm2vX1Mmcc4ykclOCYPx+gtdayJIaGi1O/jOuKdhQvITUtaWNtmiLHwzeski0ZVF6KGkUtYjUQsHBDBP83H2EXK5A4j0dJa2EyeGQG5vaMynQaEQ6vZZBz44okVhLMRDaYt6lJM3vICIrmT7B8C7XlEzwM1rg4A+OTN3QINzeTCQWJN4TZssFcdqiwWI4tMWCtEaa5lRwT0RkpVJwwJQBbtNd+6nzafODIdzZQkOjxjCIyNqm4MDrglouPCKxFs7VPQARkSn09Rm4avd7aGye+qtQ11QRkfLU4kBdU0VEqrHsg8PM3gt8CVgH/MQ590/12I+6poqIzE1dL1WZ2T+bWZ+ZvTRt+XVm9pqZvWFmt1d6Defcq8653wNuBj5Uz+MVEZHZ1fsexz7guskLzKwB+BawCzgf+ISZnW9mF5jZj6c91vvPuRF4GHikzscrIiKzqOulKufcz8xs67TFVwBvOOfeBDCze4Ddzrm/BH5lhtc5ABwws4eBH9bviEVEZDZLcY+jB/jFpJ+PAVfOtLGZXQPcBLRQocVhZp8DPgdw1llnLcRxiohIGcv+5rhz7gngiTlstxfYC94MgPU9KhGRtWspguM4sHnSz5v8ZfM2Puc4kDazVyts2gGcnmHdOuDUQhzPIqv0npbzvmp9rWqfV832s207n/U6vxZ3X/N5rXqdY3PZrtI29Ty/tsxpK+dcXR/AVuClST83Am8C24Bm4DDwvgXe595a1wOH6v07qdPvueJ7Xq77qvW1qn1eNdvP5/yZbb3Or8Xd13xeq17n2Fy2m+UcWvLzq97dcX8E/Bw418yOmdlnnHM54AvAY8CrwH7n3MsLvOuH5rl+JVrM97SQ+6r1tap9XjXbz/f80fm1fPY1n9eq1zk2l+0qbbPk55f5KSU+MzvknLtsqY9DViedX1JPi3V+qVZVqb1LfQCyqun8knpalPNLLQ4REamKWhwiIlIVBYeIiFRFwSEiIlVRcMzCzN5rZt82s38xs88v9fHI6mNmYTM7ZGZla7WJ1MrMrjGzJ/3PsGsW6nXXZHBUU+7dqay7VKmG6QRuA/Yv7lHKSlXl+eWAJBDEqwu4INZkcFBFuXd/ncq6SzX2MffpBD4KvAL0LfZByoq1j7l/fj3pnNuF9+XkzxbqANZkcDjnfgYMTFtcLPfunBsF7gF2+9sf8H/5v7m4RyorUZXn1zXAB4HfAD5rZmvyb1LmrprzyzlX8NcP4lUYXxDLvjruIipb7n2uZd1FZlH2/HLOfQHAzH4XODXpD12kGjN9ft0E7AQ6gTsXamcKjlm4OZZ1F5kP59y+pT4GWX2cc/cD9y/066pZPKFu5d5F0Pkl9bWo55eCY8IzwA4z22ZmzcDHgQNLfEyyeuj8knpa1PNrTQbHEpZ7lzVA55fU03I4v1TkUEREqrImWxwiIlI7BYeIiFRFwSEiIlVRcIiISFUUHCIiUhUFh4iIVEXBIQKYmTOzOyb9/GUz+2qd9vWH035+apbt95nZx+pxLCK1UHCIeEaAm8xsXb12YJ4AMCU4nHNX12ufIvWg4BDx5IC9wK2VNjKzr5rZD8zs52b2upl91l8eMbOfmNmzZvaime32l2/1J9f5PvAS8B2g1cyeN7O7/W2Sk17/Nv/5h83sr8rs/1Iz+6mZ/Y+ZPWZmCX/5F83sFTN7wczuWaDfiUhZqo4rMuFbwAtm9jezbHch3hwaYeA5M3sYbyKmX3POnfFbLU+b2XitoB3AJ51zTwOY2a875y6e/qJmtgtvjo4rnXNpM4tNW98E/APePAvvmtke4C+ATwO3A9uccyNm1lnb2xeZGwWHiM//0P8+8EUgU2HTB51zGSBjZgfxJtF5GPi6mX0EKODNj7DB3/7t8dCYxbXAd51zaf94pk/Wcy7wfuBxMwNoAHr9dS8Ad5vZA8ADc9iXSM0UHCJTfRN4FvhuhW2mF3hzeLNDdgGXOufGzOwtvHmeAVILdGwGvOycu6rMuuuBjwA3AH9kZhf4he9EFpzucYhM4n/L3w98psJmu80saGZxvKlfnwE6gD4/NH4J2FLh+WP+ZafpHgc+ZWYhgOmXqoDXgC4zu8pf32Rm7/NvuG92zh3Em1u6A4jM9l5FaqXgECl1B1Cpd9ULwEHgaeBrzrkTwN3AZWb2IvA7wJEKz9+Ldy/l7skLnXOP4s2hcMjMnge+PG39KPAx4K/N7DDwPHA13iWru/x9Pwf8vXNuaK5vVqRaKqsuUgV/bEfSOfe3S30sIktFLQ4REamKWhwiIlIVtThERKQqCg4REamKgkNERKqi4BARkaooOEREpCoKDhERqcr/AwF5S/rOy9LNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "label = {'kd':\"KDTree\",\n",
    "        'cg':\"Cellgrid\",\n",
    "        'ns': \"Cython NS\",\n",
    "        'bf': \"Brute Force\",\n",
    "        'oct': \"PCL Octree\"}\n",
    "\n",
    "for key in ['kd', 'cg', 'ns', 'bf', 'oct']:\n",
    "    ax.plot(Nsolvent, timing[str(key)+'.'+str(maxdist)], '-o', label = label[key])\n",
    "ax.set(xlabel='N particles', ylabel='time [s]', xscale='log', yscale='log', title='Total Time')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For very low number of particles, brute force is a good method. However, if multiple selections are required, then it is not feasible to recalculate pair distance matrix again. Other option to store the pair distances as one of the triangular matrix is also not a good option, as the memory consumption increases with the number of particles. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For lesser density of particles, while it looks like octree is a good option, but it is a non-periodic implementation. While, similar tweaks can be adapted from periodic KDtree implementation, the scaling cannot be predicted due to different traversing strategies in KDtree and Octree. In any case, it would be interesting to test an efficient periodic implementation of octree. For the range of number of particles generally considered in Molecular Dynamics, Neighbour search, definitely consumes less time as compared to periodic KDtree for fixed distance selections, also anticipated from single query benchmarks [here](https://github.com/ayushsuhane/Benchmarks_Distance/blob/master/Notebooks/OctreeComparison.ipynb).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another dimension is to test the execution time for different cutoff radius. Lets check the variation with cutoff radius for fixed number of particles as 40000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nprotein = 100\n",
    "Nsolvent = 40000\n",
    "center = box[:3]/2\n",
    "radius = 10.0\n",
    "maxdist = np.linspace(4, 10, num = 6, dtype = np.int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "    'Nprotein':Nprotein,\n",
    "    'Nsolvent':Nsolvent,\n",
    "    'maxdist':maxdist,\n",
    "    'centre':center,\n",
    "    'radius':radius\n",
    "         }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "636ced1064b2432eb9dadb411185f1f2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "64246bbf3661406787d9a620508c0c3a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=6), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "timing = timer(**params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f4c5e2f7f98>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNX5+PHPM0t29rAvguyE3QiioiwWrQqCG1q/X7V1/bVWsbVWbEXab1u12qLWrbQqVKyCCyhiFRdcQFD2fYlSUDDsEMg+y/n9cWcmM8lMkkkymWTyvF+v+5qZc7dzGXKeOefce44YY1BKKaWqyxbvDCillGpcNHAopZSKigYOpZRSUdHAoZRSKioaOJRSSkVFA4dSSqmoaOBQSikVFQ0cqskRkT0ickhE0oPSbhaRT2J4zhkiYkTkgqC0ZBF5QUROisgBEflFuX3Gi8gOESkUkWUiclpd7KtUbWngUE2VHbirPk4kIj2Bq4DccqtmAr2B04CxwL0icpFvn0zgTeABoDWwBphfR/sqVSsaOFRT9Shwj4i0rIdzPQ38Gigtl34D8H/GmOPGmO3AP4AbfesuB7YaY14zxhRjBYohItKvDvZVqlY0cKimag3wCXBPdTYWkROVLPdVst9VQIkx5t1y6a2AjsDGoOSNQJbvfVbwOmNMAfANkFWbfatzrUpVxRHvDCgVRzOAFSLyRFUbGmOirpmISDPgT8APwqzO8L3mBaXlAc2C1h8ut49/fW32VarWtMahmixjzBbgHSBijaGWZgIvGWP2hFmX73ttHpTWHDgVtL45ofzra7OvUrWmgUM1dQ8CtwCdK9tIRPIrWe6PsNt44E7fXU8HgK7AAhH5tTHmOFZn+ZCg7YcAW33vtwav890B1hOr76LG+1Z2jUpVlwYO1aQZY77GuuPoziq2y6hk+VOE3cYDA4GhvuV74DasznKAfwG/FZFWvo7rW4A5vnULgYEicoWIpGA1q20yxuyog32VqhUNHErB74H0KreKkjHmqDHmgH8BPMBxY4y/qelBrE7rvcCnwKPGmPd8+x4GrgD+CBwHRgLXBB2+NvsqVSuiEzkppZSKhtY4lFJKRUUDh1JKqaho4FBKKRUVDRxKKaWikpBPjmdmZpru3bvHOxtKKdVorF279ogxpm11tk3IwNG9e3fWrFkT72wopVSjISJ7q7ttg2+qEpHTReR5EXk93nlRSikVp8Dhm4DmkIhsKZd+kYjsFJGv/SOOGmN2G2Nuikc+lVJKVRSvGscc4KLgBBGxYw3F8ENgAHCtiAyo/6wppZSqTFz6OIwxn4lI93LJI4CvjTG7AUTkVeAyYFt1jikitwK3AnTr1q3O8qqUih2Xy8W+ffsoLi6Od1aajJSUFLp06YLT6azxMRpS53hn4Lugz/uAkSLSBmvMnWEiMt0Y81C4nY0xs4HZANnZ2TqOilKNwL59+2jWrBndu3dHROKdnYRnjOHo0aPs27ePHj161Pg4DSlwhGWMOQrcHvMTbVoAH/0e8vZBiy4wfgYMvjrmp1WqKSsuLtagUY9EhDZt2nD4cPl5vqLTkO6q2o81X4FfF19atYnIRBGZnZeXV/XGwTYtgMV3Qt53gLFeF99ppSulYkqDRv2qi3/vhhQ4VgO9RaSHiCRhDQP9djQHMMYsNsbc2qJFi+jO/NHvwVUUmuYqstKVUkqFiNftuK8AK4G+IrJPRG4yxriBO4D3ge3AAmNM/cxYlrcvcrrXWy9ZUErFR0ZGRuD9u+++S58+fdi7dy8zZ86kc+fODB06lN69e3P55ZezbZt1r86UKVMYOnQovXr1okWLFgwdOpShQ4fyxRdfxOsy6lW87qq6NkL6u8C79Zwdq08j77swKww8PggGXwWDr4F2/eo9a0qpMovW7+fR93fy/YkiOrVM5VcX9mXysEpn/a22jz76iDvvvJP333+f0047DYC7776be+65B4D58+czbtw4Nm/ezMKFCwH45JNPeOyxx3jnnXfCHtPtduNwNPiu5Kg1pKaqWqtxH8f4GeBMDU1zpMKZN0P7AbDiSXhmJPz9PFj5DOQfqrtMK6WqZdH6/Ux/czP7TxRhgP0nipj+5mYWrY+qKzSszz77jFtuuYV33nmHnj17ht1m6tSpTJgwgX//+9+VHqtLly7cd999DBs2jIULF5KTk8OFF17IGWecwXnnnceuXbsAOHjwIJdffjnZ2dmMGDGCVatW1fo66ktChUJjzGJgcXZ29i1R7ei/eyrSXVX5h2DLG7DxVXh/Oiz9LfQca9VC+l0CSWl1eyFKNUG/W7yVbd+fjLh+/bcnKPWENh0XuTzc+/omXvnq27D7DOjUnAcnZlV63pKSEiZPnswnn3xCv36VtyoMHz6cHTuqnrq9Xbt2rF+/HoCxY8fyz3/+k549e7JixQruuOMOli5dyp133sm9997LWWedxZ49e7j00kvZsmVLFUduGBIqcNTK4Ksj336b0Q7O+n/WcninFUA2LYA3b4akDOg/CYZMhe6jwWav33wr1USUDxpVpVeX0+nk7LPP5vnnn+eJJ56odNvqTrU9depUAE6cOMGqVau44oorAuvcbjcAH374ITt37gykHz9+nKKiIlJTy7V+NEAJFThEZCIwsVevXrE7Sdu+cMGDMO4B2LsCNr0K296Gjf+GZp3K+kPa62gpSkWjqprBOQ9/zP4TRRXSO7dMZf5to2p8XpvNxoIFCxg/fjx/+tOfuP/++yNuu379erKzs6s8Znp6OmAFmszMTDZs2FBhG2MMX331FUlJSTXOe7wkVB9HjW/HrQmbDXqMhsuehnt2wZUvQIdB8MVT8OwoeO5c6/2pg7HPi1JNwK8u7EuqM7RGn+q086sL+9b62GlpaSxZsoSXX36Z559/Puw2b7zxBkuXLuXaa8Pe2xNWq1at6NixY6Az3ev1snHjRgAuuOACnn766cC24YJLQ5VQgSNunKkw8Aq4bgH8cidc9AjYHLD0N/DXfvDS5VbTVmlBvHOqVKM1eVhnHrp8EJ1bpiJYNY2HLh9UZ3dVtW7dmvfee48//OEPvP229QjZrFmzArfjzps3j48//pi2bas111HAq6++ynPPPceQIUPIysoK3IH19NNPs2LFCgYPHsyAAQP4xz/+USfXUR+kum12jUl2drZpEBM5Hd5lNWVtWmDd7puUAf0nwuCp0OM87Q9RTd727dvp379/vLPR5IT7dxeRtcaYqtvh0D6O2Grbx7o7a+xv4duVVhDZ+hZsfAWadYRBV8GQa6B95W27SinVkCRUU1W99nFEw2aD7ufApL9Z/SFXzYGOQ2HVM/Ds2fDsufDF3+DUgXjnVCmlqpRQNY5GwZkCWVOspeAIbHnTqoks/S18MANOH1P2fEhyRlVHU0qpeqeBI57SM2HkrdZyJAc2zbeWhbeCMx36X2r1h5w+RvtDlFINhgaOhiKzN4z7LYy5H75bZT1kuHWRFUgyOsCgK63+kA6D4p1TpVQTp4GjobHZ4LSzreWHf4Zd71nB48vnYOVT0C7Lekp90FXQvFO8c6uUaoISqnO8xoMcNlTOFMiaDNe+Ar/cBRc/Zo2L9cEM+OsA+NdlsOEVKMmPd06VatQOHDjANddcQ8+ePTnjjDO4+OKLA4MRhuMfin3Pnj0MHDgw6vPdfPPNgSHag82ZM4c77rgj6uPVt4QKHA32rqq6kN4GRtwCN38IP18H598Lx/4Li26Hx3rDG7fA1x+Cx209NzJrIMxsab3qTIYqUcTg/7YxhilTpjBmzBi++eYb1q5dy0MPPcTBg7EZ9cHj8fDPf/6TAQMa77BECRU4mow2PWHs/XDXRvjJ+1YHes77MO8K+HMPWPT/dBpclXhiNMXzsmXLcDqd3H777YG0IUOGMHr0aB599FHOPPNMBg8ezIMPPljpcQoLC7n66qsZMGAAU6ZMYeTIkfgfRM7IyOCXv/wlQ4YMYeXKlYwZMyaw7sUXX6RPnz6MGDGCFStW1Opa6ov2cTRmItDtLGv54SOw631481bwukO380+DG2n0X6Uagv/cBwc2R16/bzV4SkLTXEXw1h2wdm74fToMgh8+XOlpt2zZwhlnnFEhfenSpeTk5PDVV19hjGHSpEl89tlnnHfeeWGP88wzz9CqVSu2bdvGli1bGDp0aGBdQUEBI0eO5C9/+UvIPrm5uTz44IOsXbuWFi1aMHbsWIYNG1ZpfhsCrXEkCkcyDJgE7uLw6/O+g3X/gqIT9ZsvpepK+aBRVXotLV26lKVLlzJs2LDAPBw5OTkRt1++fDnXXHMNAAMHDmTw4MGBdXa7PWRodb8vv/ySMWPG0LZtW5KSkgLDsTd0WuNINJGmwbU54O2fw5J7oO9FMOhq6P0DK+Ao1RBUUTNg1sDw/7dbdIUfL6nxabOysnj99dcrpBtjmD59OrfddluNj+2XkpKC3Z44z2IlVI0j4e6qqolw0+A6U2Hys3DLMsj+Cez9AuZfB4/1gcXTYO9K8NZuMhylYi7S/+3xM2p12HHjxlFSUsLs2bMDaZs2baJ58+a88MIL5Odbdy3u37+fQ4ciTxt9zjnnsGCB1d+ybds2Nm+upNnNZ+TIkXz66accPXoUl8vFa6+9VqtrqS8JVeOo8dSxiaSqaXA7D4cJf4Ddn5Q9qb72RWjZzaqFDJ5qDc6oVENT1f/tGhIRFi5cyLRp03jkkUdISUmhe/fuPP7447Rs2ZJRo6xJojIyMpg3bx7t2rULe5yf/vSn3HDDDQwYMIB+/fqRlZVFVXd4duzYkZkzZzJq1ChatmwZ0i/SkOmw6k1dST7sWGIFkN3LwHitARgHT7XmGGnWPt45VAkskYZV93g8uFwuUlJS+Oabb7jgggvYuXNng5zhT4dVV7WTnGE9iT5kqjVb4ZY3rCDy/nRrIqrTx1pBRAddVKpShYWFjB07FpfLhTGGZ555pkEGjbqggUOVadYeRv3UWg7vtO6P37TAN+hiGvQLGnTRrv91lArWrFkzmkpLh/71q/Da9oXxD8DY38B3X1q1kK0LYfMCSG8LA6+02pY7DbOeJ1FKNRkaOFTlbDY4bZS1/PARyPnACiJrnocvn4U2va1ayOCroFX3eOdWKVUPNHCo6nMkW3OE9L8Uio7Dtretpqxlf7CWrmdZtZCsKZDWOt65VUrFiD7HoWomtRWccYP14NW0LTD+QSg+AUt+YT0f8sqPrKYtV1G8c6qUqmMJFTgSenTchqxlVxj9C/jpKrjtcxh5G+xfC6/daAWRt34G//1MHzJUDVK0Q6pv2LCBd999N/B55syZPPbYY3Wer5kzZ5KWlhby0KF/OHeAP/7xj2RlZTF48GCGDh3Kl19+Wed5iESbqlTdEYGOg63lB7+3gsWmBdZMhuvnQfPO1kyGg6dC+6x451Y1Qkt2L+GJdU9woOAAHdI7cNfwu7jk9EtqfDz/kOo33HADr776KgAbN27k4MGD9OkT/kHYDRs2sGbNGi6++OIan7e6MjMz+ctf/sIjjzwSkr5y5Ureeecd1q1bR3JyMkeOHKG0tDTm+fFLqBqHakBsdug5FqY8C/fkwBXPQ/uB8MVT8OzZ8Ow5sOIJyNsf75yqRmLJ7iXM/GImuQW5GAy5BbnM/GImS3bXfJyqyoZUv/7661m0aFEg/brrruOtt95ixowZzJ8/n6FDhzJ//nzAGmJkzJgxnH766Tz55JOBff76178ycOBABg4cyOOPPw5Ykz/179+fW265haysLCZMmEBRUfgm3Z/85CfMnz+fY8eOhaTn5uaSmZlJcrI11lxmZiadOtXfjKD65LiqXwVHYMub1p1Z+9cAAj1GW7WQ/pMgpXm8c6jqUfATzI989Qg7ju2IuO2mw5so9Vb8VZ1kS2Jw28Fh9oB+rfvx6xG/jnjMJ598kv/+97/MmjWrwrpPP/2UWbNmsWjRIvLy8hg6dCg5OTnMmzePNWvW8NRTTwFWk9LSpUtZtmwZp06dom/fvhw4cIBNmzZx4403smrVKowxjBw5knnz5tGqVSt69erFmjVrGDp0KFdffTWTJk3if/7nf0LOP3PmTDIyMigsLMTj8fC73/2OjIwM8vPzyc/P59xzz6WwsJALLriAqVOncv7550e8zvJq++S41jhU/UrPhJG3wi0f+WYy/DWc+M7qB3mst9UvsvM/4C7VmQxViHBBo7L02jr//PPJycnh8OHDvPLKK1xxxRU4HOFb9y+55BKSk5PJzMykXbt2HDx4kOXLlzNlyhTS09PJyMjg8ssv5/PPPwegR48egXGpzjjjDPbs2RMxH3feeSdz587l1KlTgbSMjAzWrl3L7Nmzadu2LVOnTmXOnDl1du1V0T4OFT9tesLY6TDmPti3xqqFbHnDuhvLmQ6eYvB6rG39s72BTkiVoCqrGQBMeH0CuQW5FdI7pnfkxYterNE5Iw2p7nf99dczb948Xn31VV58MfI5/E1GYM294Xa7I24bbvtITVUALVu25Ec/+hFPP/10SLrdbmfMmDGMGTOGQYMGMXfuXG688cZKz1tXtMah4k8Eup4JlzwG9+yCa+cD3rKg4ecqgqW/rZiumoS7ht9Fij0lJC3FnsJdw++q8TEjDanurxnceOONgb4J/xzhzZo1C/n1H8no0aNZtGgRhYWFFBQUsHDhQkaPHl2jfP7iF7/g73//eyAg7dy5M2RSqQ0bNnDaaafV6Ng1oYFDNSx2pzXRlCvCTIb5B+Hh0+Bfk2HZn+DrD3VWwybiktMvYebZM+mY3hFB6JjekZlnz6zVXVX+IdU//PBDevbsSVZWFtOnT6dDhw4AtG/fnv79+/PjH/84sM/YsWPZtm1bSOd4OMOHD+fGG29kxIgRjBw5kptvvrnG08JmZmYyZcoUSkqs2Q7z8/MDQ7gPHjyYbdu2MXPmzBoduya0c1w1TJFme0ttbT2Z/t1XcGirNQw8Au36Q5czoetIa2nTU8fQagQa+rDqhYWFDBo0iHXr1lU5t0ZjosOqq8Q0fobVpxH85Lkz1Rovy9/HUXLKetDwu6+sZdsiWDfXWpfaGrqOsJYuI6wJrJLS6/86VKP14YcfctNNN3H33XcnVNCoCwkVOERkIjCxV69e8c6Kqq3qzPaW3Mwa4v30MdZnrxeO7IJ9X1kj+n73Fex6z1ondugwyFcj8QWUFl21VqIiuuCCC9i7d2+8s9EgaVOVSmyFx6w7tr770lr2rwVXobWuWceyGknXkdYT747kyo+n6lRDb6pKVNpUpVRl0lpDnwnWAuBxW30j3wXVSra9Za2zJ0Onob4ayUgroOjUuUpVoIFDNS12B3QcYi0jbrHSTh2wAsg+X1/Jl3+HL/5mrWt5WmjzVrssnf1QNXn6F6BUsw4wYJK1ALhLIHdTWfPWfz+zZj4E68HELmeUNW91yda5R1STo4FDqfIcydYDiV3PBO4AY6xbg4Obt5bPAuN7EDGzr297/63Ava2ZE4NtWlB5R7+KG7vdzqBBgzDGYLfbeeqppzj77LOjOsacOXOYMGFCVAMNzpw5k3/84x+0bdsWgIsuuoiHH344qvPGiwYOpaoiAi27WcugK6200gLYv84KJPtWw44l1tDxACktymokXUfAib3wn3vLbi3W4VNqLG/xYg7Nehx3bi6Ojh1pd/c0WkycWKtjpqamsmHDBgDef/99pk+fzqeffhqyjdvtjjhOFViBY+DAgVGPUHv33Xdzzz33RJ1nj8eD3W6Per+6ooFDqZpISrdG9e3hG0LCGDj6dWit5OsPgQh3LbqKrBqIBo5qy1u8mNwHZmCKrVEF3N9/T+4DMwBqHTz8Tp48SatWrQD45JNPeOCBB2jVqhU7duxg6dKlXHrppWzZsgWAxx57jPz8fAYOHMiaNWu47rrrSE1NZeXKlWzbto1f/OIX5Ofnk5mZyZw5c+jYsWO18vDRRx9xzz334Ha7OfPMM3n22WdJTk6me/fuTJ06lQ8++IB7772X7Oxsbr/9dg4fPozdbue1116jZ8+ePProoyxYsICSkhKmTJnC7373uzr5twmmgUOpuiACmb2tZdh1VlrRCWvo+HlXhN8n7ztY+Qx0GwkdBlvDrTRhB/70J0q2Rx5WvWjjRky5yYpMcTG5v/ktJxa8Fnaf5P796HD//ZWet6ioiKFDh1JcXExubi4ff/xxYN26devYsmULPXr0iDiC7ZVXXslTTz3FY489RnZ2Ni6Xi5///Oe89dZbtG3blvnz5/Ob3/yGF154ocK+s2bNYt48q6b6yCOPcP7553PjjTfy0Ucf0adPH66//nqeffZZpk2bBkCbNm1Yt24dACNHjuS+++5jypQpFBcX4/V6Wbp0KTk5OXz11VcYY5g0aRKfffYZ5513XqX/BtHSwKFUrKS2hF4XWA8ahhs+Rezw/nTrvTMNOp9hNW91O8saPiW1Zf3mt4ErHzSqSq+u4KaqlStXcv311wdqFSNGjKBHjx5RHW/nzp1s2bKFH/zgB4DVrBSptlG+qWrjxo306NEjMPvgDTfcwNNPPx0IHFOnTgXg1KlT7N+/nylTpgCQkmIN/rh06VKWLl0aGBMrPz+fnJwcDRxKNTqRhk+Z+CR0Pxe+XWU1b327KqjT3Tf+lj+QdDvLujU4gZ90r6pmkDNuPO7vv6+Q7ujUidNe+led5GHUqFEcOXKEw4cPA5CeXjZMjcPhwOv1Bj4XF4cfiNMYQ1ZWFitXrqyTPAULzk+kc0+fPp3bbrutzs8dTEfHVSrWBl9tBYkWXQGxXic+aaU37wQDL7fG4LrtU7jvW7j+bRh7v/Vk+5Y3YOFt8MQQ+Es/mP+/sPJp2LcWPK54X1m9anf3NCQldFh1SUmh3d3T6uwcO3bswOPx0KZNmwrr2rdvz6FDhzh69CglJSW88847gXXBQ6337duXw4cPBwKHy+Vi69at1Tp/37592bNnD19//TUAL730UtiZ/Zo1a0aXLl0CU9uWlJRQWFjIhRdeyAsvvEB+fj4A+/fv59ChQ1H8C1SP1jiUqg+Dr65eR3hyBpx+vrWANffIoe3w3Sr49kvrdfvb1jpHqtW81W0kdD3LuiU4tVXsriHO/B3gdX1Xlb+PA6xf7HPnzg17x5LT6WTGjBmMGDGCzp07069fv8C6G2+8kdtvvz3QOf76669z5513kpeXh9vtZtq0aWRlZVWZl5SUFF588UWuuuqqQOd48HzowV566SVuu+02ZsyYgdPp5LXXXmPChAls376dUaNGAdZMgfPmzaNdu3Y1+aeJSMeqUqqxOZkbGkhyN5U1b7XtVxZIuo2EVj0adPOWjlUVHzpWlVJNTfOO1pwkWVbHqPVMydqyQLJlIaydY63LaO8bLsXXT9JhMDiS4pZ1lRg0cCjV2CWlQ4/zrAWs4eUPbw/tdN++2FrnSLXmJvF3uncdkdDNWyo2GnzgEJF04BmgFPjEGPNynLOkVMNms0H7LGs58yYr7dSB0EDyxZOw/K/Wurb9ggLJSGh9eoNu3lLxF5fAISIvAJcCh4wxA4PSLwKeAOzAP40xDwOXA68bYxaLyHxAA4dS0WrWAbImWwtAaaFv9kRfX8nWoNkT09tZNZFuZ1lNXB2HhDZv6bhbTV68ahxzgKeAwM3XImIHngZ+AOwDVovI20AXYLNvM0/9ZlOpBJWUFjpkitcLh3eEdrrv8N1u6kiBTsOtznaPC1Y/D24dd6spi0vgMMZ8JiLdyyWPAL42xuwGEJFXgcuwgkgXYAOVPHciIrcCtwJ069at7jOtVCKz2aD9AGvJ/omVduqAr2nLF0i++Bt43RX3dRXB+/dbQ8y36Nrkh05pChpSH0dnIHhchn3ASOBJ4CkRuQRYHGlnY8xsYDZYt+PGMJ9KNQ3NOsCAy6wFrOatP3Ui7MCNBYfhyWEgNmjeBVqd5lu6Q8vu1mur0yC9bYPrP/EPq+52u+nfvz9z584lLS2NAwcOMG3aNFavXk3Lli1p3749jz/+OElJSSGDHUYye/Zs/vpXqx+pefPm/PWvf+Xcc8+NuP2GDRv4/vvvufjii+v0+mKhIQWOsIwxBcCP450PpZq8pDSrTyPcuFvpbWH8g9YQ8sf3wvE9kPMB5B8M3c6ZZg2d4g8qHa+AojyrD8WeBLbKhwrf9eUBVr71DfnHSshoncyoy3rSZ2SHWl1W8FhV1113Hc899xx33303U6ZM4YYbbuDVV18FrHGkDh48SNeuXas85jvvvMPf//53li9fTmZmJuvWrWPy5Ml89dVXdOgQPr8bNmxgzZo1YQNHVcO617eGkxPYDwR/I118adUmIhOBib169arLfCml/CKNu3Xhn8L3cZQWwolvrUBywhdQ/IFlz3JoOQ6O7y7b3uawAog9uSyYOJLBnsSutcdY9vJO3KXWeFH5x0pY9rI1mm5tg4ff6NGj2bRpE8uWLcPpdIY8tT1kyBCAiKPkBnvkkUd49NFHyczMBGD48OGBAQv/7//+j9WrV3PXXXdRUFBAcnIyH3zwATNmzKCoqIjly5czffp0tm/fzjfffMPu3bvp1q0b8+bN47777uOTTz6hpKSEn/3sZ4ExqepjKPVgDSlwrAZ6i0gPrIBxDfCjaA5gjFkMLM7Ozr4lBvlTSvmDQ3XvqkpKg3b9rKU8Y2DbVsjsAZ5SPn9jL0f2n7LS8fpeyxzc58ZT7vYYd6mXj1/aztbl+4GKTWCZXTMYfXWfal2a2+3mP//5DxdddBFbtmzhjDPOqNZ+4WzdurXC/tnZ2cydO5fS0lKmTp3K/PnzOfPMMzl58iRpaWn8/ve/Z82aNTz11FOANUPgtm3bWL58OampqcyePZsWLVqwevVqSkpKOOecc5gwYQI5OTn1MpR6sHjdjvsKMAbIFJF9wIPGmOdF5A7gfazbcV8wxlRvZDClVP2p7rhbVRGxmqaS0oF0SDoMDm/oNsYLGDBePJ4wHfOAx22sp+fB6mMRsV4RcJdataNKmsGCx6oaPXo0N910E88991ztry+CnTt30rFjR84880zA6v+IZNKkSaSmpgLWkOmbNm3i9ddfByAvL4+cnJx6G0o9WLzuqro2Qvq7wLs1Pa42VSnVeFVVM5h7/wryj5VUSM9olcSUn/cBT6kVKDwl1nvjBVzWbcYQsRksNTWVDetSSgzAAAAgAElEQVTXh3TaZ2VlBQromhgwYABr165l3LhxgbS1a9dWa6DDYMHDqBtj+Nvf/saFF14Yso1/uttYD6UeLKGGVTfGLDbG3NqiRYt4Z0UpVcdGXdYTR1JokeVIsjFqci/IaGc1m7U53ZrHpMNgaD8QMvtYnfHNOlpzwYsNXAVWp33ed9Z0v8YLuRvh4FY4kgMnvmXcmQMoKSpg9jN/s55dMYZNmzbx+eefVyuv9957L7/+9a85evQoYHV8z5kzh5/+9Kf07duX3NxcVq9eDViTMrnd7pCh2cO58MILefbZZ3G5rOH0d+3aRUFBQb0NpR6sIfVxKKVURP4O8GrdVSViPU9id/qawsoxxqqVeEqtbTPaldVWivMQr5uFf3+IaQ8+xiN//jMpyUl079aFx/84A07msnPnTrp07uzrVhFmzZrFVVddFTj8pEmT2L9/P2effTYiQrNmzZg3b15gJsD58+fz85//nKKiIlJTU/nwww8ZO3YsDz/8MEOHDmX69OkVsnzzzTezZ88ehg8fjjGGtm3bsmjRonobSj2YDquulIqbBjusutdTsenLXVIWbEy5vhibI7QJLOR9UoN7dkWHVQ+ifRxKqTphs4Mt1brVuDxjrCfoQ4JJiRVkSgvAc7zcDv7aT8VbjLEnW+dqYIGlKgkVOPR2XKVUzFXZDOa1+kWCayj+98V5FYdtEVtQUPEFlOD3tmp0RRceg1O51jnsSVafTlrrurneMBIqcCilGh9jDNLIfnFXSmxWjcKRHH59oBmsXGBxl0DxKaB8M5gzNJgEN4fZnVB03Oro9zefeUrLnu4PEzzqontCA4dSKm5SUlI4evQobdq0SazgUZnqNIOFq61Eagazdix3HK9VAykXOIwxHD16lJSUlFpdQkIFDu3jUKpx6dKlC/v27ePw4cPxzkojYAOTbAUWr8f36oaSk5F3OVaxmSslJYUuXbrUKid6V5VSSjVmswaGH3iyRVe4u/IRfINFc1dVQj0AqJRSTc74GRWbvZypVnqMaOBQSqnGbPDVMPFJq4aBWK8Tn4zpjIwJ1cehlFJNUl0NPFlNCVXjEJGJIjI7Ly8v3llRSqmElVCBQwc5VEqp2EuowKGUUir2NHAopZSKSqWd4yJSyZMl1iZArjGmenMzKqWUavSquqvqG2PMsMo2EJH1dZgfpZRSDVxVTVVXVOMY1dmmXuhdVUopFXuVBg5jzG4AEUkXEZvvfR8RmSQizuBtGgK9q0oppWKvup3jnwEpItIZWAr8LzAnVplSSinVcFU3cIgxphC4HHjGGHMVkBW7bCmllGqoqh04RGQUcB2wxJdmj02WlFJKNWTVDRx3AdOBhcaYrSJyOrAsdtlSSinVUFVrkENjzGdY/Rz+z7uBO2OVKaWUUg1XpTUOEZlZ1QGqs41SSqnEUVWN4+Yqnh4X4BpgZp3lqBZ06lillIq9qvo4/gE0q2TJ8G3TIOhzHEopFXuV1jiMMb+rr4wopZRqHHR0XKWUUlHRwKGUUioqGjiUUkpFpVqBwzew4UcissX3ebCI/Da2WVNKKdUQVbfG8Q+sJ8ddAMaYTVi34SqllGpiqhs40owxX5VLc9d1ZpRSSjV81Q0cR0SkJ2AARORKIDdmuVJKKdVgVWusKuBnwGygn4jsB/4L/E/McqWUUqrBqu4gh7uBC0QkHbAZY07FNltKKaUaqmoFDhFpCVwPdAccIgKAMaZBjZCrY1UppVTsVbeP412soLEZWBu0NCg6VpVSSsVedfs4Uowxv4hpTpRSSjUK1a1xvCQit4hIRxFp7V9imjOllFINUnUDRynwKLCSsmaqNbHKVDzkLV5MzrjxbO8/gJxx48lbvDjeWaqxRLmWRLkO0GtRsVXf30l1m6p+CfQyxhyJZWbiJW/xYnIfmIEpLgbA/f335D4wA4AWEyfGM2tRS5RrSZTrAL0WFVvx+E7EGFP1RiJLgcnGmMKY5KKOZWdnmzVrql8hyhk3Hvf331dc4XCQ3Ls3CAgCNhuI+BZfmgSnB6X5F5tg3YVWLl1AxFYuzZ8uILaKaSHnI+xxT773H0xhUYVLkdRUml90ERhjLRiMMdYjnf40f7rXVEwLt60xGPzvwxwnTLox3kqOU5ZevHMnuFwVvxOnk5S+fX0XVe67CEkrWwdhvhcJrKj47xg4jvUa2LfC8f1JYb7foGOf+nhZ4I+6/HfSbNy4itdY4W+y4t9o2L/bCruF28ZUvU0l5ytYviL8taSkkH7uOWGOFR3/HZt1p66PR9n320COl//pp2G/E0enTvT++KMosiFrjTHZ1dm2ujWOAmCDiCwDSvyJDe123Jpy50Z4CN7txtmhQ2jhFqlQLZceXEAarzds4VhZQRq2AA8+hu+Y5Y8bLmgAmKIiClatCh/cggvf8sExbNCLECSrShcbYrNX4/iEDxpY6fY2rSsGKwgfsCD8v3dQugn+d4QI3wcVv7/gYweft/w5w/xR+7+Toi2by4JesPKFSbjCpTppYTcpn1j9Y0e8luJiXN/tC7uu2qrxIzaux7MOWrdHq4M8RvpOIpZrdaC6gWORb0lIjo4dw9Y4HJ060fXZZ+KQo5qLVHuK9tdHvFV2Hd3+/vc45KjmKruWXu+/H4cc1Vxl13L6WwlbRDRoEb+Tjh1jds5qdY4bY+aGW2KWq3rW7u5pSEpKSJqkpNDu7mlxylHNJcq1JMp1gF6Liq14fCeV1jhEZIEx5moR2UyYOpoxZnDMclaP/B1Ih2Y9jjs3F0fHjrS7e1qj7OxLlGtJlOsAvRYVW/H4TirtHBeRjsaYXBE5Ldx6Y8zemOWsFqLtHFdKqaYums7xSpuqjDH+3pWfGmP2Bi/AT2ubUaWUUo1PdR8A/EGYtB/WZUaUUko1DlX1cfw/rJrF6SKyKWhVM2BFLDOmlFKqYarqdtx/A/8BHgLuC0o/ZYw5FrNcBRGR04HfAC2MMVfWxzmVUkpFVlUfR54xZo8x5tpyfRzVChoi8oKIHBKRLeXSLxKRnSLytYjcF2l/Xx52G2Nuqs75lFJKxV51HwCsqTnAU8C//AkiYgeexuo32QesFpG3ATtWzSbYT4wxh2KcR6WUUlGIaeAwxnwmIt3LJY8AvvZNR4uIvApcZox5CLg0lvlRSilVe9W9q6oudQa+C/q8z5cWloi0EZHngGEiMr2S7W4VkTUisubw4cN1l1ullGrgFq3fzzkPf0yP+5ZwzsMfs2j9/pieL9ZNVbVmjDkK3F6N7WYDs8F6ADDW+VKqqVq0fj+Pvr+T708U0allKr+6sC+Th0X87adibNH6/Ux/czNFLg8A+08UMf3NzQAx+17iETj2A12DPnfxpak6kih/2IlyHZA41xKPQqqmjDF4DXi8Bq+xFo/X4PWCx/fZ6zW+91jvfZ+NMXi8YfYNPmaYfcu2Jeg41mKMlRbYzmvw+Pb1Gt+xvNa+wee0jkPQcYzvONa+izd9H/g+/IpcHh59f2dCBY7VQG8R6YEVMK4BflQXBxaRicDEXr16Rb2v/mE3LI3xOozvD9xfuHh9w/m8vWE/M97eSrHLC1jXct+bmzhZVMqErI6BAiO4UCorHMoKDX9hElxweU1oIRjV/kGFlwkUaIQWqOXO8fKqvWELqfsXbmb510dCC8RAQVjxvGXnr1iIe7xl+SlfaAYXzKHBoFwh7tuvMRMBuwg2m2ALem+3CTaxFrsNCks9Yff//kT4KRbqJG91MR58xIOLvAKMATKBg8CDxpjnReRi4HGsO6leMMb8sS7PG+1YVeULKYBUp52HLh8Us0LKGIPba3B7DKUeL26PF5fH4PJ4cXm8uL2GUrf16k9zeYxvu7Jty+//t49zOFnsrnC+jGQH153VjbIpJUzolBaU/aEZYyqkl01vYcqmmqDiMUKOH0gv+0zwfpUcf/nXRyhxeytcR7LDxvBurfAGjl9W2BlT8XPIKxHSy716A8fx7eOteAz/NsH7JCq7r+CyiYT9Tvw6tUjxFXISso+9fJpNrEJQBJuNcgVhWYFo8xeUErqf3Tc5mt1mFabiSys7Tvljlj+nYPdtE9hXxCqowxTMIXkrn58w12kL5KfiviEBoNJ/g+pPanXOwx+zP0yQ6NwylRX3hZkoLIJYTORUI8aYayOkvwu8G8tzR+PR93eG/RX1wFtb2HXwVIWC2uXx4vIaXG4vbq+X0ggFurWdF5fbWNuFBIL6LWnyS9y8uGKPf1I7BAma5M43lVDw53LrxLdB0MR4FY6BP706xy93DMqlRyqgStxePF4Dgu+P14bNn18p+yP2nyfwuapXCPwRS/DnwHH92/r3C/M5aB+bzboomwiPvLcj4vfy0OWDAsfwFxr+Aqx8QRIolCMUXFXuH6Eg9l93yP5BhWywuiqkVN351YV9w/7w/dWFfWN2zgbfOR6NmjZVRarSnSp284/Pd+Ow2XDaBafdhtNuw2EXknyv1mcbSXbBYbORmmTDaZPAdtY+4tvGhsMmOB3B21R+bP86h81GkkN8eQlKD2xvpV046zO+z6s4I1hj+8OurIBacPuoOOSo5uat2hvxWq4d0S0OOaq5eBRSqnL+VpH6bGpPqMBhjFkMLM7Ozr4lmv06tUyN8Iedwor7xtdV9urFvRf1S4g/7EQqoBLpWuJRSKmqTR7WuV6/g4QKHDUV+Q+7XxxzVTOJ8oedKNcBiXUtUP+FlGp4Yto5Hi81mcgpUe6qUkqpmmgwneP1rTa34+qvKKWUqp54DDkSM8aYxcaYW1u0aBHvrCilVMJKqMChlFIq9jRwKKWUikpCBQ4RmSgis/Py8uKdFaWUSlgJFTi0j0MppWIvoQKHUkqp2NPAoZRSKioaOJRSSkUloQKHdo4rpVTsJVTg0M5xpZSKvYQKHEoppWJPA4dSSqmoaOBQSikVFQ0cSimlopJQgUPvqlJKqdhLqMChd1UppVTsJVTgUEopFXsaOJRSSkVFA4dSSqmoaOBQSikVFQ0cSimloqKBQymlVFQSKnDocxxKKRV7CRU49DkOpZSKvYQKHEoppWJPA4dSSqmoaOBQSikVFQ0cSimloqKBQymlVFQ0cCillIqKBg6llFJR0cChlFIqKho4lFJKRSWhAocOOaKUUrGXUIFDhxxRSqnYS6jAoZRSKvY0cCillIqKBg6llFJR0cDhs2T3Eia8PoHBcwcz4fUJLNm9JN5ZUkqpBskR7ww0BEt2L2HmFzMp9hQDkFuQy8wvZgJwyemXxDFnNbNk9xKeWPcEBwoO0CG9A3cNv6tRXodSqmESY0y881DnsrOzzZo1a6q9/YTXJ5BbkFsh3YaNzLRMnDYnTpuTJHtS2PdOuzPkfZItKfS9b73D5oh4jOD3EbfzHccmkSuK5YMgQIo9hZlnz9TgoZSKSETWGmOyq7Ot1jiAAwUHwqZ78XJu53Mp9ZTi8rpweVyUesveF7mLAu9dXt8633t/utu46zy/DnHgtPsCTFBgctqcfHvqW9ze0HMWe4p5Yt0TGjiUUnVCAwfQIb1D2BpHx/SO/O7s39Xq2F7jDQ0u/iAUFFz874PXhWxXLhgFB6jywWp33u6w+cgtyOXXn/2a4e2GM6z9MHq17FVpzUUppSLRwAHcNfyusM07dw2/q9bHtomNZHsyyfbkWh+rOiI1u6XYU/jqwFe8+993AWjmbMaQdkMY1m4Yw9oNY1DmIFIcKfWSR6VU46aBg7IO8EToUI4UBGeePZOLe1zMvvx9rD+03loOrmf5/uUAOGwOBrQeYAWS9lYwaZ3SOl6XoZRqwLRzPAFFc1dVXkkeGw5tYN2hdWw4tIHNRzbj8roA6N68e6BGMqzdME5rfhoiUp+XopSqJ9F0jmvgUCFKPCVsO7otUCNZf3g9eSXWoJGtU1oztO1QhrcfztB2QxnQegBOuzPOOVZK1QW9q0rVWLI9OVDDYKDVub8nbw/rDq0LNHF9/N3HgW0HZQ4KbD+k3RCaJzWP8xUopWJNaxwqakeKjrD+0HrWHbSat7Yf247HeBCEXq16MbydVSMZ3m44HdM7avOWUo1AQjVVichk4BKgOfC8MWZpVfto4Khfha5CNh/ZHKiRbDy8kQJXAQDt09oHaiTD2w+nd8ve2G32OOdYKVVegwkcIvICcClwyBgzMCj9IuAJwA780xjzcDWO1Qp4zBhzU1XbauCIL4/XQ86JnECNZO2htRwqPARAujOdIW1DbwNOc6bFOcdKqYYUOM4D8oF/+QOHiNiBXcAPgH3AauBarCDyULlD/MQYc8i331+Al40x66o6rwaOhsUYQ25BbqBGsu7QOr4+/jUGg13s9G/d32raaj+cYe2GkZmaCeiYW0pVV138rTSYwOHLTHfgnaDAMQqYaYy50Pd5OoAxpnzQ8O8vwMPAB8aYDys5z63ArQDdunU7Y+/evXV4FaqunSw9ycZDGwPBZPORzZR4SgDo2qwrbVPbhtwaDDrmllLhLNm9hLmLFjF8z4VklLYiP+k467q/zw2TJ0f1t9LQA8eVwEXGmJt9n/8XGGmMuSPC/ncCN2DVTDYYY56r6pxa42h8XB4X249tDwSSZd8uw4u3wnYp9hQu63UZrVNahy6prWmT0obmSc21M141KsYYitxFFLoLKXQVhrwWuAoC74vcRRSWFlJYWkRhcRFFJSWUlJZwfKuH7G8vxmHKbo132UrZ2P8/vPjzJ6qdj4S6HdcY8yTwZLzzoWLLaXcyuO1gBrcdzA1ZNzB47uCw2xV7ilm6ZyknSk5gqPijxyGOQCCpEFxSWtMmtU3IZx1mRVWXMQa3x0NhURGnigsoKCqgoKSIgqLCQEFeXFJMcWkpJaUuSkpKKS11Uepy4S714HZ5cLu8eNxePG6DcRuMW8Aj2LwO7F4HDuPA7nVi9zqwB713eFtg97bBbhw0w0azKvLq9CbRN2d0zP4t4hE49gNdgz538aXVmohMBCb26tWrLg6n4qiygSeXXrkUt9fNiZITHCs+xtGioxwrPha6FFmve0/u5VjxMYrcRWHPk+ZICwk0bVLahA0wrVNa0zK5ZY3vCEuk/ppdXx5g5VvfkH+shIzWyYy6rCd9RnaI+Xm9Hm9Zwet7dbt874M/+9a7XB5KSkp8BXoJxSWllJSW4nK5KS114XKFFuZel8HrNhg3GA/gsSG+xeZ1YPPasVU6950DyAh8sgOpviXkOuwesHnBYcBusDlAHII9WbA7bTgcNhxOO84kB06ng+QkJ8lJSSQnJ5GU5MDutGN32HA4bdidNj5+aTtCxVp2Rmmr2v+jV3Kl9W010FtEemAFjGuAH9XFgY0xi4HF2dnZt9TF8VT8VDXwpMPmIDM10+pIr8bfR6GrkOMlxwMB5VjxMY4WHw0JOt/nf8+WI1s4Xnwcj/FUOIYgtEppFbYmUz7wtE5pTbozHREJtEGP23NroA167reLYHLjmSjMGIPH7WXnqgN8On+H9UsZyD9Wwgf/2srR3AI69WqJx+XF7fb4CnJjfXZ5AgW7x2Vwu714fGlu33ZulweXy42r1I3b7QkEBK+7rDDH1L4J0i0uPDYXHpsbj7jx2Fy4bW68NjfG7kXsgNMgqWB3CDaHYHcKdqcdp9OOzWknyekkKdlJUpKT5OQkUpKSSU1OJjU5hbSUVNJTUklLSSMlKQlHUlAh77Bhc0idN6V+9tY2PCcrHtPRPIY3PsX4rqpXgDFAJnAQeNAY87yIXAw8jhWUXzDG/LEuz6t9HIkhXr/SvcbLyZKTgeASriYTWFd0jFOuU2GPk2RLonVqa1rs7ca531yF05sUWOeylbKq90LOOq8/NrFhExt2sQfeR1rsYkdEsHntvl/D1itum+9VEI8d4xbEY7MKXLeVbtxgfK94BK8bjAvr1Q0el/F9NnhcBq8bq/B2GTxuQ5iWwagYmwevzYPX7sEjbtw2F25x4ZZSSqUEr82D2+bCI76C3ebGI1bBHlzYG7sXu1NwOO04nDaSnE6cTgfOZAfJTuuXeWpSMikpyaQlp5CakkpacgppKWlkJKWT5kizFqf1mu5Mb9RD5+z68gAfvrQ1EMwBxGG44H+zoqoJNqjO8XjQwJEY4tUkEg3jNRSVFHMk/xhH849xrOAEJwpOcKLgJCeLTnGyKJ+0T3qT5q7YKl1iL2JHh5XYvHbsHqs92+F14vAmBb13hn1vo+YPUbql1CqgbS6r8I7wPlzaiO8uCdssYjAsHDjLV+C7AoFB7IbkpCRSkpOtgtpXWKc6UwOFdnAhnuZMI9WRaqX703zp/m0bcyEfK3Xxt9JkA0dQH8ctOTk58c6OqoVdXx5g2cs7cJeW3VnlSLIx9rp+Ff4gvF4TaOe2mj48uEsrtoFb78vatN0h6b73bi+eUk9IW7k73Da+Jhivu3Z/Pza7WM0YSfZAm7XDabV12502bE6rqcTmazaxOcHmFGyOsldxWr8wbQ4Bh9VmjsMgDsDhBbsBuzfw3ovXanoyHgwGj9eDFy9e47XS/OvKvXqNl91P22hWWnG4/VNJx7j4gT5lwcBX6Gsh33g02cDh19RrHPH4pV6+8A5u13aXVixwgwvs0MLdWvf1ukMhQcNPbEJ6i6SQgt/rqWXhbfN1SiZZ7dBW4W29dyTZrUI8qDMy5DWwvT1knbWv9fmt59biLajYqWpv7uX2P19Qq7zXtx//7S6GbP9hhWa3aG/9VA1PQt2OW18aQ7NIZYzX4PF42bXqAJ8tyMHjsgrd/GMlfPzSDo4fLKBT71bhf3kH/zov9VbsvCz/q9tdsfCvdeHt/+XtW8IFDf91dunXCrvTHqEA9xfe9sAdKvakcNv4A4Jgs8d2Ct3xVw8M2wY97oqBlezVMF15yQTmut6o+LDZJZPjnTVVjxIqcNT0dtzyzSL5x0pY9vIOgEDwMF6D12PdWeLxWHeHeD1l92R7Pb5C1WPwusvSPW6vbzvfe3fZMbz+ff3be4K28ZQVyB7feq/H96veE7qf123weiMX3B63lzXv7gUqf5reZpdyv6jtIYVtSpojULD7C+fIBbMtpHCv+Evdbv3C999tYgttN597/wryj5VUyGNG62TG3zAgmq837vz/hxrzDxO/S06/BCYnxmyZqua0qYrIhRQCjiQ7Xnftf1GHY3MIdnvZbXp2u69d2y5WIewQ3zobdnvQe4f40oLe+9JXLdod8XxTfjm8XHNMaMFevvCOp2j6OJRStadNVVEKGzQADGSN7mQVtL5f4zZ7UIFut1n3eFda+NsCv+Ttdt82vrRYDI2x5bP9EX+pd+rdss7PFyuJ9CtdqUSjgQOrUI1U2J57Ze845KjmRl3WM+wv9VGX9Yxjrmqmz8gOGiiUaoBi2yvYSIy6rCeOpNB/isZc2I69rh8ZrZMBK/hp845Sqi4lVI2jpp3jidYsor/UlVKxpJ3jSimlouoc16YqpZRSUdHAoZRSKioaOJRSSkUloQKHiEwUkdl5eXnxzopSSiWshAocxpjFxphbW7RoEe+sKKVUwkrIu6pE5DBVDcwUWSZwpA6zE0+Jci2Jch2g19JQJcq11OY6TjPGtK3OhgkZOGpDRNZU95a0hi5RriVRrgP0WhqqRLmW+rqOhGqqUkopFXsaOJRSSkVFA0dFs+OdgTqUKNeSKNcBei0NVaJcS71ch/ZxKKWUiorWOJRSSkVFA4dSSqmoaOAIIiJ2EVkvIu/EOy+1ISJ7RGSziGwQkUY9TLCItBSR10Vkh4hsF5FR8c5TTYhIX9/34V9Oisi0eOerJkTkbhHZKiJbROQVEUmJd55qSkTu8l3H1sb2fYjICyJySES2BKW1FpEPRCTH99oqFufWwBHqLmB7vDNRR8YaY4YmwL3pTwDvGWP6AUNopN+PMWan7/sYCpwBFAIL45ytqIlIZ+BOINsYMxCwA9fEN1c1IyIDgVuAEVj/ty4Vkegm84mvOcBF5dLuAz4yxvQGPvJ9rnMaOHxEpAtwCfDPeOdFWUSkBXAe8DyAMabUGHMivrmqE+OBb4wxNR3dIN4cQKqIOIA04Ps456em+gNfGmMKjTFu4FPg8jjnqdqMMZ8Bx8olXwbM9b2fC0yOxbk1cJR5HLgX8Fa1YSNggKUislZEbo13ZmqhB3AYeNHXhPhPEUmPd6bqwDXAK/HORE0YY/YDjwHfArlAnjFmaXxzVWNbgNEi0kZE0oCLga5xzlNttTfG5PreHwDax+IkGjgAEbkUOGSMWRvvvNSRc40xw4EfAj8TkfPinaEacgDDgWeNMcOAAmJU9a4vIpIETAJei3deasLXZn4ZVlDvBKSLyP/EN1c1Y4zZDjwCLAXeAzYAnrhmqg4Z61mLmDxvoYHDcg4wSUT2AK8C40RkXnyzVHO+X4UYYw5htaOPiG+OamwfsM8Y86Xv8+tYgaQx+yGwzhhzMN4ZqaELgP8aYw4bY1zAm8DZcc5TjRljnjfGnGGMOQ84DuyKd55q6aCIdATwvR6KxUk0cADGmOnGmC7GmO5YzQgfG2Ma5a8oEUkXkWb+98AErCp5o2OMOQB8JyJ9fUnjgW1xzFJduJZG2kzl8y1wloikiYhgfSeN8oYFABFp53vthtW/8e/45qjW3gZu8L2/AXgrFidxxOKgKq7aAwutv2kcwL+NMe/FN0u18nPgZV8Tz27gx3HOT435AvkPgNvinZeaMsZ8KSKvA+sAN7Cexj1cxxsi0gZwAT9rTDdfiMgrwBggU0T2AQ8CDwMLROQmrKklro7JuXXIEaWUUtHQpiqllFJR0cChlFIqKho4lFJKRUUDh1JKqaho4FBKKRUVDRyqyRCRDiLyqoh84xuO5V0R6VPFPtN8w1FUdeyrfKP3LvN9fkVENonI3eW2myki+30j5G4TkWtrcB35vtdOvltjlapXejuuahJ8D6t9Acw1xjznSxsCNDfGfF7JfnuwRoI9UsXx3wP+YIxZLiIdgOXGmAojrYrITCDfGPOYiPQG1gJtfE9hV/da8o0xGdXdXqm6pjUO1ZmrsjsAAAJmSURBVFSMBVz+oAFgjNlojPlcRMYEz8EiIk+JyI0icifWeEzLgmoS1/rmOtkiIo/40mYA5wLPi8ijWGMfdfbVKkZHypAxJgdrePVWvuPcIiKrRWSjiLzhr+mISA8RWek77x+C8tndPxeDL79PBa17x3dddhGZ48vv5vI1IKVqQgOHaioGYv26rzZjzJNYQ4aPNcaMFZFOWIPijQOGAmeKyGRjzO+BNcB1xphfYQ1i+I1v/o3KajPDgRzfmGIAbxpjzjTG+OcducmX/gTWQI+DsEakjcZQoLMxZqBv/xej3F+pCjRwKFV9ZwKf+Ab4cwMvY80XEq27RWQr8CXwx6D0gSLyuYhsBq4Dsnzp51A2vtVLUZ5rN3C6iPxNRC4CTtYgv0qF0MChmoqtWDPvheMm9G8h1lOhzjLGZAFXYDVv+c83B7jDVzP4Xbl8VNUZGfYajDHHsWa3+wS4HZ2oTNUBDRyqqfgYSA6e2EpEBvv6IPYCA0QkWURaYo346ncKaOZ7/xVwvohkiogda6TbT2uaIWPM21hNXP7RTJsBuSLixKpx+K2gbHrW4PRge4ChImITka74htIXkUzAZox5A/gtjX9YetUAaOBQTYJvUpspwAW+23G3Ag8BB4wx3wELsIafX4A14qvfbOA9EVnmm1ntPmAZsBFYa4yp7bDVvwd+ISI24AGs5qsVwI6gbe7CmpBrM9A5wnFWAP/FGnb+SazRa/Ft/4mIbADmAdNrmV+l9HZcpZRS0dEah1JKqaho4FBKKRUVDRxKKaWiooFDKaVUVDRwKKWUiooGDqWUUlHRwKGUUioq/x/33olzYrD+cwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "label = {'kd':\"KDTree\",\n",
    "        'cg':\"Cellgrid\",\n",
    "        'ns': \"Cython NS\",\n",
    "        'bf': \"Brute Force\",\n",
    "        'oct': \"PCL Octree\"}\n",
    "result_cut = defaultdict(list)\n",
    "for key in timing:\n",
    "    method = key.split('.')[0]\n",
    "    result_cut[method].append(timing[key])\n",
    "    \n",
    "for key in ['kd', 'cg', 'ns', 'bf', 'oct']:\n",
    "    ax.plot(maxdist, result_cut[key], '-o', label = label[key])\n",
    "ax.set(xlabel='Cutoff Radius', ylabel='time [s]', yscale='log', title='N = 40000')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Except cellgrid, other methods dont scale with the cutoff radius to a large extent. However, neighbour search and Octree also scales with cutoff radius but with a very small slope. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the script which results in dead kernel. I suspected it might be due to consumption of memory, but it was not the case. \n",
    "Even the above script dies if I try to run for multiple cutoff distances, and/or particles. It definitely is not due to the number of particles, because I tried with less number of runs but large number of particles. I suspect it can be due to two reasons (1) Either handling of memory on my system is not so efficient (looks like there is a lot of clutter) on my system. (2) The more likely one, the kernel might be storing the data which I can't figure out to remove after executing all the methods for one setup (one instance of number of particles)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another weird feature is the peak memory consumption. I am not sure if that is the culprit, but it could be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The memory_profiler extension is already loaded. To reload it, use:\n",
      "  %reload_ext memory_profiler\n"
     ]
    }
   ],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Methods2 = dict(KDTree=kdtree_select, NS=ns_select, OCT=oct_select, BF=bf_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 110.70 MiB, increment: 0.03 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit Methods2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "box = np.array([100,100,100,90,90,90],dtype=np.float32)\n",
    "Nprotein = 100\n",
    "Nsolvent = 1000\n",
    "center = box[:3]/2\n",
    "radius = 10.0\n",
    "maxdist = 6.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "    'Nprotein':Nprotein,\n",
    "    'Nsolvent':Nsolvent,\n",
    "    'maxdist':maxdist,\n",
    "    'centre':center,\n",
    "    'radius':radius\n",
    "         }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bbed4f93bc9d4c10b32a2882347442ed",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4491e5a135414ebdb4bb895af98d6292",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "timing = %memit -q -o timer(**params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<MemitResult : peak memory: 2862.06 MiB, increment: 1379.27 MiB>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "timing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Such a high amount of memory is consumed when the number of particles are very less. I tried with different number of particles for different methods, this looks like a background memory consumption, but according to the docs it should be only depend on the memory usage by the function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Everything below this is just for reference and to identify the problems I faced in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "Methods2 = dict(KDTree=kdtree_select, NS=ns_select, OCT=oct_select, BF=bf_select, CG=cg_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm_notebook\n",
    "\n",
    "results, struct = defaultdict(list), defaultdict(list)\n",
    "Nprotein = 100\n",
    "Nsolvent = np.logspace(3 ,5 ,num = 6, dtype = np.int)\n",
    "protein = generate_protein(center,radius,Nprotein)\n",
    "\n",
    "for N in tqdm_notebook(np.unique(Nsolvent)):\n",
    "    #Initialize\n",
    "    struct = defaultdict(list)\n",
    "    solvent = generate_solvent(center, radius, N, box)\n",
    "    print(\"Number of Points:\", N)\n",
    "    \n",
    "    \n",
    "    for key, functions in Methods.items():\n",
    "        struct[str(key)] = functions(protein, solvent, maxdist, box)\n",
    "        res = %timeit -q -o functions(protein, solvent, maxdist, box)\n",
    "        results[str(key)].append(res.average)\n",
    "        del res\n",
    "    results['N'].append(N)\n",
    "    # Testing\n",
    "    for a, b in itertools.combinations([struct[str(key)] for key in Methods], 2):\n",
    "        if (a is not None) and (b is not None):\n",
    "            np.testing.assert_equal(a.sort(axis = 0),b.sort(axis = 0))\n",
    "        else:\n",
    "            np.testing.assert_equal(a,b)\n",
    "    del struct\n",
    "    print(\"ALL TESTS PASSED\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for key in Methods:\n",
    "    ax.plot(results['N'], results[key], '-o', label = key)\n",
    "ax.set(xlabel='N particles', ylabel='time [s]', xscale='log', yscale='log', title='Total Time')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For selection purposes, brute force is better at lower particle density upto around 1000 particles. Total time for non-periodic Octree also seems promising, however, its behaviour with PBC  is unknown. Neighbour search in FATSlim scales pretty well with number of particles, where it has a significant advantage for large number of particles. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the trend can be easily extended from the single query benchmarks. But for the sake of completeness, lets check if there is any unexpected variation for selection with different cutoff radius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Methods2 = dict(KDTree=kdtree_select, NS=ns_select, OCT=oct_select, BF=bf_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "box = np.array([100,100,100,90,90,90],dtype=np.float32)\n",
    "center = box[:3]/2\n",
    "radius = 10.0\n",
    "maxdist = np.linspace(4, 10, num = 3, dtype = np.float32)\n",
    "Nprotein = 100\n",
    "Nsolvent = np.array([1000])\n",
    "protein = generate_protein(center,radius,Nprotein)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def timer(method, protein, solvent, cut, box):\n",
    "    res = %timeit -q -o method(protein, solvent, cut, box)\n",
    "    return res.average\n",
    "def memory(method, protein, solvent, cut, box):\n",
    "    res = %memit -q -o method(protein, solvent, cut, box)\n",
    "    return res.mem_usage[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "419b54c5c466424fbf9c28c83360d821",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Points: 1000\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2aaf5ee6ebb045098e3bf1459c82e93d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=3), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KDTree\n",
      "NS\n",
      "OCT\n",
      "BF\n",
      "KDTree\n",
      "NS\n",
      "OCT\n",
      "BF\n",
      "KDTree\n",
      "NS\n",
      "OCT\n",
      "BF\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm_notebook\n",
    "results2, struct2 = defaultdict(list), defaultdict(list)\n",
    "for N in tqdm_notebook(np.unique(Nsolvent)):\n",
    "    #Initialize\n",
    "    solvent = generate_solvent(center, radius, N, box)\n",
    "    print(\"Number of Points:\", N)\n",
    "    results2['N'].append(N)\n",
    "    for cut in tqdm_notebook(np.unique(maxdist)):\n",
    "        \n",
    "        for key in Methods2:\n",
    "            print(key)\n",
    "            results2[str(key)+'.'+str(N)].append(timer(Methods2[key], protein, solvent, cut, box))\n",
    "            results2['mem'+str(key)+'.'+str(N)].append(memory(Methods2[key], protein, solvent, cut, box))\n",
    "            del Methods2[key]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<MemitResult : peak memory: 5441.15 MiB, increment: 0.00 MiB>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%memit -q -o Methods2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(list,\n",
       "            {'N': [1000],\n",
       "             'KDTree.1000': [0.004731408571429451,\n",
       "              0.005371267142857375,\n",
       "              0.007642431428570749],\n",
       "             'memKDTree.1000': [113.30859375, 5314.88671875, 5291.15234375],\n",
       "             'NS.1000': [0.013609149999999447,\n",
       "              0.005241251428573637,\n",
       "              0.0028423099999989585],\n",
       "             'memNS.1000': [5313.96484375, 5290.31640625, 5440.2578125],\n",
       "             'OCT.1000': [0.0011210305714287512,\n",
       "              0.0015515371428578614,\n",
       "              0.002058074285714611],\n",
       "             'memOCT.1000': [5314.0859375, 5290.31640625, 5440.26171875],\n",
       "             'BF.1000': [0.0036874957142845233,\n",
       "              0.0042731114285730395,\n",
       "              0.0038277171428580913],\n",
       "             'memBF.1000': [5314.88671875, 5291.15234375, 5441.078125]})"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "index 1 is out of bounds for axis 0 with size 1",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-2c6414581c5a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0maxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m     \u001b[0mN\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNsolvent\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mkey\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mMethods2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m         \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmaxdist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresults2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'.'\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'-o'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1 is out of bounds for axis 0 with size 1"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8VNX5+PHPk40EyMIOCSBLQEQEhIgKVlEQUUQobtD251qpdcNdtGpR26rVby2taIvVuoOIiGBQrCguICrKvkPYAmGXQEJCtuf3x53EIZlJJpkty/PmNa+ZOffcc5+bhHnm3HPvuaKqGGOMMTUVEe4AjDHG1G2WSIwxxvjFEokxxhi/WCIxxhjjF0skxhhj/GKJxBhjjF8skdRSIjJJRN4MdxzGGFMVSyTVICIqIqtEJMKt7E8i8mol6zwkIltFJEdEMkXknZAE6yPXPqUGod1rReQHETni2u+/ikiU2/LmIvK+iOSKyHYR+VW59X/lKs8Vkdki0jwQ6xpjAs8SSfUlA2N9qSgi1wL/Dxiqqk2BNGBBEGOrTRoDdwItgTOBIcC9bsunAAVAG+DXwIsiciqA6/nfOD+7NsAx4IUArWuMCTBLJNX3V+Ax92/XlTgDmK+qWwBUdY+qTi1dKCLJIjJHRA6JyGYRuclTIyLykYjcVq5shYiMcb3uISL/c7WzQUSucqv3qohMEZF0ETkqIt+KSFfXsi9d1Va4ekxXV+cHURlVfVFVv1LVAlXdBbwFDHJttwlwOfCIquao6tfAHJwPf3CSw1xV/VJVc4BHgDEiEu/PuoHaN2PMiSyRVN8s4AhwnQ91lwDXiMh9IpImIpHllk8HMnF6OVcAfxGRCzy0Mw0YV/pGRHoCJwHprg/W/wFvA61xeksvuOqUGgs8BjQDNgN/BlDVc13L+6hqU1WtcNhNRM4RkcOVPM7x4ecAcC6wxvW6O1Ckqhvdlq8ATnW9PtX1HlecW3B6IN39XNcYEwSWSKpPcb7lPiIiMZVWVH0TuB24CPgC2CciDwCISAecb+gPqGq+qi4H/gNc46Gp94G+InKS6/2vgVmqehy4FNimqv9V1SJVXQa8B1zpvr6qfqeqRTg9g74+76zq16qaVMnj66raEJEbcA7rPesqaoqTjN1lA/Fuy7O9LPdnXWNMEFgiqQFVnYfTk/idD3XfUtWhQBJwM/CEiFyE0ws5pKpH3apvB1I8tHEUSOfnsZlxOAkBnJ7Jme69BJxE09atiT1ur4/hfNiGhIiMBp4ELlbVA67iHCChXNUE4KgPy/1Z1xgTBJZIau4PwEM4g8pVUtVCVX0XWAn0AnYDzcsdu+8I7PLSxDRgnIicDcQCn7vKdwJflOslNFXV31d/lyoSkV+4xk+8PX5RybrDgZeAkaq6ym3RRiBKRLq5lfXh50Nfa1zvS9vpAjRyrefPusaYILBEUkOquhBYDVzrrY6IXCciI1yDxBEicjHOMfxvVXUnsBh4UkRiRaQ3cCPg7dqReTi9j8eBd1S1xFX+IdBdRP6fiES7HmeIyCk+7speoEsl+/mVKzF5e3zlZd8vwOk1Xa6q35VrMxdnrOlxEWkiIoOAUcAbripvASNdSayJa59nqepRf9b18edhjKkmSyT+eRio7BqFIzi9lh3AYZwzvn7vNq4wDuiE0zt5H/ijqn7qqSHXeMgsYCjOwHpp+VFgGM5hr904h7GexvkW7otJwGuuw2JXVVW5Gh4BEoF5br2Xj9yW3wLEAftwelu/V9U1AK7nm3GSwj6c8Y1bArSuMSbAxG5sZYwxxh/WIzHGGOMXSyTGBIGIvCIi+0RktZflIiL/cF2IulJE+oU6RmMCpdYnEhHpIiIvi8jMcMdiTDW8CgyvZPnFQDfXYzzwYghiMiYogppIvH0rE5Hhrqk8NovIxMraUNUMVb0xmHEaE2iq+iVwqJIqo4DX1bEESBKRdqGJzpjA8mW+KH+8CjwPvF5a4JomZApwIc5Ffd+LyBwgEufCNXc3qOq+6m60ZcuW2qlTJ4/LcnNzadKkSXWbDDuLO7Qqi/uHH344oKqt/NxECs41QKUyXWVZ7pVEZDxOj4UmTZr079Gjh5+bNcYzf/6ug5pIVPVLEelUrngAsFlVMwBEZDowSlWfxJnuo0bc/8O1adOGZ5991mO9nJwcmjYN2YXdAWNxh1ZlcZ9//vnbQxWHa5LPqQBpaWm6dOnSUG3aNDAiUuO/62D3SDzx9E3sTG+VRaQFziSDp4vIg66EU0H5/3CDBw/22N7ChQvxtqw2s7hDKwRx7wI6uL1vj/dZDYyp1cKRSKpFVQ/iXGBmTH0yB7jN1SM/E8hW1awq1jGmVgpHIrFvYqbeE5FpwGCgpYhkAn8EogFU9V84U95cgjOt/zHg+vBEaoz/wpFIvge6iUhnnAQyFvhV5asYU7eo6rgqlitwa4jCMSaogn367zTgG+Bkce7bfaPrnhi3AfOBdcCM0nmSjDHG1D3BPmvL47cy1/085gVz27VV9ty57Hvu7xRlZRHVrh2t77qTxJEjwx2WMcbUWK0fbK8OERkJjExNTQ13KB5lz51L1iOPovn5ABTt3k3WI48CWDIxxtRZtX6KlOpQ1bmqOj4xMTHcoXi077m/lyWRUpqfz77n/h6miIwxxn/1KpHUZlpURNHu3R6XFe3ezdHPPqMkNzfEURljjP/q1aGt2qgkP5/Ds2Zx6L+veq8kQuYttyLR0TQ+I40mvziXpuedS0znzohIyGI1xpiasEQSJMWHD3Po7bf56c23KD50iNg+vWly/mCyZ7x7wuEtiY2l7R8fJbpdO3K++JKcr75k39NPs+/pp4lOSaHJub+g6bnnQkFB2To2YG+MqU0skQRY4e7dHHz1VQ7PfA89dowm551Ly9/+lri0NESExqed5jUJNDnrLNo8cD+Fu3aR89VX5Hz5FdkfzOHwtOm0jopixzsziGzZgqPzP0GPHwdswN4YE36WSAIkf8NGDr78H46kzwMREkdcQvMbbiT25O4n1EscObLKD/zolBSajR1Ls7FjKSkoIG/pUta9/TaNMraSu3hxhfqan8/eJ5+iUffuRLVsSWRSEhIZWWXMgejZ1LfeUen+tN69m03JyXV+f4wJhXqVSCo7/TcQHxDlPzRb3Xkn0W3bcPA//yH3y6+Qxo1p/ptf0/zaa4lOTg7IPkXExNBk4EByCgpIGzyYdaf0BNUK9YoPHWLrqNGulSKIbNGcqJatiGrZ0u3RgsiWLYlq2Yq81as48M/n/ToVOVCnMwcqGfnbjvv+CNbbM8ZX9SqRqOpcYG5aWtpN7uWB+IDw+KH5wAOgSmTz5rSacAfNxo0jMikpoPtUXlS7dh7P/ops2YK2Dz9C0YEDFB08QPGBAxTtP0DRgQMc37yZogMHoLCw0rY1P5/dDz7EwZdfQWJikJhoImJikOgY13u3R6MYsmfN8ng6896nnqZR165EJCQgublocbHXHlIgk5GndlQh/oILKDl6hOKjRyk5etT1nEPx0SOUHM1xLcsh+4MPvJ6ebYnEGO/qVSLxxtv1G1mPPErO5wshMhKJiHCeIyMgouLz4XdnVGgDVSISE0n9bAERsbEh2ZfWd915wgcmOAP2bR54gIThF3ldT1Upyc52Es2BA+y4zsscgUVFRCcnowUFaEEBJcfy0IJstLCAkoICtKDw52U5nk9XLj54kK1jLnfiBdbfcy8RTZoQkZBAZHw8EfHxznNCPDmfLvD4u9nzxBMUZmaihUVoURFaWOg8FxWihYVQrjx30aKycSP3drLuv5+qptSV6GgiEhIq/n5LfyRZNimvMZVpEInE2weB5ueTv2YNWlICxcUVn4uLy957+5ApOXIkZEkEfv6WXt1DOCJCZFISkUlJNEpNJSo52WPPJio5mQ4vTPEplk0XDPHcO2rRgnaPP0bxkaNs+PFHurRp43zrP3LU6QUcOUrhvr2UbN5MybFjHtsuOXKU/ZP/4QoqComORqKinIfrNdGlr6MrJBF3re+7j4j4pkQmJBDRNJ7IhPifn+PjiWjUqNL9iWpnd8A1pjINIpF4OxwUlZxM1/kf+9RGbfqQ8WXAvireejat77rT7zbaTHyA+CFDAMhrlkSrSm4Q5f3n2pbU+fMhOtqna2m8tpOcTIsbb/BhbwLzMzGmIWoQV7a3vutOpFyvoSYfmv62UZskjhxJuyceJyo5GUSISk6m3ROPVytBBaINrz/Xu+92xmJ8vCAzEL8f9/1RqNH+GNMQNYgeifvhoMLdu4muwVlbNT2kVJsFomfjbxuB+rkGsp3EkSPr7C2CjQmHBpFIIDAfEIH44DUVBernar8fY8KjXh3aEpGRIjI1Ozs73KEYY0yDUa8SSW2fRt4YY+qjepVIjDHGhJ4lEmOMMX6xRGKMMcYvlkiMMcb4xRKJMcYYv1giMcYY4xdLJMYYY/xiicQYY4xf6lUisSvbjTEm9OpVIrEr240xJvTqVSIxxhgTepZIjDHG+MUSiTHGGL80mPuRmIpmL9vFM/M3sPtwHslJcdx30cmMPj0l3GEZY+oYSyQN1Oxlu3hw1iryCosB2HU4jwdnrQKwZGKMqRY7tNUAZR8r5PEP15YlkVJ5hcU8M39DmKIyxtRV1iNpAI4VFLF0208s2nKAb7YcZPWubErUc93dh/NCG1w9JSLDgclAJPAfVX2q3PKOwGtAkqvORFWdF/JAjQkASyR1UFVjGwVFJazIPMyizQdYvOUgy3b8RGGxEh0pnN6xGXcM6cabS7ZzIKegQtut4huFclfqJRGJBKYAFwKZwPciMkdV17pVexiYoaovikhPYB7QKeTBGhMAlkjqGM9jGyvZ+dMxGkVFsGjzQb7fdohjBcWIQK/kRG44pzMDu7bkjE7NaBzj/Mo7tWhyQjulDuUe54PluxjV18ZJ/DAA2KyqGQAiMh0YBbgnEgUSXK8Tgd0hjdCYAKpXiURERgIjU1NTwx1K0Dwzf4OHsY0S/u+TjQCktm7KFf3bM7BrS87u0oLExtEe2yntwbj3bH53XmfmrshiwvTlrMrMZuLFPYiKtGG0GkgBdrq9zwTOLFdnEvCJiNwONAGGempIRMYD4wE6duwY8ECNCYR6lUhUdS4wNy0t7aZwxxIslY1hfPvQENokxPrc1ujTUyqcoTX2jJP4c/pa/vP1VtbsPsLzvzqdFk3tcFcQjANeVdX/E5GzgTdEpJeqlrhXUtWpwFSAtLQ0LyNbxoSXfd2sY9okek4UKUlx1Uoi3sRERfDYqF7835V9+HHHT4z859eszDzsd7sNzC6gg9v79q4ydzcCMwBU9RsgFmgZkuiMCTBLJHXIzkPHKCouqVAeFx3JfRedHNBtXd6/PTNvHoiIcMW/vuHdpTurXsmU+h7oJiKdRSQGGAvMKVdnBzAEQEROwUkk+0MapTEBYomkjth5tITLX1xMYbEyYUgqKUlxCE5P5MkxpwXlIsLT2icy57ZBpJ3UjPtmruSR2aspKKqYyMyJVLUIuA2YD6zDOTtrjYg8LiKXuardA9wkIiuAacB1qmqHrkydVK/GSOqrbzMO8pdv80hsHMu7N59N9zbx3HVhYHsg3rRo2ojXbxjAX+dvYOqXGazLOsILv+5H6wAcRquNSk+t3nU4j5Qln9V42hjXNSHzypU96vZ6LTDI74CNqQWsR1LLzV+zh//3ynckNRLeu2Ug3dvEhzyGqMgIHrrkFP457nTW7D7Cpf/8mh+2Hwp5HMFWemr1LtcJDaXTxsxeVn54wxjjzhJJLTb9ux38/s0f6NkugT+cGUdKUlxY4xnZJ5n3bx1IbHQkY6cu4c0l26lPR2P++vF6mzbGmBqwQ1u1kKoy5fPNPPvJRgaf3IoXft2P7xZ/He6wAOjRNoG5t53DhHeW8fDs1cxdsYudh/LIys6v8QzCgZqF2Nd2CotL2H4wl417c9i49yibXM+7s/M9tmvTxhhTOUsktUxJifLY3DW89s12fnl6Cn+9ojfRteyiwMTG0bx87Rnc/MZS/rduX1n5rsN5THxvJccKixhzentifIg7ULMQe2pn4qyV7D2az0nNG5+QNDIO5FBY7PSkRKBj88Z0ax1PVnY+OceLKrSdHOaeoDG1nSWSWuR4UTH3zFjBhyuzuOkXnXnw4lOIiJBwh+VRZISwNutohfL8ohIemrWah2atduoJNPrsY2KiIoiJjHCeXa8bRUWwLusoBeVOac4rLObBWav4YuN+RCBShMgIISJCiBQhQih7HRkhiAhvLtle4bBUfmEJT85bX/a+fbM4ureJ5/werenepind28TTtVVT4mIigYrJCIJzarUx9Y0lkloi53gRN7/xA19vPsBDl/Rg/Lldwx1SlSo75HP/8JMpKCphc8Y22ia3p6C4hIIi53Hc7XX5JFIqr7CYpdsPUVICJaoUl2jZc3GJogrFbuWlPQxPPrh1EKmtm9KkUeV/7u7Txuw6nEeK3ezLGJ9YIqkFDuQc5/r/fs/arCP835V9uLx/+3CH5JPkpLiyM5zcpSTFcctgZ76zhQt3M3hwT69tDHrqM69tfHX/BT7HUlk7fTok+dxO6bQxCxcuZPDgwT6vZ0xDVrsOvvtJREaKyNTs7Oxwh+KznYeOccWLi9m07ygvXdO/ziQRgPsuOpm46MgTyqp7KCgQbQSyHWNM9dWrHkldmLTR/cyiVvGNyCsoIiIigrd+exb9T2oW7vCqxdMMwtU9FBSINgLZjjGm+upVIqntyg/m7jt6HICJw1PrXBIp5WkG4XC0Ech2jDHV02ASSSCmvqjJ9Q4lJcqBnOPs/CmPx+auqXBmEcAbS3Zw8+D6ew8VY0z91iASSSCuVfDWRkmJcmbXFuz6KY/Mn465nvPYddh5v/twvtczk0rZBW/GmLqsQSQSz3cVLObh2avZuNe5FqL05NHSGT+UshcAvPXtDo9t3P3uigrba9m0Ee2bxXFqSiIX9WpL+6Q4UprFMfG9VWWHs9zZBW/GmLqsQSQSb9/4c44XMfXLDMR1zZ9Q9sL9CRHn4jZvnhxzGilJcbRvFkdyUhyx5c4eKvXQJUV2wZsxpt5pEImksusdFk307VqFyq5TGDfAt3tp25lFxpj6qEEkkvsuOtnvnkAg2gA7s8gYU/80iEQSiKkvrDdhjDGeNYhEAoGZ+sJ6E8YYU1G9miLFGGNM6FkiMcYY4xdLJMYYY/xiicQYY4xfLJEYY4zxiyUSY4wxfrFEYowxxi/1KpHUxTskGmNMXVevEomqzlXV8YmJieEOxRhjGox6lUiMMcaEXoOZIsVUVFhYSGZmJvn5+UHbRmJiIuvWrQta+8GSmJjI1q1bad++PdHR0eEOx5hazRJJA5aZmUl8fDydOnVCSm/KEmBHjx4lPj4+KG0H05EjRygoKCAzM5POnTuHOxxjajU7tNWA5efn06JFi6AlkbpMRGjRokVQe2vG1BeWSBo4SyLe2c/GGN9YIjFh1bRp07LX8+bNo3v37mzfvp1JkyaRkpJC37596datG2PGjGHt2rUA/PKXv6Rv376kpqaSmJhI37596du3L4sXLw7XbhjToNkYifHZ7GW7gnZjrwULFnDHHXcwf/58TjrpJADuuusu7r33XgDeeecdLrjgAlatWsX7778PwMKFC3n22Wf58MMPPbZZVFREVJT9iRsTbNYjMT6ZvWwXD85axa7DeSiw63AeD85axexlu/xu+8svv+Smm27iww8/pGvXrh7rXH311QwbNoy333670rbat2/PxIkTOf3003n//ffZtGkTF110Ef379+fcc89l48aNAOzdu5cxY8aQlpbGgAEDWLJkid/74U5EhovIBhHZLCITvdS5SkTWisgaEal8x4ypxezrmgHgsblrWLv7iNfly3YcpqC45ISyvMJi7p+5kmnf7fC4Ts/kBO4e3LHS7R4/fpzRo0ezcOFCevToUWndfv36sX79+krrALRu3Zply5YBcP755/Of//yHrl27smjRIm677TY++eQT7rjjDu6//37OOusstm3bxqWXXsrq1aurbNsXIhIJTAEuBDKB70VkjqqudavTDXgQGKSqP4lI64Bs3JgwsERifFI+iVRV7qvo6GgGDhzIyy+/zOTJkyutq6o+tXn11VcDcPjwYZYsWcLll19etqyoqAiATz/9lA0bNpSV//TTT+Tl5REXF1fdXfBkALBZVTMARGQ6MApY61bnJmCKqv4EoKr7ArFhY8LBEokB4I8jT610+aCnPmPX4bwK5SlJcbzzu7O9rnf06NFK242IiGDGjBkMGTKEv/zlLzz00ENe6y5btoy0tLRK2wNo0qQJ4CSeli1bsnz58gp1VJXvvvuOmJiYKturgRRgp9v7TODMcnW6A4jIIiASmKSqH5dvSETGA+MBOnasvHdnTLjYGInxyX0XnUxcdOQJZXHRkdx30cl+t924cWPS09N56623ePnllz3Wee+99/jkk08YN26cz+02a9aMdu3alQ3Ol5SUsGLFCgCGDh3KlClTyup6SjZBFgV0AwYD44CXRCSpfCVVnaqqaaqa1qpVqxCHaIxvLJEYn4w+PYUnx5xGSlIcgtMTeXLMaQE7a6t58+Z8/PHH/OlPf2LOnDkAPPfcc2Wn/7755pt89tlnVPfDdPr06fzrX/+iT58+nHrqqWVneE2ZMoVFixbRu3dvevbsyUsvvRSQ/XDZBXRwe9/eVeYuE5ijqoWquhXYiJNYjKlz7NCW8dno01MCljhK5eTklL3u0KEDW7duBeCyyy5j0qRJla47ePBgBg8efEJZZmbmCe+7dOnC/PnzK6zbqlUrZs6cWbOgq/Y90E1EOuMkkLHAr8rVmY3TE/mviLTEOdSVEayAjAkm65EYE2CqWgTcBswH1gEzVHWNiDwuIpe5qs0HDorIWuBz4D5VPRieiI3xT6U9EhHxfj6oqwqQpardAxeSMXWfqs4D5pUre9TttQJ3ux7G1GlVHdraoqqnV1ZBRJYFMB5jjDF1TFWHti6vYrmvdYwxxtRTlSYStwuqmohIhOt1dxG5TESi3esYY4xpmHwdbP8SiBWRFOAT4P8BrwYrKGOMMXWHr4lEVPUYMAZ4QVWvBCq/FNoYH4gI99xzT9n7Z599tuy03w0bNjB48GD69u3LKaecwvjx48MUpTGmMj4nEhE5G/g1kO4qi6ykvqmPVs6A53rBpCTneeUMv5ts1KgRs2bN4sCBAxWW3XHHHdx1110sX76cdevWcfvtt/u9PWNM4PmaSCbgzFT6vut8+C44576bhmLlDJh7B2TvBNR5nnuH38kkKiqK8ePH89xzz1VYlpWVRfv27cven3baaX5tyxgTHD5d2a6qX+KMk5S+zwDuCFZQNSUiI4GRqamp4Q6l7vloIuxZ5X155vdQfPzEssI8+OA2+OE1z+u0PQ3O+UOVm7711lvp3bs3999//wnld911FxdccAEDBw5k2LBhXH/99SQlVZiOyhgTZpX2SERkUlUN+FInVFR1rqqOT0xMDHco9U/5JFJVeTUkJCRwzTXX8I9//OOE8uuvv55169Zx5ZVXsnDhQs466yyOH/d/e8aYwKqqR/LbKq5uF5x5hCYFLCITHhc/Vfny53q5DmuVk9gBrk+vWF6qimnkS915553069eP66+//oTy5ORkbrjhBm644QZ69erF6tWr6d+/v09tGmNCo6oxkpeA+EoeTV11TH035FGILnfTp+g4pzwAmjdvzlVXXXXCNPIff/wxhYWFAOzZs4eDBw+SkhLYSSONMf6rtEeiqo+FKhBTy/W+ynle8DhkZ0JieyeJlJYHwD333MPzzz9f9v6TTz5hwoQJxMbGAvDMM8/Qtm3bgG3PGBMYNo288V3vqwKaOODEaeTbtGnDsWPHyt7/7W9/429/+1tAt2eMCTybRt4YY4xfLJEYY4zxi0+JxDVR4wIRWe1631tEHg5uaMYYY+oCX3skL+Fc2V4IoKorcU77NcYY08D5mkgaq+p35cqKAh2MMcaYusfXRHJARLoCCiAiVwBZQYuqtgrCpIXGGFPX+Xr6763AVKCHiOwCtgK/CVpUtVHppIWFec770kkLIeCnxDY0mZmZ3Hrrraxdu5aSkhIuvfRSnnnmGWJiYvjuu++499572bt3L40bN6Z///6cfvrpvPSScx3s2rVrOfnkk4mMjGT48OE89VQVV+gbYwLOpx6Jqmao6lCgFdBDVc9R1W1Bjay2WfD4z0mkVGGeU95ApGekM2zmMHq/1pthM4eRnlHJ1Cg+UlXGjBnD6NGj2bRpExs3biQnJ4c//OEP7N27lyuvvJKnn36aDRs2sGzZMoYPH84VV1zB8uXLWb58OcnJyXz++ecsX77ckogxYeJTj0REkoBrgE5AlIgAoKq1bgbggCs6Dtu+9jzPFDjlC56AVj2gdQ9o0Q2iY723t3JGUK8OD5b0jHQmLZ5EfnE+AFm5WUxaPAmAEV1G1Ljdzz77jNjY2LI5tiIjI3nuuefo3LkzIsK1117L2WefXVb/iiuuqPlOGGOCwtdDW/OAJcAqoCR44dQSxw7Bpv/BhnmweQEUHMWZn1Ir1o2Igq+fAy123ksENO/iJJZWPaD1Kc5zy26w9oNae3js6e+eZv2h9V6Xr9y/koKSghPK8ovzeXTRo8zcONPjOj2a9+CWU26pdLtr1qypMAljQkICHTt2ZPPmzVx77bU+7oExJlx8TSSxqnp3UCMJt4NbnMSx4SPY8Q1oCTRtA73GwMmXwLEDMO/eEw9vRcfByH9Az1FwcDPsWwf71zuPfeudtsoSTCSIQEm5k91KD4/V8l5J+SRSVbkxpuHwNZG8ISI3AR8CZTeEUNVDQYkqGFyHlM7LzoRl7eGChyHpJNj4kfOBf2CjU69NLzjnbid5JJ8OEW7DSJEx3g9LtTnVebgrOg4HNv2cXL58xnNs2ZmB399qemDAA5UuHzZzGFm5FU/Ua9ekHf8d/l+v6x2tYhr5nj17MnPmiT2aI0eOsGPHDoYOHcoPP/zAqFGjKm3DGBNeviaSAuAZ4A/8fHxHgS7BCCrg3M64EnAOKb3/O2dZRBR0OgfO+C10Hw7NTvLeTnUnLYxqBG17OQ+AFdO93NOjfcWyWmZCvwknjJEAxEbGMqGSlZyVAAAgAElEQVTfBL/aHTJkCBMnTuT111/nmmuuobi4mHvuuYfrrruOe++9lwEDBjBixAjOPPNMAGbNmsWgQYNo06aNX9s1xgSOr4nkHiBVVQ8EM5ig8XTGFUDjFnDHMogN0R0Vhzx64hgJQFRswO7pEUylA+qTf5zMntw9tG3Slgn9Jvg10A4gIrz//vvccsstPPHEE5SUlHDJJZfwl7/8hUaNGjF9+nTuvfde9u3bR0REBOeeey7Dhw8PxC4ZYwLE10SyGThWZa3aytuho2OHQpdEoOI9PVBI7l/rx0dKjegywu/E4UmHDh2YO3eux2Vnn302X331ldd1t23bFvB4jDHV42siyQWWi8jnnDhGUjdO/01sX3sOKbkfHvv4Qfj2385Af4uuoY/FGGMCwNcpUmYDfwYWAz+4PeqGIN8mtsYG3ekM4HsbhDfGmDrApx6Jqr4W7ECCyu2QkmZnIrXlQsD4NnDGjbDkBfjFvdAyNbzxGGNMDVTaIxGRGa7nVSKysvwjNCEGSO+r4K7VfDF4Nty1OvxJpNSgCRDZCL78a7gjMcaYGqmqR1J6buelwQ6kwWraGgbcBN887/RKWnUPd0TGGFMtlfZIVLX0CrRbVHW7+wOofO4L47tBEyAqDr54OtyRGGNMtfk62H6hh7KLAxlIg9akpdMrWf2eM7VKAxIZGUnfvn3p06cP/fr1Y/HixYBzWm9cXBx9+/YtexQU2HQsxtRGVY2R/F5EVgEnlxsf2QrUrTGS2m7gHRDTpFb3SrLnzmXTBUNYd0pPNl0whGwv135UR1xcHMuXL2fFihU8+eSTPPjgg2XLunbtWjZd/PLly4mJifF7e6EiIsNFZIOIbBaRiZXUu1xEVETSQhmfMYFUVY/kbWAkMMf1XPror6oN68ZWwdakBZz5O1jzPuxdG+5oKsieO5esRx6laPduUKVo926yHnk0IMmk1JEjR2jWrFnA2gsXEYkEpuD02nsC40Skp4d68TjjkN+GNkJjAqvSwXZVzQaygXGhCaeBO/s2+HYqfPEUXPV6SDe95y9/4fg674fV8lasQMsdWtL8fLL+8DCHZ7zrcZ1Gp/Sgye23V7rdvLw8+vbtS35+PllZWXz22Wdly7Zs2ULfvn0BGDRoEFOmTPF1d8JtALBZVTMARGQ6MAoo/w3hCeBp4L7QhmdMYPl6ZbsJhcbN4aybnQsU96z+ebLHWqB8Eqmq3Felh7YAvvnmG6655hpWr14N/Hxoqw5KAdynUsgEznSvICL9gA6qmi4iXhOJiIwHxgN07NgxCKEa4z9LJLXN2bc606Z88RRc/WbINtv2oYcqXb7pgiHOYa1yopKTOekN772nqqaRd3f22Wdz4MAB9u/f7/M6dZGIRAB/A66rqq6qTgWmAqSlpXm4s5ox4efrWVsmVOKawVm3wLq5kFV7zmdofdedSOyJtxCW2Fha33VnwLaxfv16iouLadGiRcDaDJNdQAe39+1dZaXigV7AQhHZBpwFzLEBd1NXWY+kNjrr97DkRVj4FIx7O9zRAJA4ciQA+577O0VZWUS1a0fru+4sK6+p0jESAFXltddeIzIy0u94w+x7oJuIdMZJIGOBX5UudI09tix9LyILgXtVdWmI4zQmICyR1EZxSc4hroV/gd3LIblvuCMCnGTib+Ior7i42GN5p06dysZK6hpVLRKR24D5QCTwiqquEZHHgaWqOie8ERoTWHZoq7Y662aITXJ6JabOUdV5qtpdVbuq6p9dZY96SiKqOth6I6Yus0QSYukZ6QybOYzer/Vm2MxhpGeke64YmwgDb3PuKb/rx9AGaYwx1WCJJITSM9KZtHgSWblZKEpWbhaTFk/ynkwG/M4ZfF/4ZGgDNcaYaqj1iURERovISyLyjogMC3c8/pj842Tyi/NPKMsvzmfyj5M9rxCbAANvh02fEH9kQ1BiUrUzSr2xn40xvglqIhGRV0Rkn4isLlfu0zxEAKo6W1VvAm4Grg5mvMG2J3dPtcoBGDAeGreg07bpAY8nNjaWgwcP2gemB6rKwYMHiS13yrMxpqJgn7X1KvA8UHbFmts8RBfiXPH7vYjMwTm7pfwxnBtUdZ/r9cOu9eqstk3akpWb5bHcq0bxMPAOWnz6R9j5HXQYELB42rdvT2ZmZlAvAMzPz6+TH8b5+fkkJSXRvn37cIdiTK0X1ESiql+KSKdyxR7nIVLVJ/FwAy0REeAp4CNV9Trq7D6VRJs2bVi4cKHHejk5OV6XBduFsRfyZu6blFBSVhYt0VwYe2GlMUUU9+DMqARyZ93Hyj6PhSDSwMnJyaFp06bhDqPacnJyyM7OZvv27eEOxZhaLxzXkVQ5D1E5twNDgUQRSVXVf3mqVH4qicGDB3tsbOHChXhbFmxnFp3JjOkzAMrGSn7X53f8rs/vqlx3y64xdM14lcFdYqHjWUGNM5DC+fP2R12N25hwqPWD7ar6D1Xtr6o3e0sidcVHWz8ivzifF4e+yNdjvyYmIoaD+Qd9WndXyiXQpDV8/pcgR2mMMdUTjkRS1TxE9ZKqMm39NLo160b/Nv1JbJTIkI5DmLd1HgXFVc+gWxLZCM65E7Z+AdsWhSBiY4zxTTgSSdk8RCISgzMPUb2fMmLZvmWsP7SecT3G4Qz7wOjU0WQfz+bznZ/71kjaDdC0jV1XYoypVYJ9+u804BucW/VmisiNqloElM5DtA6YoaprghlHbTBt/TTiY+IZ0XlEWdmZ7c6kTeM2zN4827dGouPgnLtg21ew9asgRWqMMdUT1ESiquNUtZ2qRqtqe1V92VVeYR6i+mzfsX18uv1Tfpn6SxpHNy4rj4yI5LKul7F492L25u71rbH+10F8O6dXYtd/GGNqgVo/2F4dIjJSRKZmZ2eHO5QTzNw4k2It5uqTK15POSp1FCVawocZH/rWWHQcdB4M2xfBY83guV6wckZgAzbGmGqoV4lEVeeq6vjExMRwh1KmsLiQdze+yzkp59AxoeKtUk9KOIl+rfsxe/Ns364wXzkD1pYeClPI3glz77BkYowJm3qVSGqj/23/HwfyDvCrU37ltc7o1NFsO7KNFftXVN3ggsehKO/EssI8p9wYY8LAEkmQTVs/jY7xHRmYPNBrnWGdhhEXFefboHt2ZvXKjTEmyCyRBNHag2tZvn85Y3uMJUK8/6ibRDfhwpMu5ONtH5NXvrdRXqKXuZ9imkBJiedlxhgTRJZIgmja+mnERcUxKnVUlXVHp44mtzCXBTsWVF5xyKPOgLu7iCgoyIEPboXiIj8iNsaY6rNEEiSH8w8zL2MeI7uMJCEmocr6/dv0J6VpStWHt3pfBSP/AYkdAHGeR78I5/8BVrwNM6+DouMB2QdjjPFFOCZtDBoRGQmMTE1NDXcozNo8i4KSAsb2GOtT/QiJYFTqKF5c/iK7c3aT3DTZe+XeVzmP8hrFw8cTYdo4uPpNiGlcsY4xxgRYveqR1JbTf4tLinln/Tuc0fYMujXr5vN6o7o6h8A+2PJBzTZ81u/hsuch43N4cwzk167raYwx9VO9SiS1xReZX7A7dzfjeoyr1nrJTZMZ0G4AH2z+gBKt4cB5v/8HV7wCmUvhtZGQe6Bm7RhjjI8skQTBtPXTaNO4Ded3OL/a647qOopdObv4Ye8PNQ/g1F/CuGmwfwP89xI4srvmbRljTBUskQRYxuEMlmQt4eqTryYqovpDUENPGkrT6Ka+T+ToTbcL4TfvOUnkleFwaKt/7RljjBeWSAJs+obpREdEM6bbmBqtHxcVx0WdLuJ/2/9HbmGuf8F0Ogeu/QCOH3GSyb51/rVnjDEeWCIJoJyCHD7Y/AEXd76YFnEtatzO6NTR5BXl8cm2T/wPKqU/XDcPUOcw1+5l/rdpjDFuLJEE0JwtczhWdKzag+zl9WnVh04Jnfw/vFWqTU+4/iOIaQqvjoTtiwPTrjHGUM8SSTinkS+9le5pLU+jV8tefrUlIoxOHc2P+35k+5HtgQmwRVe44WNIaAdvjIFNnwamXWNMg1evEkk4ryNZkrWEbUe2+d0bKTWy60giJIIPNtfwmhJPElOcw1wtU2HaWJh3v3M/k0lJdl8TY0yN1atEEk5vr3+b5rHNuajTRQFpr3Xj1gxMHsicLXMoLikOSJsANG0F134ISR3hu3879zOx+5oYY/xgiSQAduXs4oudX3B5t8uJiYwJWLujU0ez99hevs36NmBtAhCXBEX5FcvtvibGmBqwRBIA72x4hwiJ4KqTPcx/5YfzO5xPYqPEwA26u/N2kaLd18QYU02WSPyUX5TPrE2zuKDjBbRt0jagbcdExnBJ50tYsGMBx4qPBbRtr/c18VZujDFeWCLx00dbPyL7eHbABtnLG506moKSAn489mNgG/Z0XxOJcKajN8aYarBE4gdV5e31b5OalEpam7SgbOOU5qfQrVk3luQsCWzD5e9rEtcctAT2rAzsdhooERkuIhtEZLOITPSw/G4RWSsiK0VkgYicFI44jQkESyR+WLF/BesPrWdcj3GISFC2ISKM7jqa7QXb2XJ4S0DbTm/ahGEdkunduSPDunYjvfdIWPICrHw3oNtpaEQkEpgCXAz0BMaJSM9y1ZYBaaraG5gJ/DW0URoTOPUqkYT6gsS3179NfHQ8l3a5NKjbubTrpUQQEdBB9/SMdCYtnkRWbhaKkpWbxaRjG0g/qQ/MuR32rA7YtuqS9Ix0hs0cxu3bb2fYzGGkZ6TXpJkBwGZVzVDVAmA6cML9llX1c1UtHfhaAtjglKmz6lUiCeUFifuP7ed/2/7H6G6jaRwd3DsRNo9tTq+4XszdMpfCksKAtPn3H/5OfvGJpwDnF+fzXONI5/Tgd34Nxw4FZFt1hXtyBZzkunhSTZJJCrDT7X2mq8ybG4GPPC0QkfEislRElu7fv7+6cRgTEvUqkYTSzI0zKdIixp7s2610/XVW07M4mH+Qxbv8nydrxf4V7Dm2x+OyvfkHuKbrKbxWcphd710HgbwYshYr0RKeXfqsx+Q6+cfJQduuiPwGSAOe8bRcVaeqapqqprVq1SpocRjjj3p1z/ZQKSwuZMbGGZyTcg4dEzqGZJs943rSPLY5szfP5rwO59WojePFx5mybAqvrX2NCInweBfGptFNyYuM4tnmCTxbspVTpp/P0F6/YehJQ+mS2MXf3ahVjhQc4Zvd3/BV5lcs2r2IA3me7ya5J9dz0q3ELqCD2/v2rrITiMhQ4A/Aeap6vLobMaa2sERSA5/u+JQDeQeCdsqvJ5ESyaVdLuXt9W/zU/5PNIttVq31V+xfwSOLHmFr9lau6H4FvVr04qnvnjrhG3hsZCwPn/UwI7qMYOeRHSz46DY+Pbyefy77J/9c9k+6JHZhSMchXHjShfRo3gMRIT0jnck/TmZP7h7aNmnLhH4TGNFlRKB3PyBUlQ0/beDrXV/zVeZXrNi/gmItJiEmgYHJA1mStYTDxw9XWK8G1wd9D3QTkc44CWQs8Cv3CiJyOvBvYLiq7qvRDhlTS1giqYbSD82s3CwiJZLs46GdZXh06mheX/s66Rnp/Kbnb3xa53jxcaYsn8Jra16jdePW/PvCfzMweSAAsVGxXpNAh4SOXDdmBtf9dzh7s7bw2dD7WHBwJS+vfpmXVr1EStMUuiR04ds931JQUgD8PKYAVCuZBCoZeWrn3PbnsiRridPr2LWIfXnOZ/YpzU/hhl438Iv2v+C0lqcRFRFVNkZSPrlO6DehWnGoapGI3AbMByKBV1R1jYg8DixV1Tk4h7KaAu+6zvjboaqXVXunjakFLJH4qPyHTLEW8/g3jxMhESH7Bt6tWTdObXEqszfP9imRrNy/kocXPczW7K1c3u1y7k27l6YxTcuWj+gyovLYo2PhqjdoM/U8xn37FuN++yk/aRELdy7k0x2f8mXmlxVWyS/O589L/kxuYS4JjRLYmreVFvtbEB8TT0KjBOJj4omOiC6rn56RzqSvHyFfnZMIsnKzmPT1I2Xx+ar87ycrN4uHvnoIdf2Lj47n7OSz+UX7X3BOyjm0jGtZoY0RXUbAjiVMznifPRHQtgQmnDS8Rr9fVZ0HzCtX9qjb66HVbtSYWkpUNdwxBFxaWpouXbrU47KFCxcyePDgarc5bOawsrN53LVr0o5PrgjAnQyrUBr39PXT+fO3f+bdke/So3kPj3XL90IeO/sxBqYMrPnGM76AN0bDKSPhytfAdc1M79d6o1T/7ycuKs5JLDEJbM/eSqFWHNCPi4jm3I4XUFRSRHFJMYVaSFFJ0QmPYi0ue73z6E6KPbTTNLopzw95nj6t+hAVUcX3ppUznBmQC/N+LouOcy7c7P3zPGoi8oOqBucK1EpU9ndtjL/8+bu2HomPvA241mAg1i8Xd76YZ75/htmbZzNxQIULplm1fxUPL3qYjOwMLu92Ofek3UN8TLx/G+1yHgx9DP73CCyaDOfcCThjB56Sa9vGbXlrxFscLTjKF0u+ILVXKtnHszlacJSjBUc5UnCk7HnzT5vKEpO7vOICNhzaQFRElPOQqLLXcVFxZa+jI6KJlEi2HdnmMfTcwlz6t+nv234uePzEJAI/z4jcO7ATchpTn1gi8ZHXD80AT9RYlcRGiVzQ8QLSM9K5p/89REc6h4mOFx/nheUv8OqaV2kV14p/Df0Xg1IGBW7DA2+H3T/CgsegXR/oej4T+k3wOKZwZ/87ad24Na0bt2Zn7E7ObX+u5zYLjjFsy3yyoiv+GbYrKmbuL55zbsLlgxX7V9T896MKmd+77s3igc2IbEylLJH4qHfL3hU+qGoyEBsI7Zq04/Dxw/R7sx/tmrTjl6m/5ONtH5ORncGYbmO4N+1e/3sh5YnAZc/DvvUw8wb43RdlYwfVHig/uge+ewmWvsyEiONMatmc/IifL2mKLSlhwk+H4fn+kJIGfcZCr8uhcXOvTXpLapX+fvath1XvOo/DldzS2GZENqZSlkh8kJGdwec7P6dXi14czD8Y1lNd0zPSmbZ+Wtn7rNwsXljxAgnRCbw49EXOSTkneBtv1BTGvgVTB8M7v4Eb5lc9YO8ua6Uzl9eqmVBSBD1GMKJ1T1g+lckJjdkTFUnbomImHDnGiHMfc26+tfIdmHcvfPwgdL/ISSrdhkFUoxOa9jmpZe+C1TOd5LFnlTPjcefzYPBEKC6Ejx+oOEYy5FGMMd7Vq0QiIiOBkampvh0O8UWJlvDY4seIjYrln0P+6fFsn1Ca/OPkCldfA8RFxwU3iZRq0RXGTHXu+f7h3TD6BY9jHGW0BDZ8BN9MgW1fQUxTOONGOPN30Ny5wHFEy26MWPC4cwgpsT0MfebnMYlBdzgf+CumOx/+6z+E2CSnh9JnLLQ/o2z7I3JyGbFzt6udEuie67Rx7BCsm+NMRrl9EaCQ0h+GPwWnjoH4Nj/HGx0HCx5HszORxPZOErHxEWMqVa8SiarOBeampaXdFKg2Z2yYwY/7fuRPg/4U9iQC3gf39x0L4TVtJ18M5z0AXzwN7fvDGb+tWKcgF5a/zYDv/gZ5uyGhPVz4BPS7xpnLy13vqyr/sG57mvMY+hhsXegkleVvw9KXnWTUe6zTW/rsiZ97E9k74YNbYfEU2LcGSgqhRSoMfhBOu8JJiJ64Yvmihmf3GdMQ1atEEmhZOVk898NzDEweyGVda8e1YrVl0J/zJsLuZZB+Hyx8GnL3O72JgbfD0SxY+l/IP0xRfDe44hU45TKIjK663cpERkHqUOdx/CisnQMrp8PCJ8HTacjFBbB3JZx1i5M82vWtvPdkjKkRm7TRC1Xl8SWPoyiPnv1o0O43Ul0T+k0gNjL2hLKwDPpHRMDJlwAKufuc5+yd8NH98PVz0PlcuOETfuz3jHMYyt8kUl6jeDj913DtXLhzlfd6qnDRnyH5dEsixgSJ9Ui8+DDjQ77e9TUTB0wkpWllM4CHVo3PlAqGr/4Pjz2B+HZw9RvO64yFwY8jqYNzp0dPp+/aGVfGBJ0lEg8O5h3k6e+fpk+rPiGbJr46qnWmVDB5u77iaGgv0gScQXFPV6XbGVfGBJ0d2vLgqe+e4ljhMR4b+BiREZHhDqf28vZtPxy9gPL3oE/sUGFqE2NMcFiPpJzPd3zOx9s+5ta+t9I1ycuZPcZR23oBVZ39ZYwJCuuRuDlacJQ/LfkT3Zp148ZeN4Y7nNrPegHGGKxHcoK//fA3DuQfYPIFk8vmsDJVsF6AMQ2e9Uhcvsv6jpkbZ3JNz2vo1bJXuMMxxpg6wxIJkFeUx6RvJtEhvgO39L0l3OEYY0ydYoe2gBeWv8DOozt5edjLxEXFhTscY4ypU+pVj0RERorI1Oxs3++lvubAGl5f+zqXd7ucAe0GBDE6Y4ypn+pVIlHVuao6PjEx0af6hcWFPLL4EVrGtuTutLuDHJ0xxtRPDfrQ1iurX2HTT5v4x/n/ICEmIdzhGGNMnVSveiTVseXwFv698t8M7zSc8zueH+5wjDGmzmqQiaS4pJg/Lv4jjaMbM3HAxHCHY4wxdVqDTCTTN0xnxf4VPHDGA7SIaxHucIwxpk5rcInkYNFBJv84mUEpg7i0y6XhDscYY+q8BjPYnp6RzuQfJ5OVm4UgnJN8Tq25WZUxxtRlDaJHkp6RzqTFk8puUasok3+cTHpGepgjM8aYuq9BJJLJP04mvzj/hLL84nwm/zg5TBEZY0z90SASyZ5cz3fs81ZujDHGdw0ikbRt0rZa5cYYY3zXIBLJhH4TiI2MPaEsNjKWCf0mhCkiY4ypPxrEWVsjuowAKDtrq12TdkzoN6Gs3BhjTM01iEQCTjIZ0WUECxcuZPDgweEOxxhj6o0GcWjLGGNM8FgiMcYY4xdLJMYEgYgMF5ENIrJZRCrMDCoijUTkHdfyb0WkU+ijNCYw6lUiqckdEo0JNBGJBKYAFwM9gXEi0rNctRuBn1Q1FXgOeDq0URoTOPUqkVT3DonGBMkAYLOqZqhqATAdGFWuzijgNdfrmcAQscnfTB1VL8/a+uGHHw6IyHYvi1sCB0IZT4BY3KFVWdwnVbFuCrDT7X0mcKa3OqpaJCLZQIvy2xSR8cB419vjIrK66tCDIly/x3D+/TS0fT65pivWy0Siqq28LRORpaqaFsp4AsHiDq3aEreqTgWmQnhjCte2bZ9Du92arluvDm0ZU0vsAjq4vW/vKvNYR0SigETgYEiiMybALJEYE3jfA91EpLOIxABjgTnl6swBrnW9vgL4TFU1hDEaEzD18tBWFaaGO4AasrhDq8Zxu8Y8bgPmA5HAK6q6RkQeB5aq6hzgZeANEdkMHMJJNkGLKQDCtW3b5zqwXbEvQcYYY/xhh7aMMcb4xRKJMcYYvzSoRCIikSKyTEQ+DHcsvhKRJBGZKSLrRWSdiJwd7ph8ISJ3icgaEVktItNEJLbqtcJDRF4RkX3u12iISHMR+Z+IbHI9NwthPGGZXsWH7d4tImtFZKWILBCRqq6nCdi23epdLiIqIgE5PdaX7YrIVa79XiMibwdiu75sW0Q6isjnrs+slSJySYC2W+HvvdxyEZF/uOJaKSL9qmxUVRvMA7gbeBv4MNyxVCPm14Dful7HAEnhjsmHmFOArUCc6/0M4Lpwx1VJvOcC/YDVbmV/BSa6Xk8Eng5RLJHAFqCL6/e9AuhZrs4twL9cr8cC74Rou+cDjV2vfx+I7fq6bVe9eOBLYAmQFqJ97gYsA5q53rcO4e95KvB71+uewLYAbbvC33u55ZcAHwECnAV8W1WbDaZHIiLtgRHAf8Idi69EJBHnl/4ygKoWqOrh8EblsyggznWNRGNgd5jj8UpVv8Q5c8qd+xQmrwGjQxROuKZXqXK7qvq5qh5zvV2Cc31MIPiyzwBP4MxJlh/C7d4ETFHVnwBUdV8It61Agut1IgH6P+Tl793dKOB1dSwBkkSkXWVtNphEAvwduB8oCXcg1dAZ2A/819W9/Y+INAl3UFVR1V3As8AOIAvIVtVPwhtVtbVR1SzX6z1AmxBt19P0Kine6qhqEVA6vUqwt+vuRpxvrYFQ5bZdh1c6qGp6gLbp03aB7kB3EVkkIktEZHgItz0J+I2IZALzgNsDtO2qVPdvoWEkEhG5FNinqj+EO5ZqisLpgr6oqqcDuTiHWWo113jCKJxEmAw0EZHfhDeqmlOnv2/nybu4fpdpwDMh2l4E8DfgnlBsr5wonMNbg4FxwEsikhSibY8DXlXV9jiHm95w/SxqnVoZVBAMAi4TkW04XcgLROTN8Ibkk0wgU1W/db2fiZNYaruhwFZV3a+qhcAsYGCYY6quvaXdeddzoA5pVCVc06v4sl1EZCjwB+AyVT3u5zZ93XY80AtY6Po/fBYwJwAD7r7scyYwR1ULVXUrsBEnsfjLl23fiDO+iKp+A8TiTOgYbD79LbhrEIlEVR9U1faq2glncPIzVa3135BVdQ+wU0RKZ+UcAqwNY0i+2gGcJSKNXcfuhwDrwhxTdblPYXIt8EGIthuu6VWq3K6InA78GyeJBDKxVrptVc1W1Zaq2sn1f3iJK4YaTzLoy3ZdZuP0RhCRljiHujL83K6v296B838HETkFJ5HsD8C2qzIHuMZ19tZZOIemsypdIxBnAdSlB84fRV06a6svsBRYifNH3SzcMfkY92PAemA18AbQKNwxVRLrNJyxnEKcb6A34ow5LAA2AZ8CzUMYzyU433y3AH9wlT2O8+EJzgfKu8Bm4DugS4i2+ymwF1jueswJ1T6Xq7uQAJy15eM+C85htbXAKmBsCH/PPYFFOGd0LQeGBfHv/WbgZrd9nuKKa5UvP2ubIsUYY4xfGsShLWOMMcFjicQYY4xfLJEYY4zxiyUSY4wxfrFEYowxxi+WSMJMRNqKyHQR2SIiP4jIPBHpXsU6d4pIYx/avlKcGYM/d72f5prN865y9SaJyC4RWe6a5XRcDfYjx/WcLCIzq7u+MabustN/w4hRBYQAAALJSURBVMh1sd5i4DVV/ZerrA+QoKpfVbLeNpxzuw9U0f7HwJ9U9WsRaQt8raqpHupNAnJU9VkR6Qb8ALRQ56p0X/clR1Wb+lrfGFN/WI8kvM4HCkuTCICqrlDVr0RksLjdN0VEnheR60TkDpz5qz5362mME5FV4tz742lX2aPAOcDLIvIM8AmQ4up1/MJbQKq6CTgGNHO1c5OIfC8iK0TkvdKekOuK3G9c2/2TW5ydSu9z4Ir3ebdlH7r2K1JEXnXFu6p8D8kYU7dYIgmvXjjf/n2mqv/AmU76fFU9X0SScabWvgDnKvgzRGS0qj6Oc0X8r1X1PuAyYIuq9q2it9MP2KQ/T4ExS1XPUNU+ONOc3Ogqn4wzmeRpOFfJVkdfIEVVe7nW/2811zfG1CKWSOq+M4CF6kyQWAS8hXMPk+q6S0TWAN8Cf3Yr7yUiX4nIKuDXwKmu8kE4Uy2AMwVKdWQAXUTkn65puY/UIF5jTC1hiSS81gD9vSwr4sTfT7BvVfucqp4KXI5zOKx0e68Ct7l6Do+Vi6OqATaP+6DOTYL64MyZdDN16GZjxpiKLJGE12dAIxEZX1ogIr1dYxjbgZ7i3J87CdcsoC5HcabWBmfSvvNEpKWIROLcw+CLmgakqnNwDomVzi4bD2SJSDROj6TUIpwZSylX7m4b0FdEIkSkA85d4UpnUY1Q1ffg/7d3tygIRFEYht9Zh9txA67DIJgUtFjsRrOLsJkm+hPUDVjcxTGcK1hMJ4jwPmVg4MK0j7kcvsOC/6jGl/SFQfJDkSNzI2DYxn9vwBp4RsSD3EVwbc/zx9EtsO+67hBZ7zwDDmRL6DEiqpXnK2DalugsyeuunmzzfZsA43bl9W17Wk/ubr8DG+DU3g/I3RIXYAfMi98r6Ycc/5UklfhHIkkqMUgkSSUGiSSpxCCRJJUYJJKkEoNEklRikEiSSl4Cb7nAOIm+4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, axes = plt.subplots(ncols = 2,)\n",
    "\n",
    "for i in range(len(axes)):\n",
    "    ax = axes[i]\n",
    "    N = Nsolvent[i]\n",
    "    for key in Methods2:\n",
    "        ax.plot(maxdist, results2[str(key)+'.'+str(N)], '-o', label = key)\n",
    "    ax.set(xlabel='Cutoff Radius', ylabel='time [s]', yscale='log', title=' N Solvent = {}'.format(N))\n",
    "    ax.legend()\n",
    "    ax.grid(True)"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}