Skip to content

Instantly share code, notes, and snippets.

@richardjgowers
Created November 17, 2018 13:46
Show Gist options
  • Save richardjgowers/48aa5711c9fbf98b250026e09c238cf1 to your computer and use it in GitHub Desktop.
Save richardjgowers/48aa5711c9fbf98b250026e09c238cf1 to your computer and use it in GitHub Desktop.
Timing distance searching in MDAnalysis and BioPython
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Timing the performance of the new capped distance methods\n",
"\n",
"This notebook compares the performance of the new `capped` distance functions in `MDAnalysis.lib.distances` against the existing functions.\n",
"\n",
"These functions were implemented as part of Google Summer of Code 2018 by [Ayush Suhane](https://ayushsuhane.github.io) a graduate student at the University of British Columbia, Canada,\n",
"with help from [Richard J Gowers](https://github.com/richardjgowers), [Jonathan Barnoud](https://github.com/jbarnoud) and [Sébastien Buchoux](https://github.com/seb-buch). \n",
"\n",
"We will look at two cases, which are similar to finding bonds and radial distribution functions.\n",
"\n",
"This notebook was prepared using a development version of 0.19.3 of MDAnalysis. The capped distance functions are available in version 0.19.0 onwards."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.19.3-dev\n"
]
}
],
"source": [
"import numpy as np\n",
"import MDAnalysis as mda\n",
"from MDAnalysis.lib.distances import (\n",
" capped_distance, distance_array,\n",
" self_capped_distance, self_distance_array\n",
")\n",
"\n",
"from Bio import KDTree\n",
"\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib notebook\n",
"\n",
"import seaborn as sns\n",
"sns.set_context('notebook')\n",
"\n",
"print(mda.__version__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will first define a function to generate random coordinates for a given box size.\n",
"These coordinates will have a typical particle density for atomistic soft matter molecular simulations."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# typical atomistic particle density\n",
"# 4055 water atoms in 5nm cube\n",
"density = 4055 / (50. ** 3)\n",
"\n",
"def make_coords(boxsize): # boxlength in A\n",
" # makes 3d coordinates in cube of boxsize side length\n",
" n = int((boxsize ** 3) * density)\n",
" crds = np.random.random(n * 3).reshape(n, 3) * boxsize\n",
" \n",
" return crds.astype(np.float32)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### guessing bonds\n",
"\n",
"To guess which atoms are bonded (usually due to missing information in the topology) we have to find all pairwise distances within 2A. (These distances are then further processed to match a particular distance for a given pairing of atom types)\n",
"\n",
"Here we generate different box sizes from 20 to 75 A, timing how long to find all pairwise distances within 2A. We are using the `self_` variant of the distance methods as we are comparing an array of coordinates only with itself."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def kdtree_bond(coords, max_cutoff):\n",
" kdt = KDTree.KDTree(3)\n",
" kdt.set_coords(coords)\n",
" kdt.all_search(2.0)\n",
" \n",
" idx = kdt.all_get_indices()\n",
" dists = kdt.all_get_radii()\n",
" \n",
" return idx, dists"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"natoms = []\n",
"brute = []\n",
"capped = []\n",
"kdt = []\n",
"\n",
"for boxsize in np.linspace(10, 100, 26):\n",
" a = make_coords(boxsize)\n",
" \n",
" natoms.append(a.shape[0])\n",
" \n",
" t = %timeit -o -q self_capped_distance(a, max_cutoff=2)\n",
" capped.append(t)\n",
" \n",
" t = %timeit -o -q self_distance_array(a)\n",
" brute.append(t)\n",
" \n",
" t = %timeit -o -q kdtree_bond(a, max_cutoff=2)\n",
" kdt.append(t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the `capped` variant is substantially faster, whilst the brute force method is scaling $\\mathcal{O}(n^2)$ as expected."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl8VNXZwPHfyZ4QICEsCVmAsAUISQDZVFTAhQJuVRQRV1xqaalttdrqq3ZBqdpKbe0LuIO1rrhR9EXqAgouyCZbIAQSkpAA2bdJJjPP+8dMxhCyTCCZyfJ8P5/7yb13zj33mRDmmXPvuecYEUEppZRqb3y8HYBSSinVEE1QSiml2iVNUEoppdolTVBKKaXaJU1QSiml2iVNUEoppdolTVBKKaXaJU1QSiml2iVNUEqpDs0Y09cY85UxZq8xJsHb8ajWY3QkCaVUR2aMuQcoB74Ffioit3o5JNVKtAWlvMYYs9sYc8FpHvuSMeZPrVHXmTDGHDbGXFh/vZXP4dH31tz5mnufXvi3OAD4AwFAmQfPq9qYn7cDUJ2fMeYw0A+w1dk9TERGtdY5WrOu9sbT763++Zz/freJyPrTOd5dxphA4J/AhUAvIA34nYh82MyhW4A3gCpg8OmcW7VP2oJSnnKpiITWWXK8HVBHZ4zpbF8w/YAjwPlAT+B/gDeMMQObOe73OFpO3dAWVKeiCUp5Tf1LRc7te4wxO40xxcaY140xQc7XxhhjthpjSo0xrwNBjdXVVD3O18caY7Y563rT+fqfaIAx5n5jzEFn2T3GmCvP4L3+1llHoTHmxXoxNXqeBt7bfcaYnUC5McbPuZ3tPDbVGDO9gfPfYoz5oM52mjHmjTrbR4wxKQ2cbxUQB3xgjCkzxvzGeUhKE79ft/8t6hKRchF5REQOi4hdRNYAh4BxTfxeE4CbgZ/haKEnNlZWdTyaoFR7cw0wAxgEJAE3G2MCgHeBVTgu/bwJXNXSegCcdb0DvOSs699AU0nnIDAFxzf63wOvGGOiWv62ALgeuATHZahhwIOneZ7rgFlAmLOunwHjRaS7s/7DDRzzOTDFGOPjrNcfOAfAGBMPhAI76x8kIjcAmfzQAn7c+VKDv99GtKSsizGmH47f0+4mij0GfCci/wb2A8nu1K06Bk1QylPeNcYUOZd3myj3tIjkiEgB8AGQAkzC8YG6VESsIvIWjh5bTWmoHpx1+Tlft4rIauCbxioRkTed9dhF5HUcN+QnuPOGG/APETnijGkxjkRzOud52llPJY5WQyAw0hjj72x9HGzgfaQDpTh+D+cD/wdkO1sg5wMbRcTegvfS2O/3TMsCYIzxB/4FvCwi+xopMxm4ArjPuWs3jgSoOglNUMpTrhCRMOdyRRPlcuusV+D4Zt8fyJaTn4nIaOZ8DdVDI3UdaawSY8yNxpjttckVxyWk3s2cuzF1z5PhjOV0zuOqR0TSgLuBR4BjxpjXjDH9Gznuc+AC4Dzn+mc4ktP5zu2WaOz3e6ZlMcb44GgtV+NoHTbmz8BaEfnMub0bbUF1KpqgVEdwFIg2xpg6++Jasa7YhgoaYwYAz+L4kIwQkTBgF2AaKu+GuueJA3JO8zwnPbwoIq+KyLnAAOdrf27kuNoENcW5/jnuJSiPPSzp/Hd5Hkevz6tExNpIudk43se5xphcY0wu8Gu0BdWpaIJSHcFmoAZY5OwU8GNO/zLbZhyXxX7mrOvyJurqhuPD+Tg4OhpwZjfhFxpjYowxvYDfAa+f6XmMMcONMdOcXbQtQO1lv4Z8DkwFgkUkC9iI495QBLCtidPkAfHuxNMK/hcYgeOeV2VDBZwtrEeBZcBwHJcMU4CLgTBjzOl+eVHtjCYo1e6JSDXwYxw31wuBa4HVZ1jXAqAImA+swfEMTf2ye4C/4EhqecBo4MvTOa/Tq8A6IN25/KkVzhMILAFO4LiU1hdH8juFiOzH0Q17o3O7xBnHlyLSWFIDR0eEB52XH+9xM64Wc7Yk78SRbHKdvQbLjDHX1yt6A47OJPeKSG7tguNeYinaiuo0dKgj1eUZY74GlonIi214jsO04GFXpZS2oFQXZIw53xgT6bzEdxOOb9wfeTsupdTJOtuT6Eq5YziOoXFCcTx/dLWIHPVuSEqp+vQSn1JKqXap07SgnL2YxuPoRtzUDV+llFKe5wtEAd+KyCmdkhrSaRIUjuS00dtBKKWUatIU4At3CnamBHUUYOPGjcTExHg7FqWUUnVkZWUxZcoUcH5Wu6MzJSgbQExMDAMHDvRyKEoppRrh9i0Y7WaulFKqXdIEpZRSql3SBKWUUqpd6kz3oBpltVrJysrCYrF4OxSlXIKCgoiJicHf39/boSjVLnWJBJWVlUX37t0ZOHAgJ8+yoJR3iAj5+flkZWUxaNAgb4ejVLvUJS7xWSwWIiIiNDmpdsMYQ0REhLbqlWpCl2hBAZqcVLujf5OqOc9tTOedbdn4GIOPj8HXgK+PwRiDrzH4+vyw/4cyP+z3MeBbZ79rn4/Bx1mudr2h/cZ5vK+PIcDPhxsnD/To++8yCUoppTqSgvJqnlyXSkx4CLHhwdgF7CLY7I6lxm6nqkZO2W8X5z67YHPur11vaL9dcLzm3N/Y8KzB/r6dP0EZYx4GHgFGi8iueq+FAC8C43DMoHqPiKzxdIxKKeVtqzZnYLHa+d/rxzK0X3ePnVdqk1fdxCeC2D0WgotH70EZY8YCk4DMRorcA5SKyBDgUuA5Y0yop+Jrj1566SWuvvpq1/aDDz5IQkJC7ZAhbrngggtYs8aR52+77TY2bmx6yMKlS5dy7Nix0wtYKXXGLFYbKzcfZlpCX48mJ3Bcevbz9SHAz4cgf1+6BfrRI8ifniGe723qsRaUc7TxZ4B5wKeNFLsWuAlARA4YY7YAPwLerFdXGBBW71i3BuD7/Qe72ZNT0oLI3Teyfw8evnRUm9Rd6y9/+QuZmZn06dPntI5/7rnnmi2zdOlSLrzwQvr27Xta5/AWu92OMeakezs1NTX4+emVbNWxvPVdFvnl1dxxXry3Q/EqT7ag/gC8IiKHmigTB2TU2c4EYhsodzdwqN7SIUYyr6ioYM6cOYwcOZLk5GSuueYaAF5++WUmTpzIuHHjmDZtGqmpqaccO2XKFCwWC9OnT+fee+9t9Bx79uxh4sSJjB07lvnz55/UU6xua2rFihWMGDGClJQUkpKS2LdvH4sXLyYnJ4err76alJQU9uzZw3//+18mT57MmDFjGD16NK+99tpJ9d17772ce+65xMfHc//997tey87O5qqrriIpKYmkpCQee+wxAEpKSrjtttuYMGECSUlJ/OIXv8Bma3x4rtzcXKZOncq4ceMYNWoUv/nNb1yvPfLII8yfP58rrriC5ORkioqKGDhwIH/84x+ZOnUqd955Z6PHV1ZWEhUVxdGjP4xduWjRIh599NFGY1GqrdnswnMb00mO6cnEQb28HY53iUibL8Bk4BN+mCDxMJDYQLlSoE+d7X8Cv2qgXBgwsN5yLiCHDh2S+vbs2XPKPm9ZvXq1TJ8+3bVdUFAgGzZskJkzZ4rFYhERkbVr18rZZ58tIiIvvviiXHXVVa7ygJSWljZ5jrFjx8pLL70kIiKbN28WHx8f+eCDD0RE5Pzzz3et9+jRQzIzM0VExGKxSHl5uYiIDBgwQL7//vuTYqypqRERkdzcXImOjpaCggJXfddcc43YbDYpKiqSiIgI2b9/v4iIXHDBBfL444+76jl+/LiIiCxYsEBWrlwpIiI2m03mzp0rK1asaPT9VFZWut5zdXW1TJ06VT788EMREXn44YclNjbWVXdt/HfddZdbx993333yyCOPiIhIWVmZ9OnTR/Ly8pr8/bam9vS3qdqHD7/PkQH3rZE1O3K8HUqrOnTokAACDBQ3c4enrn2cDyQAh5yXX2KA/zPG3CIi6+qUywQGAMed23E0cDlQRIqAorr7OkqX3eTkZPbt28fChQu54IILmDVrFh988AE7duxg4sSJgONLQ2Fh4WnVX1JSwq5du7jhhhsAmDRpEqNHj26w7LRp07jlllu4/PLLmTVrFvHxDV9OOH78OLfeeisHDhzAz8+PgoICUlNTmTRpEgBz5szBx8eHnj17MmLECA4ePEhUVBSbNm3i448/dtXTu3dvAN5//32++eYb/vKXvwCOVmVTU6TYbDbuvfdeNm3ahIiQm5vL9u3bmTFjBgAzZ8501V3rxhtvdOv4hQsXcu655/LAAw+watUqLr744g53aVN1HiLC8g3pxPUKYUZipLfD8TqPXOITkSUi0l9EBorIQCALuKRecgLHvaY7AYwxQ3FMQviRJ2L0lPj4ePbu3ctFF13E+vXrSU5ORkS49dZb2b59O9u3b2fHjh1kZjbWj6R57ibr1atX8+ijj1JeXs7UqVP58MMPGyx31113ccEFF/D999+zfft2YmJiTrpsGBQU5Fr39fWlpqamyfOKCO+++67r/e7fv58nnnii0fJ//etfKSws5Ouvv2bnzp1cccUVJ50/NPTUfjR19zV1fGxsLOPHj+e9997jn//8JwsXLmwydqXa0paMQrZlFnHblEH4+nSML91tyesjSRhjthtj+js3nwDCjDFpwBrgDhEp9V50rS8rKwtfX1+uuOIKnnrqKY4fP86ll17KypUrycrKAhzf+L/77rvTqr9Hjx4kJiby6quvAvDNN9/w/fffn1KupqaG9PR0JkyYwP3338/FF1/Mtm3bXHUUFxe7ytbe1zHG8PHHH5OWltZsHKGhoZx99tk89dRTrn0nTpwA4LLLLmPJkiWu+04nTpzg0KHGb00WFRURFRVFUFAQ2dnZvPfee278Jtw//uc//zl33303fn5+TJ48uUV1K9Waln+eTniIP3PGNXTrvevxSoJytqR2OddTRCTHuV4uInNEZIiIDBeRln0SdQDff/89kydPJjk5mQkTJvDb3/6W8847j8WLF3PZZZeRnJxMYmJiiz+E61q5ciV///vfGTt2LM8++6zrUlxdNpuNm2++mdGjR5OcnMzRo0e58847AUdHgVtuucXVSWLJkiXcc889TJ48mbfeeoukpCS34njllVf48ssvSUxMJDk5meeffx5w9BL09fUlOTmZ0aNHM2PGDLKzsxutZ9GiRXz55ZeMGTOGu+66i+nTp7fo99Hc8eeffz5BQUH89Kc/bVG9SrWmtGNlrN+bxw2TBxIc4OvtcNqF2k4LHZ4xZiBw6NChQ6fMqLt3715GjBjhhahUR3Do0CHOOecc0tLSCAkJ8ei59W9T1br/7Z28sy2bTfdPIyI00NvhtLrDhw/XDow8SEQOu3OM1y/xKeVNDz30EFOmTOEvf/mLx5OTUrWOlVpYvTWbq8fFdMrkdLr0CcYOau3atfzud787Zf+jjz7KzJkzvRDRmbvssstO6RwSFxfH+++/32bn/MMf/sAf/vCHNqtfKXe8vOkwVrud26Z07Qdz69ME1UHNnDmzwyaixrRlIlKqvSqvquGVrzK5ZGQkg3p383Y47Ype4lNKKS96/dsjFFdaueN8bT3VpwlKKaW8pMZm5/kvDjF+YDhj48K9HU67owlKKaW85D/fHyW7qJI7zhvs7VDaJU1QSinlBSLCig3pDO7TjekJOrxWQzRBKT777DPOOuusFh2jc0wpdWY2Hcxnd04Jt0+Jx0eHNWqQJih1xp577rlmJ1DsqAnKbrdT/2H25sYaVModyzek0zs0kCvGRHs7lHar63Uz//B+yD11bLpWETkafrSk2WKbN2/m3nvvpbTUMczgE088wbp16/j888+prq6md+/evPDCCwwYMIDDhw9z1llncfPNN7NhwwYqKyv55z//yZQpU5p8DRzPSi1evBiLxUJAQABPPfWUa9ijBx98kNdee43o6GgmTJjQbMx79uzhlltuwWq1MnLkyFPmmLrnnnuYPXs2K1as4KmnniIwMBC73c4bb7zB22+/7ZpjKigoiFdffZWjR4/y4IMPYrFYqKmp4YEHHmDu3Lmu+saPH8/mzZvJycnhmmuuYckSx+81OzubRYsWceDAAQCuu+46fvvb31JSUsKvfvUrdu7cicViYerUqfz1r3/F17fhIWNyc3O57rrrKCkpwWKxMGvWLB5//HHAMcdUWloaZWVlHDx4kA0bNjBmzBgWLFjAJ598Qnx8PIsXL27w+MrKSuLj49m6dStRUVGAY6ilyMjIBp9bU13T3qMlbNh/nHsvGU6Qvw5r1Ch35+Vo7wuOOaGanw9q7X0iL8xsm2XtfQ3MgnKy/Px86devn3z55ZciIlJTUyMFBQUnzWf07LPPyrXXXisiP8yh8vLLL4uIyGeffSbR0dFisViafC0tLU0mTZokxcXFIiKya9cuiY2NFRGR999/X0aPHi2lpaVSU1Mjs2fPlnHjxjUZt84x1TZzTOl8UF3TL1/bJiP+50MpKq/2dige057ng2o/3GjhtKXNmzczcuRIzj77bMAxPUV4eDirVq3imWeeoays7JRLSAEBAcyfPx9wDGwaHBxMamoqPXr0aPS1L774goMHD3Leeee56qmpqSEvL49PP/2Ua6+91jUlxYIFC/jTn/7UaMw6x5SDzjGlWkNOUSXv78jhhskD6Bni7+1w2rWul6C8TBoYnDcjI4Nf/vKXfPvttwwaNIhNmzYxb968JutobM6n2tdEhBkzZrBy5Uq3YmhOS+aY+vbbb/nkk0+YOnUqy5Yt40c/+tEp5e666y4uu+wyVq9ejTGGYcOGtcocU40lxPrqzhEVFBTEHXfccdpzTNU/vv4cU8uXL3crJtU1vPjlIQRYcO4gb4fS7nmsk4Qx5l1jzA5jzDZjzEZjTEoDZR4xxhxzzhG13RjzjKfi85Szzz6bPXv2sHnzZsDxTTwzM5OAgAAiIyOx2+0sW7bspGOqq6td8ztt3LgRi8XC8OHDm3zt4osv5qOPPmL37t2uer799lsApk+fzhtvvEF5eTk2m40XX3yxyZh1jqmWH69zTKmGlFis/PubI8waHUVMuA5O3BxPtqBuEpFiAGPM5cALwNgGyq0UkXs8GJdH9erVi9WrV/OrX/2K8vJyfHx8ePLJJ5kzZw6jRo0iLi6O888/nw0bNriOiYiI4MCBA0ycOJGKigr+/e9/ExAQ0ORrQ4cO5ZVXXmHBggVUVlZSXV3NOeecw/jx45k9ezabN28mJSWF/v37M3Xq1CbnYwLHHFO33HILf/3rXxk3blyTc0wVFRXh4+NDbGysq3ND7RxTISEhvPrqqyxZsoSf/vSnLFmyhKSkpBbNMbVw4UJefvllfH19mTdvHvfddx9Lly7lN7/5DcnJyRhjCAwMZOnSpbXD+59i0aJFzJkzhzFjxhAbG3tac0w1dbzOMaUa8urXmZRV1XDHeTqskTu8Mh+UMeZGYJGInFVv/yNAaHMJyhgTBoTV2x0DbOxs80HV9tSrbSm4+5ryLnfnmOrIf5uqZapr7Ex5/BOG9A3lX7ed+gWvs2v380EZY54zxmQCi4GbGik21xiz0xizzhjT2LWRu4FD9ZamnxRVykN0jinVkPe2Z5NXUqXDGrWAt1pQNwDXicjMevsjgXwRsRpjLgL+BYwQkfx65bpMC8qTdI4pz9O/za5BRLhk6QZ8jOHDX0xxu9NRZ3I6LSiv9OITkVXGmBXGmIi6yUdEcuusf2yMOQIkAp/XO74IKKq7ryv+g7c2nWNKqbbxWepx9ueV8ddrkvWzqgU8conPGBNqjImts30pUOBc6paLrrOeguPh21RPxKiUUm1l+YaDRPUM4tLk/t4OpUPxVAuqG/CmMaYbYMORmC4VETHGrAUeEpEtwKPGmHHOMtXADXVbVUop1dHszCriq/QCHpg5An9fHf60JTySoEQkD2iw20rd+1Ai0ljHCaWU6pCWb0ine6AfcyfENl9YnUTTuVJKtZHM/Ao+/P4o8ybF0T1IhzVqKU1QXjBw4EASEhJISUkhISGB22+/HavVCsCyZctOGimhMY888gh9+/Z11bFgwQKqq6ubPOazzz5j3bp1ru3Dhw+fMt5ca6s7b1Ttg7wzZsygoqLC9R7GjBnDsGHDGD9+PE8//TQ2m438/HxSUlJISUlhyJAhhISEuLYfeOCBNo1Zqdby/Bfp+PoYbj1HhzU6HToWn5e89dZbJCYmYrPZmDJlCqtXr+baa6/lJz/5idt13HjjjTz55JNUVVVxwQUXsGzZMhYtWtRo+c8++4yysjIuvvji1ngLLVJdXc11112HMYb333/fNRJG7XsASE9PZ/78+aSlpfH000+zfft2V9z33HMPW7ZsabBuEcFutzc6tYZS3lBYXs0bW7K4PCWafj2Cmj9AnaLLJag/f/Nn9hXsa5O6E3olcN+E+1p0jMViwWKxEB4eDjhaRmVlZTz55JPYbDbuu+8+PvroIwBmzJjBn//851M+iAMDA5kyZQqpqak8/vjjZGZm8o9//AOAvLw8kpKSWLduHcuWLcNut7N+/Xrmzp3rmn/pgQceYO3atVRUVPD8889z7rnnAo7hjZ544gmMMQwePJjly5fTt29fXnrpJV599VXCw8PZtWsXYWFhvP3220RGRjb4HsvLy5k9ezYxMTE8++yzjSaS+Ph4XnjhBUaPHs0f//hHevbs2ejv7bnnnuOdd94hLCyMvXv3snLlSsLDw1m0aBFHjhyhsrKSG264gd/85jeA43mjX/7yl5w4cQKr1cqvf/3rk0YnV6q1rfoqg0qrTYc1OgN6ic9Lrr76alJSUoiMjGTQoEENtmpWrFjB9u3b2bp1K1u3bmXbtm2sWLHilHLFxcWsW7eOMWPGcPvtt/PWW29RVlbmqmPevHkkJyfzk5/8hBtvvJHt27dz//33A5Cfn8/kyZPZtm0bDz30EPfd50iwu3bt4v7772fdunXs3LmTxMREfv7zn7vO+e233/Lkk0+ye/duRo4cyd///vdG3+vChQuJiori+eefb7aVk5CQQEhICKmpzT9dsHHjRhYvXszWrVtJTExk/vz5/PrXv+abb77hu+++49133+XTTz/FarVy/fXX8/TTT7NlyxY2btzIH//4R7cGqFXqdFisNl7edJipw/swrF93b4fTYXW5FlRLWzhtpfYSn8Vi4aqrrmLp0qXcfffdJ5VZv349N998s+ty2C233MI777zDXXfdBThaOOvXr8fHx4fZs2dz66234uPjw2WXXcaqVau4/fbbefbZZ1m/fn2jcYSGhjJ79mzAMc/Tr3/9awA+/fRTZs6c6ZoV9s477yQ5Odl13DnnnENsbKzruLpzNNU3Y8YM1q9fz969exk5cmRLf1WNOu+881yjhpSUlPDFF1+cNDhraWkpe/fuJSIign379nHNNde4XrNarezdu5chQ4a0WjxK1Xp7axb55dU6rNEZ6nIJqr0JCgpi9uzZrFmz5pQE1dC8T3W3696/qWvRokXMmzePvn37MmLECIYNG9bo+QMDA13rdeddau7cLZmvae7cucyYMYOLLrqIjz/+uMkklZqaSkVFBQkJCY2WqVV3fia73Y6Pjw9btmzBz+/kP+sdO3YQGRnpuqelVFuy2YXnNh4iKaYnk+J7eTucDk0v8XmZ3W7n888/bzCJXHTRRbz00ktYrVasVisvv/wyF154YbN1JiYmEhERwd13383ChQtd++vPydSU6dOns3btWnJzHc9JP/vss26duzHz589nyZIlXHjhhezZs6fBMocPH2bBggXcdddd9OjRo0X1h4WFMWnSJJ544gnXvoyMDPLy8hg5ciS+vr6u+awA9uzZ47oMqlRr+nhPHodOlHPHefE6rNEZ0gTlJbX3oBITE7Hb7Tz00EOnlLnjjjtISkpizJgxjBkzhqSkJG6//Xa36r/tttvw8fFh1qxZrn1XXnklW7ZsISUlxTVPU2NGjRrFY489xkUXXURSUhI7duzgb3/7W8veZD033HDDKUlq5cqVjBkzhuHDhzNnzhyuvvpqt7rZN+S1115j+/btjB49msTERObNm0dJSQn+/v6sWbOGVatWkZSUxKhRo/j5z3/ebLd8pU7Hig0Hie0VzIxRDXcaUu7zymjmbcEYMxA4pKOZO9x2220MHz6ce++919uhqCZ0xb/NzmzL4QKuXraZ3182ipvOHujtcNqVdj8flGp7OTk5DB8+nAMHDpx0eU8p1faWb0gnLMSfOWfFeDuUTkE7SXQy/fv3d6uLtlKqdR08Xsb6vXn8fOoQQgL0o7U1dJkWVGe5lKk6D/2b7Fye25hOgK8PN+qlvVbTJRKUr6+va6w7pdoLq9V6Spd41TEdL63i7a3ZXDUuht6hgc0foNzisQRljHnXGLPDGLPNGLPROSFh/TK+xphnjDEHjTFpxpjbWuPcYWFh5OXlYbfbW6M6pc6Y3W4nLy+vyeGcVMfx8qbDWG12bp+iwxq1Jk9+fbtJRIoBjDGXAy8AY+uVuR4YAgwFIoBtxpj17vb4aEzv3r3JysrSezOqXenWrVubjyav2l55VQ2rvsrg4pH9GNS7m7fD6VQ8lqBqk5NTT6Ch5sy1wLMiYgeOG2PeBeYATzRQ1m0+Pj7ExcWdSRVKKdWgN7YcobjSqsMatQGPXgA3xjwHXAwYYEYDReKAjDrbmcAp01AaY8KAsHq7tV+nUsqjamx2nv/iEGcNCGfcgHBvh9PpeDRBichtAMaYG3C0imY2fUSj7gYebq24lFLqdKzdlUtWYSUPzW69QZDVD7zSi09EVgFTjTER9V7KBAbU2Y4DjjRQxVJgUL1lShuEqpRSDRIRVmw4SHyfblw4op+3w+mUPNKCMsaEAuEicsS5fSlQ4FzqehO43RizGkcniSuA8+rXJyJFQFG9c7RB5Eop1bDNB/PZlV3CYz8ejY+Pfv60BU9d4usGvGmM6QbYcCSmS0VEjDFrgYdEZAuwCpgIHHAe9wcRSfdQjEop5bblG9LpHRrIlWOivR1Kp+WRBCUiecCkRl6bWWfdBtzliZiUUup07cst4fP9x7nn4mEE+Tc9S7Q6fV1iJAmllGpNKzakExLgy/xJA5ovrE6bJiillGqBo8WVvL89h2vOiiUsJMDb4XRqmqCUUqoFXvzyMAIsOHeQt0Pp9DRBKaVbByWaAAAgAElEQVSUm0osVl79OpOZo6OI7RXi7XA6PU1QSinlpn9/nUlZVQ13nqeDwnqCJiillHJDdY2dF788zNmDI0iM1lHoPUETlFJKueH9HTnklli4Q1tPHqMJSimlmiEiPLshnYTI7pw/rI+3w+kyNEEppVQzPtt/nNS8Um6fEq/DqnmQJiillGrGis/TiewRxKXJ/b0dSpeiCUoppZrwfVYxm9PzufXcgQT46UemJ+lvWymlmrB8w0G6B/px3QSdldvTNEEppVQj0o6Vsfb7o8ybGEf3IH9vh9PleHRGXaWU6ggqqmtY/nk6yzccJNjfl1vO0WGNvEETlFJKOdntwttbs3ji/1I5VlrF7KQo7puRQGTPIG+H1iV5akbdCByTEQ4GqoA04E4ROV6v3EvAhcAJ5643RWSxJ2JUSnVtmw6eYPF/9rI7p4SU2DD+d/44xg0I93ZYXZqnWlACPC4inwEYY54AlgALGii7RET+4aG4lFJdXPrxMh77cB8f78kjOiyYp68bw6VJUfq8UzvgVoIyxsQCyUAYUATsEJEj7p5ERAqAz+rs+oozmDnXGBPmjKWumNOtTynV9RRVVPO3/x5g1eYMgvx9+c2M4dx6ziCdIbcdaTRBGWP8gTudSzyOy3KlQHdgiDHmELAMWCEi1e6e0BjjgyM5vd9IkV8ZY+4EDgK/FZG9DZS5G3jY3XMqpVSt6ho7q77K4On/HqDUYuXa8XH86qJh9Oke6O3QVD1NtaB2AJ/gSFBfi4it9gVjjC8wAbge2AaMasE5/w6UAQ1dxnsAOCoidmPMjcBHxpj4uud2Wgq8VG9fDLCxBXEopboQEWHdnjweW7uXw/kVTBnamwdmjSAhsoe3Q1ONaCpBXSAixxp6wZkwNgObjTFuj5xojHkSGApcKiL2BurNrrO+0hjzFI7Ek1GvXBGOS41163Y3DKVUF7Mru5g/rtnD14cKGNI3lBdvGc8Fw/ro50Y712iCaiw5ARhjggGbiFTX74nXxDGLgXHALBGpaqRMdG2SMsZcAtiA7IbKKqVUc3KLLTzxf6ms3pZFeEgAf7wikevGx+Lnq2MUdATudpJ4EnhDRL4xxswC3gLEGHOtiHzgxvGjgN8B+4FNzm8th0TkSmPMdmCmiOQALxtj+gF2oAS4TERqTuudKaW6rNoHbVdsSMdmF+44L56FU4fQQ0eD6FDc7WZ+PfCQc/0hYD5QDDwFNJugRGQ30GBbWkRS6qxf6GY8Sil1itoHbZ9cl0peSRWzkqK4f0YCsb1CvB2aOg3uJqgQEalwPnAbLyJvAxhjBrRdaEop5b7NB/P503/2sDunhOTYMP55/VjGDejl7bDUGXA3Qe03xlwPDAE+BjDG9AYq2yowpZRyx6ET5Ty2di/rnA/a/m1uCpcm9cfHRztAdHTuJqifAn8Dqvlh9IdLgHVtEZRSSjWnqKKap/+bxsrNhwn08+HeS4az4Fx90LYzcStBici3wNn19v0L+FdbBKWUUo2prrHzylcZ/E0ftO30mhpJIllEdjRXgbvllFLqTIgIH+/J47EP93HoRLk+aNsFNNWCesYYU4JjFPLPnd3AATDGRAHnAzcCocB5bRqlUqpLs9uFRa9tY83Oo/qgbRfS1IO65xpjZgM/AZ43xtj4YSw+A6wH/iEiaz0SqVKqy/rHp2ms2XmUuy8cys+mDtEHbbuIJu9BicgaYI1z4NihOEYQLwQO6AO0SilP+DT1GE+t38+VY6L5xfSh2mrqQtztJGEF9rRxLEopdZLM/Arufm07CZE9ePTK0ZqcuhhtJyul2qXKahs/eeU7RIRl88cSHKDdx7saT82oq5RSbhMRHnj3e/bmlvDCTeMZENHN2yEpL9AWlFKq3XnlqwxWb83mF9OHMjWhr7fDUV7SogRljIk1xkxqq2CUUuq7jEL+sGYP0xL6smjaUG+Ho7zIrQRljIkzxnwJ7MPRvRxjzNXGmOfaMjilVNdyrNTCT//1Hf3DgnnqmhQdT6+Lc7cFtRz4D45noKzOfR8DF7lzsDEmwhiz1hiTaozZaYxZ3dBMvMaYEGPM68aYNGPMPudzWEqpLsBqs/OzV7dRXGll2fxx9AzRuZu6OncT1ARgiXOadgEQkWKgp5vHC/C4iAwXkSTgILCkgXL3AKUiMgS4FHjOGBPq5jmUUh3YY2v38c2hAv58VRIjonT4IuV+gsrDMdWGizFmJJDpzsEiUiAin9XZ9RXQ0FxS1wLLnMccALYAP6pfyBgTZowZWHcBYtyJRSnV/ry3PZsXvjzEzWcP5PKUaG+Ho9oJd7uZP4ljRInHAD9jzHU4pnBvqBXUJGOMD3AX8H4DL8cBGXW2M4HYBsrdDTzc0nMrpdqffbkl3P/294wfGM4Ds0Z4OxzVjrg7ksQLxpgC4A7gCHAT8D8i8u5pnPPvQBnwj9M4ttZS4KV6+2KAjWdQp1LKw4orrfxk1XeEBvnxzLyx+OsYe6oOtx/UdSaj00lILsaYJ3GM6Xep835WfZk4Lv0dd27HAZ82EEsRUFSv7jMJTSnlYXa78Os3tpNVWMlrd0yib48gb4ek2hm3E5QxZgowBsf0Gi4i8qibxy8GxgGzRKSqkWJvAncCW4wxQ4HxwHXuxqiU6jie+TSN9XuP8cilIzlrYC9vh6PaIbcSlDHm78A1OC6hVdZ5Sdw8fhSOe1b7gU3O1s4hEbnSGLMdmOmcb+oJ4CVjTBpgA+4QkVJ334xSqmP4LPUYf12/nytS+nPT2QO9HY5qp9xtQV0PJNadtLAlRGQ3jjmkGnotpc56OTDndM6hlOoYjhRU8IvXtjO8X3ce+3GSXp5XjXL3juQRoLHLckop5ZbKaht3rnKMUL78hnE6QrlqkrstqAXAs8aYf+N4JspFRDa0elRKqU6ndoTyPUdLeOHms3SEctUsdxPUOBwPzJ7Hqfeg4lo7KKVU51N3hPJpCf28HY7qANxNUI/i6Bq+vi2DUUp1TrUjlE8d3odfTNcRypV73L0HVQ7opTylVIvVjlAe1TOYpdeO0RHKldvcTVAPAUuNMZHGGJ+6S1sGp5Tq2OqOUL78Bh2hXLWMu5f4XnD+vLPOPoPjHpR2w1FKNWjJh44Rypdem6IjlKsWczdBDWrTKJRSnc77O3J4/gvHCOVXjNERylXLuTtYbEbzpZRSyiE1t5T73trJWQPC+d1MHaFcnZ5GE5QxZoWI3OFcX0UjwxqJyI1tFJtSqgMqsVj5ySuOEcr/ef1YAvz0VrU6PU21oA7VWU9r60CUUh2f3S786vUdHCmo4N86Qrk6Q40mKBF5zBhznYj8W0R+78mglFIdk2OE8jwevnQk43WEcnWGmmt7L/dIFEqpDq92hPLLU/pzs45QrlpBcwlKn6hTSjXr5BHKR+sI5apVNNeLz9cYM5UmEpWIfNK6ISmlOhKL1TFCuV2EZfPHERLg9jyoSjWpub+kQOB5Gk9QAsS7cyLndO9XAQOB0SKyq4EyjwA/BWrnnfpSRBa6U79SyvNEhAfe2cWeoyU8f9NZDOytI5Sr1tNcgioXEbcSkBveBf6GY1bepqwUkXta6ZxKqTZitwvPbkzn7a1Z/GL6UKaP0BHKVevyWFtcRL4AWuXatDEmDAirtzvmjCtWSjXLZhfW7MzhH5+kceBYGReO6KsjlKs20VyC8sadzrnGmIuBXOBhEdncQJm7gYc9G5ZSXVuNzc5723N45tM00k+UM6xfKE9fN4ZZo6N0hHLVJppMUCLS3VOBOC0DFouI1RhzEfCeMWaEiOTXK7cUeKnevhiav3yolGqh6ho772zL4plPD5JZUEFCZHf+9/qxXDIqUhOTalPtqruNiOTWWf/YGHMESAQ+r1euCCiqu0+7tSrVuqpqbLz1XRb//PQg2UWVJEb3YMUN47hwRD9NTMoj2lWCMsZEi0i2cz0FR4+/VK8GpVQXY7HaeP3bIyz7/CBHiy2kxIbxpysSuWB4H/0iqDzKYwnKGPM08GMgElhvjMkXkVHGmLXAQyKyBXjUGDMOsAHVwA11W1VKqbZTWW3j1W8yWf75QY6VVnHWgHAevzqJc4f01sSkvMKTvfgWAYsa2D+zzvpNnopHKeVQXlXDK19l8OzGdE6UVTMpvhdL56YwOT5CE5PyqnZ1iU8p5TmlFisrN2fw3MZ0CiusTBnam59PG8qEQTrIq2ofNEEp1cUUV1p5edNhnv/iEMWVVi4Y3oefTxvKuAHh3g5NqZNoglKqiyiqqOaFLw7x4peHKa2q4cIR/Vg0fQhJMfWfeVeqfdAEpVQnl19WxXNfHGLlpsOUV9uYMSqSn00bQmJ0T2+HplSTNEEp1UmdKKtixYZ0Vm3OwFJjY9boKH42bQgJkT28HZpSbtEEpVQnY7XZeXnTYZauP0BFdQ2Xp0SzcOpghvT19MAwSp0ZTVBKdSJfpp3g4fd3k3asjAuG9+F/Zo9kcJ9Qb4el1GnRBKVUJ5BVWMHi/+zlw125xPUK4bkbz2L6iL76HJPq0DRBKdWBWaw2VmxI55+fpQHw64uGcft58QT5+3o5MqXOnCYopTogEWH93mP8Yc1ujhRUMmt0FL+bNYLosGBvh6ZUq9EEpVQHk368jN9/sIfP9x9naN9QXr1tImcP6e3tsJRqdZqglOogyqpq+PsnB3jhi0ME+fnyP7NHcuPkAfj7+ng7NKXahCYopdo5EeH9HTk8unYveSVVXD0uhvtmJNCne6C3Q1OqTWmCUqod25NTwiPv7+abwwWMju7J/84fx9g4HTNPdQ0eSVDGmCeBq3BMQDhaRHY1UMYXeBqYAQiwRESe80R8SrU3RRXV/PXj/bzyVQY9g/157MejueasWHx1JlvVhXiqBfUu8DdgYxNlrgeGAEOBCGCbMWa9iBxu+/CUah9sduGNLUd4/KN9FFdamT9pAL+6aBhhIQHeDk0pj/NIghKRL4DmHhq8FnhWROzAcWPMu8Ac4In6BY0xYUD9IZhjWidapbxja2YhD7+3m++zi5kwsBePXDaKkf113DzVdbWne1BxQEad7UwgtpGydwMPt3lESnnA8dIq/vzRPt76Lot+PQL529wULkvur6NAqC6vPSWollgKvFRvXwxNX0JUql2pqrHx8qbD/P2/aVhqbPzk/MH8bNoQQgM76n9LpVpXe/qfkAkMAL51btdvUbmISBFQVHeffttUHYWIsPb7XJZ8tJcjBZU6qKtSjWhPCepN4HZjzGocnSSuAM7zbkhKta6tmYUs/s9evssoJCGyO6sWTGDK0D7eDkupdslT3cyfBn4MRALrjTH5IjLKGLMWeEhEtgCrgInAAedhfxCRdE/Ep1RbO1JQwZ8/2seanUfp0z2QP181mqvHabdxpZriqV58i4BFDeyfWWfdBtzliXiU8pQSi5VnPk3jxS8O4+MDi6YN4c7zB9NN7zOpZpRUl7Aldwvl1nKGhg8lvmc8Ab5d63ED/V+iVBuw2uz8+5tMlq4/QGFFNT8eE8M9lwwjqqeONq4aZqmxsPXYVr45+g1fH/2aPQV7sIvd9bqv8SWuRxxDw4YyJHwIw8KGMSR8CDGhMfj6dM7pVTRBKdWKRIT/7j3Gox/uJf14OZPjI3hg1ggSo3t6OzTVztTYa9h1YhdfH/2ar3O/Zvux7VjtVvyMH6P7jOaOpDuYGDmR8KBwDhQdIK0wjQOFB9hXsI+PMz5GEACCfIOID4tnaNhQhoYPdSWwPsF9OnznMU1QSrWSXdnFLP7PXjan5xPfp5vOaqtOIiIcKDrgSEhHv2ZLnuPyHUBCrwTmJcxjYtRExvUbR4h/yEnHDg4b7BgozqnCWsGh4kPsL9xPWpEjcW3K2cR7B99zlekZ2JMhYUN+SFzhQxkcNpgeAR3n4W9NUEqdodxiC0+uS+XtrVmEhwTwh8tHcd2EOJ0GQ3Gk9IgrIX2T+w0FlgIA4rrHMXPQTCZGTWRC5ATCg1o2AHCIfwijeo9iVO9RJ+0vtBS6ElbtzzXpayizlrnKRHaLdCQuZ2traPhQBvUcRKBv+xsdXxOUUqepvKqG5RvSeXZDOja7cMd58SycOoQeQf7eDk15yYnKE457SLmOpJRdlg1An+A+nN3/bCZGTWRi5ESiQqPa5PzhQeGMjxzP+Mjxrn0iQm55LgeKDpyUuL4++jVWuxX44f5W/cTl7ftbmqCUaiGbXXjruyP8Zd1+jpVWMTspivtmJBDbK6T5g1WnUlpdypbcLa6ElFaUBkD3gO6M7zeeG0feyKSoSQzqOchrl3qNMUSFRhEVGsV5MT88WlpjryGzJPOkxJVakMr6jPWn3N8aEjaEYeHDmDdiHv4+nvsCZkTEYydrS8aYgcChQ4cOMXDgQO8GozolEeHLtHz+9J897MstZWxcGA/OHqnzM3nA0bKjvH3gbXLKcujm3821hPiH0M2/G6H+oa71bn7O1wO6EeIXgp9P630Pr7JVsf3Ydtdlu935u7GJjUDfQMb2HcuEqAlMiprEiF4jOmzPusqaStKL0k9pcVXZqvhi7hennWgPHz7MoEGDAAa5O0uFtqCUasbR4kre3ZbDO9uy2J9XRmyvYJ6ZN5aZoyO1A0QbsoudzTmbeS31NTZkbQAgMiSSipoKyqxl1Nhr3KonyDeo8UQW8ENCa7CMfzcsNRa25G3hq6Nfsf3YdqpsVfgaXxJ7J7Jg9AImRU0iuU9yp3lGKdgvuMH7W+XWco//vWuCUqoBZVU1fPj9Ud7Zls3m9HxE4KwB4Tx65WiuGhdNoF/H/HbcERRXFfNu2ru8nvo6R0qP0CuoFwsSFzBn2JyT7t1U26opt5afutSUU2GtoKy6zLVev8zxiuNk1GRQVl1GRU0FlTWVzcY1NHwoc4bNYVLUJMb1G0doQNcaO7GbfzePn1MTlFJONTY7X6Sd4J1t2fzf7lwsVjsDIkL4xfShXDkmmgERnv8P2pXsOrGL1/a9xkeHP6LKVsXYvmP5WcrPuHDAhQ22TgJ8AwjwDWhxD7iG2Ow2KmpOTWQV1goAUvqmEBEcccbnUS2jCUp1aSLC7pwS3tmWzXvbczhRVkXPYH+uHhfDlWNiGBsXppfx2pClxsKHhz7k9dTX2Z2/m2C/YC4ffDnXDL+G4b2GeywOXx9fugd0p3tAd4+dUzVPE5Tqko4WV/Le9hze2ZpNal4p/r6GaQl9uXJMDFMT+uglvDaWUZLBG6lv8G7au5RUlzC452B+N/F3XBp/aZe7dKYapwlKdRllVTV8tCuXd7Zlsemg477SuAHh/OmKRGYnRREW0jlucrdXNfYaNmRt4PXU19mUswk/48f0AdO5dvi1nNXvLG2pqlNoglKdms0ujvtKW7P4v915VFptxPUKYdE0x32lgb31vlJbO1F5gtUHVvPm/jfJLc+lb0hfFqYs5KqhV9EnROfCUo3zWIIyxgwDXsYxGWE+cKOIHKhX5hHgp0COc9eXIrLQUzGqzuNIQQVvfpfFW1uOkFNsoWewPz8eG82Px0YzNi5cv623MRFh67GtvL7vdT7O/Jgaew0ToyZy3/j7uCD2glZ9Nkl1Xp78K1kGPCMirxhj5gPLgWkNlFspIvd4MC7VSVTV2Ph4Tx6vf3uEL9JOADBlaB8emDWSC0f21ftKHlBuLec/6f/htdTXOFB4gO7+3Zk7fC5zhs8hvme8t8NTHYynZtTtC4wFLnLu+jfwD2NMHxE5fhr1hQFh9XbHnFmUqqPan1fK698eYfXWLAorrESHBfOL6UOZc1Ys0WE6/5InpBWm8Xrq63yQ/gHl1nISeiXwyORH+NGgH50yMrdS7vJUCyoWyHbOmouI2IwxOc799RPUXGPMxUAu8LCIbG6gvruBh9syYNW+lVfV8J+dR3nt20y2Zhbh72u4aGQ/rh0fx7lDeutU6m2spLqE9KJ09hfu58NDH7Ilbwv+Pv7MGDiDaxOuJal3kl5GVWesvV0IXgYsFhGrMeYi4D1jzAgRya9XbinwUr19McBGD8SovERE2H6kiNe/PcIHO3Ior7YxpG8oD84awZVjookIbX/TBXR0xVXFpBWlcbDoIOnF6Y6fRekcqzzmKhMdGs0vx/2SK4ZcQa+gXl6MVnU2nkpQR4BoY4yvs/XkC/R37ncRkdw66x8bY44AicDn9coVAUV19+m3tc6roLyad7Zl88a3R0jNKyXY35fZSVHMnRCrHR5aSYGlgINFB11LbTLKt/zw3TDYL5j4nvFM6j+J+J7xDA4bzOCwwUSHRuNjdO4r1fo8kqBE5JgxZjtwHfCK8+e2+vefjDHRIpLtXE/BMYdkqidiVO2LiLA5PZ9Xv85k3e48qm12kmPDePTK0VyaHEV3nXOpxUSEfEv+SYnoYLGjRVRYVegq182/G4N7DmZKzBQG9xxMfJgjGUV1i9JEpDzKk5f4fgK8bIx5CCgEbgQwxqwFHhKRLcCjxphxgA2oBm6o26pSnV9xpZW3v8viX19ncPB4OT2D/Zk3MY5rx8cyIqrjTFXtLRXWCgosBeRb8imoLCC7LJu0ojRXi6ikusRVtrt/dwaHDWZa3LSTWkT9Qvppq1S1Cx5LUCKyD5jYwP6ZddZv8lQ8qn3ZmVXEK19l8P6OHCxWOymxYTw5J5nZSVEE+Xfd7uE2u43CqkIKLAWOpdKZfBrZbmhU7p6BPRncczCXDLyEwWGDXcmoT3AfTUSqXWtvnSRUF1JZbeODHTm88nUGO7OKCfb35cox0Vw/cQCJ0T29HV6bEBHKreU/JJh6S35l/knbhZZC1+ymdfkaX3oF9XItA3oMOGk7IjiCiKAI+nXrR0RQhCYi1SFpglIel3aslFe+yuTtrVmUWmoY2jeU3182iivHRtOjA95bqr2sVptQ6ied+vusdmuD9YT6hxIRHOFKOGP6jnElG1fiCXKs9wjsofeDVKenCUp5RHWNnXV7cnnlqwy+Si/A39cwIzGK+RPjmDCoV7v7hl9lq+J4xXGOVx4/6eeJyhOnJB6LzdJgHcF+wfQK6kV4YDh9QvowLHwYvYJ70Suwl+NnUC/Cg8Jd24G+2k1eqbo0Qak2Y7MLqbmlrP3+KK99e4QTZVVEhwVz7yXDueasWPp09/wHcmVNJScqTnCs8hjHK4+71uv+PF55/KTOBLX8fPyICIogIjiC8KBwBocNJjwwnF7BjiRU29IJDwonPDBcR1BQ6gxpglKtpqiimm2ZRWzNLGRrZiHbM4sor7ZhDEwb3pf5kwZw3rA+bTbKQ2036qzSLLLLsskuyyarNIucshxX8im1lp5ynL+PP32C+9AnpA+Deg5ifOR4+ob0pXdwb/qE9HG9FhYYppfVlPIgTVDqtNjtwoFjZWzNLOS7DEdCSj9eDoCvjyEhsjs/HhvDuAHhTIzvRVTP1hkTr7S61JF8SrPJKjs1EdW/3NY7uDfRodEMCRvC5KjJJyWcPsGOpWdgz3Z3iVEppQlKuamy2sY3hwv4LqOQbc7WUWlVDQDhIf6MGxDOVWNjGBsXTnJsT0ICTu9Pyy52jlUcI7Mkk8xSx1K3RVRcVXxS+VD/UKJDoxnUcxDnRp9LdGg0Md1jiAmNoX9of4L8gs74vSulvEMTlGpUXomF/+49xn/35vHlwRNYrHZ8DAyP7MFlKf0ZGxfO2AHhDIwIaVELpG4SyijN4EjJETJKMlzJqG4ryN/Hn/6h/YkJjSExIpHo7tEnJaEeAT209aNUJ6UJSrmICLtzShxJaV8eO7McrZWY8GDmjo9jakJfxg0IJzSw6T+bKlsVhZZC12gGuRW5jhaRs1V0pPQIVbYqV/kAnwBiu8cS2yOWs/ufzYAeA4jtHsuAHgPoF9IPX5+u+6CuUl2ZJqguyC52MksyOVp+FKtN2He0jK2ZRWzNKOZEmRWDYVi/Htw8tTcTB0UQ37s7PsYHuxSwuyC90YdKa7tel1nLTjlnbRKK6xHHOf3PIa5HHHE94hjQfQB9Q/pqElJKnUITVCdXbi1nf+F+UgtSSS1MZX/BfvYX7m/42Z3e0K23YzULyMqFt5sYCdHX+Dqe43E+RDq69+hTHiztFdSLPsF96Netn/aAU0q1iCaoTkJEyC3PZW/BXlciSi1M5UjpDzOadPPrTgixVBeNp7K0L+H+/TlrYDhj4noyIqo7fr6O1pVd7AiCzW7Djh0RwSY2DIbwoHAdzUAp5RGaoDogm91GRkkGewv2sq9gn+tnbQ83gyG2eywJvRKYNegySop7801qENsOgZ+PDxeO6Mfci2KZMrTtnklSSqkzpQmqHSiuKuZI6RHHOG02K1W2Kqrt1VTbnIu9mipbFfmV+ewt2MuBwgOuUav9ffwZGj6UC+MuJKFXAgm9EhgWPoyMEzW89k0myz/LpsRSw8CIEO7/URxXjY3xyggOSnUodhvUVEGNBWxWQEDsIHLyutid21Lntfr7G1q34xgDuJF6TzqHzRGDzQp2axPr1WCraYX1Rur39Yf7Mz36z6AJygPsYud4xXGOlh8luyz7h2d8nD+LqoqarwTHMz/Dew3nqqFXuZJRXPeBlFqEgvJq8suq+T69nIe/28aOI0UE+Pnwo8RI5o6PY1J8+xvvTrVjtR+mTS5y6nr9D+/6H+ANlaF+PY0lAbvzA7PamTyqwOZMIjXVzvXa/dXO/XXX65VxHV933VlWbN787Z8ZHz/w8QffAPBtaN251K77B0Ngj1P3+/g5j3Ou+7fOw/Yt4bEEZYwZBrwMRAD5wI0icqBeGV/gaWAGju8XS0TkOU/F2FIV1gpyK3I5VnGMoqoiSqpKKK4qpriqmMKqQnLKcjhafpS8ijxq7DWu4wyGyG6RxHWP46IBFxHXPY7YHrFEBPbGWuOLxWqoqDJUVhnKLFBuMZRahOIKG/n5VXyXUc3H5dXkl2VQVJnm+P9cx7B+oTw0eyRXjokmvFuAh38rbmr0A7D+/rofVo2UqfvttNF66ix2m3PdVmdbftg+pYycesxJ5eqfq36ZhuKoLVNbt407Sz8AAAy6SURBVA3sNc7F5lyc21Jvu8HXG9hvtzVw3jrvqank05H5+IFfkOPD1S/QsfgGgl+Ac38gBPVopEwD6z5+4OMLGDA+YMzJ68bHuV27Tr39zR1j6h1f7xgf35MTh2+AM3nUW699vRN9EfVkC2oZ8IyIvGKMmQ8sB6bVK3M9MAQYiiORbTPGrBeRw54I0GarocRSQH7lCQoq8ymyFJFXeZzc8jxOVJ6gpLqMsupySqwlHLecoLTm1HHdAAJMAN18utHTN5wIn74MDBhKqHSnm4QSaAsloDoU6wmw5FRTVlXD1uoiNlmPU1VtxRcbftjwxf7DT2PDDzuhATA8wIdxgb70CDB0i/QhNMDQLcCHUH9DsL9jOzzYF1O9Eb6yn/ohVfeD1vWBZnV8oP1/e+caJEdVxfHff2Y2a5INCSEkkBeUQOQZFcPToFJItJAIVYgYAbGEVMEHkCpBLURCeJiS4KOEwg+CpCq8kYfIQwlFBREBAyQgKI8IgTw0JMQVQkFIdo8f7p1sb2d2M7uzPdM7Ob+qrul7T9/uc/pu39O3u/ecjs1dg1tyfWu5o6tdt9/E/rbuO12Xcig7JOWBJy6FYvf1QiksKibKqV8ltiu1QmFYol2ha73cRgptksfttqi7HpXkPcpSA3OlAbbbduW6CoN1pXbJAbxcLrakHEl0OqUhXfX+LwtNQ10clKSxwMHAsbHqVuBaSbua2brEpqcAvzGzTmCdpHuBk4H5qf2NAkalDjOxFh0XPHAFP1t/e0VZyYwxHR2M6OxkRGcnEzuNaVu2sPuWDnbbsoVxHR3s3NHJyM4OdurspLW/4+/2JjsGbIrLdkkOPMXEgKju5UJL4g6tVKFcvhttSQyCxa5BdOtvofvxtsoK3XXYZtBLD2o9DZA9bZcqQw/7KOuY0qfb4FxMDdiVtqlQLrfrdZvyHbHjONVSrxnUJGC1WXiwa2YdktbE+qSDmgy8mSi/FbdJcz4wZyAV3HfSocxcu4ShtNKmoYygleFqZaSG0VYYSqlUoFAsUSoWKRaLqDSEQssQiqVWiqUWKIiNMj4QlAohYOqQlhItxRJDhrTQUiqhrQN8L3ek3e6Ak3fQpa5Bf5s730oDqQ+GjuMMbgbrRxK/BBak6iYCj/d3h4dPncHhU2fUopPjOI4zgNTLQa0EJkgqxtlTERgf65O8BewBLInl9IwKADNrB7p9+uZfqDmO4zQXdQkDYGZvA8uAWbFqFrA09f4J4E5gtqSCpF2BE4G76qGj4ziOky/qGafmbOBcSa8C58Yykh6UNC1usxB4HXgNeAq4zMxer6OOjuM4Tk6o2zsoM3sZOKxC/XGJ9Q7gnHrp5DiO4+QXj/TpOI7j5BJ3UI7jOE4uGayfmVeiCLBq1apG6+E4juOkSIzNVYf6kKUDuQ1SJE2nhv+DchzHcerCUWb2l2o2bCYH1QocAvwb6E8o4vI/+h5FSCg7mGkWW5rFDmgeW5rFDnBb6k0R2B1YYmZVBWxrmkd80eCqvHIlEv/ou6pewWmzollsaRY7oHlsaRY7wG1pEP/qy8b+kYTjOI6TS9xBOY7jOLnEHZTjOI6TS9xBddEOzCUVhHaQ0iy2NIsd0Dy2NIsd4Lbknqb5is9xHMdpLnwG5TiO4+QSd1CO4zhOLnEH5TiO4+QSd1CApCmSnpT0avzdp9E6pZG0QtLLkpbF5Uux/nBJz0fdH5Y0NtGmX7IB1vtqSW9IMkkHJup7POdZyDK2pWLfRFnu+kfSLjEP2yuSXpB0d0wQmom+WdmyHTss1pX75KBEu5mxv5ZLul3SsFplA2TPvfE8LZX0uKRPxfpBd60MGGa2wy/Ao8Bpcf004NFG61RBxxXAgak6AcuB6bF8MfDbWmQZ6D0dmJTWv7dznoUsY1u26Zs89w8wGvhCojwfuCELfbO0pSc74roBbRXatAH/AfaJ5euBS2qRDWC/jEysnwA8N1ivlQE7J41WoNELMJbwaWYxlouxvGujdUvpuc0gSIg9+GKiPAbYWIusHvr3ds6zkGXdF5X6ZjD1D3AS8EgW+tbTlrIdcb0nB3UycH+iPA14qRZZRrZ8C3hmsF8rtS5NE4uvBiYBqy1k88XMOiStifXrGqrZttwsSYSYgxcBk4E3y0IzWy+pIGl0f2VmtqEOdvR2zpWBrB792K1vzKydQdA/kgqELNb3ZaFvvWxJ2VFmsaQS8BBwqYV4nd30Ad4i/I1Qg2zAkHQ9MIPw9/xlmvNaqRp/BzV4OMrMPkm4IxVwbYP1cboYzH1zDbCRwaVzJdJ2TDazacDngP2BHzdKsb5gZmeZ2WTCDej8RuvTaNxBwUpggqQiQPwdH+tzg5mtjL+bgOuAzxLu4vYobyNpTNjENtQgqwe9nfMsZJnSQ99AzvtH0tXAPsApZtaZkb6Z21LBjmSfvEt4X1SxTwgzo5U1ygYcM1sIHE1IndE010pf2eEdlJm9DSwDZsWqWcBSM8vNVFfScEkj47qAbxB0fhYYqpCsEeBs4I643l9Z5vR2zrOQZWlLL30DOe4fSVcCnwFOtK7cPFnom6ktleyQtLOkoXG9BHyNrj75I3BI4qu1pD79lQ2EHW2SJiXKM4ENQNNcK/2i0S/B8rAA+wJPA6/G3080WqeUfh8HlgIvAC8BdwK7R9mRwN+B14BFwLhEu37JBlj3XxHuArcQvoIqv3Tu8ZxnIcvKlt76Jq/9AxxA+IjgFcJAtQy4Jyt9s7KlJzuAI2J/PA/8gzCDaku0OyG2WR77a3itsgGwZRzwVDxPywhf2R08WK+VgVo8Fp/jOI6TS3b4R3yO4zhOPnEH5TiO4+QSd1CO4zhOLnEH5TiO4+QSd1CO4zhOLnEH5Tg1IGmBpCsadGxJulHSfyX9rRE6OE6WuINymgqF1BdrJQ1P1J0laXED1cqK6cCxwEQzO7SvjeO5+uLAq+U4A4M7KKcZKQHfbbQSfaUceqYP7AGsMLP3s9DHcRqNOyinGZkPXCBpVFogaU+FZHalRN1iSWfF9W9LekLSLyS1S3pd0pGxfqWktyWdkdrtGEmLJL0n6TFJydhz+0bZBoXEel9PyBZI+rVC0r33CbHX0vqOl3RfbL9c0uxYfyYhQsIRkjZKmluh7V6SHpX0jqT1km4unxNJCwnx5P4Q238/1n9V0kvR9sWS9kvsb4WkCxUSAb4v6QZJ4yQ9FG1/RNLOcduPSbopHrtd0hJJ46roO8fZijsopxl5BlgMXNDP9ocRQuXsAtwC3EaIVL43IbnbtZLaEtufClxOyHW0DLgZQpw+QmifWwg5eGYB10k6INH2m8CVwAhCqo40txLCK40nxJT7iaRjzOwGQjy4J82szczmVGgrYF5sux8hncKlAGZ2OiEA6szY/ipJU+LxzifkDXqQ4MCGJPZ5EuGx4hRgJiGVxUXR9gJwXtzuDGBkPOYuUdcPKujoOD3iDsppVi4BzlVMAd5H3jCzGy3ky7mdMMheZmabzOxh4COCsyrzgJn92UKw0h8RZjWTgOMJj+BuNLMtZvYccBfB0ZT5vZk9YWadZvZhUom4j+nAD8zsQzNbRpg1nV6NEWa23MwWRb3XAT8HPt9Lk1OiLYvMbDNwNTCUEEuvzDVmttbMVgOPA0+b2dJo+z3Ap+N2mwmOaW8z6zCzZy1EFnecqnEH5TQlZvYicD/ww340X5tY/yDuL12XnEFtTVNgZhsJUajHE94RHRYfcbVLaifMtnar1LYC44ENZvZeou5NYEI1RkgaK+k2SaslvQvcRJjp9Ha8ZHLBzqhf8njp89DTeVkI/Am4TdIaSVdJaqlGb8cp4w7KaWbmALPpPsCWPygYlqhLOoz+kEyT0AaMBtYQBvfHzGxUYmkzs3MSbXuL1rwGGC1pRKJuMrC6Sr3mxf1PNbOdCI8n1cux19A9d5OibdUer2vHZpvNbK6Z7U+YgR1PSGPuOFXjDsppWsxsOeER3XmJunWEAfc0SUVJ3wH2qvFQx0maHt/VXE547LWSMIObIul0SS1xOST54cF29F8J/BWYFz86mAqcSXzHVQUjCFlm2yVNAC5MydcS0oWUuQP4iqRj4mzne8CmqEOfkHS0pIPil4nvEh75dfR1P86OjTsop9m5DBieqptNGKzfIeQU6vMAnOIWwmxtAyF53qkA8dHcDEISwzWEHFI/BVr7sO9ZwJ6x/T3AHDNbVGXbucDBwP+AB4C7U/J5wMXx8eMFZvYKYZZ1DbCe8BHETDP7qA/6ltkN+B3BOf0TeIzwiNFxqsbzQTmO4zi5xGdQjuM4Ti5xB+U4juPkEndQjuM4Ti5xB+U4juPkEndQjuM4Ti5xB+U4juPkEndQjuM4Ti5xB+U4juPkkv8DTuqK36jX92EAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"ax.plot(natoms, [t.average for t in brute], label='self_distance_array')\n",
"ax.plot(natoms, [t.average for t in capped], label='capped_distance_array')\n",
"ax.plot(natoms, [t.average for t in kdt], label='BioPython KDTree')\n",
"\n",
"ax.legend()\n",
"\n",
"ax.set_title('Finding all pairs within 2 $\\AA$')\n",
"ax.set_xlabel('Number of atoms')\n",
"ax.set_ylabel('Time (s)')\n",
"\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### RDF test\n",
"\n",
"Here we test finding all pairwise distances up to 12.0 A between two sets of coordinates. This is a proxy for a radial distribution function, where we want to know the probability distribution of pairwise distances between two chemical species.\n",
"\n",
"\n",
"The particle density is now reduced by a factor of 20, and boxsizes vary from 50 to 175 A cubes."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"def KDTree_rdf(coords1, coords2, max_cutoff):\n",
" kdt = KDTree.KDTree(3)\n",
" kdt.set_coords(coords1)\n",
" \n",
" idx, dists = [], []\n",
" for pos in coords2:\n",
" kdt.search(pos, max_cutoff)\n",
" idx.append(kdt.get_indices())\n",
" dists.append(kdt.get_radii())\n",
" \n",
" return idx, dists"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"rdf_cutoff = 12.0\n",
"\n",
"natoms_rdf = []\n",
"brute_rdf = []\n",
"capped_rdf = []\n",
"kdt_rdf = []\n",
"\n",
"for boxsize in np.linspace(50, 175, 26):\n",
" c = make_coords(boxsize)\n",
" s1, s2 = c[::20], c[5::20]\n",
" \n",
" t_brute = %timeit -o -q distance_array(s1, s2)\n",
" t_capped = %timeit -o -q capped_distance(s1, s2, max_cutoff=rdf_cutoff)\n",
" t_kdt = %timeit -o -q KDTree_rdf(s1, s2, max_cutoff=rdf_cutoff)\n",
" \n",
" brute_rdf.append(t_brute)\n",
" capped_rdf.append(t_capped)\n",
" kdt_rdf.append(t_kdt)\n",
" \n",
" natoms_rdf.append(s1.shape[0])"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1a20735940>"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEcCAYAAAA/aDgKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4FNX6wPHvSe9ACKEl9BJ6DD1U6U1+KlgQRBAb6MUCWK6oeFEuIgr3KkqxIE1UUC9NUHqkS5EWAgEiKYSS3jfZPb8/dokhJJAA2U15P88zz+7OnJl5d7OZd8+cmXOU1hohhBCiOOxsHYAQQoiyR5KHEEKIYpPkIYQQotgkeQghhCg2SR5CCCGKTZKHEEKIYpPkIUqMUmqxUuo9y/OeSqmoYqy7XSn1lOX5SKXUryUVpxCi+CR5iDtmOdAnKKWcS2L7WuvlWut+RYgjN1lVNEqpaUqpZXew/r1KqW1KqSSlVES+Zb5KqW+VUjGW5buUUh1vsi2llPpAKRVnmWYppdQt9l9fKWVSSn12u+9BWJckD3FHlFL1gG6ABobaNBhxJ9KAr4ApBSzzAA4AbQFv4BtgvVLKo5BtPQPcD7QBWgNDgGdvsf/RQALwaEn9CBF3lyQPcadGA3uBxcATt7sRpVRfpdQpyy/bTwGVZ9kYpdTvludKKTVHKXXZUvaoUqqlUuoZYCTwqlIqVSm11lL+daXUWaVUilLqpFLqgfzbVUrNttScziulBuZZ7q2U+tryiztBKfVznmVDlFJHlFKJSqndSqnWeZa9ppSKtuwzTCnVu5D3nHtqLv/7tLzWSqmJSqlzSqmrSqkPlVI3/M8qpQYA/wQesbz3Py3zayml1iil4pVS4Uqppwv7/LXW+7XWS4FzBSw7p7X+WGt9UWtt1FovBJyApoVs7gngI611lNY6GvgIGFPYvi1GA1OBbOC+W5QVpYAkD3GnRgPLLVN/pVT14m5AKeUDrMZ88PABzgJdCineD+gONAEqA48AcZYD2nJgltbaQ2t97QB0FnPNqBLwLrBMKVUzz/Y6AmGW/c4CvsxzimUp4Aa0AHyBOZZ4gzD/Sn8WqAosANYopZyVUk2BF4D2WmtPoD8QUdzPJI8HgHZAEPB/wJP5C2itNwIzgO8s772NZdG3QBRQCxgOzCgskRWHUioQc/IIL6RIC+DPPK//tMwrbHvdAD9gJfA95u+UKOUkeYjbppTqCtQFvtdaH8R8oH7sNjY1CDiptV6ltc4G5gKxhZTNBjyBAEBprUO11hcL27DW+getdYzW2qS1/g44A3TIU+QvrfUirbUR8+mYmkB1S4IZCDyntU7QWmdrrXdY1nkaWKC13mf5Jf4NkAV0AoyAM9BcKeWotY7QWp+9jc/kmg+01vFa6wuYP5cRRVlJKeUPdAVe01pnaq2PAF8Aj99BLCilvDAn1Xe11kmFFPMA8i5LAjxu0u7xBPCL1joBWAEMVEr53kmcouRJ8hB34gngV631VcvrFdzeqataQOS1F9rcW2dkQQW11luBT4F5wCWl1ELLAa1ASqnReU4vJQItMdcyrslNUlrrdMtTD8AfiLcc0PKrC0y6tk3Ldv2BWlrrcOAlYBpwWSm1UilV6xbv/2byfg5/Yf6siqIW5vhT8q1f+3YDUUq5AmuBvVrrf9+kaCqQ92/iBaTqAnphtWzzIcy1RrTWe4AL3N6PEGFFkjzEbbH80z8M9FBKxSqlYoGXgTZKqTY3X/sGFzEffK9tW+V9nZ/W+r9a67aYT4U04e9G3usOTkqpusAizKeRqmqtKwPHydOechORgLdSqnIhy97XWlfOM7lprb+1xLdCa32tVqaBDwrZRxrm02LX1CigTN7PoQ4QU8i28h+YYyzxe+ZbP7qQ9W/K0oj9s2X9WzV+n8DcWH5NG8u8gjyAObl8lud7VBs5dVXqSfIQt+t+zKdomgOBlqkZEELx//HXAy2UUg8qpRyAiRR8IEUp1V4p1VEp5Yj54JtpiQPgEtAgT3F3zAfVK5Z1x2KuedyS5VTYL5gPalWUUo5Kqe6WxYuA5yxxKKWUu1JqsFLKUynVVCnVy3KwzQQy8sSX3xHgQaWUm1KqETCugDJTLPv3B14EvitkW5eAetca1LXWkcBu4N9KKRdLg/44LL/w81NK2SmlXABH80vlopRysixzBFZZ3storbWpsM/NYgnwilKqtqXWNQnzBRUFeQJz+1Er/v4edQEClVKtbrEfYUOSPMTtegL4Wmt9QWsde23CfEpppCUJFInltNdDwEwgDmgM7CqkuBfmg3cC5tMwccBsy7IvMbc1JCqlftZan8R8pc8ezAfXVjfZbkEex9zGcgq4jPl0FFrrPzC3e3xqiSOcv68mcra8j6uYT4n5Yr4SqiBzAIMltm8o+MD+P+Ag5kSz3vIeC/KD5TFOKXXI8nwEUA9zLeQn4B2t9W+FrN8dc3LYgLmGkgFcuzEzGPPltv2ARMsVXamWhm6UUt2UUql5trUA8+mtY5hreust866jlKoN9Abm5v0OWdrPNnIHV++JkqdkMCghSiellAYaW9pRhChVpOYhhBCi2CR5CCGEKDY5bSWEEKLYpOYhhBCi2Ip8RUxpZ7k0sj3mewYKuzRSCCHE9ewx96xwQGudVdSVyk3ywJw4QmwdhBBClFHdgN9vWcqiPCWPiwAhISH4+fnZOhYhhCgToqKi6NatG1iOoUVVnpKHEcDPz4969erZOBQhhChzinW6XxrMhRBCFJskDyGEEMVWnk5bFcpkMhEVFUVaWpqtQxE25u7ujp+fH3Z28rtJiDtRIZLH1atXUUrRtGlTOWhUYCaTiejoaK5evYqvr4w1JMSdqBBH0sTERKpXry6Jo4Kzs7OjevXqJCUVNgCeEKKoKsTR1Gg04ujoaOswRCng6OhITk6OrcMQ4q5KzEwkx2Td73WFSB4AhQ+fLCoS+R6I8iQpK4n/Hvov/Vf3Z925dVbdd4Vo8xBCiPIkKSuJpSeXsix0GWnZafSv15/WPq2tGkOFqXmUZkopUlNTCQwMJCMjo9ByiYmJzJo1y4qRCSFKk2RDMp8d+YyBqwey4OgCgmsFs3roamb3mE2Dyg1uvYG7SGoepciRI0duuvxa8nj11VetFNHdk5OTg4ODwy3nCWEtEUkRbL6wmX51+1HHq46tw7mpVEMqy0KXseTkElIMKfSu05vxbcbT1LupzWKqcP+57649wcmY5BLZdvNaXrxzX4tblvvxxx/55z//ibe3N4MGDcqdr5QiJSUFNzc3XnjhBbZu3YqzszMeHh7s2rWL559/nsTERAIDA3Fzc2P37t189NFHrFy5kpycHFxcXPj8888JDAzM3d7777/PTz/9RFxcHB9++CHDhg0DYM+ePUyZMoWUlBQAPvzwQ/r160dYWBgvvfQSV69exWAw8NJLLzF27NhC38uWLVuYOnUqmZmZ5OTk8Oabb/Loo48C0LNnT4KDg9m3bx8uLi7MmzePdu3a8cILL7B582ZGjRpF48aNC1z/wIEDjB07luPHj+fuq02bNnz++ecEBwcX/48jRD6zDswiJDqE/xz6Dx1rdGRYk2H0rtMbJ3snW4eWKy07jeWhy/nmxDckG5Lp6d+TCW0m0KxqM1uHVvGSh61dvnyZp59+mt27d9O0adMCT0P9+eefbN68mVOnTmFnZ0dCQgJA7sE3bw1l9OjRTJo0CYDNmzfz3HPPsXfv3tzlXl5eHDhwgF27dvHwww8zbNgw4uPjeeCBB/jxxx8JDg7GaDSSnJxMTk4Ojz32GMuXLycgIICUlBTatWtH586dCQgIKPD9BAUF8fvvv2Nvb8+lS5do27Yt/fv3p0qVKgAcP36cTZs24eDgQEREBHFxcTRr1oxp06YBkJCQUOD67du3x8PDgx07dtCjRw9CQkKws7OTxCHuioikCEKiQ3gs4DG8Xbz58cyPvLrzVSo7V2Zow6EMazKMBpWsexoor/TsdFacWsE3J74hMSuRHn49GB84nhZVb/3j1FoqXPIoSs2gJO3du5egoCCaNjVXN5955hlee+2168o0aNAAo9HIuHHj6NWrF0OGDCl0ewcPHmTGjBnEx8djZ2fH6dOnr1t+rRbQqVMnYmJiyMzMZM+ePTRv3jz3QGxvb0+VKlU4efIkoaGhuesAZGVlERoaWmjyuHLlCk8++SRnzpzBwcGB+Ph4wsLC6NSpEwCPPfbYdaemXFxcePjhh4u0/sSJE/nss8/o0aMH8+bN4/nnn7/l5ytEUaw4tQJHO0eebv00Pq4+PN36afbG7GXVmVWsCF3BkpNLCPINYliTYfSr2w8XBxerxJWenc7KsJUsPr6YhKwEutbuyoQ2E2hVrZVV9l8cFS552FpRhv2tVKkSJ06cYPv27WzZsoXXXnuNQ4cO3VDOYDAwfPhwdu7cSVBQEDExMdSuXfu6Mi4u5i+9vb09YG5nKCwGrTU+Pj63bHvJa/z48QwdOpQff/wRpRRNmjQhMzMzd7mHh8d15d3d3a+7XPZm6z/00EO88cYbHD58mG3btvHVV18VOS4hCpNsSObn8J8ZWH8gPq4+ANgpO4JrBxNcO5irGVdZc3YNq0+v5s3f32TmvpkMaTiEYY2HlVgbQ0ZOBt+Hfc9Xx78iPjOeLrW6MD5wPG2qtSmR/d0NcrWVlXXu3JnDhw9z5swZAL744osbyly5coWMjAwGDBjAzJkzqVSpEufOncPLy4v09PTcm9yutRP4+/sD8NlnnxUphuDgYE6ePMmePXsA802UCQkJNG3aFDc3N5YuXZpb9tSpUyQnF95GlJiYSL169VBK8dtvvxEeHl60D6II6zs6OvLkk08ydOhQRo4ciZubW7G2LURBfjrzExk5GYxsNrLA5T6uPjzZ8knWPbCOr/p/RTe/bqw6vYrha4fz2PrHWH16NenZ6XcllsycTJaeXMrA1QOZ/cdsmlRpwtKBS5nfd36pThwgycPqfH19WbhwIffddx/BwcEFXm0UGRlJnz59aNOmDa1bt2bgwIF06tQJb29vRo4cSatWrQgODsbLy4t//etftG/fnu7du+Pu7l6kGLy9vfnxxx955ZVXaN26NW3btuXgwYM4ODiwdu1aVq5cSevWrWnRogUTJkzAYDAUuq2ZM2cyefJkOnfuzKpVq2jdunjXmt9q/aeeeoro6GjGjx9frO0KURCjyci3p74lyDeI5lWb37SsUor2NdrzQfcP2PrQVl5t/yrp2elM2zONe7+/l3f3vMuJqyeKdDYhvyxjFstDlzPox0HMOjCLhpUbsnjAYhb1W0Sgb+Dtvj2rUrfzxksjpVQ94Pz58+dvGAwqNDSUZs1sf3WCKL5ly5bx7bffsn79+ru2Tfk+VFxb/trCS9tfYk7POfSp26fY62ut+fPKn/xw+gd+jfiVTGMmAd4BDGs8jMENBuPp5HnT9Q1GA6vPrOaLo19wOeMybau35fnA52lfo/3tvqU7FhERQf369QHqa60jirqe1do8lFJNgG+AqkAcMFprfSZfGV/ga8AfcAK2AhO11tIZUQXUv39/zp49y5o1a2wdiignloYupZZ7Le71v/e21ldKEegbSKBvIK91eI0N5zaw+sxq3t/3Ph/98RH96/VneJPhtKnW5rq2PYPRwE9nfmLRsUVcSr9EkG8QM7rNoEONDmW2yxxrNpjPB+ZprZcppUYBC4Be+cr8EwjVWg9WSjliHoz9QeB7K8Yp8rl8+TL9+vW7Yf6DDz7I22+/XWL73bRpU4ltW1Q8oXGhHLx0kMntJmNvZ3/H2/Ny8uLRgEd5pOkjnIw7yQ+nf+CX87/wv7P/o1HlRjzY+EEG1h/ItshtLDy6kNi0WAKrBTK9y3Q61exUZpPGNVZJHpYaRRDQ1zLrW+BTpVQ1rfWVPEU14KmUsgOcMdc+oq0Royicr69vsa7AEqI0Wha6DFcHVx5o/MBd3a5SihY+LWjh04JX27/KL+d/YfWZ1cw6MItZB8z3cbX2ac20ztMIrhVc5pPGNdaqefgD0VprI4DW2qiUirHMz5s8pgOrgYuAO/Cp1npX/o0ppSoDlfPN9iuJwIUQZd/VjKv8cv4XhjUehpeTV4ntx83RjWFNhjGsyTDC4sPYfGEzrX1a07V213KTNK4pbfd5PAQcBXoDnsAvSqnhWutV+cq9BLxj7eCEEGXTD2E/kG3K5rFmj1ltn029m9q076mSZq1LdSOB2kopewDLYy3L/Lz+ASzXWpu01knA/4CCWrbmAvXzTd1KKHYhRBlmMBr4Luw7utXuRv1K9W0dTrlhleShtb4MHAFGWGaNAA7na+8AOA8MAFBKOQF9gOP5yqC1TtRaR+SdgKiSil8IUXZtjNhIXGYco5qNsnUo5Yo1bxJ8DviHUuo05hrGcwBKqQ1KqXaWMi8B3ZRSxzAnm9PAIivGWCFt376ddu3a3bpgHj179mTdOvPIZU899RQhISE3LT937lwuX7582zEKcTu01iw7uYwGlRrQuVZnW4dTrlitzUNrfQroWMD8QXmen+XvK7JEGVFQFyv5zZ07lz59+uDr62uFiIQwO3T5EKHxobzV6a1y12Bta6Wtwbzk/fI6xB4rmW3XaAUDZ96yWEFjafz666/s2LEDg8GAj48PX331FXXr1iUiIoJ27doxZswYdu7cSUZGBp999hndunW76TKADRs28P7775OZmYmTkxNz5szJ7e126tSprFy5ktq1a9OhQ4dbxnzy5EnGjh1LdnY2zZs3v67zw549ezJ58mSGDBnCwoULmTNnDs7OzphMJr7//ntWr15NTEwMw4cPx8XFhRUrVnDx4sWbjgPSvn179uzZQ0xMDA8//DAzZ5o/1+joaCZOnJjbN9iIESN44403SE5O5pVXXuHo0aNkZmZy77338vHHH+d2CCkqpuWhy/Fy8uK+hvfZOpRyp+IlDxsrbCyNoKAgZs+eDZh/yb/22musXLkSgLi4OFq3bs3s2bPZsWMHI0aM4OzZszddFhUVxfTp09m0aRNeXl6cOHGCgQMHcuHCBdauXcuaNWs4cuQIrq6u3H///beM+/HHH2fixIk88cQT7N27ly5duhRYbsqUKRw/fhx/f3+ysrIwGo28+eabLFq0iFWrVtGyZUsAatasedNxQC5cuMDOnTtJSUmhYcOGjBs3jsaNGzNq1CgGDRrE6tWrAbh69SoAr7zyCj169OCLL77AZDIxcuRIvvrqK55++uk7+GuJsiw6NZotF7YwpsUYXB1cbR1OuVPxkkcRagYlqbCxNJYuXcq8efNITU3N7TX3GicnJ0aNMjf29ejRA1dXV8LCwvDy8ip02e+//87Zs2fp3r177nZycnK4dOkS27Zt45FHHsntLn3cuHG89957hcacnJzM8ePHefzxxwHz2CCtWhU8vkCvXr0YO3Ys//d//8fgwYNp0KDgAXVuNQ7IQw89hJ2dHZUqVaJZs2acPXuWmjVrsnv3bn777bfc7fj4mLvUXrNmDfv37+ejjz4CID09HT8/ufWnIlt5aiUKxYiAEbcuLIqt4iUPGyuoI8q//vqLl19+mQMHDlC/fn12797NY48Vfj261rrQ87fXlmmtGTBgAEuWLClSDLdS1PPFP/74IwcOHGDr1q3ce++9zJ8/n4EDB95Q7lbjgFwbhwTMCTZ/Qs1Pa83PP/9caLISFUt6djqrz6ymT90+1HCvYetwyiXpkt3KChpL48KFCzg5OVGjRg1MJhPz58+/bh2DwcCKFSsACAkJITMzM3ckwsKW9evXj40bN3LixInc7Rw4cACA3r178/3335OWlobRaOTrr7++acxeXl60bNkydz/79+/n2LEb241ycnI4d+4cHTp04PXXX6dfv34cPnw4dxtJSUm5ZW9nHBAPDw+Cg4OZM2dO7rxrp62GDh3KzJkzMRqNufPPnz9/y22K8mnN2TWkGFLk8twSJMnDygoaSyMrK4uHHnqIFi1a0KtXr2vdI+eqWrUqZ86coWPHjkyYMIFvv/0WJyenmy5r3Lgxy5YtY9y4cbRp04ZmzZqxYMECAIYMGcKQIUMIDAykV69e3HPPPbeMe8mSJXzyyScEBQWxaNGi3NNLeRmNRsaMGUOrVq1o06YNFy9e5NlnnwVg4sSJjB07lsDAQE6ePHnb44AsW7aMXbt20bJlS9q0acOXX34JmK/msre3p02bNrRq1YoBAwYQHS3dolVEJm1ieehyWlZtWeoHVCrLZDyPUu7aFVXXfmEXdZkoXFn+PohbC4kKYcKWCczsNpPBDQbbOpxS73bH85CahxCiXFkeupxqrtXoV/fGYQTE3SPJo5SrV69eoTWLmy27HRs2bCAwMPCGacOGDXdtH0KUpHOJ59gVs4tHAx7F0d7R1uGUa3K1lcg1aNAgBg0adOuCQpRSy0OX42TnxPAmw20dSrknNQ8hRLmQlJXEmrNrGNxgMN4u3rYOp9yT5CGEKBdWn1lNpjGTkc1G2jqUCkGShxCizMsx5fDtqW/pUKNDuR6AqTSR5CGEKPO2XNhCbFqs3BRoRZI8bKBevXoEBAQQGBhIQEAATz/9NNnZ2QDMnz//ujuoi0opRWpqKmDuNyooKIipU6cC5l5qGzRoQGBgII0aNaJv376sX78egE2bNuVeVVWjRg18fX1zX//000936R0LUbKWnVyGn4cf3f2637qwuDu01uViAuoB+vz58zq/kydP3jDPlurWrauPHTumtdY6JydHd+7cWa9cufKOtgnolJQUHRkZqQMCAvTs2bNzl/Xo0UOvXbs29/W2bdt09erV9apVq67bxjvvvKMnTZpU6D6ys7PvKMbSorR9H8SdOXblmG65uKVeemKprUMpk86fP68BDdTTxTjmSs3DxjIzM8nMzMztinzatGlMnjwZMHf3MXnyZFq2bEnLli2ZPHlybt9NBQkPD6dHjx5MmjSJSZMmFVquZ8+eTJs2LXeMjJsZNWoUL774Iv3798/tCXjPnj307NmTdu3a0a5dO3755Zfc8uvWrSM4OJi2bdsSHByc25+WECVlWegy3B3dub/RrYcWEHdPhbvP44P9H3Aq/lSJbDvAO4DXOrxWpLLXBkY6e/Ys/fr1o1+/G++GXbhwIUeOHOHQoUMADBw4kIULFzJ+/PgCt9mnTx9mzZrFk08+ecv9d+zYkVdeeaVIse7du5dt27bh5uZGfHw8EyZMYOPGjVSvXp3o6Gg6duzIyZMniY2NZcaMGfz66694eHhw9OhRhg4dSkRERJH2I0RxXUm/wqaITTza9FE8nDxsHU6FIjUPG1m1ahVHjhzhypUrZGZmMnfu3BvKbN68mTFjxuDk5ISTkxNjx45l8+bNhW5z8ODBfPnll7kjFN6MLkafZg899BBubm4A/P7775w/f57+/fsTGBjI4MGDUUpx7tw5Nm7cSHh4OF27diUwMJDRo0djMBiIi4sr8r6EKI7vwr7DaDLyWEDhQxiIklHhah5FrRlYi4uLC0OGDGHdunW89NJL1y3TBYzbcbNxNT799FOmTJnCgAED2LhxI56enoWWPXDgQO6ofrdybdCoazEFBQWxdevWG8rt2LGDIUOG8NVXXxVpu0LciSxjFj+c/oEe/j3w9/K3dTgVjtQ8bMxkMrFjxw6aNGlyw7K+ffuyePFisrOzyc7O5ptvvqFPnz6Fbkspxeeff07Lli0ZMGBAoTWQkJAQpk2bxmuvFT+RdunShZMnT7Jz587cefv37wdgwIABrF+/ntDQUMCcaKTNQ5SUDec2EJ8ZL5fn2kiFq3mUFtfaPAwGAy1btuTtt9++ocwzzzxDeHh47ngb/fv3v+WY3Eop5s+fz7PPPptbAwHzeBpTp04lLS2NunXrsmjRIoYMGVLsuH18fPjf//7Hq6++SmJiItnZ2TRo0IB169bRtGlTFi9ezBNPPEFWVhYGg4Hu3bvTvn37Yu9HiJvRWrMsdBmNqzSmQ40Otg6nQpLxPESFI9+Hsu9A7AGe3PQk7wa/y4ONH7R1OGWajOchhKgwlp5cShXnKgyqL71A24okDyFEmRKZEsn2yO0MbzIcFwcXW4dTYUnyEEKUKStCV2Cv7Hk04FFbh1KhVZjkUV7adsSdke9B2ZZqSOWn8J/oV68fvm6+tg6nQqsQycPFxYW4uDg5cFRwWmvi4uJwcZFTHWXV/87+j7TsNLk8txSoEJfq+vn5ERUVxZUrV2wdirAxFxcX/Pz8bB2GuA0mbWJ56HLaVGtDq2qtbB1OhVchkoejo+O1S9GEEGXUzqidRKZEMjFooq1DEVSQ01ZCiLJvWegyqrtVp3ed3rYORSDJQwhRBpxOOM2+i/sYETACRztHW4cjkOQhhCgDVoSuwMXeheFNhts6FGEhyUMIUaolZCaw7tw67mt4H5WcK9k6HGEhyUMIUaqtOr2KLGMWI5uNtHUoIg9JHkKIUivblM3KUysJrhVMw8oNbR2OyEOShxCiVIpIimDuwblczrgstY5SyGr3eSilmgDfAFWBOGC01vpMAeUeBt4CFKCBPlrrS9aKUwhhGzmmHA5fPsyOyB3siNpBRHIEAF1rd6Vr7a62DU7cwJo3Cc4H5mmtlymlRgELgF55Cyil2gHTgF5a61ilVCUgy4oxCiGsKMWQwq7oXWyP2k5IVAjJhmQc7BzoUKMDIwJG0NO/J7U8atk6TFEAqyQPpZQvEAT0tcz6FvhUKVVNa523z5CXgdla61gArXWSNeITQlhPZHIk26O2syNyBwcvHSRH51DFuQo9/XvS078nwbWCcXd0t3WY4hasVfPwB6K11kYArbVRKRVjmZ83eTQHziuldgIewI/A+zpfj4ZKqcpA5Xz7kA6LhCiFjCYjR68eZXukOWGcTToLQMNKDRndYjQ9/XvS2qc19nb2No5UFEdp69vKAWiNuYbiBGwELgBL8pV7CXjHuqEJIYoqLTuNXdG72BG1g5CoEBKyEnBQDrSt0ZbhTYbTw78H/p7+tg5T3AFrJY9IoLZSyt5S67AHalnm5/UXsEprnQVkKaX+B3TgxuQxF1icb54fEHLXIxdCFElMaoy5dhG1gwOxB8g2ZePl5EU3v2709OtJl9pd8HTytHWY4i6xSvLQWl9WSh0BRgDLLI+H87V3AKwABimllloRbnaaAAAgAElEQVRi6w2sKmB7iUBi3nlKqZIIXQhRiBxTDkcuH2Fn9E5CokIITwwHoJ5XPUY2G0kPvx4E+gbiYFfaTnCIu8Gaf9XngG+UUm8DCcBoAKXUBuBtrfUfwEqgHXASMAGbgC+tGKMQ4ibiM+P5Pfp3dkbtZHf0blKyU8yno6q35f5299PDrwf1KtWzdZjCCqyWPLTWp4COBcwflOe5CXjFMgkhbMykTYTGhxISFUJIVAjHrh5Do/Fx9aFP3T508+tG55qd8XDysHWowsqkPimEuE5adhp7YvawM2onIdEhXM24ikLR0qcl4wPH092vO828m2GnpIOK0iA6MYNFO8/xWMc6NKluvTYlSR5CVHBaayKSI8zJIiqEg5cPkmPKwdPRk+DawXT3606XWl2o6lrV1qGKPM5eSWX+9rP8dDgagGY1PSV5CCFKVo4ph/0X97Mzemfu8K5gvvfi8WaP082vG4G+gTLwUil0IiaJz7adZcPxizjZ2zGqU12e7t6A2pVdrRqHJA8hKhCtNdsjtzPn0BzOJ53H2d6ZDjU6MLr5aLr5daO2R21bhygK8UdEPPO2hbMt7Aqezg6M79GQJ7vWx8fD2SbxSPIQooI4duUYHx38iIOXDlLPqx4f9viQHn49cHWw7i9WUXRaa3aeucq8beHsPx+Pt7sTU/o3ZVSnulRytW2tUJKHEOVcZHIk/z38XzZGbMTbxZupHafyYJMH5ZRUKWYyaTadiGXe9nCORydTw8uFt4c059EO/rg5lY7DdumIQghx1yVmJrLg6AJWhq3E0c6RZ1s/y9iWY6XTwVIs22hizZEYPtseztkradSr6sYHw1rxwD1+ODmUrqvbJHkIUc5k5mSyPHQ5Xx77krScNB5o9AATAifg6+Zr69BEITKzjfzwRyTzd5wjOjGDgBqefDLiHga1qom9XensPUOShxDlhEmbWHduHZ8c/oTYtFh6+PXgpaCXaFSlka1DE4VIzcph2d6/+CLkPFdTswiqU5np97fg3qa+pb7LJUkeQpQDu2N28/EfHxOWEEaLqi2Y0XUG7Wu0t3VYohAJaQa+3h3B4l3nSc7MoVtjHyb0vIdODbxLfdK4RpKHEGVYWHwYHx/8mN0xu6ntUZtZ3WfRv15/ufu7lLqUnMminedYsf8C6QYj/VtUZ0LPRrTxzz88UelXpOShlPIH2mAegCkR+FNrnb87dSGElcSmxfLJ4U9Ye3Ytnk6eTGk3hUcDHsXJ3snWoZUZm07E8unWcDycHajk6mie3Bz/fm6ZKueZ5+nieFttEH/FpTF/xzlWH4zCqDVD29RifM+GVr0j/G4rNHkopRyBZy1TAyAcSAE8gUZKqfOYxyVfqLU2WCFWISq8FEMKXxz7guWhy9FaM6bFGMa1Gkcl50q2Dq1MScrI5s2fjuHsYI+Lox1nr6SSlJFNUkY2WTmmQtdTCjydHa5LMpVdnfAqJOHYKcV3By6w5s8YHOzseKidH892b0idqm5WfLcl42Y1jz+BrZiTx75rQ8gCWAZz6gCMBA4DLUoySCEEbPlrC9P2TCMxK5H7GtzHC/e8QC2PWrYOq0ya89tp4tIMrH2hKy1rX594M7ONJGdkk2hJJknp5sdrr5Ovzc/IJjHdQGxSMkkZOSRnZGMw3ph43JzseapbA57qWh9fLxdrvcUSd7Pk0VNrfbmgBZZEsgfYo5SqViKRCSFyLQ9dzgf7P6BF1RYs7LuQZlWb2TqkMutkTDJL9kQwsmOdGxIHgIujPS6O9sU+0Gutycw2kZhhyE06qVk5BNWpQhX38nc6sdDkUVjiAFBKuQJGrbWhgNEAhRB3iUmbmHtwLl+f+Jpe/r2Y2X2mdCdyB7TWvLPmOJVcHZncr+ld3bZSClcne1ydXKlZqfz/jYp0SYZSarZSqoPl+WAgHkhUSt1XksEJUZEZjAbeCHmDr098zSNNH+Hjnh9L4rhDPx2O5kBEAq8PDKCyW/mrDVhTUS/VHQm8bXn+NjAKSALmAGtLIC4hKrQUQwovb3uZfbH7eDHoRca1HFdmrv8vrZIzs5mx4RSB/pV5qK2/rcMp84qaPNy01ulKqapAA631agClVN2SC02IiulS2iUmbJnAucRzzOg6g/saSgX/bjA3kmfx1Zh22JXSLj/KkqImj9NKqZFAI+A3AKWUD5BRUoEJURGFJ4Tz3ObnSM1OZV6feQTXCrZ1SOXCqdhkluz5i8c61KG1X9m7Ia80KmrymAD8BzAA4yzz+gO/lkRQQlREB2IP8OK2F3Gxd2HxgMUEeAfYOqRyQWvN2z+fwMvFgSn9724jeUVWpOShtT4ABOebtxxYXhJBCVHRbIzYyD9D/om/pz+f9/lc7t+4i34+Es3+iHj+/WAraSS/iwq92kop1aYoGyhqOSFEwZacWMKUHVNo5dOKJQOXSOK4i1IsjeRt/CvzSDtpJL+bblbzmKeUSgaWAju01jHXFiilagI9gNGAB9C9RKMUohwyaROz/5jN0pNL6Vu3L//u9m+c7W0zHnV5NXfzGa6mZvHlE9JIfrfd7CbBrkqpIcBzwJdKKSN/922lgM3Ap1rrDVaJVIhyJMuYxZu/v8mmiE2MbDaSKe2mYG9nb+uwypVTscks3h3Bo+2lkbwk3LTNQ2u9Dlhn6SSxMeZedROAM1rrHCvEJ0S5k5SVxIvbXuTgpYNMajuJJ1o8Ifdw3GVaa97+3wk8XRx4VRrJS0RRG8yzgZMlHIsQ5V5sWizjN48nIjmCD7p9wKAGg2wdUrm05s8Y9p+PZ8YDrcplv1KlgQwGJYSVhMWHMWHzBNJz0lnQZwEdanawdUjlUkpmNu+tD6W1XyUeaS+N5CVFhhsTwgr2XdzHmI1jQME3A7+RxFGC/mNpJJ/+fy1va+AmUTSSPIQoYUcuH+G5zc9Rw70Gywctp0mVJrYOqdwKi03h690RPNrev0wO7VqWFCt5KKX8lVKdSioYIcqbpKwkXt35KjXcarB4wGJquNewdUjllrmR/DieLg5M6S9355e0onbJXkcptQs4hfkSXZRSw5VSX5RkcEKUZVprpu6aypWMK8zuMVuGii1ha/6MYd/5eKb0b4q3NJKXuKLWPBYA6zHf45Ftmfcb0LckghKiPFgeupztkduZ1HYSLXxkpOaSlJqVw/vrQ2lVuxKPtq9j63AqhKJebdUBGKy1NimlNIDWOkkpJT+lhCjA8avH+ejgR9zrfy8jm420dTjl3n82n+ZyShYLHm8rjeRWUtSaxyXM3bHnUko1By7c9YiEKONSDClM3jGZaq7VmN5lutwAWMLOXErh613mRvJ76lSxdTgVRlGTx2zMd5qPBRyUUiOA74APSiwyIcogrTXv7H6H2LRYZnWfJe0cJezaneTuzg68OkAaya2pqHeYf6WUigeeASKBJ4C3tNY/l2RwQpQ134d9z29//cbLbV8m0DfQ1uGUe2uPXmTPuTim399SGsmtrMiX6mqtf9ZaD9Jat9BaDyhu4lBKNVFK7VFKnbY8Nr5J2aZKqXSl1Ozi7EMIWzoVf4pZB2bRtXZXxrQYY+twyj1zI/lJWtb24rEO0khubUXunkQp1Q24B3MX7Lm01jOKuIn5wDyt9TKl1CjMV3D1KmA/9pZlUqsRZUZadhpTdkyhsnNl3u/6PnZK7r8taZ9sOcOl5Cw+HyWN5LZQpOShlPoEeBgI4fpxy3UR1/cFgvj70t5vgU+VUtW01lfyFX8dWIc5SXkgRCmntWb63ulcSLnAF/2+wNvF29YhlXtnLqXw5e/nebidH0HSSG4TRa15jARa5h0Qqpj8gWittRFAa21USsVY5ucmD6VUa8xjo98LvFXYxpRSlTF3D5+X323GJsQd+Tn8Z9afW8/zgc/TvkZ7W4dT7mmteWfNCdyc7HlNGsltpqjJIxLIKslALGOGLALGWpLLzYq/BLxTkvEIURThCeHM2DeDjjU78nSrp20dToWw/thFdp+NY/r/taCqh4y8aCtFTR7jgEVKqW8x3/ORS2u9swjrRwK1lVL2lsRgD9SyzL+mJtAQ2GBJHJUBpZTy0lo/k297c4HF+eb5YT6tJoRVpGenM3nHZNwc3ZjZbaaMBGgFaVk5vLculBa1vHisY11bh1OhFTV5tAUGYh6rPH+bxy0vc9BaX1ZKHQFGAMssj4fztndorS8APtdeK6WmAR5a68kFbC8RSMw7T27EEtY2c/9MziWdY0HfBfi4+tx6BXHH/rv1DLHJmcwbGSSN5DZW1EtCZgD3aa19tNb+eabiXB/3HPAPpdRp4B+W1yilNiil2hUvbCFsa+3ZtfwU/hNPt36azrU62zqcCiH8cgpfhpznobZ+tK0rjeS2VtSaRxpQlNNThdJanwI6FjC/wHE4tdbT7mR/QpSUiKQIpu+dTpBvEOPbjLd1OBXCdY3kA6WRvDQoas3jbWCuUqqGUsou71SSwQlR2mQZs5i8YzLO9s580P0DHOxkJGdr2HAsll3hcUzu3xQfaSQvFYr6zf/K8vhsnnkKc5uHtBKKCuPDAx8SlhDGvN7zZGAnK0nLyuG99SdpXtOLkdJIXmoUNXnUL9EohCgDNkVs4ruw7xjbYizd/brbOpwK45Ot4VxMyuTTx+6RRvJSpKgdI/5V0oEIUZpFJkcybfc0WldrzT+C/mHrcCqMo1GJfPn7OYa39aNtXblzvzQpNHkopRZeu79CKbWUQroi0VqPLqHYhCgVDEYDk3dORinFrO6zcLRztHVI5dqFuHTWHYth/dGLnIhJprKbI69LI3mpc7Oax/k8z8NLOhAhSqOkrCTe3fMuJ+NOMvfeudT2qG3rkMqlqIR0Nhy7yLqjFzkalQRAoH9lpg5uxtA2taSRvBQqNHlorf+tlBqhtf5Wa/2uNYMSojTY8tcWpu+dTlJWEi+3fZnedXrbOqRy5WJSBhuOxbLuaAyHL5jv+W3tV4k3BgYwqFVN/L3dbByhuJlbtXkswNwDrhAVRlxGHP/e/282RWyimXcz5vedT4C3nDa5Gy4nZ7Lh2EXWH7vIgYgEAJrX9OLVAU0Z3Komdau62zhCUVS3Sh5yaYOoMLTWbDi/gZn7Z5KWncbEeyYypuUYaeO4Q1dTs/jleCzr/oxhf0Q8WkNADU8m9W3C4NY1aVBNRl4oi26VPOyVUvdykySitd56d0MSwvoupV3ivb3vsT1qO619WvOvLv+iYeWGtg6rzIpPM7DxeCzrj8Ww52wcJg0Nq7kzsVdjhrSuSePqnrYOUdyhWyUPZ+BLCk8eGmhwVyMSwoq01vwc/jMfHvgQg8nA5HaTGdVslPSQexsS0w38euIS645dZFf4VYwmTX0fd56/txGDW9ekaXVP6cC0HLlV8kjTWktyEOVSTGoM03ZPY8/FPbSr3o53g9+ljpeMhX07/nckmtdWHyUz20Qdbzee6d6AIa1r0rymlySMcko65hEVjkmb+D7se+YcnAPA1I5TeajpQzLu+G3QWvPfLeHM2XyaDvW8mTqkGa1qV5KEUQFIg7moUP5K/ot3dr/DwUsHCa4VzDud36GWRy1bh1UmZeUYeX31MX46HM2DQbX594OtcHaQ030VxU2Th9ZaWrVEuWA0GVkWuoxPDn+Ck50T/wr+F/c3ul9+Id+m+DQDzy79gwMRCUzu14Tn720kn2UFI6etRLl2MfUiv/71K2vPriUsIYyefj15q/Nb+Lr52jq0Miv8cirjvjmQ21nhkNZSc6uIJHmIcudawvj1r185euUoAAHeAXzQ7QMG1h8ov5DvwO7wqzy37CCO9nasfKYTQXVkRL+KSpKHKBdyE0bErxy9+nfCeDHoRfrW7UtdLxkH4k59fyCSf/50jPo+7nw1pr10H1LBSfIQZVZMagy//fXbdQmjmXczXgx6kX51+8llt3eJyaSZtSmM+TvO0q2xD/NGBuHlInfdV3SSPESZUlANQxJGyckwGHn5uyNsPBHLyI51eHdoCxzs5ZJmIclDlBFnEs6w8OhCNkVsQqMlYVjB5eRMnlryB8eik5g6uBnjutaX9iKRS5KHKNVC40JZcHQBWy5swc3BjbEtxzKs8TBJGCUs9GIy4xYfICE9m4WPt6Nv8+q2DkmUMpI8RKl07MoxFhxdwI6oHXg6evJs62d5vPnjVHKuZOvQyr1tpy7zwopDeLg48MNznWlZWz5zcSNJHqJUOXz5MAv+XMCumF1Ucq7EC4EvMKLZCLycvGwdWoXwze4I3l17gmY1vfjyifbUqORi65BEKSXJQ9ic1po/Lv3Bgj8XsC92H94u3rwU9BKPBjyKu6MMDmQNOUYT760PZfHuCPo0q85/Hg3E3VkOD6Jw8u0QNqO1Zs/FPSz4cwGHLh/Cx9WHKe2mMLzJcNwc5R4Ca0nNyuEfKw6xLewKT3WtzxuDmmFvJw3j4uYkeYi7xmA0kGxIJi07jdTsVNIM5sfU7FRSDal/z7c8nks8R2h8KL5uvrzR4Q0ebPwgLg5ymsSaohMzGLf4AGcup/L+Ay0Z2VFuphRFI8lD3BW//fUbr+98HYPJcNNyDnYOeDp64uHkQWXnyrzV6S3ub3Q/TvZOVopUXPNnZCJPLfmDTIORr8e0p3uTarYOSZQhkjzEHQuLD+PN39+kSZUm/F+j/8Pd0R0PRw88nDzMj44euDuZ50mSsD2tNd//Eck7a07g4+HM8qc60kSGhRXFJMlD3JHEzERe3PYino6e/LfXf6nmJr9eS7OIq2m88eMx9pyLo1MDbz4ZEUQ1T2dbhyXKIEke4rblmHKYvGMyV9KvsHjAYkkcpVi20cSikHP8Z/MZnBzs+PeDrXiknT920jAubpMkD3HbPj74Mfti9zG9y3RaVWtl63BEIY5GJfLa6mOEXkxmYMsavDu0Bb5ecmGCuDOSPMRtWXN2DUtPLmVks5Hc3+h+W4cjCpBuyGHOb6f58vfz+Hg4M39UWwa0rGHrsEQ5IclDFNuJqyd4d/e7dKjRgUntJtk6HFGAkDNX+OdPx4iMz+CxjnV4bUAAlVylG3Vx90jyEMVyNeMqL257ER9XH2b3mI2jnRyQSpOENAPT15/kx0PRNPBx57tnOtGxQVVbhyXKIUkeosiyjdm8sv0VkrKSWDpoKVVcZAjS0kJrzZo/Y/jX2pMkZWTzj16NeP7eRrg42ts6NFFOSfIQRTZz/0wOXz7Mh90/JMA7wNbhCIuohHSm/nyc7WFXCPSvzPJhrQioIR1JipJlteShlGoCfANUBeKA0VrrM/nKvAU8CuRYpn9qrTdZK0ZRuB9O/8D3p7/nyZZPMqD+AFuHIwCjSbNkTwQfbgoD4J37mjO6cz3pl0pYhTVrHvOBeVrrZUqpUcACoFe+MvuBj7TW6UqpNsAOpVRNrXWGFeMU+Ry+fJgZ+2bQpXYXJt4z0dbhCCAsNoXXVh/lSGQiPZtW4737W+JXRTqTFNZjleShlPIFgoC+llnfAp8qpappra9cK5evlnEUUJhrKlHWiFPcKDYtlpe3vUwt91p80O0D7O3kHLotZWYb+WxbOJ9tP4uXqyP/eTSQoW1qyfCwwuqsVfPwB6K11kYArbVRKRVjmX+lkHVGA2e11jckDqVUZaByvtl+dzHeCkdrTUZOBsmGZFINqaRkp5BiSOHzI5+TkZPBF/2+kFH8bEhrzc4zV3l37QnOXUnjwXtqM3VIc7zdpa8wYRulssFcKdUDmM7fNZX8XgLesV5E5cvRK0eZd2QeCZkJpBhSSMlOIdWQitGc269jp+z4uOfHNKrSyAaRCoADEfF8uCmM/efjqePtxpInO0gPuMLmrJU8IoHaSil7S63DHqhlmX8dpVRnYBnwf1rrsEK2NxdYnG+eHxBy90IunyJTInlhyws42DnQvGpzGlRugKejJ55Of08eTh54OXrh4eRBdbfqVHevbuuwK6Tj0UnM/jWM7WFXqObpzPT/a8Ej7evg5GBn69CEsE7y0FpfVkodAUZgTgwjgMN52zsAlFLtge+A4VrrQzfZXiKQmG/dux53eZNsSOb5Lc9j1EaWDlhKXS8Z+Kc0Cr+cypzfTrP+2EUquzny+sAAnuhcD1cnaW8SpYc1T1s9B3yjlHobSMDcpoFSagPwttb6D+AzwBVYkCcZPK61PmbFOMulbFM2k7dPJjIlkoV9F0riKIWiEtL5z+YzrD4UhaujPRN7N+apbvXxcpG7+EXpY7XkobU+BXQsYP6gPM/bWyueikRrzcx9M9lzcQ//Cv4X7WvIx1yaXE7J5LNtZ1m+7y+UUjzZpT7jezakqoeMsyFKr1LZYC7urmWhy3Jv8Hug8QO2DkdYJKVnM3/nWRbvisBgNPFwO38m9m5EzUqutg5NiFuS5FHO7YjcwYcHPqR3nd68GPSircMRQFpWDl/vOs+CnedIzcphaJtavNynCfV83G0dmhBFJsmjHAuLD2PKzikEeAcwo+sM7JRcpWNLmdlGVuy7wLxt4cSlGejbvDqT+jWRfqhEmSTJo5y6kn6FF7a+gKeTJ5/2/hQ3R+m6wla01vxwMIq5v50mJimT4IZVmdy/KUF1pFdiUXZJ8iiHMnIymLh1IklZSSwesBhfN19bh1ShLQo5x4wNpwj0r8yHD7WhSyMfW4ckxB2T5FHOmLSJN39/kxNxJ5h771yaV21u65AqtB2nrzDzl1MMblWTTx+7R+5HEuWGJI9yIjIlkp1RO/ntr984eOkgk9pOoled/J0WC2s6fzWNf6w4RJPqnnz4UGtJHKJckeRRRuWYcjhy+Qg7o3ayI2oH55LOAVC/Un1ebvsyT7R4wsYRVmwpmdk8veQP7O0Ui0a3w81J/tVE+SLf6DJm78W9/Hj6R36P+Z0UQwoOdg60r96eh5s+TPfa3fH38rd1iBWeyaR5+bs/OX81jaXjOuDvLRcriPJHkkcZcuTyEcb/Nh4vZy961+lND78edK7VGXdHuT+gNJm7+TSbQy8x7b7mBDeUxnFRPknyKCPiM+OZvGMyNdxr8N193+HlJPcGlEa/HLvIf7eG83A7P54IrmfrcIQoMZI8ygCjycgbIW+QkJnA0kFLJXGUUqEXk5n0w5/cU6cy0+9vKQ3kolyTW47LgIXHFrI7Zjevd3xdLr0tpeLTDDy95A88XRxYMKotzg7Sfboo36TmUcrtidnD50c+Z0iDIQxvPNzW4YgC5BhNvLDiEJdTsvj+2c74ernYOiQhSpzUPEqxS2mXeD3kdRpUasBbnd6S0yCl1PsbQtl9No5/P9CKQP/Ktg5HCKuQmkcplW3K5tWdr5KRk8HX/b+WvqlKqe//iOTrXRE82aU+w9r62TocIaxGkkcp9cmhTzh0+RAzu82kQeUGtg5HFODQhQSm/nScLo2q8s9BAbYORwirktNWpdDWC1v5+sTXPNzkYQY3GGzrcEQBLiVn8tzSg9So5MKnI4JwsJd/JVGxSM2jFNBaE5MWQ1h8GGHxYSw9uZRm3s14tcOrtg5NFCAz28izSw+SmpXD0nEdqeLuZOuQhLA6SR42oLVm84XNHLx0kFPxpzgdf5qU7BQAFIqm3k35qOdHONvLGNaljdaaqT8f50hkIvNHBdG0hqetQxLCJiR5WNmF5Au8u+dd9sfux9XBlSZVmjCw/kCaejelqXdTGlduLI3jpdji3RGsOhjFxN6NGdCypq3DEcJmJHlYSY4ph2UnlzHvyDwc7Bx4p/M7PNDoAezt5GaysmJX+FXeWx9K3+bVeal3Y1uHIwQYsyH+HFwOhZqtwdt6F9dI8rCCsPgw3tn9DifiTnCv/7282fFNqrtXt3VYohguxKXz/IpDNKzmzpxHArGzk3tuhBXlTRJXTpmny6cgLhxM2eYyA2ZCp/FWC0mSRwnKzMlk4dGFfH38a7ycvZjdYzb96vaTm/3KmLSsHJ5Z+gcmk2bh4+3wcJZ/G1FCjNkQdzZPggiFK2HXJwkUVKkHvs2g6QCo1gx8A8CniVVDlf+CErInZg/v7X2PCykXGNpwKFPaTaGyi9x9bAtpWTlsOXWZ2KQMMgwmMnOMZBiMZGabp4xsIxnZpr9fG4yWMqbc5VprFo/tQD0f6f5e3AV3kiSqNgYn27eLSvK4y+Iy4pj9x2zWnVtHHc86LOq3iE41O9k6rArHZNLsOx/P6kNRbDh2kXSDMXeZo73CxdEeV0d7XJ3scXGwx8XJHldHO6q6O+FS2bzMPM8eF0c72tfzpnuTajZ8R6LM0frv001lOEkURpLHXWI0Gfk5/Gc+Pvgx6TnpPNfmOZ5q9ZRcbmtlF+LSWX0oitWHoohKyMDD2YH7WtdiWFs/mtfywsXBTm7oq4hMRkiPh7QrkH7V/JgWB5mJYDRYphzzoynbfNA3GiyP2ZZ5eV4bDWDKuXHd6+ZnAzpPEGUzSRRGkscd0lqzM2on/zn8H84knKFt9ba83elt6VLEilKzcthw7CKrDkax/3w8SkGXhj5M7teU/i1q4OokV7TdNdmZkJUMmUmWKREyLa+zUsDBGRzdwNEVnNzNz3Mf3SzLLJPdHSRxkxEyEiDtap6EcLXw1xkJXH8gz0uBvRPYO1omJ7DL8zz/fCd3sK9SSPk869hZnleuU6aTRGEkedyBQ5cOMffQXA5fPkwdzzp82P1D+tfrLw3iVmAyafaei2PVoSh+ORZLRraR+j7uTOnflAfuqU2tyq62DrH00hpSL0NCRL4kkJQvMVyb8swzZt29OBxcLQnF3fLomud5nqQDBSSDeNCmgrfrWgXcq4GbD1RrCvW6/P3a/dpkee1a2XygF8UmyeM2pGenM233NH6J+IVqrtV4q9NbPND4ARzt5EtYkrJyjETGp7PmSAyrD0UTnZiBp7MD999Tm+FtaxNUp4ok7ry0hpRYuGI5z37t8s4rp8zJoiD2zuBS6fqpcp3rXzt7gUvlPPO8zI9OHuZTNdlpYEg3P2Zn/P3ckA7Z6WBIu/4xOyPPvHRIvXR9WfTfB/uqDaFOx8KTgVtVsJfDmqTJfcwAAA4mSURBVDXIp1xMsWmx/GPrPzidcJoJgRMY02IMrg7yK/d25BhNxKcbiEu1TGlZuY/xaQauphqISzU/j0s1kJKVA4BS0LWRD68OMJ+WcnEshaelTEbz1TSXjkHsMbhyGpw9zAc5j+qWyffvR1fv2z+NozUkR1saZf+/vbsNjqs6Dzj+f/Zd2tW7VjK2iAnGNhTwBNdMmxYopCVJialJM07iCU2bD52hH2gy09BOO5mhZIZxJslAB6ctpU2nHWBIQpOUZhKmhKFh0hBchhgCJMh2ABvbaPVmaVfSSlrtPv1wjrQrWSu8a2mllZ7fzJ29957ru+ce6+5zz73nntNbfCg70AtTo8XtGtrcffYrPwrJy90PcUPb/KAQXo6BrDqWYR9mrbPgUYFXB1/lzmfuJDuT5dAHDnFDzw2rnaW6MJnL86uBMY6lxjjWn+Foaozj/WOcGBqnsMht6GBAaI9H6IhH6EhE2NXWSkfCLSebotywI8lFLWsoYE+PQ+oX0PdzFyj6XoHUazCTdemBsPuhzk2420Uzk+fuQ4I+mHRBvOvc4DI7HwzD4LH5tYiBXpjOFPfV2Okeyu7a74JEcqf7jCdd5DVmGVjwOA8FLfDE8Se49/C9dDZ08tDND7G9zbqnWCg77YNEf4ZjqTEfJDKcHJ6YCxKhgHBJZ5wrLmpi766L6GqO0RmPuGCRiNKZiNAcC6/NN7hV3S2Vvld8oHjVzQ8dZ+5hbKwFNu2CPZ+BTVe7qXMnhCLFfUxlXBAZS8F4f3F+LAVjA+4z9ZpLK8yUz0+i2wWG9x3wQcIHinjniheFMRY83sVzp5/j/p/dz+vDr7O7azf33XgfHQ0bs1qemcxxeiTL6bPZuc9TJcsDmeLD1FBAuDQZ58rNLdx2zRa2dzWxvTvBJR1xIqEaNpXNzxQfAuey7qp/ZmqRz2yZ9ZNuSp9xgWJ8oLjv1q0uOFy9vxgoWnqWvroX8c8ImqHzsqXzXii4ZxOlgWUm61rtJHdCY/vylJExVbDgsYjJmUmePfUsj/c+zuG+w2xJbOHg9Qe55b23EJD1+45ALl/g9NksJ4YnODE0zomhCU4MTfhAMUF6cv5VcCQYYHNrjC1tDdy0M0lPWyOXdSXY0Z1ga0ec8IW+TzEz5a7SZ5uBTqX9cnp+q6CpdMm6BZ+5ieq+OxCCUMw1PQ3F3IPY7R8qBonuK11LnZUUCLgA0djubkMZs4ZY8CjRO9zLo798lKdOPMV4bpxkQ5K79tzFJy//JJHg+hjwZ2J6hreHs5wsDRCDadLDKaZGU7TrCJ2M0ilpuoNp9kSztEVmaG7J09Seo1FyxCRHVKcIFqaQmUkYzELfpHsxKhhxTS5DDe7h6zmfMZ8ec9vmJoo/9POCQ+b8moWGGootfqL+ir6lp2S5pbg+3Dg/IISiPi+zy35dMGotdox5F3aGAC+mXuTBlx/k+XeepyHUwAe3fpC92/Zybfe1dddl+tRMntToFG8PpRlInWGk/xTjZ98hN5pCxlI05oZJigsQ10majwbStJIhSAEWtDTWQBiJdS4SDFrnB4HZz1DU1RZyWXd7Jedv+czeLpoY8ut8Wn7aNe+MNUO0CRKb3C2Z2eVoSUCYW27yyz5gWBt9Y1ZFzYKHiOwA/h3Xjm8I+LSqHluwTRB4APgw7gnkl1T1X1Y6b2+OvskbI2/w2d2fZf+O/bREW5Znx/mce6kpPw2BIEjAtaqZnVd1D0QLM6B5P19wP7RjKQqZPqZH3mFm9Aya7kPGUgQnUjAzRaGg5BUKquQLkFfIq5JXoZFp3k+GgJzblCkXiTHdkETjSSLNVxNu6UbiSd/KJ+mbkrp5ibVY6xxjzKJqWfN4EPh7VX1ERG4H/gn4wIJtPgVcBmzHBZkjIvK0qr61khnbt20f+7btI7zUVWwhX3wTNzviP8+6/nImhorTWD+MD6BjKSR79oLyFQBiQFobGNBW+rWNfi5mQqMISiggNEQCNEQCxHwHfrFwkKlojGzrJhLtm2nu3EKwuXuuCWg4mlhYwTDGmIrVJHiISBewG7jZr3oM+JqIJFW1pPkKnwD+WVULwICI/CewH/jKgv21AgufVvZUm7/B//s2wSMPu8aWqgTyUwTzkwTzWUIzY4RzGSL5pR+8ZohzlgSD2kJ/oYUB7WFAWxmkhRxBBCVIgSCFufkCAoEQ0UiEaCRMNBIhFo0SiUbJNybR+CakqZvGeBNNsTCJWIieaGjufYe4jSthjFkltfr1uRg4rap5AFXNi8gZv740eLwHOFGyfNJvs9DngLuXK3O/ONlPe1/f3PIUYSY1QpZ2MtpDhkbS2kiGRsYDCSaCTUwGm8hFWtDGDsLxdprijbQ0hImGA0SCAcLBAA3BAJcGhXg0RFMsRFMsTFMsRHMsRCIaprUxvDbfjjbGmHdRr5eufwf824J1PcCPq9nZNXvvYOjGzxAICIJ7wzkgQiAgRIIBIqEA0ZALCmvy5TVjjKmxWgWPt4EtIhL0tY4gsNmvL3US2Aq84JcX1kQAUNURYF7PbhfSIV67f8PZGGPM+anJG2+q2g+8BBzwqw4ARxY87wB4HPhTEQmISBK4Dfh2LfJojDHm/NXydek7gDtF5Chwp19GRH4gInv8Ng8DbwDHgOeBL6rqGzXMozHGmPNQs2ceqvo68BuLrL+lZD4P/Fmt8mSMMaY667ejJmOMMSvGgocxxpiKWfAwxhhTsXp9z2MxQYBTp06tdj6MMaZulPxmVvTGsqguMg5oHRKR66jyJUFjjDFcr6r/e74br6fgEQWuBd4B8otsMvsG+vWAVU/OZeVTnpXN0qx8lrbWyycIXAS8oKrnMYiOs25uW/mDLhs1S95AP7XSvfTWIyuf8qxslmbls7Q6KZ9fVfoP7IG5McaYilnwMMYYUzELHsYYYyq2kYLHCHAPC3rjNXOsfMqzslmalc/S1mX5rJvWVsYYY2pnI9U8jDHGLBMLHsYYYyq2IYKHiOwQkZ+KyFH/uX2187TSRKTDj5XSKyI/F5Hv+AG2EJHfFJGXfXk8JSJdJf+uqrR6JSJ3i4iKyFV+2coGEJGYiPyjiBwTkVdE5CG/vuy5VG1aPRKRvSJyRERe8ufXH/r1G6d8VHXdT8AzwO1+/nbgmdXOUw2OuR24sWT5K8DXAQGOA9f59V8A/tXPV5VWrxOwG3gSN9TxVVY288rmAeB+is9Fu/1n2XOp2rR6m/z/91ngKr+8C8jgLsY3TPmsegZq8B/dhWvlEPTLQb+cXO281bgcPgY8jevC5dWS9Z3AmJ+vKq0eJyAK/BR4L/CWDx5WNi7/CX+OJBasL3suVZu22sdaZfkIMAT8tl++ATi60cpn3XRPsoSLgdPqRilEVfMicsavXziG+rokIgHcCI3/BbwHd6UNgKoO+jHj26tNU9XhWh3LMvoi8IiqvlnSfYSVjbMN9+N4t4jcBIzhalNZyp9LUmVa3Z2Dqqoi8nHgCREZB5qAj7D0b826K58N8czDcAj3A/C11c7IWiAi78fVFv5htfOyRoWAS4EjqroH+CvgO7gayYYnIiHgr4F9qroVuBX4JhusfDZCzeNtYIuIBH1EDwKb/fp1T0S+CmwHblXVgoicBLaWpHfiLqaGq02r1bEso98BLgdmax09wH/j7vNv9LIBV4uaAR4DUNXDIjKIq3mUO5ekyrR69D5gs6r+BEBVf+JrIJNsoPJZ9zUPVe0HXgIO+FUHcFdUa7Y6uFxE5F7g14HbtNjV8otAg7jxTwDuAL51gWl1RVW/pKqbVfUSVb0E1032h3CNCjZ02YC77Qb8D3AzuJZAuPvyRylzLi11nq3Dc/AU0CMiOwFE5ApgE3CMjVQ+q/3QpRYT7irzMO6P/zCwc7XzVINjvhJQoBf3h/kS8F2f9lvAK7g/9h/iW9JcSFo9T/gH5lY288rkUuBH/ph+Bvy+X1/2XKo2rR4n4FO+bF72020brXysexJjjDEVW/e3rYwxxiw/Cx7GGGMqZsHDGGNMxSx4GGOMqZgFD2OMMRWz4GGMMaZiFjyMWQEi8paI/N5q58OYlWLBwxhjTMUseBjj+drC5/3gPqMi8k0RiZXZdpuIPCMiQyIyKCKPikirT3sY19Pu90RkTET+0q//AxF5TURGRORHvluL0u++y3/3uIh8XUS6ReRJEcmIyNMi0ua3jYnII/67R0TkBRHpXvkSMqbIgocx830c+DBunI9dwJ+U2U6Ag7gO7K7AdZ/9twCq+kfASVxnlAlV/bLvH+ox4HO4MRx+gAsukZJ9fgzXn9QOXE+tTwJ/gxsfJAD8ud/uj4EW/50duL60shd22MZUxoKHMfM9oKpn1PWI+z1cD6rnUNXjqvpDVZ1S14Hdfbjeesv5BPB9/29ywFeBBlyfWLMOqWpKVU8DPwYOq+oRdZ1afhe4xm+XwwWNy1Q1r6ovqmr6Ao7ZmIpZ8DBmvr6S+QnKjNEgIl0i8g0ROS0iaeARXA2hnM3MHzCqgB8uoGSbVMl8dpHl2bw8jOtC/hsickZEviwi4aUPy5jlZcHDmOocxPVavEtVm3HjTktJ+sIeR88wf8wPwY88V+kXq2pOVe9R1V/D1Vz2Ap+udD/GXAgLHsZUpwk3OuOIiGwB7lqQnsJ1az7rW8BHROR3fS3hL4Ap4LlKv1hEbhKRq/2gQWncbax8FcdgTNUseBhTnXuA3cAo8H3cMK2lDgJf8K2hPq+qvbjaySFgEPdA/FZVna7iuzcB/4ELHL8EnsXdNjOmZmw8D2OMMRWzmocxxpiKWfAwxhhTMQsexhhjKmbBwxhjTMUseBhjjKmYBQ9jjDEVs+BhjDGmYhY8jDHGVMyChzHGmIr9P3AR58DZrmeWAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"ax.errorbar(natoms_rdf, [t.average for t in brute_rdf], label='distance_array')\n",
"ax.errorbar(natoms_rdf, [t.average for t in capped_rdf], label='capped_distance')\n",
"ax.errorbar(natoms_rdf, [t.average for t in kdt_rdf], label='Bio KDTree')\n",
"\n",
"ax.set_title('All distances up to 12.0 A')\n",
"ax.set_ylabel('Time (s)')\n",
"ax.set_xlabel('n atoms')\n",
"\n",
"ax.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Conclusion\n",
"\n",
"The new `capped` functions in `lib.distances` are much faster!\n",
"They also impement quite complicated algorithms with minimal difficulty to a user;\n",
"all that is required is to define a maximum \"interesting\" distance,\n",
"which is often known. Ie for bonds we weren't interested in distances larger than 2A, for an RDF often 12A is enough.\n",
"\n",
"These benchmarks represent only the first iteration of this functions,\n",
"and it is anticipated that further optimisations can be made to these."
]
}
],
"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.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment