{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Evaluating convergence\n", "\n", "Here we evaluate the convergence of a trajectory using the clustering ensemble similarity method and the dimensionality reduction ensemble similarity methods.\n", "\n", "**Last executed:** Sep 25, 2020 with MDAnalysis 1.0.0\n", "\n", "**Last updated:** September 2020\n", "\n", "**Minimum version of MDAnalysis:** 1.0.0\n", "\n", "**Packages required:**\n", " \n", "* MDAnalysis (Michaud-Agrawal *et al.*, 2011, Gowers *et al.*, 2016)\n", "* MDAnalysisTests\n", "* [scikit-learn](https://scikit-learn.org/stable/)\n", " \n", "**Optional packages for visualisation:**\n", "\n", "* [matplotlib](https://matplotlib.org)\n", "\n", "\n", "
\n", " \n", "**Note**\n", "\n", "The metrics and methods in the `encore` module are from (Tiberti *et al.*, 2015). Please cite them when using the ``MDAnalysis.analysis.encore`` module in published work.\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import MDAnalysis as mda\n", "from MDAnalysis.tests.datafiles import PSF, DCD\n", "from MDAnalysis.analysis import encore\n", "from MDAnalysis.analysis.encore.clustering import ClusteringMethod as clm\n", "from MDAnalysis.analysis.encore.dimensionality_reduction import DimensionalityReductionMethod as drm\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The test files we will be working with here feature adenylate kinase (AdK), a phosophotransferase enzyme. (Beckstein *et al.*, 2009)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "u = mda.Universe(PSF, DCD)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Evaluating convergence with similarity measures" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The convergence of the trajectory is evaluated by the similarity of the conformation ensembles in windows of the trajectory. The trajectory is divided into windows that increase by `window_size` frames. For example, if your trajectory had 13 frames and you specified a `window_size=3`, your windows would be:\n", "\n", " - Window 1: ---\n", " - Window 2: ------\n", " - Window 3: ---------\n", " - Window 4: -------------\n", " \n", "Where `-` represents 1 frame.\n", "\n", "These are compared using either the similarity of their clusters (`ces_convergence`) or their reduced dimension coordinates (`dres_convergence`). The rate at which the similarity values drop to 0 is indicative of how much the trajectory keeps on resampling the same regions of the conformational space, and therefore is the rate of convergence. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using default arguments with clustering ensemble similarity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See [clustering_ensemble_similarity.ipynb](clustering_ensemble_similarity.ipynb#Calculating-clustering-similarity-with-default-settings) for an introduction to comparing trajectories via clustering. See the [API documentation for ces_convergence](https://docs.mdanalysis.org/stable/documentation_pages/analysis/encore/similarity.html#MDAnalysis.analysis.encore.similarity.ces_convergence) for more information." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "ces_conv = encore.ces_convergence(u, # universe\n", " 10, # window size\n", " select='name CA') # default" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output is an array of similarity values, with the shape (`number_of_windows`, `number_of_clustering_methods`)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.48194205],\n", " [0.40284672],\n", " [0.31699026],\n", " [0.25220447],\n", " [0.19829817],\n", " [0.14642725],\n", " [0.09911411],\n", " [0.05667391],\n", " [0. ]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ces_conv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This can be easily plotted as a line." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Jensen-Shannon divergence')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZfr/8fedRgkdAiq996JGmlQVlbWAHXtdRLFg23W/v3V33dUtrt1FimBZLNg7dqWD0nsxhCooIEgvCdy/P2bQbAxkQmZykszndV25zJw5c+ajV5x7znmecz/m7oiISPxKCDqAiIgES4VARCTOqRCIiMQ5FQIRkTinQiAiEudUCERE4lxMC4GZnWlmy8wsw8zuzeP5Xma2zczmhn/+FMs8IiLya0mxOrCZJQJDgT7AOmCGmb3n7otz7TrJ3c+O9Lg1atTwBg0aRC+oiEgcmDVr1mZ3T8vruZgVAqAjkOHumQBmNhboB+QuBAXSoEEDZs6cGYV4IiLxw8xWH+65WF4aqg2szfF4XXhbbl3MbJ6ZfWRmrWOYR0RE8hDLMwLLY1vufhazgfruvtPMfgO8AzT91YHMBgIDAerVqxftnCIicS2WZwTrgLo5HtcB1ufcwd23u/vO8O/jgGQzq5H7QO4+0t3T3T09LS3PS1wiInKUYlkIZgBNzayhmaUAA4D3cu5gZseYmYV/7xjO82MMM4mISC4xuzTk7tlmdgvwCZAIPOvui8xsUPj54cCFwE1mlg3sAQa42qGKiBQpK2mfu+np6a5ZQyIiBWNms9w9Pa/ndGexiEici5tCsGXXfu5/fxF7sw4EHUVEpFiJm0IwOWMzz09dxSUjp7Nxx96g44iIFBtxUwjObX8cw684keXf7+C8oVNZ9v2OoCOJiBQLcVMIAM5ofQyv3diFrAMHuWDYVCYs3xR0JBGRwMVVIQBoW6cy795yMnWrlee652cwZvph22+IiMSFuCsEAMdWLsfrg7rQs1ka972zkL99sJgDB0vWNFoRkWiJy0IAUKFMEs9clc41XRswevJKbhwzk137soOOJSJS5OK2EAAkJhh/Obc195/bmi+XbuSi4dPYsG1P0LFERIpUXBeCQ67u2oDR15zE6h930X/oFBZ+ty3oSCIiRUaFIKx385q8cVNXEs24aPg0Plv8Q9CRRESKhApBDi2PrcQ7g0+maa0KDBwzk1GTMilpvZhERApKhSCXmpXK8urALpzR6hge+HAJ9727kOwDB4OOJSISMyoEeSiXksjTl5/AoJ6NeXH6Gq59fgbb92YFHUtEJCZUCA4jIcG4t28L/nl+W6at+JELh01l7ZbdQccSEYk6FYJ8DOhYjxeu68iGbXs57+kpzFmzNehIIiJRpUIQgZOb1ODtm7tSLiWRASOnM27BhqAjiYhEjQpBhJrUrMg7N59Mm9qVufml2Qz9KkMzikSkVFAhKIDqFcrw0g2dOLf9cfz7k2X87o357M/WjCIRKdlitnh9aVU2OZEnBnSgQY1UnvziW9Zu3c3wK06kSvmUoKOJiBwVnREcBTPjzj7NeOyS9sxe/RPnPz2VVZt3BR1LROSoqBAUwnnH1+HFGzqxdfd+znt6Ct+s3BJ0JBGRAlMhKKSODavx9s0nU7V8CleM+pq356wLOpKISIGoEERBgxqpvHVzV06oX4U7Xp3Ho58t14wiESkxVAiipEr5FP57XScuPLEOT37xLbePncverANBxxIRyZdmDUVRSlIC/76wHY3SUnno42Ws/2kPI648keoVygQdTUTksHRGEGVmxs29mjD0shNY8N02znt6KhkbdwYdS0TksFQIYuSsdscydmBndu/P5rynpzAlY3PQkURE8qRCEEPH16vK2zefzLGVy3L1s9/w6ow1QUcSEfkVFYIYq1utPG/c1JUujavz+zcX8M+PlnLwoGYUiUjxoUJQBCqVTea5a07i8k71GD5hBYNfns2e/ZpRJCLFgwpBEUlKTOCB/m247+xWfLzoewaMnMbG7XuDjiUiEttCYGZnmtkyM8sws3uPsN9JZnbAzC6MZZ6gmRnXd2vIyCvTWf7DTvoPncLS77cHHUtE4ly+hcDMmpnZF2a2MPy4nZn9MYLXJQJDgb5AK+BSM2t1mP3+BXxS0PAlVZ9WtXh9UBcOuHPhsGl8tWxj0JFEJI5FckbwDPAHIAvA3ecDAyJ4XUcgw90z3X0/MBbol8d+twJvAnH1adimdmXeHdyN+tXLc/3zM3htxtqgI4lInIqkEJR3929ybcuO4HW1gZyfbuvC235mZrWB84DhERyv1Dmmclleu7ELJzepwe/enM/IiSuCjiQicSiSQrDZzBoDDhC+jh/Jor2Wx7bc8yYfB37v7kecQmNmA81sppnN3LRpUwRvXXKklkli9NUncVa7Y/n7uKX86+OlalgnIkUqkl5Dg4GRQAsz+w5YCVwRwevWAXVzPK4DrM+1Tzow1swAagC/MbNsd38n507uPjKcgfT09FL3KZmSlMCTA46ncrlkho1fwU+7s3igfxsSE/KqpSIi0ZVvIXD3TOA0M0sFEtx9R4THngE0NbOGwHeExhUuy3Xshod+N7PngQ9yF4F4kZhgPNi/DVXLJzP0qxVs35PFo5e0p0xSYtDRRKSUi2TW0N/NrIq773L3HWZW1cweyO917p4N3EJoNtAS4DV3X2Rmg8xsUOGjlz5mxj1ntOCPZ7XkwwUbuOGFmezaF8lwjIjI0bP8rkeb2Rx3Pz7XttnufkJMkx1Genq6z5w5M4i3LlKvzVzLvW/Op33dKjx3zUlUKZ8SdCQRKcHMbJa7p+f1XCSDxYlm9nNDfTMrB6jBfoxdnF6XYVecyKL127l4xDS+36a7kEUkNiIpBC8CX5jZ9WZ2HfAZ8EJsYwnAGa2P4flrT+K7rXu4cPhUVm3eFXQkESmF8i0E7v4Q8CDQEmgN/C28TYpA18Y1eGVgZ3bvP8CFw6exeL1aUohIdEXUa8jdP3L3u939LnePm1YQxUW7OlV47cYuJCcal4ycxoxVW4KOJCKlSCSzhs43s2/NbJuZbTezHWamr6VFrEnNCrxxU1fSKpbhytFf89XSuOrIISIxFMkZwUPAue5e2d0ruXtFd68U62Dya7WrlOP1G7vQtGZFfvvfmbw797ugI4lIKRBJIfjB3ZfEPIlEpHqFMrz8206kN6jKkFfn8t9pq4KOJCIlXCQtJmaa2avAO8C+Qxvd/a2YpZIjqlg2meev7citr8zhT+8uYuuuLG47tQnhVh0iIgUSSSGoBOwGTs+xzQEVggCVTU5k2OUncO9bC3js8+Vs3b2fP53digT1JxKRAoqk19C1RRFECi4pMYGHLmhH5XLJjJ68km17snjownYkJ2oFUhGJXMxWKJOikZBg/PGsltxzRnPenvMdg8bMYm/WEbt6i4j8j1iuUCZFxMwY3LsJD/Rvw5fLNnLV6G/Yvjcr6FgiUkLEcoUyKWJXdK7PkwOOZ87arQwYMZ1NO/bl/yIRiXuxXKFMAnBO++MYdfVJrNy8i4tHTGPtlt1BRxKRYi6SQjAYGMEvK5QNAW6KaSoplJ7N0njxho78uHMfFw2fxrc/RLqWkIjEo0iazmW6+2lAGtDC3bu5+6qYJ5NCObF+NV4b1IWD7lw0Yhpz1/4UdCQRKaYiWZjmzjw2bwNmufvcmKQ6gnhZmCZa1vy4mytGf83mnfsYeWU63ZrWCDqSiASgsAvTpAODgNrhn4FAL+AZM/tdtEJKbNSrXp43BnWhXrXyXPf8DD5aoOEdEflfkRSC6sAJ4RbUdxEqDGlAD+CaGGaTKKlZqSyvDuxC2zqVGfzybMZ+syboSCJSjERSCOoB+3M8zgLqu/secvQekuKtcvlkxlzfke5N07j3rQUMn7Ai6EgiUkxE0mvoZWC6mb0bfnwO8IqZpQKLY5ZMoq58ShLPXJXOXa/P458fLWXr7v3ce2YLNasTiXNHLAQW+oR4HhgHdAMMGOTuh0ZrL49pOom6lKQEHr+kA5XLJTFiQibbdmfx4HltSVSzOpG4dcRC4O5uZu+4+4nArCLKJDGWmGD8rV8bqpVP4ckvM9i2J4vHB3SgTFJi0NFEJACRjBFMN7OTYp5EipSZcefpzbnv7FZ8tPB7rn9+Jrv2qXOISDyKpBD0JlQMVpjZfDNbYGbzYx1Misb13RryyEXtmZb5I5eN+pqtu/bn/yIRKVUiGSzuG/MUEqgLTqxDpXLJDH55NhePmMaY6ztxTOWyQccSkSISSYuJ1UBd4JTw77sjeZ2ULH1a1eK/13Vkw7a9XDBsKis37wo6kogUkUgWpvkz8HtCaxIAJAMvxjKUBKNzo+qMHdiZvVkHuGj4VBZ+ty3oSCJSBCL5Zn8ecC6wC8Dd1wMVYxlKgtOmdmVeH9SFMkmJXDpyulpSiMSBSArBfg91pju0HkFqbCNJ0BqlVeD1QV1omJbKTS/N5q7X5rFDK56JlFqRFILXzGwEUMXMfgt8Tmj5SinFjqtSjjdv6sptpzTh7Tnr6PvEJL5ZuSXoWCISA5EMFj8MvAG8CTQH/uTuT8U6mAQvOTGBO09vzuuDupKYYFwychr/+ngp+7MPBh1NRKIoksHiO4Al7n6Pu9/t7p9FenAzO9PMlplZhpndm8fz/cL3Jsw1s5lm1q2A+aUInFi/KuNu684l6XUZNn4F/YdOYblWPRMpNSK5NFQJ+MTMJpnZYDOrFcmBzSwRGEroPoRWwKVm1irXbl8A7d29A3AdMCry6FKUUssk8c8L2vHMVen8sH0vZz81mWcnr+TgwSMvbCQixV8kl4bud/fWhNYuPg6YYGafR3DsjkBGeKnL/cBYoF+uY+/0X5ZISyU8IC3FV59Wtfh4SA+6N6nBXz9YzFXPfsP32/YGHUtECqEgN4ZtBL4HfgRqRrB/bWBtjsfrwtv+h5mdZ2ZLgQ8JnRVIMZdWsQyjrk7n7+e1ZdbqrZzx+ETen7c+6FgicpQiGSO4yczGE7qMUwP4rbu3i+DYefU1/tU3fnd/291bAP2Bvx0mw8DwGMLMTZs2RfDWEmtmxmWd6jHu9u40rJHKra/MYcjYOWzbo2mmIiVNJGcE9YEh7t7a3f/s7pEuRrOOUGuKQ+oAh/3a6O4TgcZm9qvV1d19pLunu3t6WlpahG8vRaFhjVTeGNSFO05rxvvzN9D38YlMW/Fj0LFEpAAOWwjMrFL414eANWZWLedPBMeeATQ1s4ZmlgIMAN7L9R5NwovfYGYnACmELj1JCZKUmMDtpzXlzZu6UiY5kctGTefBDxezL/tA0NFEJAJH6j76MnA2oQVpnP+91ONAoyMd2N2zzewW4BMgEXjW3ReZ2aDw88OBC4CrzCwL2ANckmPwWEqYDnWr8OFt3fj7uCU8M2klk77dzOMDOtDimEr5v1hEAmMl7XM3PT3dZ86cmf+OEqivlm7knjfms31PFvec0ZzruzUkQcthigTGzGa5e3qezx2uEIQv1RyWu8+OQrYCUyEoOX7cuY8/vLWATxf/QOdG1Xjk4g7UrlIu6FgiceloC8FX4V/LAunAPEKXh9oBX7t7IHcBqxCULO7O6zPXcf/7i0gIr5Xcr8NxhIeGRKSIHKkQHHaw2N17u3tvYDVwQnjWzonA8UBGbKJKaWNmXHxSXT66vQfNalVkyKtzufWVOWzbrWmmIsVFJNNHW7j7gkMP3H0h0CF2kaQ0qle9PK/d2IV7zmjOxwu/54zHJzL5281BxxIRIisES8xslJn1MrOeZvYMsCTWwaT0SUwwBvduwts3n0xqmUSuGP0197+/iL1ZmmYqEqRICsG1wCLgdmAIsDi8TeSotK1TmQ9u7c41XRvw3JRVnPPUZBat17KYIkHR9FEJ1ITlm7jn9Xls3b2fO/s0Z2CPRiRqmqlI1B3VYLFIUejZLI1PhvSgT6ta/OvjpVw6cjprt+wOOpZIXFEhkMBVTU1h6GUn8MhF7Vm8YTt9n5jEG7PWUdLOVkVKKhUCKRbMjAtOrMNHt3en1bGVuPv1edz80my27tofdDSRUu9IvYYAMLNmwD2EupD+vL+7nxLDXBKn6lYrzysDO/PMpEwe+XQZs1Zv5aEL29GreSRLYIjI0ch3sNjM5gHDCTWf+3men7vPim20vGmwOH4sWr+NO16dy/IfdnJVl/r8oW9LyqUkBh1LpEQ60mBxvmcEQLa7D4tyJpF8tT6uMu/d0o1/f7KM0ZNXMiVjM49fcjxt61QOOppIqRLJGMH7ZnazmR1bwPUIRAqtbHIi953dipdu6MSufQc47+kpPPXFt2QdOBh0NJFSI5JLQyvz2OzufsT1CGJFl4bi17bdWfzx3YW8P289TWtW4M/ntKZb018taCcieTiq7qPFlQqBfLroex74cAlrtuymT6ta3HdWK+pVLx90LJFirVCFwMySgZuAHuFN44ER7h5I+0gVAgHYl32A0ZNX8p8vM8g+4NzQvSGDezchtUwkw14i8aewhWAUkAy8EN50JXDA3W+IasoIqRBITj9s38u/PlrKW3O+o1alMtzbtwX9O9TWegciuRS2EMxz9/b5bSsqKgSSl9lrtnL/e4uYt24bJ9Srwp/PaU37ulWCjiVSbBS219ABM2uc42CNyHE/gUhxcEK9qrx988n8+8J2rNmyh35Dp3DP6/PYuGNv0NFEir1ILqjeA3xlZpmElqqsj9pQSzGUkGBclF6XM9scw3++zODZKSv5aOH33HZqE67p2pCUJHVUEclLRLOGzKwM0JxQIVjq7vtiHexwdGlIIpW5aScPfriEL5ZupGGNVO47uyWntKgVdCyRQESjDfWJQBugPXCJmV0VrXAisdIorQKjrzmJ5689CTO47vmZXPPcN6zYtDPoaCLFSiSDxWOAxsBcfhkbcHe/LcbZ8qQzAjka+7MP8t9pq3ji82/Zk3WAa7o24LbTmlKpbHLQ0USKRGFnDS0BWnkxufNMhUAKY/POfTz8yTJenbmW6qkp3HNGcy46sS4JWhVNSrnCXhpaCBwT3UgiwahRoQz/vKAd7w3uRoPqqfz+zQX0GzqFmau2BB1NJDCRFIIawGIz+8TM3jv0E+tgIrHUtk5lXh/UhScGdGDTjn1cOHwat4+dw4Zte4KOJlLkIpk++pdYhxAJgpnRr0Nt+rSqxbDxKxgxMZNPF/3A4N6NuaF7I8oma+0DiQ9qOicStnbLbh78cAkfL/qeOlXL8cezWnJG62PUrkJKhUKNEZjZ+Wb2rZltM7PtZrbDzLZHP6ZIsOpWK8/wK0/k5Rs6kZqSxKAXZ3P5qK9Z+r3+3KV0i2TWUAZwjrsvKZpIR6YzAikK2QcO8vI3a3jk0+Xs2JvFFZ3rc2efZlQpnxJ0NJGjUthZQz8UlyIgUlSSEhO4qksDxt/diys61+fF6avp9fB4xkxbRbZWR5NSJpJCMNPMXjWzS8OXic43s/MjObiZnWlmy8wsw8zuzeP5y81sfvhnqpkF0tFU5HCqpqbw135tGHd7d1oeU4n73l3E2U9NZuqKzUFHE4maSApBJWA3cDpwTvjn7PxeZGaJwFCgL9AKuNTMWuXabSXQ093bAX8DRkYeXaTotDimEi//thPDLj+BHXuzueyZr7npxVms3bI76GgihZbv9FF3P9pOox2BDHfPBDCzsUA/YHGOY0/Nsf90oM5RvpdIzJkZfdseS+8WNXlmYiZPj1/Bl0s3cmOPRgzq1ZjyKVodTUqmfP9yzawscD3QGih7aLu7X5fPS2sDa3M8Xgd0OsL+1wMf5ZdHJGhlkxO59dSmXHBiHf750VKe/DKD12et496+LTi3/XGabiolTiSXhsYQajFxBjCB0Lf2HRG8Lq//G/KcomRmvQkVgt8f5vmBZjbTzGZu2rQpgrcWib3jqpTjyUuP5/VBXaiWmsLtY+dywbCpzFq9NehoIgUSSSFo4u73Abvc/QXgLKBtBK9bB9TN8bgOsD73TmbWDhgF9HP3H/M6kLuPdPd0d09PS0uL4K1Fis5JDarx3i3d+NcFbVm7dQ8XDJvK4Jdna/xASoxICkFW+J8/mVkboDLQIILXzQCamllDM0sBBgD/06PIzOoBbwFXuvvyiFOLFDOJCcYlJ9Vj/N29uO3Upnyx5AdOfWQC/xi3hG17svI/gEiAIikEI82sKnAfoQ/yxcBD+b3I3bOBW4BPgCXAa+6+yMwGmdmg8G5/AqoDT5vZXDPTnWJSoqWWSeLOPs0Yf3dvzu1wHCMnZdLr31/xwtRVZOn+Aymm1GtIJIYWfreNBz9cwrTMH2mUlsof+rbktJY1NaAsRa6wC9OUAS4gdDno51lG7v7XKGaMmAqBlDTuzhdLNvL3j5aQuWkXnRtV449ntaJN7cpBR5M4UtgWE+8Smv+fDezK8SMiETAzTmtVi0+G9OCv/Vqz/IednPOfydz12jy+37Y36HgiEZ0RLHT3NkWUJ186I5CSbvveLIZ+lcFzk1eRkAADuzfixp6NSS2jG9Ikdgp7RjDVzCKZLioiEahUNpk/9G3JF3f15LSWtXjyywx6PTyesd+s4cDBkjVmJ6XDYc8IzGwBoRvAkoCmQCawj9CNYh7uD1TkdEYgpc3sNVt54IPFzF7zEy2Oqcj/O6sl3ZvqfhmJrqMaLDaz+kc6qLuvjkK2AlMhkNLI3Rm34Hv++fES1m7ZQ6/mafzfb1rSrFbFoKNJKXG0haA8kOXuWeHHzYHfAKvd/a1Yhc2PCoGUZvuyD/DC1FU89WUGu/ZlM6BjPe44rRlpFcsEHU1KuKMdI/iY8B3EZtYEmAY0Agab2T+iHVJEoExSIgN7NGbCPb25qksDXpuxlt4Pj2foVxnszToQdDwppY44RuDubcO//w2o5u6Dw+0iZh16rqjpjEDiSeamnfzjo6V8tvgHjqtclnvObE6/9rVJSNANaVIwR3tGkLNCnAJ8BuDu+wHdKy9SBBqlVeCZq9J55bedqVYhhTtenUf/p6fwzcotQUeTUuRIhWC+mT1sZncATYBPAcysSpEkE5GfdWlcnfcGd+ORi9qzcfs+Lh4xjRvHzGTlZt3bKYV3pELwW2AzoXGC0939UE/dVsDDMc4lIrkkJBgXnFiHr+7uxV19mjHp282c/tgE/vr+Yn7avT/oeFKCFajpnJmd4O6zY5gnXxojEAnZuGMvj322nFdnrKVi2WRuPaUJV3VpQEpSJPeJSrwp7J3FOY2KQh4RiYKaFcvyj/PbMe727rSrU5kHPlxCn8cm8NGCDZS0rsISrIIWAk1VEClmWhxTiTHXd+L5a0+iTFICN700m4tHTGPu2p+CjiYlREELwf0xSSEihdareU3G3dadv5/XlpWbd9F/6BRuHzuHdVu1ZKYcWURjBGZWG6jP/65HMDGGuQ5LYwQi+du5L5vh41fwzKRMHLiqc30G9WpMjQq6QzleFXZhmn8BlxBaovLQrY3u7udGNWWEVAhEIrf+pz088uly3p6zjjJJiVzdtQE39mhE1dSUoKNJEStsIVgGtHP3fbEIV1AqBCIFl7lpJ0988S3vzVtP+eREruvWkBu6NaJy+eSgo0kRKeysoUxAfy0iJVijtAo8MeB4PhnSg57N03jqywy6PfQlT3z+LTv2ZgUdTwIWyRnBm0B74AtC6xEA4O63xTZa3nRGIFJ4i9dv57HPl/PZ4h+oUj6ZgT0acXWXBlolrRQr7KWhq/Pa7u4vRCFbgakQiETPgnXbePSzZXy1bBPVU1MY1LMxV3SuT7mUxKCjSZQVqhCED1AOqOfuy6IdrqBUCESib/aarTz22XImfbuZtIpluLlXYy7tWI+yySoIpUWhxgjM7BxgLqH1CTCzDmb2XnQjikiQTqhXlTHXd+K1G7vQOC2V+99fTK9/j2fM9NXsz1az4dIuksHivwAdgZ8A3H0u0DCGmUQkIB0bVmPswC68fEMnalctx33vLKT3w+MZ+80asg6oIJRWkRSCbHfflmubGpmIlGJdm9TgjUFdeOG6jtSoWIZ731rAqY9M4M1Z68hWQSh1IikEC83sMiDRzJqa2VPA1BjnEpGAmRk9m6Xxzs1dGX11OhXLJnHX6/M4/fGJvDv3Ow4c1PfB0iKSQnAr0JrQ1NFXgO3AkFiGEpHiw8w4tWUtPri1G8OvOJGUxARuHzuXvk9MZNyCDRxUQSjxCroeQSKQ6u7bYxfpyDRrSCRYBw864xZu4PHPvyVj405aHluJO05rSp9WtTBTg+LiqrCzhl42s0pmlgosApaZ2T3RDikiJUNCgnF2u+P4ZEgPHr+kA3uzDjBwzCzO/c8Uvlq6UWshlECRXBpqFT4D6A+MA+oBV8Y0lYgUe4kJRv/ja/PZHT3494Xt+GnPfq59fgbnD5vK5G83qyCUIJEUgmQzSyZUCN519yw0a0hEwpISE7govS5f3NmLv5/Xlh+27eWK0V9zycjpTM/8Meh4EoFICsEIYBWQCkw0s/qEBozzZWZnmtkyM8sws3vzeL6FmU0zs31mdndBgotI8ZKSlMBlnerx1T29+Gu/1qzavIsBI6dz+ajpzFq9Jeh4cgQFGiz++UVmSe6enc8+icByoA+wDpgBXOrui3PsU5PQgjf9ga3u/nB+763BYpGSYW/WAV76eg3Dxmeweed+ejZL444+zehQt0rQ0eLSkQaL8201aGZlgAuABrn2/2s+L+0IZLh7Zvg4Y4F+hBa4AcDdNwIbzeys/HKISMlSNjmR67s15NKOdRkzbTXDJ6yg/9ApnNayJkNOa0ab2pWDjihhkVwaepfQB3g2sCvHT35qA2tzPF4X3iYicaR8ShI39mzMpN+fwj1nNGfGqq2c/dRkBo2ZxcLvcjctkCBE0ny8jrufeRTHzmtC8VENMpvZQGAgQL169Y7mECISsAplkhjcuwlXdqnPs5NXMnrSSj5e9D0nN6nOwB6N6dG0hu5DCEgkZwRTzaztURx7HVA3x+M6wPqjOA7uPtLd0909PS0t7WgOISLFRKWyyQw5rRlT/nAKf+jbgoyNO7n62W/o+8Qk3p6zTs3tAhBJIegGzArP/plvZgvMbH4Er5sBNDWzhmaWAgwA1L5aRIBQQbixZ2Mm/e4U/n1hOw66c8er8+j50FeMmpTJzn1HnI8iURTJCmX18+U7IikAAA00SURBVNru7qvzPbjZb4DHgUTgWXd/0MwGhV8/3MyOAWYClYCDwE5+uYEtT5o1JFI6uTvjl21ixMQVTM/cQsWySVzeqT7XntyAWpXKBh2vxIvGCmXdgKbu/pyZpQEV3H1llHNGRIVApPSbv+4nRkzM5KMFG0J3MHeozcAejWhaq2LQ0Uqswq5Z/GcgHWju7s3M7DjgdXc/OfpR86dCIBI/1vy4m1GTM3lt5lr2Zh3klBY1GdijEZ0aVtPAcgEVthDMBY4HZrv78eFt8929XdSTRkCFQCT+bNm1nzHTVvPCtFVs2bWf9nWrcGOPRpzR+hgSE1QQIlGo7qPAfg9VCw8fLDWa4URE8lMtNYXbT2vK1HtP4YH+bfhp935ufmk2pzwynjHTVrFn/4GgI5ZokZwR3A00JdQq4h/A9cDL7v5k7OP9ms4IROTAQefTRd8zYmImc9f+RLXUFK7sXJ+rutSneoUyQccrlqIxWNwHOD388BN3/zyK+QpEhUBEDnF3ZqzaysiJK/h8yUbKJidw4Yl1uKFbIxrU0MWLnI6qEJjZDn65Ezj3Rbi9wArg/7n7F9EKGgkVAhHJS8bGHTwzcSVvz/mOrIMH6dvmGAb2aKwmd2GFPiPI44CJQBvgJXdvU8h8BaJCICJHsnH7Xp6buooXp69mx95sOjasxo09GtG7eU0S4nhgOeqFIMeBb3T3EUd9gKOgQiAikdi5L5ux36zh2ckrWb9tL01qVmBg90b0O/44yiQlBh2vyMWsEARBhUBECiLrwEE+nL+BERMzWbJhOzUrluHakxtyWad6VC6XHHS8IqNCICJxz92ZnLGZERMymZyxmdSURC7tWI/rujXkuCrlgo4XcyoEIiI5LPxuG89MyuSD+Rsw4Jz2xzGwRyNaHlsp6Ggxo0IgIpKHdVt38+zkVYydsYbd+w/QvWkNbuzRmJObVC91LSxUCEREjmDb7ixe/Ho1z01Zxead++jSqDp/69+GJjUrBB0talQIREQisDfrAK/OWMsjny5jT9YBBvZoxC29m1IupeTPMipsryERkbhQNjmRq7s24Iu7enFO++MY+tUK+jw2gS+W/BB0tJhSIRARySWtYhkevbgDYwd2plxyIte/MJPf/ncm67buDjpaTKgQiIgcRudG1fnwtu7c27cFk7/dTJ9HJzJs/Ar2Z5eudZVVCEREjiAlKYFBPRvz2Z096N60Bv/6eClnPTmJ6Zk/Bh0talQIREQiUKdqeUZelc7oq9PZk3WAASOnc+erc9m0Y1/Q0QpNhUBEpABObVmLz+7oyS29m/D+/PWc+sh4xkxfzYGDJWsGZk4qBCIiBVQuJZG7z2jOR7f3oE3tytz3zkLOf3oKC9ZtCzraUVEhEBE5Sk1qVuClGzrxxIAOfPfTXs4dOpk/vbuQbXuygo5WICoEIiKFYGb061CbL+/uydVdGvDi9NWc+sgE3pnzHSXlhl0VAhGRKKhUNpm/nNuadwd3o3aVsgx5dS6XPfM1GRt3BB0tXyoEIiJR1LZOZd66+WQe6N+GReu30feJSTz08VL27D8QdLTDUiEQEYmyxATjis71+fLuXpzbvjZPj1/BaY9O4PPFxbNVhQqBiEiM1KhQhkcubs+rAztTPiWRG/47kxteKH6tKlQIRERirFOj6oy7PdSqYkrGZk57dAJPj88oNq0qVAhERIpAcmKoVcXnd/WkZ7M0Hvp4Gb95chLTVgTfqkKFQESkCNWuUo4RV4ZaVezNOsClz0znjoBbVagQiIgEIGerig/mr+eUR8YzZtqqQFpVqBCIiATkUKuKj4f0oG3tytz37iLOe3oK89f9VKQ5YloIzOxMM1tmZhlmdm8ez5uZPRl+fr6ZnRDLPCIixVHjtF9aVWzYtpd+Q6dw3ztF16oiZoXAzBKBoUBfoBVwqZm1yrVbX6Bp+GcgMCxWeUREirNDrSq+uCvUquKlr1dz6iPjeXvOupi3qojlGUFHIMPdM919PzAW6Jdrn37Afz1kOlDFzI6NYSYRkWLtUKuK927pRu2q5bnj1Xlc+sz0mLaqiGUhqA2szfF4XXhbQfcREYk7bWpX5q2buvLgeW1YvH47fZ+YxKhJmTF5r6SYHDXE8tiW+/wmkn0ws4GELh1Rr169wicTESkBEhOMyzvV54zWx/CPcUupXz01Ju8Ty0KwDqib43EdYP1R7IO7jwRGAqSnp5eMvq4iIlFyqFVFrMTy0tAMoKmZNTSzFGAA8F6ufd4DrgrPHuoMbHP3DTHMJCIiucTsjMDds83sFuATIBF41t0Xmdmg8PPDgXHAb4AMYDdwbazyiIhI3mJ5aQh3H0fowz7ntuE5fndgcCwziIjIkenOYhGROKdCICIS51QIRETinAqBiEicUyEQEYlzFutmRtFmZpuA1Uf58hrA5ijGiZbimguKbzblKhjlKpjSmKu+u6fl9USJKwSFYWYz3T096By5FddcUHyzKVfBKFfBxFsuXRoSEYlzKgQiInEu3grByKADHEZxzQXFN5tyFYxyFUxc5YqrMQIREfm1eDsjEBGRXOKmEJjZmWa2zMwyzOzeoPMAmNmzZrbRzBYGnSUnM6trZl+Z2RIzW2RmtwedCcDMyprZN2Y2L5zr/qAz5WRmiWY2x8w+CDrLIWa2yswWmNlcM5sZdJ5DzKyKmb1hZkvDf2ddikGm5uH/Tod+tpvZkKBzAZjZHeG/+YVm9oqZlY3q8ePh0pCZJQLLgT6EFsOZAVzq7osDztUD2Elo3eY2QWbJKbxu9LHuPtvMKgKzgP7F4L+XAanuvtPMkoHJwO3h9a4DZ2Z3AulAJXc/O+g8ECoEQLq7F6s58Wb2AjDJ3UeF1ysp7+4/BZ3rkPBnxndAJ3c/2vuWopWlNqG/9VbuvsfMXgPGufvz0XqPeDkj6AhkuHumu+8HxgL9As6Eu08EtgSdIzd33+Dus8O/7wCWUAzWkvaQneGHyeGfYvFNxszqAGcBo4LOUtyZWSWgBzAawN33F6ciEHYqsCLoIpBDElDOzJKA8uSxkmNhxEshqA2szfF4HcXgg60kMLMGwPHA18EmCQlffpkLbAQ+c/dikQt4HPgdcDDoILk48KmZzQqv/V0cNAI2Ac+FL6WNMrPYLMZ79AYArwQdAsDdvwMeBtYAGwit5PhpNN8jXgqB5bGtWHyTLM7MrALwJjDE3bcHnQfA3Q+4ewdC61t3NLPAL6mZ2dnARnefFXSWPJzs7icAfYHB4cuRQUsCTgCGufvxwC6gWIzbAYQvVZ0LvB50FgAzq0roCkZD4Dgg1cyuiOZ7xEshWAfUzfG4DlE+tSptwtfg3wRecve3gs6TW/hSwnjgzICjAJwMnBu+Hj8WOMXMXgw2Uoi7rw//cyPwNqHLpEFbB6zLcTb3BqHCUFz0BWa7+w9BBwk7DVjp7pvcPQt4C+gazTeIl0IwA2hqZg3D1X4A8F7AmYqt8KDsaGCJuz8adJ5DzCzNzKqEfy9H6H+QpcGmAnf/g7vXcfcGhP62vnT3qH5jOxpmlhoe7Cd86eV0IPAZau7+PbDWzJqHN50KBDoRIZdLKSaXhcLWAJ3NrHz4/81TCY3bRU1M1ywuLtw928xuAT4BEoFn3X1RwLEws1eAXkANM1sH/NndRwebCgh9w70SWBC+Hg/wf+E1qIN0LPBCeEZHAvCauxebqZrFUC3g7dBnB0nAy+7+cbCRfnYr8FL4i1kmcG3AeQAws/KEZhfeGHSWQ9z9azN7A5gNZANziPIdxnExfVRERA4vXi4NiYjIYagQiIjEORUCEZE4p0IgIhLnVAhEROKcCoHENTN7LGeHSTP7xMxG5Xj8iJn9qaAda83seTO7MJpZRWJFhUDi3VTCd2maWQJQA2id4/muwCfu/s8AsokUCRUCiXdT+OV2/daE7rzdYWZVzawM0BJob2b/gZ+/6T9pZlPNLPPQt34L+Y+ZLTazD4Gah97AzE4NN1dbEF6DooyZdTSzt8LP9zOzPWaWEl5zIbMI//1FVAgkvoV78WSbWT1CBWEaoU6rXQitLTAf2J/rZccC3YCzgUNnCucBzYG2wG/55SyjLPA8cIm7tyV0h+9NhO4SPT782u6ECtBJQCeKSadXiR8qBCK/nBUcKgTTcjyemsf+77j7wfBCPbXC23oAr4S7o64Hvgxvb06oYdjy8OMXgB7ung1kmFlLQo3gHg0fozswKdr/giJHokIg8ss4QVtC38ynEzoj6EqoSOS2L8fvOVuc59WvJa8W6IdMItTpMgv4nNBZRjdgYqTBRaJBhUAk9GF/NrAl/I1+C1CFUDGYFuExJgIDwgvnHAv0Dm9fCjQwsybhx1cCE3K8Zggwzd03AdWBFkDgDRElvsRF91GRfCwgNFvo5VzbKrj75nD3zvy8DZwSft1ywh/27r7XzK4FXg8vMzgDGB5+zdeELi0dOgOYT2iBG3WClCKl7qMiInFOl4ZEROKcCoGISJxTIRARiXMqBCIicU6FQEQkzqkQiIjEORUCEZE4p0IgIhLn/j+pjS1VqW70PAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ces_fig, ces_ax = plt.subplots()\n", "plt.plot(ces_conv)\n", "ces_ax.set_xlabel('Window')\n", "ces_ax.set_ylabel('Jensen-Shannon divergence')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparing different clustering methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You may want to try different clustering methods, or use different parameters within the methods. `encore.ces_convergence` allows you to pass a list of `clustering_methods` to be applied, much like [normal clustering ensemble similarity methods](clustering_ensemble_similarity.ipynb#Calculating-clustering-similarity-with-multiple-methods).\n", "\n", "
\n", " \n", "**Note**\n", "\n", "To use the other ENCORE methods available, you need to install [scikit-learn](https://scikit-learn.org/stable/).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The KMeans clustering algorithm separates samples into $n$ groups of equal variance, with centroids that minimise the inertia. You must choose how many clusters to partition. [(See the scikit-learn user guide for more information.)](https://scikit-learn.org/stable/modules/clustering.html#k-means)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "km1 = clm.KMeans(12, # no. clusters\n", " init = 'k-means++', # default\n", " algorithm=\"auto\") # default\n", "\n", "km2 = clm.KMeans(6, # no. clusters\n", " init = 'k-means++', # default\n", " algorithm=\"auto\") # default\n", "\n", "km3 = clm.KMeans(3, # no. clusters\n", " init = 'k-means++', # default\n", " algorithm=\"auto\") # default" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we pass a list of clustering methods to `encore.ces_convergence`, the similarity values get saved in `ces_conv2` in order." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(9, 3)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ces_conv2 = encore.ces_convergence(u, # universe\n", " 10, # window size\n", " select='name CA',\n", " clustering_method=[km1, km2, km3]\n", " )\n", "ces_conv2.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the number of clusters partitioned by KMeans has an effect on the resulting rate of convergence." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1yVZRvA8d/NUARBcSsgoLgRcOPeqJWmpVlqpWWOssyWVm+ZlWVlqZlpaWllVqZpy733QsE9QEFxgrKRfb9/HEQMxMM4HJDr+/nw0fOcZ1zHt5fr3M99P9eltNYIIYQovSzMHYAQQgjzkkQghBClnCQCIYQo5SQRCCFEKSeJQAghSjkrcweQV1WqVNFubm7mDkMIIUoUf3//CK111ZzeK3GJwM3NjQMHDpg7DCGEKFGUUqF3e09uDQkhRCkniUAIIUo5SQRCCFHKlbg5AiFEyZSSkkJYWBiJiYnmDuW+ZmNjg7OzM9bW1kYfI4lACFEkwsLCsLe3x83NDaWUucO5L2mtuX79OmFhYbi7uxt9nNwaEkIUicTERCpXrixJwISUUlSuXDnPoy5JBEKIIiNJwPTy829cahLBlehEpvx9jJS0dHOHIoQQxUqpSQQBF6JYuDOELzeeMXcoQggzeeaZZ6hWrRqenp53bH/99ddp2LAhXl5eDBgwgKioKKPPOXz4cJYtW5bnWEJCQliyZEmejzOFUpMIenvW4LGWzszZHMSBkBvmDkcIYQbDhw9nzZo12bb37NmTo0ePcvjwYerXr8/HH39s8ljykwjS0tJMEkupSQQA7/ZtgrOjLS//FkBsYoq5wxFCFLFOnTpRqVKlbNv9/PywsjIsovT19SUsLCzH4z/99FOaNm2Kt7c3kyZNyva+m5sbERERABw4cIAuXboAsHXrVnx8fPDx8aFZs2bExsYyadIktm/fjo+PDzNmzCAtLY3XX3+dVq1a4eXlxTfffAPAli1b6Nq1K0OGDKFp06bEx8fz4IMP4u3tjaenJ7/99luB/11K1fLR8mWtmDHYh0HzdvHeX8f5/DFvc4ckRKk05e9jHL8UU6jnbFzLgcl9mxT4PN9//z2DBw/Otn316tWsXLmSvXv3Ymtry40bxt9ZmD59OnPmzKF9+/bExcVhY2PDtGnTmD59Ov/88w8A3377LRUqVGD//v0kJSXRvn17/Pz8ANi3bx9Hjx7F3d2d5cuXU6tWLf79918AoqOjC/yZS9WIAKCFqyPjutVj+cEw/j182dzhCCGKkalTp2JlZcXQoUOzvbdhwwZGjBiBra0tQI4ji7tp3749r7zyCl9++SVRUVGZo4+s1q1bx48//oiPjw9t2rTh+vXrnDljmNNs3bp15nMBTZs2ZcOGDUycOJHt27dToUKF/HzUO5SqEcEtL3bzYNvpcN5acYTmrhWpWaGcuUMSolQpjG/uhe2HH37gn3/+YePGjTkuwdRa33NpppWVFenphpWJWdfyT5o0iQcffJBVq1bh6+vLhg0bcjz/7Nmz6dWr1x3bt2zZgp2dXebr+vXr4+/vz6pVq3jzzTfx8/Pj3XffzdNn/S+TjgiUUr2VUqeUUkFKqWw31JRSXZRS0UqpgIyfgn0aI1lbWjBjsA8paem89nsg6em6KC4rhCim1qxZwyeffMJff/2V+Y3/v/z8/Pj+++9JSEgAyPHWkJubG/7+/gAsX748c3twcDBNmzZl4sSJtGzZkpMnT2Jvb09sbGzmPr169WLu3LmkpBjmL0+fPk18fHy2a1y6dAlbW1uGDRvGa6+9xsGDB/P/wTOYLBEopSyBOUAfoDHwhFKqcQ67btda+2T8vG+qeP7LvYod7z7UmJ1B1/l+57miuqwQwoyeeOIJ2rZty6lTp3B2dua7774DYNy4ccTGxtKzZ098fHwYM2ZMtmN79+5Nv379aNmyJT4+PkyfPj3bPpMnT2b8+PF07NgRS0vLzO0zZ87E09MTb29vypUrR58+ffDy8sLKygpvb29mzJjByJEjady4Mc2bN8fT05PRo0eTmpqa7RpHjhyhdevW+Pj4MHXqVP73v/8V+N9FaW2ab8NKqbbAe1rrXhmv3wTQWn+cZZ8uwGta64eMPW/Lli11vhrTXD4MG6fAwIVg40BGLIz6yZ+tp8L568X2NKzhkPfzCiGMcuLECRo1amTuMEqFnP6tlVL+WuuWOe1vyltDTsCFLK/DMrb9V1ulVKBSarVSKscbh0qpUUqpA0qpA+Hh4fmLJjkegjfD3y9BRvJTSvHJo15UsLVm/C8BJKaYZo2uEEIUZ6ZMBDnNqvx3+HEQcNVaewOzgZU5nUhr/a3WuqXWumXVqjm23Lw317bQ/R04tgIOfJe5uZJdGT4b6MWpq7F8tvZU/s4thBAlmCkTQRjgkuW1M3Ap6w5a6xitdVzG31cB1kqpKiaLqN148OgJa96Ey4GZm7s0qMbTbV35bsc5dpyJMNnlhRCiODJlItgP1FNKuSulygCPA39l3UEpVUNlrMdSSrXOiOe6ySKysIAB34BdVVj6NCTefqDlzQca4VGtPK/+HkBkfLLJQhBCiOLGZIlAa50KjAPWAieApVrrY0qpMUqpW1PyA4GjSqlA4EvgcW2q2etb7CrDwO8h6jz89WLmfIGNtSUzB/twIz6Zt1YcwdRhCCFEcWHS5wi01qu01vW11nW11lMzts3TWs/L+PtXWusmWmtvrbWv1nqXKePJVNsXur8Lx1fC/gWZmz2dKvCqXwNWH73C8oMXiyQUIYQwt1JXYiJTu5egnh+sfQsuBWRufq5jHdq4V2Lyn0c5fz3BjAEKIQpbVFQUAwcOpGHDhjRq1Ijdu3cbfWz58uXzdc2VK1dy/PjxfB1bVEpvIrCwgP7zDPMFvw+HREPhJksLxReDfbCwUExYGkCqNLIR4r4xfvx4evfuzcmTJwkMDCyS5xrykwhyepDMlEpvIoCM+YKF2eYLnCqW48P+nviHRjJ3S7CZgxRCFIaYmBi2bdvGs88+C0CZMmWoWLFitv2uXr3KgAED8Pb2xtvbm1277rxjvWXLFh566PYzsOPGjWPRokWAoaZQ48aN8fLy4rXXXmPXrl389ddfvP766/j4+BAcHExwcDC9e/emRYsWdOzYkZMnTwKGXgmvvPIKXbt2ZeLEiTmWrjaVUll07g6120CPybD+XcN8QevnAHjYx4lNJ68xc+MZOtavio9L9v9ghBD5tHoSXDlSuOes0RT6TLvr22fPnqVq1aqMGDGCwMBAWrRowaxZs+4o6Abw0ksv0blzZ1asWEFaWhpxcXFGXf7GjRusWLGCkydPopQiKiqKihUr0q9fPx566CEGDhwIQPfu3Zk3bx716tVj7969PP/882zatAkw1BfasGEDlpaW9O3bN1vpalMp3SOCW9q+CPV6ZcwXHMrc/P7DntRwsGHCbwHEJxXtUE0IUbhSU1M5ePAgY8eO5dChQ9jZ2TFtWvbEsWnTJsaOHQuApaWl0WWeHRwcsLGxYeTIkfzxxx85Fq+Li4tj165dDBo0CB8fH0aPHs3ly7fL4Q8aNCizRpExpasLi4wIIOP5gnkwr4NhvmD0NrCpQIVy1nz+mDdPzN/Dh/+e4ONHmpo7UiHuD7l8czcVZ2dnnJ2dadOmDQADBw7MMRHcS9ZS03C73LSVlRX79u1j48aN/Prrr3z11VeZ3/RvSU9Pp2LFigQEBJCTrKOTnEpXN2zYMM/xGkNGBLfYVjLMF0SHwZ/jMucLfOtUZnSnuvyy7zzrj181c5BCiPyqUaMGLi4unDplKCWzceNGGjfOXhC5e/fuzJ07FzD0CI6JubOTmqurK8ePHycpKYno6Gg2btwIGL7tR0dH88ADDzBz5szMX/ZZy007ODjg7u7O77//DhgKXwYGBpKTnEpXm4okgqxqt4Huk+HEX7BvfubmV3rWp3FNByYuP8y12MRcTiCEKM5mz57N0KFD8fLyIiAggLfeeivbPrNmzWLz5s00bdqUFi1acOzYsTved3Fx4bHHHsPLy4uhQ4fSrFkzAGJjY3nooYfw8vKic+fOzJgxA4DHH3+czz77jGbNmhEcHMzPP//Md999h7e3N02aNOHPP//MMdacSlebisnKUJtKvstQGys9HX59AoI3wTNrwak5AEHXYnnwyx20q1uZ74e3umenIiHEnaQMddEpTmWoSyYLC+g/F+yqGeYLbkYB4FHNnrcfbMTmU+Es3hNq3hiFEKIQSSLIiW0lGLQQYi7CX7fnC570daVz/ap8+O8Jgq6Zbk2vEEIUJUkEd+PSOmO+4G/Y9y1gaGTz2SAv7Mpa8fJvASSnylPHQoiSTxJBbtq9CPX7wNq34aKhQXQ1exumPdKUoxdjmLHhtJkDFEKIgpNEkBuloP/XYF/jjvkCvyY1eKK1C/O2BrP3rOnaJwghRFGQRHAvtpUM/QtiLsKfL2TOF/zvwca4VrLllaWBxCSmmDlIIYTIP0kExnBpDT3eg5P/wN5vALAra8WMwT5ciUnk3ZVHzRqeEOLeEhMTad26deb6/cmTJ+fpeClDLaDtOMN8wbr/wUV/AJrVdmR893qsDLjEnwHSyEaI4qxs2bJs2rSJwMBAAgICWLNmDXv27DH5de+LMtRKqfpKqY1KqaMZr72UUv8zfWjFzF3mC57vUpfmtSvyv5VHuRh107wxCiHuSimV+a0+JSWFlJSUHB8MlTLUOZsPvA58A6C1PqyUWgJ8aLKoiqtb9YgW9jbMFwxejJWlBTMG+/DArO28ujSAJSN9sbCQp46FyM0n+z7h5I3CrZ3TsFJDJraemOs+aWlptGjRgqCgIF544YXMAnRZSRnqnNlqrff9Z1vprcns0gp6TMmYL5gHgGtlOyb3a8KeszeYv/2smQMUQtyNpaUlAQEBhIWFsW/fPo4ezT6/J2WocxahlKoLaACl1EDgcu6H3OfavgChO2HdO+DcGpxbMKiFM5tPXmP6ulN0qFeFJrWM+49HiNLoXt/cTa1ixYp06dKFNWvW4OnpmadjS2sZ6hcw3BZqqJS6CLwMjDVJNCWFUvDwHLCvCcuGw81IlFJ8NKApjrZlePnXABJT0swdpRAii/DwcKKiDHN7N2/evOsvVilDnQOt9VmtdQ+gKtBQa91Bax1isohKisx6RJdgpeH5Ake7Mkwf5M2Za3FMW226/9GEEHl3+fJlunbtipeXF61ataJnz553TPreImWoc9pBqY+AT7XWURmvHYFXtdZmWTlk8jLUebV7jqHFZa+Poe3zAEz5+xgLd4awaEQrujSoZuYAhSgepAx10TFFGeo+t5IAgNY6EnigQFHeT3yfhwYPwvp3IczwfMHE3g2pX708ry87zI34ZDMHKIQQuTMmEVgqpcreeqGUKgeUzWX/0kUp6J8xX/D7cLgZiY21JTMHNyM6IYVJyw9T0pr/CCFKF2MSwWJgo1LqWaXUM8B64AfThlXClHOEQYsg9nLmfEHjWg683qsB645fZemBC+aOUIhiQb4UmV5+/o2NmSz+FJgKNAKaAB9kbBNZObeAnu/DqX9hz9cAPNvBnXZ1KzPl7+OERMSbOUAhzMvGxobr169LMjAhrTXXr1/P88Nn0rO4MGkNvw6FM2sN/Y6dW3I5+ia9ZmyjTtXyLBvTFitLKe8kSqeUlBTCwsIy190L07CxscHZ2Rlra+s7tuc2WWzMqqFHgE+AaoDK+NFaa4d7BaSU6g3MAiyBBVrraXfZrxWwBxistV6W2zmLdSIAuBkJ33QyPH43eivYVuKfw5cYt+QQ47vXY0LP+uaOUAhRChV01dCnQD+tdQWttYPW2t7IJGAJzAH6AI2BJ5RSje+y3yfAWiNiKf7KOcLARYb5goz+BQ951eKRZk58tTkI/9BIc0cohBB3MCYRXNVan8jHuVsDQRkPpCUDvwIP57Dfi8By4Fo+rlE8ObcAvw/g1CrDcwbAlIebULOCDa8sDSAuqfSWahJCFD/GJIIDSqnflFJPKKUeufVjxHFOQNblMmEZ2zIppZyAAcC83E6klBqllDqglDoQHh5uxKWLgTZjoOFDsGEyXNiPvY01Xzzmw4UbCXzwd/FuUiGEKF2MSQQOQALgB/TN+Mn+XHZ2OdVi/u+ExExgotY618I8WutvtdYttdYtq1atasSliwGl4OGvwKEWLBsBCTdo7V6JsV3q8tuBC6w5esXcEQohBGBE9VGt9Yh8njsMcMny2hm49J99WgK/ZjSHqAI8oJRK1VqvzOc1i5dbzxd81wtWPg9P/ML47vXZdjqCSX8cplntilR3MF2NcSGEMIYpO5TtB+oppdyVUmWAx4G/su6gtXbXWrtprd2AZcDz900SuMWpBfh9CKdXw+6vKGNlwczHfUhMSeO13wNJTy9Zy3eFEPcfY24NzQfeBFLA0KEMwy/1XGmtU4FxGFYDnQCWaq2PKaXGKKXG5D/kEqjN6Iz5gvfgwj7qVi3P/x5szPYzEfywO8TMwQkhSjuTdijTWq/SWtfXWtfVWk/N2DZPa51tclhrPfxezxCUWLf6FzjUgt8N8wVD29SmW8NqfLz6JKevmq4XqRBC3IsxiUA6lBWGchUN8wVxV2HlWBTwyaNe2Je1YvyvAcQmppg7QiFEKSUdyoqSUwvoNRVOr4Fds6lqX5bPBnlx6koMfWZtZ+/Z6+aOUAhRCkmHsqLWehQ06muYLzi/l24Nq/P7mHZYWSgen7+Hqf8elzaXQogiZUytoVdy2BwN+Gutc+7AbELFvtaQMW5GGeoRpafBmO1gW4mE5FQ+WnWCxXvO06C6PV8M9qZJrQrmjlQIcZ8oaK2hlsAYDE8FOwGjgC7AfKXUG4UVZKlya74g/hqsGAPp6diWseLD/k1ZNKIVkQnJ9J+zkzmbg0iT5aVCCBMzJhFUBpprrV/VWr+KITFUBToBw00Y2/3NqTn4TTWUrN49O3NzlwbVWPtyJ/ya1OCztad47Jvd0stACGFSxiSC2kDWxrspgKvW+iaQZJKoSovWz0GjfrBhChxemrnZ0a4MXz3RjFmP+3DmaiwPfLmdn/eGSkMPIYRJ3LPEBLAE2KOU+jPjdV/gF6WUHSDV0wriVj2ihBvwx3MQcRq6vAUWFiileNjHidbulXj998O8veIoG45f5ZNHvagmZSmEEIUo18liZSgC5IyhKU0HDIXkdmitzTZbe19MFv9XajL8+woc+gkaPwz950EZ28y309M1P+0J5ePVJyhnbcnUAU15oGlNMwYshChpCtqhzF9r3cIkkeXDfZkIwNDmcvdXsO4dqOUDj/8CDnf+sg8Oj+OV3wIIDIumv08tpjzsSYVy1nc5oRBC3FbQRDAHWKS13m+K4PIqv4kgJS2FhNQEE0RUcNYW1thaZ4wATq2GZc+CTQV44hdDUsgiJS2dOZuDmL0piGr2ZZk+yJv2HlXMELUQoiQpaCI4DjQAQoB4bvcs9irkOI2S30SwNmQtr219zQQRFZyVsuKl5i8xvMlwlFJw5QgseRxu3oAB30DjftmOCbwQxYSlAZwNj2d4Ozcm9WmIjbWlGaIXQpQEBU0Erjlt11qHFkJseZbfRBAaE8qOiztMEFHB7b+yn43nN9LLrRfvt3vfMDqIvQq/DYWw/dD9XejwimFyOYubyWl8suYki3aFULeqHV885oO3S0UzfQohRHFWoESQcYIOQD2t9UKlVFWgvNb6XCHHaZT7cY5Aa82iY4uYeXAmdSrUYVbXWdR2qA0pifDnC3B0GXg/AX1ngVXZbMfvOBPB68sCuRabxIvdPHihqwfWlsasDBZClBYFerJYKTUZmIihJwGANbC48MITSilGeI5gbo+5hN8M5/F/H2db2DawtoFHF0DXtyHwF/ihH8RHZDu+Q70qrHm5E329ajJzwxkGzt1FcHicGT6JEKIkMuZr4wCgH4b5AbTWlwB7UwZVWrWr1Y7fHvoNp/JOjNs4jnmB80hHQ+c3YOBCuBwA87vCtRPZjq1QzpqZjzdjzpDmhN5I4MEvt/PDrhDpgCaEuCdjEkGyNtw/utWPwM60IZVuTuWd+LHPjzxY50HmBMxh/ObxxCbHgucjMHwVpCbBgp5wZn2Oxz/oVZN1L3fCt05lJv91jKe+38fl6JtF/CmEECWJMYlgqVLqG6CiUuo5YAOG9pXCRMpZleOjDh8xqfUkdoTtYMi/QzgbdRacW8Bzm6GSGyx5DPbMMzx/8B/VHGxYOLwVUwd44h8aSa8Z2/gz4KKUqBBC5MjYyeKegB+GpaNrtdY5fx0tAvfjZHFuDlw5wKtbXyUxNZGpHabSw7UHJMXBitFw8h9o+Qz0+RQsc36wLCQinleWBnDwfBQPetXkw4c9cbQrU8SfQghhbgVdPjoB+F1rHWaK4PKqtCUCgCvxV3h1y6scjjjMyKYjGeczDksUbJwCO2eCe2d47Aco55jj8alp6Xyz7Swz1p+mkl0ZPh3oRZcG1Yr4UwghzKmg/QgcgLVKqe1KqReUUtULNzxxLzXsarCw90IerfcoC44s4IWNLxCdEgs9p8DDX0PoLljQA64H53i8laUFL3T1YOUL7aloa83whft5e8UREpJTi/iTCCGKI6NuDQEopbyAwcCjQFhG+8oiVxpHBFktO72Mj/Z+RDXbaszqOosGlRoYEsGvQ0Gnw+DF4N7xrscnpqTx+bpTLNhxDtdKtnz+mA8tXHMeSQgh7h8FHRHccg24AlzHUI1UmMHA+gNZ2HshKekpDFs1jFVnV4FrO3huE5SvDj/1B/8f7nq8jbUlbz/YmCUjfUlJ0wyat4vpa0+RnJpehJ9CCFGcGPNA2Vil1BZgI1AFeM5cdYaEgXdVb3576DcaV27MxO0T+XT/p6RWdIGR6w3zBX+/BGvfNvREvou2dSuz5uWOPNLcma82BzHg652cvhpbhJ9CCFFcGDMicAVe1lo30VpP1lpLM5pioEq5KizotYChjYby0/GfGLV+FNd1KgxZCq1HG0pa/zoEku7+y93exprpg7yZN6wFl6MTeWj2DhZsPysPoQlRytx1jkAp5aC1jlFKVcrpfa31DZNGdhelfY4gJ38H/82U3VNwtHFkRpcZeFbxhP0LYNUbULUhDPkVKtbO9RzhsUm8+cdhNpy4hm+dSkwf5I2zo22uxwghSo78zhEsyfjTHziQ8ad/lteimOhbty8/9vkRCyx4evXTrDizAlqNhGHLIDoM5neDC7m3k6hqX5b5T7Xk00e9OBIWTZ+Z21nmHyYPoQlRChi9aqi4kBHB3UUmRvLGtjfYc3kPgxsMZmKriVjfOGd4CjnmEjw8B7wG3fM8F24k8OrSQPaF3KB7w2pMHdCUGhWkT7IQJVm+HihTSjXP7aRa64NGXLg3MAuwBBZoraf95/2HgQ+AdCAVw1xErk0DJBHkLjU9lS8PfcnCowvxqerD510+pxpW8NuTELoDOr0BXd4Ei9ynh9LSNQt3nmP6ulNYW1rwzkONGdTC2dA4RwhR4uQ3EWzO+KsN0BIIxFBiwgvYq7XucI+LWgKngZ5AGLAfeCLrZLNSqjwQr7XWGc8pLNVaN8ztvJIIjLMmZA3v7nwXO2s7vujyBc0qNYF/J8ChxdC4P/SfC2XuPQdwLiKeicsOsy/kBp3qV+XjR5riVLFcEXwCIURhytccgda6q9a6KxAKNNdat8xoYt8MCDLiuq2BIK31Wa11MvAr8PB/rhGnb2ciOzIqnIqC6+3Wm58f+BlbK1ueWfMMvwb9ge47G3p+AMf/hEUPQOyVe57HvYodv47y5b2+jdl/7ga9Zmxjyd7zMncgxH3EmOWjDbXWR2690FofBXxy2f8WJ+BCltdhGdvuoJQaoJQ6CfwLPGPEeYWR6jnW45eHfqGdUzum7p3KO7veJcl3NDy+BMJPw7dd4XLgPc9jYaEY3t6dtS93oqlTBd5acYRh3+3lwo2EIvgUQghTMyYRnFBKLVBKdVFKdVZKzQeyd0bJLqebydm+RmqtV2TcDuqPYb4g+4mUGqWUOqCUOhAeHm7EpcUtDmUcmN1tNmO8x/Bn8J88tfopLjs3g2fXgrKA73vDib+NOlftyrb8PLINH/b3JOB8FL1mbuPH3dL8RoiSzphEMAI4BowHXgaOZ2y7lzDAJctrZ+DS3XbWWm8D6iqlquTw3rcZt6ZaVq1a1YhLi6wslAUv+LzAl12/5HzMeQb/M5h9+qahLEW1xvDbMNgxI8feBtnOZaEY5uvKulc608LVkXf/PMYT8/cQEhFfBJ9ECGEKJls+qpSywjBZ3B24iGGyeIjW+liWfTyA4IzJ4ubA34CzziUomSwumHPR53h588uExoQyocUEnqo3CPXXi3B0GXgPgb4zwaqsUefSWvP7gTA++Pc4KWnpvN6rIcPbuWFpISuLhChuCqvoXJ5orVOBccBaDLeSlmqtjymlxiilxmTs9ihwVCkVAMwBBueWBETBuVdwZ8mDS+jq0pXpB6Yzcfd7JPT7Erq8BYFL4MeHIT7CqHMppXislQvrJ3SmXd0qfPDPcR77ZjfB4XEm/hRCiMIkD5SVUlprvjv6HV8e/BIPRw9mdZmFy/n9sHKsoYrpkKVQLdeVvNnO98fBi0z5+xhJqem80rM+IzvWkdGBEMWEWUYEonhTSjGy6Ujm9ZjH1firDP53MDscq8HwVZCaCN/1hDMb8nS+R1s4s+GVzobnDVaf5JG5uzgjFU2FKPaMaVVZH3gdQxVSq1vbtdbdTBtazmREUPguxF5gwuYJnI48zbhm4xhZuzcWvw6Bq8egfm9oMQI8uoOFpVHn01rz9+HLTP7zKPFJaYzvUY/RnepgZSnfO4Qwl4L2LA4E5mEoNpdZ4F5r7V+YQRpLEoFp3Ey9yXu73mPVuVV0c+nG1NZvUX7PPDj0E8SHQwUXaP4UNHsSHGoadc6IuCQm/3mMf49cxtPJgc8GetOopoOJP4kQIicFTQT+GU8UFwuSCExHa83iE4v5/MDn1HaozcyuM6lj5wynVoH/Iji7GZSlYZTQcgTU7WbUKGH1kcu88+dRom+m8EJXD57v4kEZKxkdCFGUCpoI3sPQpnIFkHRru/QjuH/tv7Kf17a+RlJaEp93/pz2Tu0Nb9w4a2iDGfBzxiihdsYoYdg9Rwk34pOZ8vcx/kFbsWsAACAASURBVAy4RMMa9kwf5I2nU4Ui+DRCCCh4IjiXw2atta5TGMHllSSConEl/grjNo4jOCqYKe2n0K9uv9tvpibDqX/hwEI4t9UwSmjQxzCXULdrrqOEdceu8PbKo9yIT2Zs57q82N2DslbGzT0IIfKvQImguJFEUHTikuN4ecvL7L28l/HNx/Os57PZy1BfD4aDP8ChnyEhwjBKaJExl2BfI8fzRiek8P4/x1l+MIz61cvz2UBvvF0qFsEnEqL0KuiIwBoYC3TK2LQF+EZrnVKYQRpLEkHRSklL4X87/8eqc6sY3GAwb7Z+E8ucvvGnJsPJf8B/IZzbBhZWGaOE4VCnW479DzafvMabfxzhWmwiz3Wqw4Qe9bGxltGBEKZQ0ESwALAGfsjY9CSQprUeWahRGkkSQdFL1+nM9J/JwmML6V67O9M6TsPGKpeOZdeDDZPLAT9DwnVDv+TmT2eMEqrfsWtMYgof/XuCX/dfoE5VOz4b6EUL1xzbZAshCqDAy0e11t732lZUJBGYz+Lji/l0/6f4VPNhdrfZVCh7j8ne1CTDKOHAQgjZnjFKeCBjlND1jlHC9jPhTFp+hEvRN3mmvTuv+TWgXBkZHQhRWAqaCA4Cg7TWwRmv6wDLtNa5trI0FUkE5rUmZA1vbX8LF3sX5vWYR83yxj1TQEQQHFwEAUsMowRHt4xRwjAoXw2AuKRUpq0+weI953GrbMsnj3rRpk5lk30WIUqTgiaC7sBC4CyGHgOuwAit9eZcDzQRSQTmt//KfsZvGk85q3J83eNrGlRqYPzBqUmG/gf+i26PEho+aFhx5N4ZLCzYFRzBxOWHuXDjJk+3deWN3g2xK2t1z1MLIe6uwKuGlFJlgQYYEsFJrXXSPQ4xGUkExcOZyDOM2TCGhJQEZnadSZuabfJ+kogzGXMJS+DmDXB0hxZPg89QEspU4tM1p/hhdwjOjuX45BEv2nlka1UhhDBSYSSCdoAbd9Ya+rGwAswLSQTFx5X4K4zdMJaQmBA+6vARfdz75O9EKYm35xJCd4CFtWGU0HIE+5Unbyw/yrmIeIa0qc2bfRpib2NduB9EiFKgoLeGfgLqAgHcrjWktdYvFWqURpJEULzEJMfw0qaX8L/qz2stX+PpJk8X7IThpw3PJQT8DDcjoVIdUryf4qvIVszeG00NBxumPepFp/rSqU6IvChoIjgBNC4uDWMkERQ/SWlJvLX9LdaFruPJxk/yWsvXsFAFrCWUkpgxl7AQQneChTWRrn58eMWXPyLd6d6oBq/0bEDjWlLETghjFDQR/A68pLW+bIrg8koSQfGUrtP5dP+n/HziZ3q59eKjDh9RxrJM4Zw8/JShxlHgErgZSVS52rxx8ynWJTbmQa+aTOhRH49q5QvnWkLcpwqaCDYDPsA+7iw61++uB5mQJILiS2vNomOL+ML/C1rVaMXMrjNxKFOI39hTEuHEX7BtOjriNHudnmbUhZ7EpSgGNHNmfPd61K5sW3jXE+I+UtBE0Dmn7VrrrYUQW55JIij+/jn7D+/sfAf3Cu7M7T6X6nbV731QXiQnwJpJcPAHUmq1Yl7lN5l9KJn0dM3gVi6M6+ZBzQrlCveaQpRwUnROFLndl3YzYcsE7MvYM7f7XDwcPQr/IkeXw1/jwcKCyJ4z+fxCPX7bfwGlFE/6ujK2S12qlC9b+NcVogQqUM9ipdQjSqkzSqlopVSMUipWKRVT+GGK+0nbWm1Z1HsRqempPLXmKfyvmqChneejMGYbVKqD498j+LDMj2wa78vD3rVYuPMcnT7dzKdrThKdYJb6iEKUGMbcGgoC+mqtTxRNSLmTEUHJcjHuImPWj+FS3CWmdZpGT9eehX+R1GTYOAV2fwU1msLARQTrGszccIa/Ay9hb2PFcx3r8EwHd8rLE8qilCroHMFOrXV7k0SWD5IISp6oxCjGbRrH4fDDTGo9iSGNhpjmQqfXwooxhjIWD30B3o9z4nIMX6w/zfrjV3G0tWZsl7o86esmBe1EqVPQRDALqAGs5M5VQ38UZpDGkkRQMt1Mvckb295gy4UtPOv5LOObj8/e5KYwRF+EP54zPHvgPQQe+AzKlifwQhSfrz/NttPhVLMvy7huHgxu5SLd0USpUdBEsDCHzVpr/UxhBJdXkghKrtT0VD7a+xG/n/6dvnX6MqXdFKwtTVAuIj0Ntn4KWz+Byh4waKHhlhGw9+x1Pl93mn0hN3CqWI7x3evxSHMnrCwL+ACcEMWcrBoSxYbWmvlH5jP70Gza1mzLjK4zsLO2M83Fzm2H5SMNpSp6TYVWI0EptNZsPxPB5+tOERgWjXsVO17uUY++XrWwsDDBKEWIYqCgIwIb4FmgCZDZlkpGBKIgVgat5L1d71HfsT5f9/iaKuVMVFk0PsIwbxC0Hhr1g36zoZyhP7LWmg0nrvH5ulOcvBJLg+r2TOhZn15NqpvmtpUQZlSg5aPATxjmCHoBWwFnILbwwhOlUX+P/szuNpuQmBCGrRrGuehzprmQXRUYshR6fgCnVsG8jnBhPwBKKXo2rs6qlzoy+4lmpKSlM2axP/2+2smWU9coaaNlIfLLmBHBIa11M6XUYa21V0Yz+7Va625FE+KdZERwfzkacZQXNr5Auk7nq+5f4V3VhB1Qw/xh2QiIDoPu70C78Xe0y0xNS2fFoYvM2niGsMibtHJz5FW/BvhKlzRxHyjoiODW0zhRSilPoAKG3gTGXLi3UuqUUipIKTUph/eHKqUOZ/zsUkqZpQ+yMB/PKp4s7rMY+zL2jFw7ks3nTdj4zrkFjNkOjfrChvfg54EQF575tpWlBYNaurDp1S580N+T8zcSePzbPQxbsJeD5yNNF5cQZmbMiGAksBzwwtCysjzwrtZ63j2OswROAz2BMGA/8ITW+niWfdoBJ7TWkUqpPsB7WutcW13JiOD+dP3mdV7Y+AInbpzgf77/Y1D9Qaa7mNaGzmhrJoFNBXjkW6jTJdtuiSlpLN4TytwtwVyPT6Z7w2q84lefJrUqmC42IUzELKuGlFJtMfxi75Xx+k0ArfXHd9nfETiqtXbK7bySCO5fCSkJvLr1VXZc3MEY7zE87/28aSdtrx6D30dAxGno+Cp0eRMssz95HJ+UyqJdIXyzNZiYxFQebFqTCT3r4VHN3nSxCVHICrpqqCzwKNlbVb5/j+MGAr211iMzXj8JtNFaj7vL/q8BDW/t/5/3RgGjAGrXrt0iNDQ015hFyZWSnsIHuz9gRdAKHqn3CO/4voOVhQnLQiTHw+o34NBiqN0WHl0AFZxz3DX6ZgoLtp/l+x3nuJmSRv9mTrzcvb6UvhYlQkETwRogGvDndqtKtNaf3+O4QUCv/ySC1lrrF3PYtyvwNdBBa309t/PKiOD+p7VmTsAcvjn8DR2dOjK983RsrU38y/bw7/DPy2BpDQ9/DQ0fuOuu1+OS+GbbWX7YFUJauuaxVi68KKWvRTFX0ERwVGvtmY+LGnVrSCnlBawA+mitT9/rvJIISo+lp5Yyde9UGldqzJwec6hkU8m0F7webFhVdDkQ2oyBnu+D1d3LWF+NSWTO5iB+2XcepRTD2rgyunMdqjvY3PUYIcyloIngW2C21vpIHi9qhWGyuDtwEcNk8RCt9bEs+9QGNgFPaa13GXNeSQSly6bzm3hj2xtUt63OvB7zcHFwMe0FU5Ng/WTYOxdqesPAhVC5bq6HXLiRwOxNZ1h+8CKWSvFIcydGdapDnarSPlMUH/lKBEqpI4DGMC9QDziLoeicwlBryMuICz8AzAQsge+11lOVUmMwnGCeUmoBhvmHWzf9U+8W6C2SCEqfgGsBjNs0Dktlydfdv6ZJlSamv+jJVfDn85CWAg/NBK97r2IKvR7P/O1nWXogjJS0dHo3qcGYznXxdqlo+niFuIf8JgLX3E6qtTbLjK0kgtLpXPQ5xqwfQ2RSJJ93/pyOzh1Nf9HoMEOtovO7odkw6PMplLl3XaTw2CQW7TrHj7tDiU1Mpb1HZcZ0rksHjypSukKYTX4TgS2QorVOyXjdAHgACDVXCWqQRFCahSeE8/zG5zkdeZr6jvXxqOiBR0UP6jnWw6OiBzXtahb+L9q0VNg6DbZNhyr1DZVMqxs3IolNTOGXfedZsP0c12KT8HRyYEznuvTxrImlFLcTRSy/iWAb8KzW+oxSygPYB/wMNAb2aa3fNFXAuZFEULrFJcfx/dHvOX7jOEGRQVxNuJr5np21HXUr1qVeRUNi8HA0JIrKNpULniDOboE/RkFiNPT+GFqMACPPmZSaxoqDF/l221nORsTjVtmW5zrV4dHmzthYSz8EUTTyPUegtW6a8fcPgEpa6xeUUmUA/1vvFTVJBCKrmOQYgqOCORN5hqCoIIKigjgTeYaopKjMfRzLOmYmhVsjiLoV6+JQxiFvF4u7BitGQ/AmaNwf+s7KrGRqjLR0zbpjV5i3NZjAsGiqlC/Lsx3cGepbGwcbE/RlECKL/CaCw7cmhJVSO4HPtNYrM14Haq3NUhdIEoG4F6011xOvGxJDZEZyiDpDUGQQCakJmftVt62Oh6PHHSOIOhXqUM4ql+cB0tNh1yzY+AFUcIKBiww1jPIY3+7g68zdGsz2MxHYl7ViqK8rz7R3o5osPRUmkt9EsBi4gmHp5yTAXWudoJSqCGyVRCBKGq01l+MvZ44abo0gzkadJTk9GQCFwsXeJTMx3EoSrhVcsbbI8q39wj5Y9izEXoLuk6HtuDsqmRrr6MVo5m0NZtWRy1hZWPBoC2dGd6qDWxUTNesRpVZ+E0E5YDxQE8PSz8CM7e2Aulrrn0wUb64kEYjClpqeyoXYC4ZbTBkjh6CoIEJjQknThofprSyscHNwMySGjNtM9cpVw2njx1ic/Ac8esKAeYb+B/kQEhHPt9vPssw/jNS0dPp41mRM57o0dZYCd6JwFFrROaVUc631wUKLLB8kEYiikpyWzLnoc5kjh6BIwy2mi3EXM/exsbShrrUDHhGh1MOaDm49qeveA1xag32NPF/zWmwiC3eGsHh3KLFJqXTwqMLYLnVpV7cQJrxFqVaYieCg1rp5oUWWD5IIhLklpCQQHBV8x9xD0I0ThGdMUHskJ+MXn0Avi4rUqdXakBScW0GNpoZaRkaISUxhyd7zfLfjHOGxSTR1qsDYLnXp1aSGLD0V+VKYieCQ1rpZoUWWD5IIRHF1LeEaG8+tZW3QSg5GnUYDHqmaXrHR9IpPwF1bQ61m4NIKnDMSRPlquZ4zMSWNFYcu8s3WYEKuJ+BexY5RnerwSHMnylrJ0lNhvMJMBP1vrRwyF0kEoiS4lnCN9aHrWReyjkPXDqHR1Leyp1diGn5XQ3BLTjTs6Oh2Oyk4t4Lqnjn2REhL16w9doW5W4I5cjGaqvYZS0/b1MZelp4KIxQ4ESilnABX7uxHsK3QIswDSQSipLkaf5UN5zewNmQth64dAqBBeWd6lamBX1wcrhcPQ9wVw87WtlCr+Z2jhiwT0FprdgVfZ+6WYHYERWBvY8UwX1dGtHejmr0sPRV3V9Dqo58Ag4Hj3O5HoLXW/Qo1SiNJIhAl2ZX4K5kjhYDwAAAaVWqIX/U2+CkHakcEG5amXjkM6amGgxzdb48YXFpDtSZgacWRsIylp0cvY21pwcAWzozqKEtPRc4KmghOAV5a6yRTBJdXkgjE/eJK/BXWhaxjbehaDocfBqBRpUb0cuuFn3NnXGKuGZJC2H7Dn/HXDAda24FT84zk0JpQ28bM2x/Ncv8wUtPTeaCpYempp5MsPRW3FTQRrAYGaa3jTBFcXkkiEPejy3GXWRe6jnUh6zgcYUgKjSs3NiQFVz+cyztBVChc2A9h+zJGDUcg4zkHKtXlZo0WbE1wZ/65yhxKqkX7etUY27kubWXpqaDgiWA54A1sxNCPAACt9UuFGaSxJBGI+93FuIusD1nPutB1HIkw9IPyrOyJn5sffm5+OJV3MuyYHA+XDt05akiIMLxlaUtAel32ptTlRqVmNOvSn95etSljlfenn8X9oaCJ4OmctmutfyiE2PJMEoEoTcJiw1gfup61IWs5dt3Q3K9plaaZI4Wa5Wve3llriDyXOWpIP78Prh3DQqcRkF6XKdYv07FNG4a0caVGBZlYLm0KY9VQOaC21vpUYQeXV5IIRGl1IfZCZlI4fv04AF5VvfBz9aOXWy9q2OXwJHNSHOkn/iFt1RukpSTxXvKT/K674te4Bk+1dcO3TiW5bVRKFHRE0BeYDpTRWrsrpXyA92XVkBDmcyHmAmtD17IuZB0nbpwAwLuqN73cetHTtWf2pBB9EVaOgXPbOOXYmdFRTxFysxz1q5fnybZuDGjmRPmy2Z9fEPePgiYCf6AbsOXWU8VZexUUNUkEQtzpfMx51oWuY23IWk7eOAmAT1WfzKRQ3a66Ycf0dNgzBza+j7ZxZFuT9/ks2ImjF2MoX9aKgS2cGebrike18mb8NMJUCpoI9mqt22QtL5G1V0FRk0QgxN2FxoQalqSGrOVUpOFObvNqzXmj9Rs0qZzRYvPKEUMv5vCT6NajCWg4gR/3X+Xfw5dJTkunvUdlnmrrRveG1bCylMnl+0VBE8F3GFYMTQIeBV4CrLXWYwo7UGNIIhDCOCHRIawLXcfSU0uJS4ljdrfZtKrRyvBmyk3Y8B7snQdVG8Gj84koX5/f9l/g5z2hXIpOpFYFG4b6uvJ4Kxcqly9r1s8iCq6gicAWeBvwAxSwFvhAa51Y2IEaQxKBEHlzNf4qo9aP4mLcRb7o8gWdnDvdfjNoA6x8Hm5GQvd3wfcFUjVsOHGNn/aEsDPoOmUsLXjIqyZPtnXFx6WiTC6XUIVZdM4SsNNaxxRWcHkliUCIvItMjGT0+tGciTzDxx0/prd779tvxkfAXy/BqX/BvTP0n2towwkEXYvlp92hLD94kbikVLycK/Ckryt9vWthYy3VT0uSgo4IlgBjMNQZ8gcqAF9orT8r7ECNIYlAiPyJTY5l3MZxHLp2iPfavccj9R65/abWcPBHWDMJLMtA35nQZEDm23FJqaw4GMaPu0M5cy2OirbWDG7lwrA2rrhUsjXDpxF5VdBEEKC19lFKDQVaABMBf5ksFqLkuZl6kwlbJrDz4k5eb/k6TzV56s4drgcbJpIvHQTvIdDnE7BxyHxba83us9f5aXco645fJV1rujesxpNt3ejoUQULaZpTbBU0ERwDfIAlwFda661KqUBpXi9EyZSSlsLE7RNZH7qesd5jGes99s77/mkpsPVT2D4dKrjAI/Ohdpts57kcfZMle8/zy77zRMQl417FjmG+rgxs4UyFctIjobgpaCJ4CcMoIBB4EKgNLNZadyzsQI0hiUCIgktNT2XK7imsDFrJsEbDeKPVG9kngc/vgT9GQfQF6PgadH4jx1abSalprDl6hR93h+IfGkk5a0v6N3PiqbauNKrpkG1/YR6FNlmc5YRWWuvUAkeWD5IIhCgc6Tqdz/Z/xuITixngMYDJbSdjafGfCeDEGFj9BgT+Ak4t4ZFvoXLdu57z6MVoftodysqAiySlptParRJPtnWlt2cNrOWZBLMq6IigLIbnB9y4s0PZ+4UYo9EkEQhReLTWfB34NfMC5+Hn6se0jtOwzuFbP0f/gH8mGG4b9f4Ymj8FuSwjjUpI5vcDYfy0J5TzNxKoZl+WJ1rXZkib2lR3kIJ35lDQRLAGiMawYuhWhzK01p8bceHewCzAEligtZ72n/cbAguB5sDbWuvp9zqnJAIhCt8Px35g+oHpdHDqwBddvqCcVbnsO2WpV0TDh6Dvl2BXOdfzpqdrtp4O54fdIWw9HY6lUvTyrMFTvq60dpeCd0WpoIngqNbaMx8XtQROAz2BMGA/8ITW+niWfaph6IXcH4iURCCE+Sw7vYz3d79P8+rN+arbV5Qvk0PNoSz1iihXCfp/DR7djTp/6PV4Fu8J5bf9F4hJTKVhDXtGtHfjYR8neSahCOSWCIy5abdLKZWfAnOtgSCt9VmtdTLwK/Bw1h201te01vuBlHycXwhRiAbWH8gnnT4h8FogI9eNJDIxMvtOFhbQ7kV4bhOUqwiLH4HVEw0lK+7BtbIdbz/YmL1v9eCTR5uilGLi8iN0+GQzszeeITI+2QSfShjDmETQAfBXSp1SSh1WSh1RSh024jgn4EKW12EZ2/JMKTVKKXVAKXUgPDw8P6cQQhihj3sfZnWbRVBUECPWjOBawrWcd6zRFEZtgTZjDPWKvu0KV44adY1yZSwZ3Ko2q17qwM8j2+Dp5MDn60/TdtpG3ll5lJCI+EL7PMI4xtwacs1pu9Y69B7HDQJ6aa1HZrx+EmittX4xh33fA+Lk1pAQxcP+K/sZt3EcjjaOzPebj4u9y913vqNe0WTwfd4wcsiD01djWbD9LCsPXSIlPR2/xtV5rmMdWrg6yjxCISnQraGMX/guQLeMvycYcxyGEUDW/3qcgUtGHCeEMLNWNVqxwG8BcSlxDF89nOCo4Lvv7NEDxu4Cj56w7m34qT/E5O3/6vWr2/PpQG92TOrKC1082HP2BgPn7WbA17tYdeQyael5X+YujGfMiGAy0BJooLWur5SqBfyutW5/j+OsMEwWdwcuYpgsHqK1PpbDvu8hIwIhip0zkWcYtX4UqempzOs573ZPg5xkq1c0C5r0z9d1E5JTWeYfxnc7zhF6PQGXSuV4tr07g1q6YCed1PKlwLWGgGbAwbw2plFKPQDMxLB89Hut9VSl1BgArfU8pVQN4ADgAKQDcUDj3KqbSiIQomidjznPc+ueIzo5mjnd59CieovcD8har8hnKPSedke9orxIS9esP36F+dvP4R8aiYONFUN9XRnezk2eR8ijgiaCfVrr1kqpg1rr5kopO2C3FJ0TovS4En+FUetHcTnuMjO6zqCDU4fcDzCyXlFe+IdGsmD7WdYeu4KlhaKftxPPdXKnYQ0pY2GMgiaC14B6GJ4H+Bh4Fliitf6ysAM1hiQCIczjRuINxqwfw5moM3zS8RP83PzufVDWekWdXjf85PTkch6EXo/n+x3nWHogjJspaXSsV4XnOtahY70qMrGciwLXGlJK9cTQoQxgrdZ6QyHGlyeSCIQwn5jkGMZtHEdgeCDvtX2PAfUG3PugPNYrMlZUQjI/7z3Pol0hhMcm0bCGPSM71qGfdy3KWEldo//KVyJQSsUCt978b5pNBIIxlIXYWFiBGkMSgRDmlZCSwIQtE9h1aRcTW01kWONhxh2YtV5Rn2nQ7Mlc6xUZKyk1jb8CLrFg+zlOXY2lmn1Zhrd3Y2hrVyrYSjnsW0xRfdQS8AR+zk/5iYKQRCCE+SWnJTNx20Q2nN/A8z7PM8ZrjHG3Zf5br6j1c1ClAdjXKHBS0Fqz7UwE87edZUdQBLZlLHmspQvPdnCXLmqYIBFkOfForfU3+T5BPkgiEKJ4SE1PZfKuyfwV/BdPNX6K11q+ZlwyyFqvKC2jrETZClC1viEpZP2zoiv8tzS2EY5fimHB9rP8FXiJdK3p41mTkR3daVbbMc/nul+YLBGYgyQCIYqPdJ3OtH3T+OXkLzxa71He8X0ne0+Du4mPgKtHIfw0RJyC8FMQcRrirt7ex7IsVKkHVepD1Qa3/6zsAVZl73mJK9GJLNoVws97Q4lNTKWlqyPPdapDj0bVsSxlbTUlEQghTEZrzexDs5l/ZD693XrzUcePsLYowL35m5HZk0P4KYg6T+a0pbIAR7f/jCAyEkUOzyzEJaWydP8FvttxjotRN3GrbMuzHeswsLkz5cqUjsqnkgiEECa38OhCvvD/gk7Onfi88+fYWBXyA18pNyHizO3EEHHKkDCuB0F6lgLG9jWzjyCqNIDy1UhN16w5doX5284SGBaNo601w3xdeaqtG1Xt7z3CKMkkEQghisTSU0v5cM+HtKjegtndZufc06CwpaVCZEj2EUTEaUiOu72fTYXMEYSu0oDT6bVYdLoMvwUprCytGODjxMiO7tSrbm/6mM1AEoEQosisOruKt3a8RaNKjZjbYy4VbSqaJxCtDcXvbo0csv4Zf7ucfbqlDVetnfFPqMbptFqUrdmIjn0ex6uus3niNhFJBEKIIrXlwhZe3fIqtR1q823Pb6lqW9XcId0p4cbtkUPGbaa0a6ewjDG0UAnXDmyuOYqOgydQ07EIRjVFQBKBEKLI7b28lxc3vUhlm8rM95uPs30J+IadHM/NkP1E/P0uLrGBnNK1Odr0DR7oN6TETyoXtFWlEELkWZuabVjgt4CY5BieXvM0Z6PPmjukeytjR7n6XXB5ZSsRfeZT2TqFR4+O4+C0nmzYto30+7QvgiQCIYTJeFX14vte35OWnsbw1cM5fv24uUMyjlJUafMYVSYFcr7Fm/joE3TZ+DCrPx1G4OlcmvSUUJIIhBAm1aBSA37o8wM2VjY8u/ZZDl49aO6QjGdVltp9J1HulcOEuD1Gr8RVuP/cgZVzJnEpIsrc0RUaSQRCCJNzdXDlxz4/UqVcFUavH83OizvNHVKeWNhXxWPENyQ/t50Ix2b0D59LyuzW/PXrPBKSUu59gmJOEoEQokjUsKvBot6LcKvgxrhN41gfut7cIeWZrZMndV5eRXj/X7AuW45+JydyalpHNm1aU6LnD2TVkBCiSMUkx/D8huc5EnGEHrV70LZWW3xr+paMVUVZpaUSsmEeFfd8RkUdxaay3any8Ad4Nc6lr7MZyfJRIUSxkpCSwBf+X7D5/Gau3bwGgFN5J3xr+uJby5c2NdrgaFMyKoWm34zm9PL3cQ/6gXSt2FLlCbweewen6sXr2QlJBEKIYklrzbmYc+y5tIc9l/ew/8p+4lIMZSEaVWpEm5pt8K3pS/PqzSn3//buPbiK8ozj+PeXC5AQuRgSTAjIRQYRFKEIAcW2oh20VOpYbzPa0T9K7ajVtqOtdUrvre20tnScKeOgi2GfKwAACwtJREFUVaeKFxRLW6eoQ1WqECFc5BKkXEQCCTk0KpeAScjTP3YTDiGA4DnZA/t8Zs7k7J7ds79kTs6z++7u++bkRZz22PbXbWbrc/dx7q5XqbPeLD/nTiZd9226d+sSdTTAC4Fz7hTR3NLMuv+tY0lNUBhW1q2kqaWJ3KxcLiy+kPKScsaXjGdE4QhysnKijtuhxLpF7Jt/LwMPVLGeQeycMINJV1xDVsTdXnshcM6dkhqaGlhRt4KKmgqW1Cyhqr4KgILcAi4666KgKamknEE9B2XWwPVmbH79SQoW/ZzilgSLcyfQ/Su/5IILPhdZJC8EzrnTQv2Bet6pfYclO5ZQUVNB9d5qAIrziikvLW87YijOL444aaDlkwaqXvoNA6tm0cWaWNT7GoZd/zP6lfbr9CxeCJxzp6Vte7a1HS1U1FTw0SfBTV5Deg5pO78w9qyxnNEl2q6lG+p3sOnZ+zmv9m/sIZ/KQdMZf/33KcjvvPMeXgicc6e9Fmthw4cb2k48V+6s5MDBA2Qrm5F9RrYVhlFFo+iSHc0J3MTGSurn3cuwfZVspZRtY+9n4pU3k5Wd/lu6vBA452Kn8WAjqxKr2k48r9m1hhZrIS8njzF9xzChJLh/YWjvoWSpE++tNWPT2/PounAGZQe3sTLnArKn/Irzx05K62a9EDjnYm93426W1S5ra0Zq7Q21d9febUcL5aXl9CvonPZ7a25k9d9nMmDVTHrYXhb3mMLZ1/2asgGD0rI9LwTOOdfOzn07qaitaDvx3HpjW1lBGUN6DaEov4ji/GKK84qDn+GjV9deKb1Caf/ueqqe/REjq+fQRA6V/W9l9A0PcMYZPVO2DfBC4Jxzx2RmbPl4C4trFrO0dinb926nrqGO+gP1Ryybm5VLcX4xRXlFhxWIovwi+ub3bZufn5t/QhkSW9ezY+59jNrzBrUUsmXU9xh39e1kZ6dmQJzICoGkKcBMIBuYbWYPtntd4etXAQ3ArWZ2zD5qvRA45zpL08EmEvsT1DXUUddQR2J/gp0NO0k0HJpX11BHQ3PDEesW5BZ0WCCSH4V5heRm5R623sZ3FqBXHmBI8395L3soTZf/gpETpnzm3yWSQiApG9gAXAFUA0uBm8xsXdIyVwF3ERSC8cBMMxt/rPf1QuCcyzT7mvZ1WCCSC0eiIUGzNR+2nhBndjvziCOL4m5F7Fm7jGFVz3PuwXo25V1MybW/pd/g804647EKQTrv0R4HbDSzzWGIZ4BpQPIQRdOAJy2oRksk9ZJUYmY1aczlnHMp1T23O4N7DmZwz8FHXabFWvjwwIdHPbKo3VfL6l2rD2+OKssH8sm19ylaeC0XvzWGGbc8lfL86SwE/YBtSdPVBHv9x1umH3BYIZA0HZgOMGDAgJQHdc65dMtSFoV5hRTmFTKc4UddrvFgI7v27zrsyGJrYiPvv/cahT3OTku2dBaCjk6rt2+H+jTLYGaPAI9A0DT02aM551xm6pLdhdKCUkoLSg9/4fM/Sds203kXRTXQP2m6DNhxEss455xLo3QWgqXAUEmDJHUBbgTmt1tmPvB1BcqBj/38gHPOda60NQ2ZWbOkO4EFBJePPmZmayXdHr4+C3iZ4IqhjQSXj96WrjzOOec6ltaRHczsZYIv++R5s5KeG3BHOjM455w7tk7sack551wm8kLgnHMx54XAOedizguBc87F3CnX+6ikBLD1JFfvA+xKYZxUydRckLnZPNeJ8Vwn5nTMdbaZFXX0wilXCD4LScuO1ulSlDI1F2RuNs91YjzXiYlbLm8acs65mPNC4JxzMRe3QvBI1AGOIlNzQeZm81wnxnOdmFjlitU5Auecc0eK2xGBc865drwQOOdczMWmEEiaIuk9SRsl/SDqPACSHpNUJ2lN1FmSSeov6d+SqiStlXR31JkAJHWT9I6kVWGun0adKZmkbEkrJP0j6iytJL0vabWklZIyZrDvcFjauZLWh5+zCRmQaVj4d2p97JZ0T9S5ACR9J/zMr5E0R1K3lL5/HM4RSMoGNgBXEAyGsxS4yczWHXPF9Oe6FNhLMG7zyCizJJNUApSY2XJJZwCVwFcz4O8loLuZ7ZWUC/wHuNvMlkSZq5Wk7wJjgR5mNjXqPBAUAmCsmWXUzVGSngAWmdnscLySfDP7KOpcrcLvjO3AeDM72RtYU5WlH8Fn/Twz2y/pOeBlM3s8VduIyxHBOGCjmW02s0bgGWBaxJkwszeB+uMu2MnMrMbMlofP9wBVBGNJR8oCe8PJ3PCREXsyksqALwOzo86S6ST1AC4FHgUws8ZMKgKhycCmqItAkhwgT1IOkE+KR3KMSyHoB2xLmq4mA77YTgWSBgKjgYpokwTC5peVQB3wqpllRC7gj8B9QEvUQdox4BVJlZKmRx0mNBhIAH8Jm9JmS+oedah2bgTmRB0CwMy2A78DPgBqCEZyfCWV24hLIVAH8zJiTzKTSSoAXgDuMbPdUecBMLODZnYhwfjW4yRF3qQmaSpQZ2aVUWfpwMVmNga4ErgjbI6MWg4wBvizmY0G9gEZcd4OIGyquhp4PuosAJJ6E7RgDAJKge6Sbk7lNuJSCKqB/knTZaT40Op0E7bBvwA8ZWYvRp2nvbAp4XVgSsRRAC4Grg7b458BLpP012gjBcxsR/izDphH0EwatWqgOulobi5BYcgUVwLLzWxn1EFClwNbzCxhZk3Ai8DEVG4gLoVgKTBU0qCw2t8IzI84U8YKT8o+ClSZ2UNR52klqUhSr/B5HsE/yPpoU4GZ3W9mZWY2kOCztdDMUrrHdjIkdQ9P9hM2vXwJiPwKNTOrBbZJGhbOmgxEeiFCOzeRIc1CoQ+Ackn54f/mZILzdimT1jGLM4WZNUu6E1gAZAOPmdnaiGMhaQ7wBaCPpGrgx2b2aLSpgGAP9xZgddgeD/DDcAzqKJUAT4RXdGQBz5lZxlyqmYH6AvOC7w5ygKfN7F/RRmpzF/BUuGO2Gbgt4jwASMonuLrwm1FnaWVmFZLmAsuBZmAFKe5qIhaXjzrnnDu6uDQNOeecOwovBM45F3NeCJxzLua8EDjnXMx5IXDOuZjzQuBiTdIfknuYlLRA0uyk6d9LmnGiPdZKelzS11KZ1bl08ULg4u5twrs0JWUBfYARSa9PBBaY2YMRZHOuU3ghcHH3Fodu1x9BcOftHkm9JXUFhgOjJD0MbXv6f5L0tqTNrXv9CjwsaZ2kfwLFrRuQNDnsXG11OAZFV0njJL0Yvj5N0n5JXcIxFzZ34u/vnBcCF29hXzzNkgYQFITFBD2tTiAYW+BdoLHdaiXAJcBUoPVI4RpgGHA+8A0OHWV0Ax4HbjCz8wnu8P0WwV2io8N1JxEUoIuA8WRIT68uPrwQOHfoqKC1ECxOmn67g+VfMrOWcKCevuG8S4E5Ye+oO4CF4fxhBB2GbQinnwAuNbNmYKOk4QQdwT0UvsckYFGqf0HnjsULgXOHzhOcT7BnvoTgiGAiQZFo75Ok58ldnHfUX0tHXaC3WkTQ02UT8BrBUcYlwJufNrhzqeCFwLngy34qUB/u0dcDvQiKweJP+R5vAjeGA+eUAF8M568HBko6J5y+BXgjaZ17gMVmlgAKgXOByDtEdPESi95HnTuO1QRXCz3dbl6Bme0Ke+88nnnAZeF6Gwi/7M3sgKTbgOfDYQaXArPCdSoImpZajwDeJRjgxnuCdJ3Kex91zrmY86Yh55yLOS8EzjkXc14InHMu5rwQOOdczHkhcM65mPNC4JxzMeeFwDnnYu7/pUcHt4WACq8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "labels = ['12 clusters', '6 clusters', '3 clusters']\n", "\n", "ces_fig2, ces_ax2 = plt.subplots()\n", "for data, label in zip(ces_conv2.T, labels):\n", " plt.plot(data, label=label)\n", "ces_ax2.set_xlabel('Window')\n", "ces_ax2.set_ylabel('Jensen-Shannon divergence')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using default arguments with dimension reduction ensemble similarity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See [dimension_reduction_ensemble_similarity.ipynb](dimension_reduction_ensemble_similarity.ipynb#Calculating-dimension-reduction-similarity-with-default-settings) for an introduction on comparing trajectories via dimensionality reduction. We now use the `dres_convergence` function ([API docs](https://docs.mdanalysis.org/stable/documentation_pages/analysis/encore/similarity.html#MDAnalysis.analysis.encore.similarity.dres_convergence))." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "dres_conv = encore.dres_convergence(u, # universe\n", " 10, # window size\n", " select='name CA') # default" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Much like `ces_convergence`, the output is an array of similarity values." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.53259849],\n", " [0.40801451],\n", " [0.3137567 ],\n", " [0.25160251],\n", " [0.19061233],\n", " [0.12228127],\n", " [0.05960165],\n", " [0.03835858],\n", " [0. ]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dres_conv" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Jensen-Shannon divergence')" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUZdr/8c+VkBAIEFqkhSYgSC8RxIKI4IIFLOjqou7aENuK6+6qv2ftu4+7PpZV14JiX+uChVUROxZQSZDepUik1yBICbl+f8yAWTaESTKTM8l8369XXsk5M3PmC2Kuucu5b3N3REQkcSUFHUBERIKlQiAikuBUCEREEpwKgYhIglMhEBFJcNWCDlBaDRs29FatWgUdQ0SkUsnNzd3g7pnFPVbpCkGrVq3IyckJOoaISKViZisO9pi6hkREEpwKgYhIglMhEBFJcCoEIiIJToVARCTBqRCIiCQ4FQIRkQSXMIVg3bad3PHvuewuKAw6iohIXEmYQpCzfDPPfLmc2ybMQXswiIj8LGEKwSldmnBV/za8/M1Knp960BvsREQSTsIUAoDfn9yegUc24s635/Hlkg1BxxERiQsJVQiSkoy/n9edNpnpXPXidJZv2B50JBGRwCVUIQCoVb0aYy86iiSDy57PIX/nnqAjiYgEKuEKAUCLBjV5dEQvlm/YzuhXZrC3UIPHIpK4ErIQAPRt04Dbh3bi4wXruGfSgqDjiIgEptLtRxBNFxzdkgVr8hkzeSntG9XmrJ5ZQUcSEalwCdsi2Oe20ztx9OH1uen12Xz7/eag44iIVLiYFgIzG2xmC81siZndVMzj/c1sq5nNCH/dGss8xUlJTuLREb1oVKc6V7yQy5qtOys6gohIoGJWCMwsGXgEGAJ0BM43s47FPPVzd+8e/rozVnlKUj89lad+fRTbdxUw8oUcdu7ZG0QMEZFAxLJF0BtY4u5L3X038AowLIbvVy5HNKrNg+f1YPYPW/njuFlahkJEEkYsC0EzYGWR47zwuQP1NbOZZjbRzDoVdyEzG2lmOWaWs379+lhkBWBgx0b8/uT2TJi5ikc//S5m7yMiEk9iWQismHMHfsyeDrR0927Aw8CbxV3I3Z9w92x3z87MzIxyzP90Vf82DO3WlHvfX8gH89bG9L1EROJBLAtBHtC8yHEWsKroE9w9391/DP/8LpBiZg1jmOmQzIx7hnelS7MMRr/yLQvXbAsyjohIzMWyEEwD2plZazNLBc4DJhR9gpk1NjML/9w7nGdjDDNFJC0lmScuzCa9ejUue34am7bvDjqSiEjMxKwQuHsBcA0wCZgPvObuc81slJmNCj9tODDHzGYCDwHneZyM0jbOSGPMhb1Ym7+Lq17MZc9ebWgjIlWTxcnv3YhlZ2d7Tk5Ohb3fG9/mcf2rM7ng6Bb8+YwuFfa+IiLRZGa57p5d3GMJvcREJM7skcWCNdtCy1A0rsOFR7cMOpKISFQl/BITkfjjLzowoMNh3DFhLlO+04Y2IlK1qBBEIDnJePC87rRqGNrQ5vuNO4KOJCISNSoEEaqdlsLYi7Jxh8uen8Y2bWgjIlWECkEptGqYzqMjevLd+u1c/+oMCrWhjYhUASoEpXRs24bcelpHPpy/jnvfXxh0HBGRctOsoTK4qG9LFqzZxqOffkf7xrUZ1r24JZRERCoHtQjKwMy4Y2gnereuzx/HzWLmyi1BRxIRKTMVgjJKrZbEYyN60rBWdUa+kMPafG1oIyKVkwpBOTSoVZ2xv85m284CRr6Qqw1tRKRSUiEopyOb1OH+c7szc+UWbn59tja0EZFKR4UgCgZ3bswNg47gjW9/YMxnS4OOIyJSKioEUXLNgLac2rUJf3tvAR8v0IY2IlJ5qBBEiZlx7/BudGpah9++PIPFa7WhjYhUDocsBGZ2hJl9ZGZzwsddzexPsY9W+dRIDW1ok5aSzGXP57BZG9qISCUQSYvgSeBmYA+Au88itNuYFKNp3RqMubAXq7fs5OqXpmtDGxGJe5EUgpru/s0B5wpiEaaq6NWyHv97VhemfLeRP789L+g4IiIlimSJiQ1m1gZwADMbDqyOaaoqYHivLBaszmfsF8to37gOv+rTIuhIIiLFiqQQXA08AXQwsx+AZcAFMU1VRdx8ypEsXvcjt741hzaZ6fQ5vEHQkURE/sshu4bcfam7DwQygQ7ufpy7L495siogOcl46PwetGhQkytfnM7KTdrQRkTiTySzhv7XzOq6+3Z332Zm9czszxURrirIqBHa0KZgbyGXP5/D9l0aXhGR+BLJYPEQd9+/vKa7bwZOiV2kqufwzFr841c9WbR2mza0EZG4E0khSDaz6vsOzKwGUL2E50sx+h2Ryf+c2pH3563lgQ8XBR1HRGS/SAaL/wl8ZGbPEJo5dAnwXExTVVGXHNuKhWvyefjjJRzRqDand2sadCQRkUMXAne/x8xmAycBBtzl7pNinqwKMjPuOqMzS9dv5w/jZtK6YTqdm2UEHUtEElxEaw25+0R3/72736AiUD7VqyXz2AW9qF8zlcufz2HdNm1oIyLBimTW0FlmttjMtppZvpltM7P8ighXVWXWrs6Tv85my449XPFCLrsKtKGNiAQnkhbBPcBQd89w9zruXtvd68Q6WFXXqWkG953bjW+/38L/e32ONrQRkcBEUgjWuvv8slzczAab2UIzW2JmN5XwvKPMbG94+YqEcUqXJlx3UjvGT8/j1rfmslfTSkUkAJHMGsoxs1eBN4Fd+066++slvcjMkoFHgEFAHjDNzCa4+7xinvc3ICHHHkYPbMfOgr2MmbyUzTt2c/+53Umtpm0iRKTiRFII6gA7gJOLnHOgxEIA9AaWuPtSADN7BRgGHLgc57XAeOCoSAJXNWbGzUOOpH7NVO6euICtP+3h8Qt6kV49kv80IiLlF8n00YvLeO1mwMoix3lAn6JPMLNmwJnAAEooBGY2EhgJ0KJF1VzF84oT2lAvPZWbxs9ixNiveeY3R1EvPTXoWCKSAGK5Q5kVc+7ATvC/Aze6e4nTZtz9CXfPdvfszMzMCN66cjo3uzmPXdCLeavzOWfMVFZv/SnoSCKSAGK5Q1ke0LzIcRaw6oDnZAOvmNlyYDjwqJmdEcG1q6xfdGrMcxf3Zs3WnQx/bCrfrf8x6EgiUsXFcoeyaUA7M2ttZqmEiseEok9w99bu3srdWwHjgKvc/c0Irl2l9W3TgFdGHs2ugr2c8/hUZudtDTqSiFRhkRSCMu1Q5u4FwDWEZgPNB15z97lmNsrMRpUjc0Lo3CyDf406hhopyZz3xFSmLNkQdCQRqaLsUDcymdnhhHYoOwbYTHiHsqA2p8nOzvacnJwg3joQa7bu5KKnv2b5hh08dH53BnduEnQkEamEzCzX3bOLe0w7lMW5xhlpvHZFXzo3q8NVL07nlW++DzqSiFQxh5w+ama/O+AYYCuQ6+4zYpRLiqhbM5V/XtaHq16czk2vz2bTjt1ceUKbff8tRETKJZIxgmxgFKH7ApoRms/fH3jSzP4Yu2hSVM3Uajx5UTbDujflnvcW8pd35munMxGJikhuX20A9HT3HwHM7DZCM3z6AbmEFqWTCpCSnMQD53anbo0Uxn6xjM079vDXs7uQkqwlKUSk7CIpBC2A3UWO9wAt3f0nM9t1kNdIjCQlGbcP7UT99Oo88OEitv60m3/8qidpKclBRxORSiqSQvAS8JWZvRU+Ph142czS+e91g6QCmBnXDWxH/Vqp3PrWHC566hvG/iabOmkpQUcTkUqoxD4FC41GPgtcDmwhNEg8yt3vdPft7j4i9hHlYC48uiUPndeDb1du5pdjvtJuZyJSJiW2CNzdzexNd+9FaDxA4szp3ZqSUSOFK17I5ZzHp/LCJX1o0aBm0LFEpBKJZJTxKzNLyCWiK4t+R2Ty4uV92PrTHs5+fArzV2snURGJXCSF4ERCxeA7M5tlZrPNbFasg0np9GxRj39d0ZdkM84dM5VpyzcFHUlEKolICsEQ4HBCewacDpwW/i5xpl2j2oy7si+Ztapz4VNf8/GCtUFHEpFKIJIlJlYQWk56QPjnHZG8ToKRVa8m/xrVl3aH1eby53N549u8oCOJSJyLZGOa24AbCe1JAJAC/DOWoaR8GtSqzkuX96FP6/pc/+pMnv5iWdCRRCSORfLJ/kxgKLAdwN1XAbVjGUrKr3ZaCk//5igGd2rMnW/P495JCznUSrMikpgiKQS7PfQbZN9+BOmxjSTRkpaSzCMjenLeUc35xydL+J8357BX6xOJyAEiubP4NTMbA9Q1s8uBSwhtXymVQHKScfdZXaiXnspjn37Hlh27eeCX3aleTUtSiEjIIQuBu99rZoOAfKA9cKu7fxDzZBI1ZsaNgztQv2Yqf3l3Pvk/5fD4hb2oVT2SzwEiUtVFsh/B9cC/9Mu/8ru83+HUS0/lxvGzGPHkVzxzcW/qp6cGHUtEAhbJGEEdYJKZfW5mV5tZo1iHktgZ3iuLMRf0YsGabZzz+BR+2PJT0JFEJGCR3Edwh7t3Aq4GmgKTzezDmCeTmBnYsRHPX9Kbdfm7GP7YFJas2xZ0JBEJUGluDFsHrAE2AofFJo5UlD6HN+CVK45mz17nnMenMmPllqAjiUhAIrmh7Eoz+xT4CGgIXO7uXWMdTGKvU9MMxl/Zl1pp1fjVk1/xxeINQUcSkQBE0iJoCYx2907ufpu7azOaKqRlg3TGjzqGFvVrcvGz3/Du7NVBRxKRCnbQQmBmdcI/3gN8b2b1i35VTDypCIfVSePVkX3pllWXq1+azotfrwg6kohUoJJaBC+Fv+cCOeHvuUWOpQrJqJnCC5f2of8RmfzPG3P4x8eLtSSFSII46H0E7n5a+HvriosjQaqRmswTF2Xzx3GzuPf9RWzcvps/ndqR5CQLOpqIxNBBC4GZ9Szphe4+PfpxJGgpyUncd0436tZM4ZkvlzNt+SbuGtaZHi3qBR1NRGLEDtb8N7NPwj+mAdnATMCArsDX7n5chSQ8QHZ2tufkqGcq1tydt2et5s/vzGPdtl2cd1QL/viL9tTTncgilZKZ5bp7dnGPHXSMwN1PdPcTgRVAT3fPDm9i3wNYEuEbDzazhWa2xMxuKubxYeHtL2eYWY6ZBVJc5L+ZGad3a8pHN/Tn0mNb81rOSgbc9ymvTvueQq1gKlKlHLRFsP8JZjPcvfuhzhXzumRgETAIyAOmAecXnX5qZrWA7e7uZtYVeM3dO5R0XbUIgrFgTT63vDmHacs307NFXe46ozOdmmYEHUtEIlSmFkER881srJn1N7MTzOxJYH4Er+sNLHH3pe6+G3gFGFb0Ce7+o/9cidIJ73kg8adD4zq8dkVf7j2nGys27uD0h7/g9glzyd+5J+hoIlJOkRSCi4G5wHXAaGBe+NyhNANWFjnOC5/7D2Z2ppktAN4htNeBxCkzY3ivLD6+oT8j+rTkuanLOem+ybz57Q+aaipSiUWy6NxOd3/A3c8Mfz3g7jsjuHZxcw7/67eFu78R7g46A7ir2AuZjQyPIeSsX78+greWWMqomcJdZ3TmrauPpWlGGqNfncH5T37F4rVavE6kMirNonOllQc0L3KcBaw62JPd/TOgjZk1LOaxJ8KD1dmZmZnRTypl0jWrLq9fdSx/ObMz81dvY8iDn3P3xPls31UQdDQRKYVYFoJpQDsza21mqcB5wISiTzCztmZm4Z97AqmEVjeVSiI5yRjRpyUf33ACZ/ZoxpjJSxl4/2Qmzl6t7iKRSiJmhcDdC4BrgEmEBpdfc/e5ZjbKzEaFn3Y2MMfMZgCPAL90/faolBrUqs7/ndONcaP6klEjhStfnM6vn5nGsg3bg44mIocQyfTRI4A/EFqFdP+dyO4+ILbRiqfpo/GvYG8hz09dwf0fLGJ3QSGj+rfhqv5tSEtJDjqaSMIqafpoJIVgJvA4ocXm9u477+650QwZKRWCymNd/k7+8u583pqxiub1a3DH0E4M6KCdTkWCUN5CkBu+ozguqBBUPlOWbOCWt+bw3frtDOrYiNtO70hWvZpBxxJJKOW9oezfZnaVmTXRfgRSFse0bcjE6/px4+AOfLF4AwPvn8wjnyxhd0Fh0NFEhMhaBMuKOe3ufnhsIpVMLYLK7YctP3HXv+fx3tw1HJ6Zzl3DOnNs2/+aMSwiUVaurqF4o0JQNXyycB23T5jLio07OK1rE245rSON6qQFHUukyipX15CZpZjZb81sXPjrGjNLiX5MSSQntj+MSaP7MXpgO96ft5aT7pvM2M+XUrBX3UUiFS2SrqGxQArwXPjUhcBed78sxtmKpRZB1bNi43ZumzCXTxeup0Pj2tx1RmeOaqVhKJFoKvf0UXfvdqhzFUWFoGpydybNXcud/57Lqq07ObtnFjef0oGGtaoHHU2kSijvrKG9ZtamyMUOp8j9BCLRYGYM7tyYD284gSv7t2HCzB8YcO+nvPDVCvZqIxyRmIqkEPwB+MTMPjWzycDHwA2xjSWJqmZqNW4c3IGJ1x1Pp6YZ3PLmHM589EtmrtwSdDSRKiuiWUNmVh1oT2hp6QXuvivWwQ5GXUOJw92ZMHMVf3lnPut/3MWverfgD79oT92a2jdZpLRK6hqqVtzJYvQCWoWf383McPfno5RPpFhmxrDuzRjQ4TAe+GAxz01dzsQ5a7hpSAeG98wiKam4LS9EpLQiGSx+AWgDzODnsQF399/GOFux1CJIXPNW5XPLW3PIXRHaN/nOYZ3p3Ez7JotEoryzhuYDHeNleWgVgsRWWOiMn57H395bwMbtu9VdJBKh8s4amgM0jm4kkbJJSjLOyW7ORzf05zfHtOKVaSs58d5Peenr7zW7SKSMImkRfAJ0B74B9g8Su/vQ2EYrnloEUtT81fncNmEu3yzbRNesDO4Y2okeLeoFHUsk7pS3a+iE4s67++QoZCs1FQI50L7ZRf/77nzW5u/i3Ows/jhYN6OJFKVF5yQh/LirgIc/WsxTXyyjZmoyN5zcnhF9WlAtOZZbc4tUDuVddO4sM1tsZlvNLN/MtplZfvRjipRPrerVuPmUI3lvdD+6ZtXltglzOe3hL/hm2aago4nEtUg+Kt0DDHX3DHev4+613b1OrIOJlFXbw2rxwqW9eWxET/J/2sO5Y6Zy/aszWJe/M+hoInEpkkKw1t3nxzyJSBSZGUO6NOGjG/pz7YC2vDNrNSfe+ylPfraUPVrqWuQ/RDJY/CCh6aNv8p+zhl6PbbTiaYxAymL5hu3c8e+5fLJwPW0Pq8UdQztpZzRJKOW9j6AOsAM4GTg9/HVa9OKJxF6rhuk8c3Fvxl6Uza6CvYwY+zVXvzidVVt+CjqaSOA0a0gSzs49e3nis6U88skSksy4ZkBbLju+NdWrJQcdTSRmynsfQRpwKdAJ2L+prLtfEs2QkVIhkGhZuWkHf35nHpPmrqVVg5rcNrQTJ7Y/LOhYIjFR3q6hFwiNEfwCmAxkAduiF08kGM3r12TMhdk8d0lvksy4+JlpXP58Dis37Qg6mkiFiqQQtHX3W4Dt7v4ccCrQJbaxRCrOCUdk8t7oftw4uANfLtnAwPsn88AHi9i5RxvxSWKIpBDsCX/fYmadgQxCexOIVBmp1ZK4sn8bPrrhBE7u1JgHP1rMwPsn8/7cNVS2cTSR0oqkEDxhZvWAW4AJwDxCN5mJVDlNMmrw8Pk9ePnyo6mZmszIF3L5zTPTWLZhe9DRRGImprOGzGww8CCQDIx1978e8PgI4Mbw4Y/Ale4+s6RrarBYKsqevYU8P3UFf/9gEbsKCrns+NZcM6AtNVMj3dhPJH6Ud9ZQdeBsft6qEgB3v/MQr0sGFgGDgDxgGnC+u88r8pxjgPnuvtnMhgC3u3ufkq6rQiAVbd22nfx14gJen/4DTTLS+NOpHTmlS2PMtFWmVB7lnTX0FjAMKAC2F/k6lN7AEndf6u67gVfC19nP3ae4++bw4VeEZiSJxJXDaqdx/7ndGTeqL/VqpnL1S9MZMfZrFq/V5DmpGiJp42a5++AyXLsZsLLIcR5Q0qf9S4GJxT1gZiOBkQAtWrQoQxSR8stuVZ9/X3scL329gv+btJAhD37Ob45pxXUD21E7LSXoeCJlFkmLYIqZlWW6aHHt5mL7oczsREKF4MbiHnf3J9w9292zMzMzyxBFJDqSk4wL+7bik9/355zsLJ76chkD7pvMG9/maXaRVFoHLQRmNtvMZgHHAdPNbKGZzSpy/lDygOZFjrOAVcW8T1dgLDDM3TeWLr5IMBrUqs7dZ3XljauOpWlGGte/OpNzx0xVd5FUSgcdLDazliW90N1XlHhhs2qEBotPAn4gNFj8K3efW+Q5LYCPgYvcfUokgTVYLPGmsND5V+5K/jpxAdt37eX6QUdw+fGttTOaxJWyDhavB1a5+4rwL/004Cyg16GKAIC7FwDXAJOA+cBr7j7XzEaZ2ajw024FGgCPmtkMM9NveKl0kpKMXx7Vgg9+dwInHXkYf3tvAWc/rtaBVB4ltQg+Ay5198Vm1hb4BngR6Ah84+43V1zMn6lFIPHM3Xl71mpufWuOWgcSV8raIqjn7ovDP/8aeNndrwWGoP0IRIplZpzeralaB1KplFQIijYVBgAfAITvCdBefyIlaFirOo+O6MnD5/fg+43bOfWhL3js0+8o0DaZEodKKgSzzOxeM7seaAu8D2BmdSskmUglp9aBVBYlFYLLgQ2ElpY42d33LdLeEbg3xrlEqgy1DiTelWrROTPr6e7TY5jnkDRYLJXZhh93ccubc5g4Zw3dsjK495xutGtUO+hYkgDKu9ZQUWOjkEckYf1H62DTDk596Ase/XSJWgcSqNIWAi23KFJOB44d3PPeQs5+bIrGDiQwpS0Ed8QkhUgCUutA4kVEhcDMmoX3DthkZv3MrF+Mc4kkBLUOJB5EsjHN34BfEtqict9u3u7uQ2OcrVgaLJaq6sC7kkcPasfI4w/XXckSFSUNFkeyH8EZQHt33xXdWCJS1L7WQd82DbjlzTnc895CJs1Zo5lFEnORfNRYCmjXDZEKorEDqWiRtAh2ADPM7CNgf6vA3X8bs1QiCU6tA6lIkYwR/Lq48+7+XEwSHYLGCCTRaOxAoqGkMYKI7iw2sxpAC3dfGO1wpaVCIIlKdyVLeZTrzmIzOx2YAbwXPu5uZhOiG1FEDkVjBxIrkbQtbwd6A1sA3H0G0DqGmUTkIHTfgcRCJIWgwN23HnAu8pXqRCTq1DqQaIqkEMwxs18ByWbWzsweBiLaaF5EYketA4mWSArBtUAnQlNHXwbygdGxDCUikVPrQMqrtPsRJAPp7p4fu0gl06whkYPTzCI5mPLOGnrJzOqYWTowF1hoZn+IdkgRKT+1DqQsIuka6hhuAZwBvAu0AC6MaSoRKbPixg5Oe/gLcldsCjqaxKlICkGKmaUQKgRvufseNGtIJO7tax08fkFPtv60h7Mfm8pN42exefvuoKNJnImkEIwBlgPpwGdm1pLQgLGIxDkzY3DnJnz4uxO4/PjW/Cs3jwH3fcpr01ZSWKjPcxJSqsHi/S8yq+buBTHIc0gaLBYpuwVr8vnTG3PIWbGZ7Jb1uOuMzhzZpE7QsaQClGutITOrDpwNtKLIaqXufmcUM0ZMhUCkfAoLnXHT87j73fnk7yzg4mNaMXrQEdSqHslixFJZlWvWEPAWMAwoALYX+RKRSigpyTg3uzkf39Cfc7OzGPvFMgbeN5mJs1dTlh4CqfwiKQRZ7v5Ld7/H3e/b9xXJxc1ssJktNLMlZnZTMY93MLOpZrbLzH5f6vQiUmb10lO5+6yujL/yGOqlp3Lli9P5zTPTWLFRn/MSTSSFYIqZdSnthcM3nz0CDAE6AuebWccDnrYJ+C1wb2mvLyLR0atlPf59zbHcclpHcpZvYtADn/Hgh4vZVbD30C+WKiGSQnAckBv+ZD/LzGab2awIXtcbWOLuS919N/AKoS6m/dx9nbtPA/aUOrmIRE215CQuPa41H93Qn0EdG/HAh4sY/PfP+Xzx+qCjSQWIpBAMAdoBJwOnA6eFvx9KM2BlkeO88LlSM7ORZpZjZjnr1+sfpkisNM5I45Ff9eT5S3rj7lz41Ddc89J01ubvDDqaxNAhC4G7rwCaAwPCP++I5HWAFXe50sXbn+EJd8929+zMzMyyXEJESqHfEZm8N7of1w88gvfnreWk+ybz9BfLtFRFFRXJWkO3ATcCN4dPpQD/jODaeYQKyD5ZwKrSBhSRYKSlJHPdwHa8P7ofPVvW48635zH0H18y/fvNQUeTKIvkk/2ZwFDCU0bdfRUQyXKG04B2ZtbazFKB8wBtcSlSybRqmM5zFx/FoyN6smn7bs5+bAo3vz6bLTu0VEVVEckdJLvd3c3MAcKrkB6SuxeY2TXAJCAZeNrd55rZqPDjj5tZYyAHqAMUmtlofl7kTkTihJlxSpcm9Dsik79/sIhnpixn0tw13DykA8N7ZWFWXE+wVBaR3Fn8e0KDxYOAu4FLgZfc/aHYx/tvurNYJHjzV+fzpzfnkLtiM0e1qsefz+hC+8ba9yCelWuJifAFBhGaNQQwyd0/jGK+UlEhEIkPhYXOuNw87p4YWqri0uNac91J7UjXUhVxqUyFwMy28fMsnwPbfTuB74D/cfePohU0EioEIvFl0/bd/G3iAl7NWUmTjDRuO70jv+jUWN1FcaZMaw25e213rxP+ql30C2gMXAE8GKPMIlJJ1E9P5W/DuzL+yr5k1Ehh1D+nc8mz0/h+446go0mEIpk19F/cfa+7zwQejnIeEamkerWsz9vXHsefTj2Sb5ZtYtADk3n4Iy1VURmUqRDs4+5johVERCq/aslJXHb84Xx4wwkMPLIR932wiCF//5wvl2wIOpqUoFyFQESkOE0yavDIiJ48e/FR7HVnxNiv+e3L37JOS1XEJRUCEYmZ/u0PY9Lofowe2I735q7hpPsm8+yXWqoi3qgQiEhMpaUkM3rgEUwa3Y/uLepy+7/nMeyRL/lWS1XEjTLtWRwkTR8VqbzcnXdnr+HOt+eybtsujmxch27NM+iaVZeuWRm0b1Sbasn6fBoL5b6hLJ6oEIhUftt27uHZL5fz9eTNFbYAAAsbSURBVLJNzMrbQv7OAgCqV0uic7MMumZl0C1cHFo1SCcpSfcklJcKgYjELXdn+cYdzMrbwsyVW5mVt4U5q7ayc09oHKF2WjW6ZoVaDd3C35tkpOmGtVIqqRDoXnARCZSZ0bphOq0bpjOse2jvqoK9hSxe92OoOOSFisOTny2loDD0wbVhrer7i0LX5qHWQ/301CD/GJWaCoGIxJ1qyUkc2aQORzapwy+PCp3buWcv81bnM2vlFmblbWVm3hY+XriOfZ0aWfVq0K35z62Gzs0yqKV1jyKivyURqRTSUpLp2aIePVvU239u2849zP5hK7PCrYYZ32/hnVmrATCDtpm1Ql1K4QHpI5vUpnq15KD+CHFLhUBEKq3aaSkc06Yhx7RpuP/chh93MTvcYpiVt5XJi9YxfnoeACnJRofGdX4ejG6eQbvDapOc4IPRGiwWkSrN3Vm1dSezVv483jA7byvbdoVmKtVISaZzszr7p7B2b16XFvVrVrnBaA0Wi0jCMjOa1a1Bs7o1GNKlCRDaS2HZxu37ZyrNzNvCP79awa6C0EylNpnpDO/VnDN7NKNxRlqQ8SuEWgQiIsCevYUsXLON6d9vZsKMVeSs2EySwfHtMhneK4tBHRuRllJ5xxd0H4GISCkt27Cd8bl5vD49j1Vbd1I7rRqnd2vK8F5Z9Ghet9J1HakQiIiUUWGhM3XpRsbl5jFxzmp27ink8Mx0hvfK4qweWZWm60iFQEQkCrbt3MPE2WsYl5vHN8s3YQbHtW3I8F5Z/KJT47juOlIhEBGJsuUbtvP69DzGT/+BH7b8RO3q1Tgt3HXUs0X8dR2pEIiIxEhhofPVsnDX0ew1/LRnL60bhrqOzuzRjKZ1awQdEVAhEBGpED/uKuDd2atDXUfL/rPr6OSOjamRGlzXkQqBiEgF+37jDsZPz2P89DzyNv9ErerVOK1rE4b3yqJXy3oV3nWkQiAiEpDCQufrZZv2zzrasXsvrRrUDHUd9cyiWQV1HakQiIjEge27Cpg4Zw3jclfy1dJQ19ExbRowvFcWgzs1iWnXkQqBiEicWbnp566jlZtCXUendmnC8OwssmPQdRRYITCzwcCDQDIw1t3/esDjFn78FGAH8Bt3n17SNVUIRKQqKSx0pi0PdR29MzvUddSyQU3O7pnFWT2bkVWvZlTeJ5BCYGbJwCJgEJAHTAPOd/d5RZ5zCnAtoULQB3jQ3fuUdF0VAhGpqrbvKuC9OaEb1qYu3QgU6Trq3JiaqWVfJzSoQtAXuN3dfxE+vhnA3e8u8pwxwKfu/nL4eCHQ391XH+y6KgQikghWbtrBG9/+wLjcPL7ftIP01GSuH3QElx1/eJmuF9Qy1M2AlUWO8wh96j/Uc5oB/1EIzGwkMBKgRYsWUQ8qIhJvmtevyW9Pase1A9oybflmxufm0SQjNjOMYlkIihvpOLD5EclzcPcngCcg1CIofzQRkcrBzOjduj69W9eP2XskxezKoU/3zYscZwGryvAcERGJoVgWgmlAOzNrbWapwHnAhAOeMwG4yEKOBraWND4gIiLRF7OuIXcvMLNrgEmEpo8+7e5zzWxU+PHHgXcJzRhaQmj66MWxyiMiIsWL6Z7F7v4uoV/2Rc89XuRnB66OZQYRESlZLLuGRESkElAhEBFJcCoEIiIJToVARCTBVbrVR81sPbCijC9vCGyIYpxoiddcEL/ZlKt0lKt0qmKulu6eWdwDla4QlIeZ5RxsrY0gxWsuiN9sylU6ylU6iZZLXUMiIglOhUBEJMElWiF4IugABxGvuSB+sylX6ShX6SRUroQaIxARkf+WaC0CERE5gAqBiEiCS5hCYGaDzWyhmS0xs5uCzgNgZk+b2TozmxN0lqLMrLmZfWJm881srpldF3QmADNLM7NvzGxmONcdQWcqysySzexbM3s76Cz7mNlyM5ttZjPMLG72eDWzumY2zswWhP+d9Y2DTO3Df0/7vvLNbHTQuQDM7Prwv/k5ZvaymaVF9fqJMEZgZsnAImAQoc1wpgHnu/u8gHP1A34Ennf3zkFmKcrMmgBN3H26mdUGcoEz4uDvy4B0d//RzFKAL4Dr3P2rIHPtY2a/A7KBOu5+WtB5IFQIgGx3j6ubo8zsOeBzdx8b3q+kprtvCTrXPuHfGT8Afdy9rDewRitLM0L/1ju6+09m9hrwrrs/G633SJQWQW9gibsvdffdwCvAsIAz4e6fAZuCznEgd1/t7tPDP28D5hPaSzpQHvJj+DAl/BUXn2TMLAs4FRgbdJZ4Z2Z1gH7AUwDuvjueikDYScB3QReBIqoBNcysGlCTKO/kmCiFoBmwsshxHnHwi60yMLNWQA/g62CThIS7X2YA64AP3D0ucgF/B/4IFAYd5AAOvG9muWY2MugwYYcD64Fnwl1pY80sPehQBzgPeDnoEADu/gNwL/A9sJrQTo7vR/M9EqUQWDHn4uKTZDwzs1rAeGC0u+cHnQfA3fe6e3dC+1v3NrPAu9TM7DRgnbvnBp2lGMe6e09gCHB1uDsyaNWAnsBj7t4D2A7ExbgdQLiraijwr6CzAJhZPUI9GK2BpkC6mV0QzfdIlEKQBzQvcpxFlJtWVU24D3488KK7vx50ngOFuxI+BQYHHAXgWGBouD/+FWCAmf0z2Egh7r4q/H0d8AahbtKg5QF5RVpz4wgVhngxBJju7muDDhI2EFjm7uvdfQ/wOnBMNN8gUQrBNKCdmbUOV/vzgAkBZ4pb4UHZp4D57n5/0Hn2MbNMM6sb/rkGof9BFgSbCtz9ZnfPcvdWhP5tfezuUf3EVhZmlh4e7Cfc9XIyEPgMNXdfA6w0s/bhUycBgU5EOMD5xEm3UNj3wNFmVjP8/+ZJhMbtoiamexbHC3cvMLNrgElAMvC0u88NOBZm9jLQH2hoZnnAbe7+VLCpgNAn3AuB2eH+eID/F96DOkhNgOfCMzqSgNfcPW6masahRsAbod8dVANecvf3go2037XAi+EPZkuBiwPOA4CZ1SQ0u/CKoLPs4+5fm9k4YDpQAHxLlJeaSIjpoyIicnCJ0jUkIiIHoUIgIpLgVAhERBKcCoGISIJTIRARSXAqBJLQzOyBoitMmtkkMxtb5Pg+M7u1tCvWmtmzZjY8mllFYkWFQBLdFMJ3aZpZEtAQ6FTk8WOASe7+1wCyiVQIFQJJdF/y8+36nQjdebvNzOqZWXXgSKCbmf0D9n/Sf8jMppjZ0n2f+i3kH2Y2z8zeAQ7b9wZmdlJ4cbXZ4T0oqptZbzN7Pfz4MDP7ycxSw3suLK3AP7+ICoEktvBaPAVm1oJQQZhKaKXVvoT2FpgF7D7gZU2A44DTgH0thTOB9kAX4HJ+bmWkAc8Cv3T3LoTu8L2S0F2iPcKvPZ5QAToK6EOcrPQqiUOFQOTnVsG+QjC1yPGUYp7/prsXhjfqaRQ+1w94Obw66irg4/D59oQWDFsUPn4O6OfuBcASMzuS0EJw94evcTzwebT/gCIlUSEQ+XmcoAuhT+ZfEWoRHEOoSBxoV5Gfiy5xXtx6LcUtgb7P54RWutwDfEiolXEc8FmkwUWiQYVAJPTL/jRgU/gT/SagLqFiMDXCa3wGnBfeOKcJcGL4/AKglZm1DR9fCEwu8prRwFR3Xw80ADoAgS+IKIklIVYfFTmE2YRmC710wLla7r4hvHrnobwBDAi/bhHhX/buvtPMLgb+Fd5mcBrwePg1XxPqWtrXAphFaIMbrQQpFUqrj4qIJDh1DYmIJDgVAhGRBKdCICKS4FQIREQSnAqBiEiCUyEQEUlwKgQiIgnu/wNMtefINggeHQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dres_fig, dres_ax = plt.subplots()\n", "plt.plot(dres_conv)\n", "dres_ax.set_xlabel('Window')\n", "dres_ax.set_ylabel('Jensen-Shannon divergence')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparing different dimensionality reduction methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, you may want to compare the performance of different methods.\n", "\n", "Principal component analysis uses singular value decomposition to project data onto a lower dimensional space. [(See the scikit-learn user guide for more information.)](https://scikit-learn.org/stable/modules/decomposition.html#pca)\n", "\n", "The method provided by MDAnalysis.encore accepts any of the keyword arguments of [sklearn.decomposition.PCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html) *except* `n_components`. Instead, use `dimension` to specify how many components to keep." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "pc1 = drm.PrincipalComponentAnalysis(dimension=1,\n", " svd_solver='auto')\n", "pc2 = drm.PrincipalComponentAnalysis(dimension=2,\n", " svd_solver='auto')\n", "pc3 = drm.PrincipalComponentAnalysis(dimension=3,\n", " svd_solver='auto')" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(9, 3)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dres_conv2 = encore.dres_convergence(u, # universe\n", " 10, # window size\n", " select='name CA',\n", " dimensionality_reduction_method=[pc1, pc2, pc3]\n", " )\n", "dres_conv2.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, the size of the subspace you choose to include in your similarity comparison, affects the apparent rate of convergence over the trajectory." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1xV9f/A8deHLSiggKLgnoiKe+89U9Nypbmz0rTUym+Wo2znNsss9640c+fee6Io4gZRwYWDzef3x7n6I0O4cM/lXuDzfDzuQ++553zO+9qj+z6fLaSUKIqiKDmXjaUDUBRFUSxLJQJFUZQcTiUCRVGUHE4lAkVRlBxOJQJFUZQczs7SAaSXp6enLFasmKXDUBRFyVKOHTsWKaX0SumzLJcIihUrxtGjRy0dhqIoSpYihLj2ss9U05CiKEoOpxKBoihKDmfWRCCEaC2EuCCECBFCfJzC56OFECcNr0AhRKIQIp85Y1IURVH+zWx9BEIIW2AW0AIIBY4IIdZKKc89O0dK+R3wneH8DsD7Usp75opJURQlNfHx8YSGhhITE2PpUDLMyckJX19f7O3tjb7GnJ3FNYEQKeVlACHEcqAjcO4l5/cAlpkxHkVRlFSFhoaSJ08eihUrhhDC0uGkm5SSu3fvEhoaSvHixY2+zpxNQz7AjWTvQw3H/kMI4Qy0Bv4wYzyKoiipiomJwcPDI0smAQAhBB4eHumu0ZgzEaT0L/mypU47APte1iwkhBgshDgqhDgaERGhW4CKoigvyqpJ4JmMxG/OpqFQoHCy977AzZec251UmoWklHOAOQDVq1fP0LrZlx5cYtPVTbg5uOHm6Iarg6v2p6Or9ncHN+xtjW9TUxRFyS7MmQiOAKWFEMWBMLQf+54vniSEcAMaAW+YMRYuPrjIz6d+Rr60UgLOds7/ShLP/u7q6PqfBJL87852zln+KUJRFOvQv39/1q1bR/78+QkMDASgb9++7Nq1C1dXV6Kjo6lduzZfffUVPj4ptranm9kSgZQyQQgxFNgM2AK/SSnPCiGGGD7/yXBqZ2CLlPKJuWIBaF2sNS2KtOBx/GMexj4kKi6Kh7EPtVfcQ6Jio3gYp72Pio0iKi6Kyw8uPz8WnxT/0rLthN3/1yxeSBJuDm4v/czVwRU7myw3uVtRFDPq27cvQ4cOpU+fPv86/t1339G1a1eklEydOpUmTZoQGBiIg4ODyfc066+QlHIDsOGFYz+98H4+MN+ccTxja2P7/Mc4PaSUxCTGPE8cUXFRzxNH8gTy7LOIpxFcenCJqNgoHsU/SrVsF3sXPJw8qFagGvV86lG7YO10x6coSvbRsGFDrl69+tLPhRC8//77rF69mo0bN9KxY0eT75lzHkfjnkDUTfAsne5LhRDksstFLrtceLt4p+vahKQEHsU9+nctJFkCiYqNIvxJOFuvb2V1yGpshA2VPCtRz6ce9X3qU96jPDZCTQBXlMw24e+znLsZpWuZ5Qu5Mq6Dvy5lVa1alfPnz6tEkC7Bm+H3fpCvJJRtA2VaQ5HaYOYOYjsbO/I65SWvU95Uz0tISiAwMpC9YXvZF7aPH0/+yKyTs8jrmJfahWpT36c+dQvVxTOXp1njVRQla9Bzv/mckwiK1oW238OFjXB4DhyYCU5uUKo5lGkDpZtDrtR/rM3JzsaOyvkrUzl/ZYZWGcq9mHscuHmAfWH72HdzHxuvbATAL58f9XzqUa9QPQLyB2Bvo0Y6KYo56PXkbi4nTpygWbNmupSVcxJBHm+oOUh7xT6CSzsgeJNWUwj8A4QtFKkDZVtricGzlEXDzeeUj3Yl2tGuRDuSZBIX7l1g38197A3by/zA+cw9MxcXexdqedd63oxUKHchi8asKIr5SSmZMWMG4eHhtG7dWpcyhZ7Vi8xQvXp1qet+BEmJEHYcgjfChU1w56x2PJObkNLjUdwjDocfZu9NrRkp/Ek4AMXdilOvkJYUqhWohpOdk4UjVZSsJSgoCD8/P4vG0KNHD3bu3ElkZCQFChRgwoQJ7Nmz5/nw0adPnz4fPurr65tiGSl9DyHEMSll9ZTOV4ngRfevabWE4E1wdQ8kxhmakFpoiaFUM4s2Ib1ISsmVqCtaE1LYPo7cOkJcUhyOto5U965OvUL1qOdTj+KuxdVcB0VJgzUkAj2oRKCnF5uQnkZqTUhF60KZVlbRhPSi6IRojt0+9rxv4crDKwAUcimk9S341KOWdy1yO+S2cKSKYn1UIsgiMjURJJeUCGHHtM7m4E1wx7CIqkcprfmobBsoXBtsravbJexx2PPawsHwgzxNeIqd0Dqmn3U6l81XVg1RVRRUIsgyLJYIXvS8CWkjXNkDSfFW3YQEEJ8Yz8mIk+wL28f+m/sJuhcEgIeTx/OkUKdQnTSHuipKdqUSQRZhNYkgudhHcGm71tl8cTM8vZusCam19rKyJiSAyOhI9t/cz96wvRy4eYAHsQ8QCCp4VqBuobrU96lPgFeA6ltQcgyVCLIIq0wEySUlQuhRQ79C1mlCSkxK5Nzdc+y9uZf9Yfs5HXmaJJmEXz4/BlUaRLMizVTzkZLtqUSQRVh9InjR/ataE9KFjXB1r6EJyR1Kt9ASgxU2IQE8jH3I1mtbmXd2HteirlHcrTgDKgygbYm2ahKbkm2pRJBFZLlEkFxMFFze8d8mpCJ1oFh9bb6Cbw1wtJ4RPYlJifxz7R9+OfMLwfeDKeRSiH4V+tGpVCc1T0HJdqwhEdy4cYM+ffpw69YtbGxsGDx4MMOHD0/XUtQqEWQVz5uQNsLFrXA7EJBaYihYSUsORepoySF3fktHi5SS3aG7mXNmDqcjTuPh5EEf/z50K9sNF3sXS4enKLqwhkQQHh5OeHg4VatW5dGjR1SrVo01a9bw7bff0r59+38tRT179uwUl6JObyKwrobqnMTGForU0l7Nx0PMQ7hxBK4fgOsH4ehvcPBH7dx8Jf8/KRStC/lKQCZ34AohaFS4EQ19G3L09lF+Of0LU45N4dczv9LTrye9yvXC3ck9U2NSlOyoYMGCFCxYEIA8efLg5+dHWFjYv87ReylqlQishZObtvBd6eba+4RYCD/1/4nhwno4uVj7zMVLSwrPkoN3QKZ1PgshqOFdgxreNTgTcYa5Z+by06mfWHB2Aa+XeZ0+/n3I72z5GoyimGzjx3DrjL5leleENl8bffrVq1c5ceIEtWrVYsmSJf/5XK+lqFUisFZ2jlC4pvaqNxySkiAy+P8Tw/UDEPS3dq69C/hW///EkEn9DBW9KjKt6TRC7ofwa+CvLA5azNLzS+lUqhP9KvSjcJ7CaReiKEqKHj9+TJcuXZg6dSqurq4pnqNX075KBFmFjQ3kL6e9qvfTjkXd/Hdi2PUN/+1nMNQczNjPUCpvKb5q8BXvVH6HeYHzWBOyhj8v/kmb4m0YWHEgJd1Lmu3eimI26Xhy11t8fDxdunShV69evPrqqy89T6+lqFUiyMpcC0GFLtoLtH6G0CNaYrh24OX9DEXqgEdJ3fsZCucpzGd1PmNIwBAWnF3AquBVrLu8jmZFmjGo4iD8Pa17fXdFsQZSSgYMGICfnx8ffPDBS8/RcylqNWooO0uIS9bPYHhF39c+y4R+hgcxD1hyfglLgpbwKO4RdQvVZWDFgVQvUF3NVlaskjWMGtq7dy8NGjSgYsWK2Nhokzi//PJLVq5cafRS1FY1fFQI0RqYBtgCc6WU/6lrCSEaA1MBeyBSStkotTJVIjBBUhLcvfj/zUnX9sODa9pnZuxneBz3mJXBK1l4diF3Y+5S2asygyoNooFPA5UQFKtiDYlAD1aTCIQQtkAw0AIIBY4APaSU55Kd4w7sB1pLKa8LIfJLKe+kVq5KBDqLumnoYzD0M9wOBJlk6GcIgIAeENBNG9VkopiEGFaHrGZe4DzCn4RTNm9ZBlYaSIsiLbC1sdXhyyiKaVQi0JkQog4wXkrZyvB+DICU8qtk57wDFJJSjjW2XJUIzCwmCkIPa4nh4j8QflKrLVR6HWoM0Ia/mSg+KZ4Nlzcw98xcrkZdpZhrMfpX6E/7Eu2xt6Kd4JScJ6cmAnOuIuYD3Ej2PtRwLLkyQF4hxE4hxDEhRJ+UChJCDBZCHBVCHI2IiDBTuAoATq5Qqjk0HQtv7YJBO6BCZzi1DH6qD7+2hNMrtXkOGWRvY0/HUh1Z03ENPzT6gVx2ufhs/2e0Xd2WJUFLiEmI0fELKYqSFnMmgpQaf1+sftgB1YB2QCvgUyFEmf9cJOUcKWV1KWV1Ly8v/SNVXs6nKnScBSPPQ6uv4Ekk/DkIJvvB1vHavgwZZGtjS8tiLVnRfgWzm8+mkEshvj78Na3+aMXcM3N5FPdIv++hKMpLmTMRhALJZxT5AjdTOGeTlPKJlDIS2A0EmDEmJaNy5YU678DQo9B7jdapvG8aTAuAJa9D8BZt/aQMEEJQ36c+C9osYH7r+fh5+DHt+DRa/d6KGSdmcD/mvs5fRlGU5Mw5j+AIUFoIURwIA7oDPV845y9gphDCDnAAagFTzBiTYiobGyjZRHs9DIVjC+D4Alj6GrgXher9oUpvcPHIUPHVClSjWoFqnLt7jrln5vLL6V9YdG4RXUp3oa9/Xwq4FND5CymKYrYagZQyARgKbAaCgJVSyrNCiCFCiCGGc4KATcBp4DDaENNAc8Wk6MzNF5p+AiMCoes8cC8CW8fB5HLw52C4cRgyOBihvEd5JjeezJqOa2hRtAXLzi+j9Z+tGb9/PNejruv8RRTFesTExFCzZk0CAgLw9/dn3LhxAPTt25fixYsTEBBAmTJl6NOnz38Wo8soNaFM0ded83D0Vzi5DOIeaaOMagyEiq+BQ8aXqw57HMa8wHmsvriaBJnAe1Xeo3+F/moegqIraxg1JKXkyZMn5M6dm/j4eOrXr8+0adP46aefzLYMdZo1AiFEGSHENiFEoOF9JSGE0cM9lRwmfzlo+53Wudx+ilYj+Hs4/FAONnwIERcyVKxPbh/G1h7L5q6baVG0BVOPT+WjPR8RnRCt8xdQFMsSQpA7tzaZMz4+nvj4+P888Dxbhtrb25uNGzeafE9j+gh+AUYDPwNIKU8LIZYCX5h8dyX7csyt9RdU66c1ER2ZC8fmweGfoVgDrZZQrh2kc96AZy5Pvmv4HeXylWP68elcfXiV6U2n4+3ibaYvouRU3xz+hvP3zutaZrl85fio5kdpnpeYmEi1atUICQnh3XffpVatWsyePfs/5+m1DLUxfQTOUsrDLxxLMOmuSs4hhLb5Tpdf4IMgbROe+9dg1ZswpQLs+BIepq+dUwjBwIoDmdlsJjce3aDbum4cv33cLOEriiXY2tpy8uRJQkNDOXz4MIGBKXedZuYy1JFCiJIY5gAIIboC4brcXclZXDyh/vtQ9z0I2arVEnZ9C7u/h3JtofoAKNHY6FVRG/o2ZEnbJby34z0GbBnAJ7U+oWuZrmb9CkrOYcyTu7m5u7vTuHFjNm3alOLnei1DbUyN4F20ZqFyQogwYATwtsl3VnIuG1so0wp6rYL3TkDdoXB1HyzqBDOrw4Ef/3+V1DSUcC/B0nZLqVWwFhMOTGDSwUnEJ8Wb+QsoivlERETw4MEDAKKjo9m6dSvlypX71zlSSqZPn67bMtRpJgIp5WUpZXPACygnpawvpbxq8p0VBSBfcWgxUWs26jwHcuWDzWPgBz/4ayjcPJlmEa4OrsxqOot+/v1YfmE5g7cM5l7MvUwIXlH0Fx4eTpMmTahUqRI1atSgRYsWtG/fHoDRo0c/Hz565MgRduzY8Z8RQxmR5vBRIcSXwLdSygeG93mBkelZKE5PGR0+Gh2XyPIj13mzTjFsbNSQQ6sWfgqO/ApnVkH8U/Cpri14598Z7HOleum6y+sYv388Hk4eTG86nbL5ymZS0Ep2YA3DR/VgjkXn2jxLAgBSyvtAW5OitIC/T99kwt/n+PSvQN06WBQzKRgAr0zXagltvoXYKFjztra+0ZaxcO/ySy9tX6I9C1ovIEEm0Htjb7Zc3ZKJgStK1mRMIrAVQjg+eyOEyAU4pnK+VXqtmi9vNy7JkkPXGb/2rEoGWUEud6j1Frx7GN78G4o31PoPpleFk0tfepm/pz8r2q+gTN4yjNw1khknZpAkkzIxcEXJWowZNbQY2CaEmIc2cqg/sMCsUZmBEIIPW5UlMUkyZ/dlbGwEn7Uvr2amZgVCaEmgeEOICofVb8HaYZDHG0o2TfESz1ye/NbqN744+AVzTs8h+H4wX9X/itwO+uy6pmRfUsos/buQkYdcYzqLvwUmAX6AP/C54ViWI4RgTJty9K9XnHn7rvLlhiBVM8hqXAtCt0XgWRZW9IFbZ156qoOtAxPqTmBMzTHsCd3DGxveUOsUKalycnLi7t27WfZ3QUrJ3bt3cXJyStd1OXKtISkl49eeZcGBawxpVJKPWpfN0k8AOdLDMJjbXPv7wK3g9uKeR/92KPwQI3eNRErJd42+o26hupkQpJLVxMfHExoaSkxM1t0cycnJCV9fX+zt/z1r36StKoUQrwLfAPnRNpsRgJRSuuoSdTrpteiclJJP/wpk8cHrDG1SipEty6hkkNXcCoR5bcCtMPTfmOa+yjce3WD4juFcenCJkdVG0rt8b/XfXMkxTB019C3wipTSTUrpKqXMY6kkoCchBBNfqUCPmoWZuSOEadsuWjokJb28K8DrCyHyAqzoDQlxqZ5eOE9hFrdZTJPCTfju6HeM3TeW2MSMb7mpKNmFMYngtmHfgGzHxkYwqVNFXqvmy9StF5mhkkHWU7IJvDIDruyCv99Lc/8DZ3tnJjeezDuV32HtpbX029SPO0/vZFKwimKdjBk1dFQIsQJYAzx/fJJS/mm2qDKRjY3g6y6VSJSSH/4JxtZW8E7jUpYOS0mPyj213dJ2TNKaiZp+kurpNsKGtwPepox7GcbsHUP3dd2Z2mQqlbwqZVLAimJdjKkRuAJPgZZAB8OrvTmDymy2NoLvugbQsXIhvt10gTm7L1k6JCW9Go6GKm/A7m/h+EKjLmlWtBmL2y7GwdaBvpv68lfIX2YOUlGsU5o1Aillv8wIxNJsbQQ/vBZAYpLkyw3nsRGCgQ1KWDosxVhCQPup2jyDv0dAnkJQunmal5XJW4bl7ZYzatcoxu4by/l75xlZfSR2NubczltRrIvaoSwZO1sbpnarTNuK3nyxPoj5+65YOiQlPWzt4fUFUKC8tt9B+CmjLnN3cuenFj/xht8bLA5azNtb3+Zh7EMzB6so1sOYpqFfgDFAPGg7lAHdzRmUJdnZ2jCtexVa+Rdg/N/nWHTgqqVDUtLDMQ/0XAVO7rDkdXhww6jL7Gzs+KjmR0ysO5Fjt4/RY30PQu6HmDlYRbEOZt2hTAjRWghxQQgRIoT4OIXPGwshHgohThpenxlTrrnZ29owo0dVmvsV4NO/zrL0kJqNmqW4FoQ3fof4aFjSFaIfpH2NQefSnfmt1W88jX9Krw292H59uxkDVRTrYEwiyNAOZUIIW2AW0AYoD/QQQpRP4dQ9UsrKhtdE40M3Lwc7G2b1qkLTcvn53+ozrDxi3JOlYiXy+0H3xXD3Eqx4AxKMny9QOX9llrdfTnG34gzfMZyfTv2UZZccUBRjmHOHsppAiGFjmzhgOWDaDsuZzNHOlh97VaVRGS8++vM0vx8LtXRISnoUbwidfoSre+Cvd9OcY5Cct4s381vPp32J9sw6OYuRu0byNP6pGYNVFMsx5w5lPkDyx+hQw7EX1RFCnBJCbBRC+KdUkBBisBDiqBDiaEREhBG31o+TvS0/965G/VKejP79FKtPqGSQpVR6HZp+qm1ysy19FU4nOye+rP8lo6qPYtv1bfTe2Juwx2FmClRRLMeYUUMfCCE+AN4CBhneDxBCVE7r0hSOvfhIdhwoKqUMAGagTVr770VSzpFSVpdSVvfy8korZN052dsyp3d16pTwYOTKU6w9dTPTY1BM0GAkVOsLeyfD0d/SdakQgjf93+THZj8S/jic7uu6c+TWEfPEqSgWYkzTUHVgCNrTvA8wGGgM/CKE+DCV60KBwsne+wL/+gWVUkZJKR8b/r4BsBdCeBodfSbK5WDL3DerU6NYPt5fcZL1p9PsJlGshRDQ9gco3RLWj4Tg9O9aVs+nHkvbLSWvU14GbRnEsvPLVL+Bkm0Ykwg8gKpSypFSypFoicELaAj0TeW6I0BpIURxIYQD2pDTtclPEEJ4C8Pyj0KImoZ47qb7W2QSZwc7futbg6pF3Hlv+Qk2BapkkGXY2kHXeeBdEVb1hZsn0l1EMbdiLG27lPo+9fny0JdMODCB+MR4/WNVlExmTCIoAiRf1jEerTknmmRrD71ISpkADAU2A0HASinlWSHEECHEEMNpXYFAIcQpYDrQXVr5Y5aLox3z+tUkwNeNoUtPsOXsLUuHpBjLMbc2x8DZQ5tjcP9auovI7ZCbaU2mMajiIP64+AcDtgwgMjrSDMEqSuYxZj+CT4HOwLOFWDqgPdn/AMyRUvYya4Qv0Gs/AlM9iomn96+HOXvzIT+9UY1mfgUsHZJirIgL8GsLyF0A+m8G53wZKmbTlU18uu9T3BzdmNZ0Gv4eKY51UBSrkOH9CAzNNvOBQcAD4CEwREo5UUr5JLOTgDXJ42TPgv418SvoytuLj7PzglrKOMvwKgvdl8H9q+meY5Bc6+KtWdhmITbChjc3vsmGyxv0jVNRMkmqicDQTLNGSnlMSjlNSjlVSmn5x3Er4ZbLnkX9a1G6QG4GLzrGnouZO7RVMUGxetBpNlzbB6uHQFJShorx8/BjWbtl+Hv489Gej5h+fDpJMmNlKYqlGNNHcFAIUcPskWRRbs72LB5Qi5JeuRm44Cj7Q1R7cZZRsSs0nwBn/4Rt4zNcjEcuD+a2nEuX0l345cwvjNo1iuiEaP3iVBQzMyYRNEFLBpeEEKeFEGeEEKfNHVhWktfFgSUDa1Hc04X+C45w8LLVDnxSXlRvOFQfAPumweFfMlyMva094+qMY1T1UWy9tpX+m/qrTmQlyzCms7hoSsellOkfcqEDa+ksTknk41h6zDlI6P1oFvSvSc3iGeuEVDJZYoLWV3BxM3RbAuXamlTc9uvb+XjPx7g5ujGz6UzK5iurU6CKknEmbV5v+MEvDDQ1/P2pMdflRJ65HVk6qDaF3J3oO+8wR6/es3RIijFs7aDrr1AwAH7vD6HHTCquaZGmLGi9gCSZRJ+NfdgdulunQBXFPIxZYmIc8BHangQA9sBicwaVlXnlcWTZoNp4uzrRd94Rjl+/b+mQFGM4uEDPlZA7Pyx9He6ZtimRn4cfS9supahrUYZtH8aic4vUTGTFahnzZN8ZeAV4AiClvAnkMWdQWV1+VyeWDqqNR24H3vz1MKduGL8evmJBufPDG3+ATNT2MXhqWo2ugEsB5reeT5PCTfj2yLd8cfAL4pPUTGTF+hiTCOIMw0if7UfgYt6QsgdvNyeWDaqNu4s9vX89xJlQtfVhluBZWptj8OAGLOuubW5jAmd7ZyY3nkz/Cv1ZGbySd7e+S1RclE7BKoo+jEkEK4UQPwPuQohBwFa07SuVNBRyz8WyQbXJ42TPG78e4uxNlQyyhKJ14NWf4cYhWP1WhucYPGMjbHi/2vtMrDuRI7eO0HtDb248UhsdKdbDmM7i74HfgT+AssBnUsoZ5g4su/DN68zywbVxcbDljbmHCApXT4NZgn9naPkFnPsL/vlUlyI7l+7MnJZziIyOpNf6Xpy4k/6F7xTFHIzpLH4fCJJSjpZSjpJS/pMJcWUrhfM5s2xwbRztbOk19xAXbj2ydEiKMeoMhZpvwYGZcPAnXYqs4V2DJW2X4OroyoDNA1h3eZ0u5SqKKYxpGnIFNgsh9ggh3hVCqNXVMqCohwvLBtfGzkbQa+5BQu6oZGD1hIDWX0HZdrDpYwj6W5dii7kVY0nbJQR4BTBmzxhmnpipRhQpFmVM09AEKaU/2t7FhYBdQoitZo8sGyruqSUDIQQ9fjnEpYjHlg5JSYuNLXSZCz7V4I+BcEOf3cncHN2Y02IOnUp14ufTP/Ph7g+JSYjRpWxFSa/0TAy7A9xC2zgmv3nCyf5KeuVm2aBaSCnpMecgVyKfWDokJS0OztBzBeQpCMu6wd1LuhRrb2vPxLoTGVF1BJuublJ7GygWY0wfwdtCiJ3ANsATGCSlrGTuwLKzUvnzsHRQbRKStGRw7a5KBlbPxdMwx0Bqcwye6PODLYRgQMUBTGk8heB7wfRa34uL9y/qUraiGMuYGkFRYISU0l9KOU5Kec7cQeUEZQrkYcnAWsQmJNJjzkECw9TQUqvnUVKrGUTd1OYYxD3VrejmRZszv8184pPi6b2xN3tC9+hWtqKk5aWJQAjhavjrt8B1IUS+5K/MCS978yvoyuKBtUiUkk6z9vHjzhASk1SnoVUrXBNe/QVCj8KfgyApUbei/T38WdpuKYXzFGbo9qEsCVqiW9mKkprUagRLDX8eA44a/jyW7L2iA/9Cbmwe0ZCW/gX4dtMFesw5yI17+j1pKmZQ/hVtNNH5dbD5f1pzkU68XbxZ0HoBDX0b8vXhr5l0cBIJSQm6la8oKXlpIpBStjf8WVxKWcLw57NXCWMKF0K0FkJcEEKECCE+TuW8GkKIRCFE1/R/hazP3dmBWT2r8sNrAZwLj6LttD2sPhGqhhRas9pvQ+134NBPcPBHXYt2tndmauOpvFn+TZZfWM7Q7UN5HKdGmCnm89L9CIQQVVO7UEp5PNWChbAFgoEWQChwBOjxYh+D4bx/gBjgNynl76mVa837Eejhxr2nvL/iJEev3ad9pYJM6lQRN2d7S4elpCQpCVa9qc0veG0++HfS/Ra/B//OpIOTKOZWjJnNZuKT20f3eyg5Q0b3I/jB8JoFHALmoK0xdAiYbsR9awIhUsrLUso4YDnQMYXzhqEtX6F2f0ebhbzirTqMblWWTYG3aD1tt9r+0lrZ2MCrc8C3Bvw5GK4f1P0WXct0ZXaL2dx+epue63ty8s5J3e+hKKk1DTWRUjYBrgFVpZTVpZTVgCpAiBFl+wDJV9YKNRx7Tgjhg7bMdarz94UQg4UQR4UQR/t1vTsAACAASURBVCMisv8G8bY2gneblOLPd+qSy96WnnMPMWn9OWIT9OuYVHRinwt6LAc3X20kUcQF3W9Ru2BtFrddjIu9CwM2D2DjlY2630PJ2YwZPlpOSnnm2RspZSBQ2YjrRArHXmyHmgp8JKVM9RdOSjnHkIiqe3l5GXHr7KGSrzvr3qvPG7WL8MueK3ScuY/zt9SidVbHxQPe+B1s7GF+e4gI1v0WJdxKsKTtEip4VuDD3R8y+9Rs1Yek6MaYRBAkhJgrhGgshGgkhPgFCDLiulC0LS6f8QVuvnBOdWC5EOIq0BX4UQihf0NrFubsYMcXnSryW9/qRD6O5ZWZ+/h17xWS1DBT65KvBLxpWItofjuzJIO8Tnn5peUvvFLyFX48+SNj9o4hNjFW9/soOY8xm9c7AW8DDQ2HdgOzpZSpLowihLBD6yxuBoShdRb3lFKefcn584F1Ob2zODWRj2P5+I/TbA26Q/1Snnz/WgDebk6WDktJ7s55WNBBW7DuzXXgVUb3W0gpmXtmLtNPTKeyV2WmNpmKRy4P3e+jZC+mbl4fI6WcIqXsbHhNSSsJGK5LAIYCm9FqECullGeFEEOEEEPS+yUU8MztyC99qvNl54ocu3afVlN3s+FMuKXDUpLLX06rGcgkWGCeZiIhBIMqDeKHRj8QdC+IXht6cemBPusfKTlTmjUCa5OTawTJXY54zPsrTnIq9CFdqvoy/pXy5HFSw0ytxp3zWiIQNmarGQAERgYybPswYhJi+KHRD9T1qWuW+yhZn0k1AsU6lfDKze9v1+W9pqVYfSKUNtP2cOSqaZutKzrKX05LAM9qBpHmWUiugmcFlrZdSqHchXhn2zusOL/CLPdRsjeVCLIwe1sbPmhZllVD6iAEdPv5AN9vvkB8oml77Co6Sd5MNL+d2ZJBwdwFWdhmIfV86vHFoS/45vA3JOq4BpKS/RnTWVwGGI22Cqnds+NSyqbmDS1lqmkoZY9jE5iw9iyrjoVSydeNKd0qU9Irt6XDUgDuBBk6kG2h7zrwLG2W2yQmJfL90e9ZHLSYhr4N+bbht7jYu5jlXkrWk1rTkDGJ4BTahK9jwPPHDCnlMT2DNJZKBKnbFBjOx3+eISY+kU/aleeNWkUQIqUpHUqmuhOkzTGwsTNrMgBYeWElXx76khLuJZjVdBYFcxc0272UrMPURHDMMKPYKqhEkLbbUTGMWnWKPRcjaVouP990qYRXHkdLh6VkYjLYf3M/o3aOwsHWgRlNZ1DRq6LZ7qVkDaZ2Fv8thHhHCFFQ7UeQNRRwdWJBv5qM71CevSGRtJ66m63nbls6LCW/n5YAkhK0hBBpzEotGVO3UF0WtV2Ek50T/Tb3Y3/YfrPdS8n6jKkRXEnhsDR2KWq9qRpB+gTffsTw5ScJCo+iZ60ijG3nh7ODXdoXKuZz+5zWZ2BjB33Xg2cps93qXsw9Bm0ZxM3HN1nUZhGl8prvXop1M6lpyNqoRJB+sQmJTP4nmDm7L1PMw4Wp3SoTUNjd0mHlbM+Sga29NszUjMng1pNb9FjfA0dbR5a2W0o+J1Whz4lMahoSQtgLId4TQvxueA0VQqiZS1mIo50tY9r4sXRgbWLjE3l19n5mbLtIghpmajkFymtDSxPjtXkGd803M9jbxZvpTaYTGR3J8O3D1fpEyn8Y00cwG6gG/Gh4VTMcU7KYOiU92DiiIe0rFeSHf4LpNucg1++qbTEt5nkyiNPmGZgxGVT0qsik+pM4GXGS8fvHq5VLlX8xJhHUkFK+KaXcbnj1A2qYOzDFPNxy2TOtexWmda9M8O1HtJm2m1VHb6gfBkspUF5rGsqEZNCqWCuGVRnGusvrmHN6jtnuo2Q9xiSCRCFEyWdvhBAlSDafQMmaOlb2YdOIhlTwcWP076d5Z8lx7j+Js3RYOVMm1gwGVRxE+xLtmXlyJpuubjLbfZSsxZhEMBrYIYTYKYTYBWwHRpo3LCUz+LjnYumg2oxpU46tQbdpNXU3u4Oz/w5wVqmAf7JkYL4+AyEEE+pOoEr+KozdO5YzEWfSvkjJ9owaNSSEcATKou06dl5KabHeJjVqyDzO3nzIiOUnuXjnMX3rFuPjNuVwsre1dFg5z+2zhtFEjtqcA4+SaV+TAfdi7tFzfU9iEmJY1m6Zmn2cA+ix+mg1oAIQAHQTQvTRKzjFOvgXcuPvYfXpW7cY8/dfpcOMvZy9+dDSYeU8Bfyhz1pIjDVrzSCfUz5mNZtFbGIsQ7cP5Un8E7PcR8kajBk+ugj4HqiP1klcA22LSSWbcbK3Zfwr/izsX5OH0fF0mrWP7zdfIDpOdQllKu8KWjJIiDFrMijpXpLvG33PpQeX+Gj3R2rF0hzMmJnFQUB5aSXDSlTTUOa4/ySOievOsfpEGD7uufisQ3lali+gFrDLTLcCtWYiOyezNhMtP7+cSYcm0ad8H0bXGG2WeyiWZ2rTUCDgrW9IirXL6+LAlG6VWTG4Nnmc7Hhr0TH6zjvClUjVhJBpvCtoHcgJMVpCMFPNoHu57vTy68XCcwtZFbzKLPdQrJsxNYIdQGXgMPC8k1hK+Yp5Q0uZqhFkvoTEJBYeuMbkf4KJS0jirUYleKdxKXI5qM7kTPGsZmCfS6sZ5NN/ma+EpASGbR/GwZsHmd1iNrUL1tb9HoplmboMdaOUjkspd+kQW7qpRGA5d6Ji+Grj+efNReM6lKeFai7KHLfOwIJXzJoMHsc9pvfG3tx+cpvF7RZTws0i60oqZmKxReeEEK2BaYAtMFdK+fULn3cEPgeSgARghJRyb2plqkRgeYcu3+Wzv85y4fYjmpT1YlwHf4p5qp2wzC4TkkHY4zB6ru+Ji70LS9suxd1JLU6YXZi66NyrQoiLQoiHQogoIcQjIUSUEdfZArOANkB5oIcQovwLp20DAqSUlYH+wNy0ylUsr1YJD9a9V5+x7fw4cvU+LafsZvIWNbrI7LwrwptrIT4a5neAe5d1v4VPbh+mNZnG7Se3GbFzBPGJ8brfQ7E+xnQWfwu8IqV0k1K6SinzSCldjbiuJhAipbwspYwDlgMdk58gpXycbDSSC2AVI5OUtNnb2jCwQQm2j2xE24reTN8eQospu9hy9pZat8icnieDJ2ZLBpXzV+bzep9z7PYxJhyYoP575gDGJILbUsqgDJTtA9xI9j7UcOxfhBCdhRDngfVotYL/EEIMFkIcFUIcjYhQSyBYk/yuTkztXoXlg2vj7GDL4EXH6D//CNfuqtFFZuNdURtN9DwZpLR3lGnalmjL2wFv89elv/gt8Dfdy1esizGJ4KgQYoUQooehmehVIcSrRlyXUg/ifx4tpJSrpZTlgE5o/QX/vUjKOVLK6lLK6l5eXkbcWslstUt4sP69Boxt58fhK/doMWU3k/8JJiZeNReZhXdFbdJZ/BNt0pkZksHbAW/Tplgbph6fytZrW3UvX7EexiQCV+Ap0BLoYHi1N+K6UKBwsve+wM2XnSyl3A2UFEJ4GlG2YoWeNxeNakybCt5M33aR5pN3qf2SzaVgJbMmAyEEE+tNpJJXJcbsGcPZu2d1LV+xHmYbNSSEsAOCgWZAGHAE6CmlPJvsnFLAJSmlFEJUBf4GfFObxaxGDWUdBy7d5bO/Arl45zFNy+VnXIfyFPVQo4t0F34aFr4CDrm1JqN8xXUtPjI6kp7re5KYlMjSdksp4FJA1/KVzGHqPAInYADgDzg9Oy6lTLE9/4Vr2wJT0YaP/ialnCSEGGK4/ichxEdAHyAeiAZGq+Gj2Ut8YhLz911l6tZg4pMkbzcqyduNS6qVTfUWfgoWdtSSQd91kLeYrsUH3w+m94beFHUtyvzW83G2d9a1fMX8TE0Eq4DzQE9gItALCJJSDtc7UGOoRJA13Y6KYdL6INaeuknhfLkY196f5uXVk6Wuwk9p8wwc85glGewO3c2w7cNo7NuYKU2mYCOMXbxYsQamrjVUSkr5KfBESrkAaAdU1DNAJfsr4OrE9B5VWDqoFo52tgxceJQB84+oPZP1VDBAG1oa+0jrM7h/VdfiG/o2ZHT10Wy/sZ1px6fpWrZiWcYkgmczSh4IISoAbkAxs0WkZGt1S3qycXgD/te2HAcv36X5lF1M3apGF+nGzMmgl18vupXtxm+Bv7H64mpdy1Ysx5hEMEcIkRf4FFgLnEObZKYoGWJva8PghiXZNrIxrfy9mbr1Ii2m7GJbkBpdpIt/JYMOcP+abkULIfio5kfUKViHiQcmcuTWEd3KVizHrGsNmYPqI8h+9odE8tnas4TceUxzv/yM6+BP4XyqM9JkN09qHciOrtD3b137DKLionhjwxvci7nHkrZLKOpaVLeyFfMwtbPYEeiC1hxk9+y4lHKijjEaTSWC7CkuIYl5+64wbdtFEpMk7zQuxVuNSqjRRaZ6lgzsc0HvNZC/nG5F34i6Qc8NPXF3dGdx28W4ObrpVraiP1M7i/9CWyMoAXiS7KUounGws+GtRiXZNrIRLcoXYMrWYFpO2c3286q5yCSFKkO/jSAlzGsDYcd0K7qwa2GmNZlG2OMwRu4cSXySWqAuqzKmRhAopayQSfGkSdUIcoZ9IZF89lcglyKe0NyvAOM6lFfNRaa4dxkWdoKnd6HHcijeQLei115ayyd7P6FL6S6MqzNO7U9hpUytEewXQqjhokqmqlfKk43DG/Jxm3LsvxRJ88m7mLb1ohpdlFH5SkD/zeDmC4u7wIWNuhX9SslXGFRxEH9c/IOF5xbqVq6SeV5aIxBCnEFbJM4OKA1cRtuqUgBSSlkps4JMTtUIcp7wh9F8sT6I9afDKerhzPgO/jQpl9/SYWVNT+9piSD8FHSaDQHddCk2SSYxatcotl7byvSm02lcuLEu5Sr6yVBnsRAi1WEAUkr9xqSlg0oEOdfei5F8tjaQyxFPaFzWi0/a+lG6QB5Lh5X1xD6CZT3g6h5o+z3UHKRLsdEJ0fTb1I/LDy+zqM0iyuYrq0u5ij4ymgicgXgpZbzhfVmgLXBNSvmnuYJNi0oEOVtcQhLz919hxrYQnsYn0rNmEUY0L41HbkdLh5a1xMfA7/3hwnpoOhYajAId2vYjnkbQY30PAJa1W4aXs1o23lpktI9gE4YZxIZVQg8AJYB3hRBf6R2kohjDwU6bjLZzdGN61SrC0sPXafz9TubsvkRsguo/MJq9E7y+ECp1h+1fwJax2sgiE3k5ezGz2Uyi4qJ4b/t7RCdE6xCsYm6pJYK8UsqLhr+/CSyTUg5D24PYmP0IFMVsPHI7MrFjBTaPaED1onn5csN5mk/exYYz4WprRWPZ2mn9BDUHw4GZsHYYJJmeTMvlK8c3Db7h7N2zfLL3E5Jkkg7BKuaUWiJI/n9TU+AfAMP+w+q/rGIVSuXPw7x+NVnYvybO9na8s+Q4r/98gFM3Hlg6tKzBxgbafAsNP4QTi+D3fpAQa3KxTYo0YWT1kfxz7R9mnZylQ6CKOdml8tlpIcT3aJvKlAK2AAgh3DMjMEVJj4ZlvKhb0oNVx0L5YcsFOs7aR+cqPoxuVZZC7rksHZ51EwKafgK53GHz/7TO5G6LwcG0TYT6lO/DlYdXmHN6DsVci9GhZAedAlb0llpncS5gOFAQbVOZU4bjdYGSUspFmRZlMqqzWEnLo5h4Zu+8xNy9VxDA4IYlGNKoJC6OqT33KAAcXwR/vwe+NaDnSi05mCA+MZ4hW4dw4s4J5racS9UCVXUKVEkvk9YaeqGgqlLK47pFlgEqESjGCr3/lG82XeDvUzfxyuPI6JZl6VLNF1sbNfM1Vef+gt8HgFc56P0n5DZtzsbD2If02tCLqNgolrRbQuE8hdO+SNGdqTOLk5urQzyKkil88zozo0cV/ni7Lr55c/HhH6dpP2Mv+0MiLR2adSvfEXqugHuX4LfW8OC6ScW5Oboxq9ksEmUiQ7cN5VHcI50CVfSS3kSgHqWULKda0bz8+XZdZvSoQlR0PD3nHmLggqNcinhs6dCsV6lm2mqlTyO1ZBARbFJxRV2LMrXJVK5HXWfUrlEkJCXoFKiih/QmgglmiUJRzEwIQYeAQmwb2YgPW5fl4OW7tJqym/Frz3L/SZylw7NORWpB3w2QGA/zWmtLWpughncNPq3zKftv7ufrw1+rYb5WxKhEIITwMXQS3xNCNBRCNDTyutZCiAtCiBAhxMcpfN5LCHHa8NovhAhIZ/yKki5O9ra807gUO0Y15vUahVl44CqNv9/Jr3uvEJegRkX/h3cF6L8J7F1gQQe4us+k4l4t/Sr9/Pux4sIKlp5fqlOQiqmMWYb6G6Ab2haVz2abSCnlK2lcZwsEAy2AUOAI0ENKeS7ZOXWBICnlfSFEG2C8lLJWauWqzmJFT+dvRTFpfRB7LkZSzMOZMW39aFm+gFpK+UUPw2BRJ62/4PVFUKZlhotKTErk/Z3vsyt0F+PrjKdTqU7q3zsTmLpD2QWgkpQyXbNMhBB10H7YWxnejwGQUqa4PIVhX+RAKaVPauWqRKDoTUrJzuAIJq0PIuTOY2qXyMfYduWp4KN23PqXJ5HayqW3A6Hzz1Cxa4aLehr/lGHbh3H41mFaFWvFp7U/VTucmZmpo4YuA/YZuK8PcCPZ+1DDsZcZAKS4SLoQYrAQ4qgQ4mhEREQGQlGUlxNC0KRsfjYNb8DnHf0Jvv2YDjP3MmrVKW5HxVg6POvh4glv/g2Fa8EfA+HobxkuytnemTkt5jC86nC2XdtGl7VdOHLriI7BKulhTI3gDyAA2Ia2HwEAUsr30rjuNaCVlHKg4X1voKZhvaIXz20C/AjUl1LeTa1cVSNQzO1hdDw/7ghh3r6r2NoIhjQqyaCGxXF2UBPSAIiPhpVvwsXN0Hw81H/fpOLORp7loz0fcT3qOgMqDuCdgHewt83Is6eSGlNrBGuBz4H9wLFkr7SEAslnjvgCN1MIrhLa/ISOaSUBRckMbrnsGdPWj38+aEiTcl5M2RpM0+938cexUJKS1EgX7HNB9yVQoStsHa+9TBgB5O/pz8r2K3m19KvMPTOX3ht7c/XhVb2iVYxg1Mxiw3ITRaSUF4wuWAg7tM7iZmjrFR0BekopzyY7pwiwHegjpdxvTLmqRqBktiNX7/H5unOcDn1IRR83xrbzo1YJD0uHZXlJibBhlNZEVL2/tsmNja1JRW69tpXxB8YTlxjHxzU/pnOpzqojWScm1QiEEB2Ak2j7EyCEqCyEWJvWdVLKBGAosBkIAlZKKc8KIYYIIYYYTvsM8AB+FEKcFEKoX3jF6tQolo8179RjSrcAIh/H0m3OQYYsOsbVyCeWDs2ybGyh3WSo/4GWDP4cpM05MEHzos35o8MfVPKqxLj94/hg5wc8iFEryZqbMX0Ex9CWod4ppaxiOHZGSmmRDe1VjUCxpOi4RObuuczsXZeIT0zizTrFGNasNG65cnib9t6psHUclG4Jry0AB2eTikuSSSw8u5BpJ6aRzzEfkxpMonbB2joFmzOZ2keQIKV8+MIx1VCq5Ei5HGwZ1qw0O0c15tUqvvy67wqNv9vBgv1XiU/MwRPS6o+A9lPh4j/aENOYF38y0sdG2NC3Ql+Wtl2Ki4MLg7cMZvLRycQlqlng5mBMIggUQvQEbIUQpYUQM9A6jhUlx8rv6sQ3XSuxblh9ynm7Mm7tWVpN3c2GM+E5t0O5ej/o+iuEHtZmIT8xfXE/Pw8/VrRfwWtlXmPe2Xm8seENLj+4rEOwSnLGNA05A58ALdEWndsMfC6ltMgAa9U0pFgbKSXbgu7w1cYgLkU8oZx3HoY3K00rf29scuKS1xf/gRW9wb0w9F4Nbr66FLvj+g7G7R9HdEI0o2uM5rUyr6mO5HTQcz8CW8BFShmlV3DppRKBYq0SkyTrTt9k2raLXDYkhBHNS9OyfA5MCNf2w9Ju4OQGff4Cj5K6FBvxNIKx+8ay/+Z+Gvs2ZkK9CeRzyqdL2dmdqUtMLAWGoK0zdAxwAyZLKb/TO1BjqESgWLvEJMnaU2FM3xbClcgn+BV0NSSEHLaG0c2TWn+BEFrNwFuf8SVJMoklQUuYcmwKbo5ufFHvC+r51NOl7OzM1M7i8oYaQCdgA1AE6K1jfIqSrdjaCDpX8eWf9xsy+fUAouMSeGvRMdpN38uWs7dyzvLLhSprK5faOsK8dnD9oC7F2ggbepfvzbJ2y3B3dGfI1iF8c/gbYhPTtRyakowxicBeCGGPlgj+klLGo0YNKUqa7GxteLWqL1s/aMQPrwXwJC6BwYuO0X7GXv45dztnJATP0loyyO0FCztByFbdii6bryzL2i2jZ7meLA5aTI/1Pbh4/6Ju5eckxiSCn4GrgAuwWwhRFLBYH4GiZDV2tjZ0qebLtg8a8V3XSjyKSWDQwqO8MnMf24JyQEJwLwz9NoFnKVjaHc6u0a1oJzsnxtQaw6xms7gbfZfu67qzJGhJ9v831Vm6OoufXySEnWHmcKZTfQRKVhefmMTqE2HM2H6RG/eiqeTrxojmpWlSNn/27kOIfqB1IIcehg7ToGofXYuPjI7ks32fsSdsD/V96vN5vc/xzOWp6z2yMlM7ix2BLkAx4Pnyi1LKiTrGaDSVCJTsIj4xiT+PhzJjewih96MJ8HVjRPMyNC7rlX0TQtxTWNlbayJq+QXU/c9ixCaRUrL8wnJ+OPoDLvYufF7vcxr6GrWhYrZnaiLYBDxEGzH0bIcypJQ/6BmksVQiULKb+MQk/jimJYSwB9FULuzOiOalaVQmmyaEhDhtXaJza6BCF2g+QWs+0lHI/RA+2vMRwfeD6V62OyOrj8TJzknXe2Q1piaCQCllBbNElgEqESjZVVxCEn8cD2WmISFUKeLOiOZlaFjaM/slhKRE2PUN7Jumva87DOqNAMfcut0iNjGWacensejcIkq6leSbht9QNl9Z3crPakxNBHOAGVLKM+YILr1UIlCyu7iEJFYdu8Gs7SHcfBhDVUNCaJAdE8KDG9p+BoG/Q25vaPYZBPQAG2PGsRhnf9h+Ptn3CQ9jHzKi6gjeKP8GNkK/8rMKUxPBOaAUcAVthzKBtnl9Jb0DNYZKBEpOEZuQyKqjoczaEUL4wxiqFc3L+83LUK+UR/ZLCDcOw6YxEHYUCgZAq6+gmH6TxO7F3GPc/nHsvLGTOgXr8EX9L8jvnF+38rMCUxNB0ZSOSymv6RBbuqlEoOQ0sQmJrDwayo+GhFCjWF5GNC9D3ZLZLCEkJUHgH1oNISoU/F6BFhMhX3FdipdSsip4Fd8d+Q4nOycm1J1A0yJNdSk7KzB5rSEhRH2gtJRynhDCC8gtpbyic5xGUYlAyaliExJZceQGs3aEcDsqlprF8jGiRWnqlsxmQyTjnsKBmbB3CiQlQO23ocEocHLVpfjLDy/z8e6PCboXxGtlXmNU9VE425u2f0JWYGqNYBxQHSgrpSwjhCgErJJSWmRxD5UIlJwuJl5LCD/u1BJCreL5GNG8DHVKZrPtM6NuwrbP4dRScPaEpmO1uQcmbocJEJ8Yz4yTM5gfOJ+irkX5puE3lPcor0PQ1svURHASqAIcT7ZD2WnVR6AolhUTn8iyw9eZvfMSdx7FUruElhBqZ7f9lMOOw+b/wfUDkN8fWn8JJRrrUvSh8EP8b+//uBdzj2FVhtHXv2+27Ug2NREcllLWFEIcl1JWFUK4AAdUIlAU6xATn8jSQ9eZvesSEY9iqVPCg/dblKFm8Wy0PLOUcO4v+OdTeHAdyrTRJqR5ljK56AcxD5hwYAJbr2+llnctvqj/Bd4u3joEbV1MXX10pRDiZ8BdCDEI2AbMNfLGrYUQF4QQIUKIj1P4vJwQ4oAQIlYIMcqYMhVF+Tcne1v61y/Ong+b8Gn78ly885jXfz5Ar7kHOXL1nqXD04cQ4N8J3j0CzcfD1b3wYy1tpFH0fZOKdndyZ3LjyUysO5HTkafpsrYLO2/s1CPqLMPYzuIWaDuUAWyWUqa5hKBhE5tgoAUQChwBekgpzyU7Jz9QFG1l0/tSyu/TKlfVCBQlddFxiSw5dI2fdl0i8nEc9Ut50qdOURqXzY+DXTZp9nh8B7Z/AScWaZvfNP6ftlWmrb1JxV6LusboXaM5f+88Q6sMZVDFQdlmZFaGmoaEEI/4/+WmX/yXiAEuAZ9IKbe95Po6wHgpZSvD+zEAUsqvUjh3PPBYJQJF0U90XCKLD17j592XiXwci7uzPW0rFqRzFR+qFcmbPXZNuxWo9R9c2QWeZaDVl1C6hUlFxiTEMP7AeNZfXk+Loi34ot4X2WJUkW5bVSYr0BaoACx52fITQoiuQGsp5UDD+95ALSnl0BTOHU8qiUAIMRgYDFCkSJFq165ZZAqDomRJ8YlJ7A2JZM2JMDafvUVMfBK+eXPRsXIhOlfxoVT+PJYO0TRSwoWNsGUs3LsEJZtBq0mQ38+EIiULzi5gyvEplHIvxbQm0/DNo8/ey5aieyJIVvBbUsqfX/LZa0CrFxJBTSnlf5YbVDUCRckcj2MT2HL2FmtO3mTvxQiSJFTwcaVTZR9eCShEftcsvDBbQhwc+UVbwyj2MVTrC03+By4Zn2exL2wfo3ePxlbY8n2j76lVsJZ+8WYysyWCNG6qmoYUxYrdeRTD36fCWXMijDNhD7ERUK+UJ50q+9Cqgje5He3SLsQaPbkLu76GI7+CQ25o9CHUHAx2Dhkq7lrUNd7b/p7Wf1BjND3L9cyS/QaWSgR2aJ3FzYAwtM7inlLKsymcOx6VCBTFYkLuPOavk2GsPhFG6P1onOxtaFHem85VCtGgtBf2tlmwk/nOea25KOQfyFcCWnwO5dppI5DS6XHcY8bsHcPOGzvpXKozY2uPxcE2Y4nFUiySCAw3bgtMBWyB36SUk4QQQwCklD8JIbyBo4ArkAQ8BspLKV+6FaZKBIpiPlJKjl+/z+oTYaw7Hc6Dp/Hkc3GgQ6WCdKziQ5XC7lnvafjiVq1Dkj3zMwAAD7VJREFUOfICFGugdSgXTP80qCSZxI8nf+Tn0z9TyasSUxtPxcvZywwBm4fFEoE5qESgKJkjLiGJXcERrDkZxtZzt4lNSKKohzOdKvvQqYoPxT1dLB2i8RIT4Ng82PGlNu+gam9oMhbyFEh3UVuubmHsvrHksc/DlCZTqORlkbm16aYSgaIoJomKiWdT4C3WnAjjwOW7SAkBhd3pXLkQ7QMK4Znb0dIhGif6Puz+Hg79DHaO0OADqP0u2Kevk/zCvQsM3zGcO0/vMK7OODqW6mimgPWjEoGiKLq59TCGtafCWH3iJkHhUdjaCBqW9qRTFR9alC+As0MW6GS+ewm2fAoX1oN7EW27TP//a+/Og6us7z2Ov7/ZyAIhYU9CWEQ2EQkWQVBQUFtaKakz9l7sqHOd6fXq1A31OrU617bOWG3BhVumjlctdkQdqtLgFS/KQGVfFBJ2EVEhCyRICIFsJznf+8fvicSYsJ6T58Dzfc0w5CzPOd/DhPN5nt/veb6/m85o/qCyrpKHP36YDQc2cOvwW3lozEMkxMXuZ7cgMMZExWcHqvlHYQkFm0soraojNSmeqSP6kD86h6sGdSch1ieZ937s5g8OboPcK11Du5wfnPbmoXCI2Z/MZv7O+YzLGsesSbPISM6IYsFnz4LAGBNV4bCy4avDFBS6SebqukZ6dO7E9FHZ/Gx0NiNzusbuJHO4CTa/DsuehOMVcNkMmPKYO1I4TQs/X8iT656kV2ov5kyZw5DMIVEs+OxYEBhjOkxdqIl/flbOws0lLN9VQUNTmIt6pnGTN8mc2y1G2zXUHYVVz8LauaBht3byxAfdqaenoaiiiJnLZ3IsdIynrn6K6/tfH+WCz4wFgTHGF1U1IRZvK2Ph5hI2fOk6oY7pn0n+6BymjcwiMy0Gz8WvKobVL8Cnr7kV0kb+HCY+BD1PvZdfXlPOzOUz2XJoC3eNuou7R90dM+sbWBAYY3xXXFnDoqJSFm4q4fPyYyTGC9cM6cn0vByuH94r9iaZqw/Amv+GT16FUK1rgz3pP6H3iJNuVt9Uz+/X/p5FXyxicu5k/jDxD6Ql+n+qrQWBMSZmqCo7yo5SUFjKosJSDhx1k8w3XNKb/LwYvJL5+CE3XLThf6ChGoZNg0kPQ/bodjdRVebvnM+sT2YxIH0Ac6bMoV/66c85RIMFgTEmJp2YZC7lg23uSuZMr112fl4OY/rHULvsmsPu+oP1f4G6Krj4BtfHKHdsu5usK1vHwx8/TFjDzJo0iwk5Ezqw4O+yIDDGxLyGxjArP6+goLCUj3YcpDbURHbXZH46KpvpedlckpUeG2ce1R11XU7XzoWab2DgNS4QBlzd5tP3V+/nvmX3sbdqLw/+4EFuv+R2Xz6HBYEx5rxyvL6RpTsPUlBYyordFTSGlYt7dSbfC4X+3f0fc6fhuJs/WD0HjpdDvwluyGjQlO9dmFYTquHx1Y/z0dcfMe2iaTwx/gmSEzq25bcFgTHmvFV5vIHF28ooKCz99syjvNwM8vOyufGyLHp18XkNhVAtbPobrHoeqkshZ4ybVB7yo+8EQljDvLTlJeYWzmVE9xE8P/l5+qT16bAyLQiMMReE0iO1vFdUSkFhKTvKjhInMGFQD6bnZTP10j6kJ5/bmsXnpLEeCt9w1yIc2Qd9LnOBMGwaxJ2Y/F62bxmPrnyUlIQUnpv8HKN7tT/pHEkWBMaYC86e8moWFZZSUFTK19/UkJQQx5ShvcjPy2bysF4kJ8b7U1hTCLYsgJWz3dKZPYe7IaMRN0Gcq2lP5R7uW34fZcfLeGzcY9w85Oaol2VBYIy5YKkqRcVVFBSW8F5RGYeO1dOlUwI/urQP00dlM8GvnkfhJtj2LqycBRW7oPvF7sK0kT+H+ESq6qt4ZMUjrCldw4yhM3hk7CMkxkXviMaCwBgTCE1hZd3ebygoLOGDbQe8nkdJTLvMTTL7srBOOAy73oMVf4IDWyGjv2tdMeoXNMbF8cKmF5i3fR5jeo9h9rWz6ZbcLSplWBAYYwLH9TyqYFFRCUt3ltPQGCa3Wwr5o3LIz8tmcO8uHVuQKuxeAiv+CCWfQnoOXPUAXH4b7+1bym/X/JbuKd2ZM2UOw7oNi/jbWxAYYwKtui7Eh9sPUlBUyqrPKwgrDOvThfy8HH46Kou+mR3YCE8VvljmjhD2rYXOvWHCvWwfOJ77V/2aqvoqnrzqSaYOnBrRt7UgMMYYT0V1PYu3llFQWMKmfUcAuGJAJtPzcrhxZBbdOqoRnip8tcodIXy5AlK7c+iKO5h5fCeF32zllyN/yT159xAfF5lJbwsCY4xpw/7DrhFeQWEJuw8eIyFOuHpwD4ZnpZORkkhGaiIZqUlkpCSSmZbk3ZdEUkKEJ5/3rXdHCHs+oiE5g6cGX847x/YwMWciz0x6hi5J5z6M5VsQiMhU4AUgHnhZVZ9u9bh4j/8EqAH+TVU3new1LQiMMdGw64BrhPfB1jJKjtQSamr/uzE1KZ7M1CS6piSSmZZIRkqSFxqJJ+5PTToRJKmJZKQknvrspZJNsGIW+tn7LMjsztMZnenbpS9zrpvLwK4Dz+nz+RIEIhIP7AZuAIqBjcAtqrqjxXN+AtyLC4JxwAuqOu5kr2tBYIyJNlWlpqGJypoGjtSE3J/aBiprQlTVuL/d/Q0cqQ1RWdNAVU2II7UhmsLtf6d26ZRAxneCI4lMLyS6Nv+cmkhW3Rfkbp3LzpKlPNSrB6GEJJ4Z/zsmDZ5+1p/pZEEQzQbgY4E9qrrXK+ItIB/Y0eI5+cDf1KXROhHJEJEsVS2LYl3GGHNSIkJapwTSOiXQN/P0t1NVqusbOXL8RHAcaRUmzQFSWROiuLLWhUhtiO/vk/+CQXINd9S/w/tZX3HP6t/wL+sW8Phtr0fyowLRDYIcYH+L28W4vf5TPScH+E4QiMidwJ0A/fr529PbGGPaIyKkJyeSnpxIP07/TKRwWDla58Ki0jvKcAFyCZU11zL+8C66Vv2RrJ6DolJ3NIOgras2Wmfe6TwHVX0JeAnc0NC5l2aMMbEjLk68uYQkBtBWZ9UhwNkPC53y/aP2ym7vPrfF7b5A6Vk8xxhjTBRFMwg2AoNFZKCIJAEzgEWtnrMIuF2cK4Eqmx8wxpiOFbWhIVVtFJF7gCW400dfVdXtInKX9/iLwGLcGUN7cKeP3hGteowxxrQtmnMEqOpi3Jd9y/tebPGzAr+KZg3GGGNOzoferMYYY2KJBYExxgScBYExxgScBYExxgTcedd9VEQqgK/PcvMewKEIlhMpsVoXxG5tVteZsbrOzIVYV39V7dnWA+ddEJwLEfmkvaZLforVuiB2a7O6zozVdWaCVpcNDRljTMBZEBhjTMAFLQhe8ruAdsRqXRC7tVldZ8bqOjOBqitQcwTGGGO+L2hHBMYYY1qxIDDGmIALTBCIyFQR+UxE9ojIr/2uB0BEXhWRchHZ5nctLYlIrogsF5GdIrJdRO73uyYAEUkWkQ0iUuTV9Tu/a2pJROJFZLOI/K/ftTQTka9EZKuIFIpIzCz27S1L+7aI7PJ+z8bHQE1DvX+n5j9HReQBv+sCEJGZ3u/8NhF5U0SSI/r6QZgjEJF4YDdwA24xnI3ALaq646QbRr+uScAx3LrNl/pZS0sikgVkqeomEekCfAr8LAb+vQRIU9VjIpIIrALuV9V1ftbVTEQeBMYA6ao6ze96wAUBMEZVY+riKBF5DVipqi9765WkquoRv+tq5n1nlADjVPVsL2CNVC05uN/1S1S1VkQWAItVdV6k3iMoRwRjgT2quldVG4C3gHyfa0JVVwCH/a6jNVUtU9VN3s/VwE7cWtK+UueYdzPR+xMTezIi0he4EXjZ71pinYikA5OAVwBUtSGWQsBzHfCF3yHQQgKQIiIJQCoRXskxKEGQA+xvcbuYGPhiOx+IyABgNLDe30ocb/ilECgHPlLVmKgLeB54BAj7XUgrCnwoIp+KyJ1+F+O5CKgA/uoNpb0sIm0t1OunGcCbfhcBoKolwCxgH1CGW8nxw0i+R1CCQNq4Lyb2JGOZiHQG3gEeUNWjftcDoKpNqpqHW996rIj4PqQmItOAclX91O9a2nCVql4O/Bj4lTcc6bcE4HLgL6o6GjgOxMS8HYA3VDUd+LvftQCISCZuBGMgkA2kicitkXyPoARBMZDb4nZfInxodaHxxuDfAear6rt+19OaN5TwT2Cqz6UAXAVM98bj3wKmiMjr/pbkqGqp93c5sBA3TOq3YqC4xdHc27hgiBU/Bjap6kG/C/FcD3ypqhWqGgLeBSZE8g2CEgQbgcEiMtBL+xnAIp9rilnepOwrwE5VfdbvepqJSE8RyfB+TsH9B9nlb1Wgqo+qal9VHYD73VqmqhHdYzsbIpLmTfbjDb38EPD9DDVVPQDsF5Gh3l3XAb6eiNDKLcTIsJBnH3CliKR6/zevw83bRUxU1yyOFaraKCL3AEuAeOBVVd3uc1mIyJvAtUAPESkGnlDVV/ytCnB7uLcBW73xeIDfeGtQ+ykLeM07oyMOWKCqMXOqZgzqDSx03x0kAG+o6v/5W9K37gXmeztme4E7fK4HABFJxZ1d+B9+19JMVdeLyNvAJqAR2EyEW00E4vRRY4wx7QvK0JAxxph2WBAYY0zAWRAYY0zAWRAYY0zAWRAYY0zAWRCYQBOR51p2mBSRJSLycovbs0Xkv860Y62IzBORmyNZqzHRYkFggm4N3lWaIhIH9ABGtHh8ArBEVZ/2oTZjOoQFgQm61Zy4XH8E7srbahHJFJFOwHBglIj8Gb7d058jImtEZG/zXr84fxaRHSLyPtCr+Q1E5DqvudpWbw2KTiIyVkTe9R7PF5FaEUny1lzY24Gf3xgLAhNsXi+eRhHphwuEtbhOq+NxawtsARpabZYFXA1MA5qPFG4ChgIjgX/nxFFGMjAP+FdVHYm7wvdu3FWio71tJ+IC6ApgHDHS6dUEhwWBMSeOCpqDYG2L22vaeP4/VDXsLdTT27tvEvCm1x21FFjm3T8U1zBst3f7NWCSqjYCe0RkOK4R3LPea0wEVkb6AxpzMhYExpyYJxiJ2zNfhzsimIALidbqW/zcssV5W/1a2mqB3mwlrtNlCFiKO8q4GlhxuoUbEwkWBMa4L/tpwGFvj/4wkIELg7Wn+RorgBnewjlZwGTv/l3AABG52Lt9G/Bxi20eANaqagXQHRgG+N4Q0QRLILqPGnMKW3FnC73R6r7OqnrI6955KguBKd52u/G+7FW1TkTuAP7uLTO4EXjR22Y9bmip+QhgC26BG+sEaTqUdR81xpiAs6EhY4wJOAsCY4wJOAsCY4wJOAsCY4wJOAsCY4wJOAsCY4wJOAsCY4wJuP8HAA6uzp224rYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "labels = ['1D', '2D', '3D']\n", "\n", "dres_fig2, dres_ax2 = plt.subplots()\n", "for data, label in zip(dres_conv2.T, labels):\n", " plt.plot(data, label=label)\n", "dres_ax2.set_xlabel('Window')\n", "dres_ax2.set_ylabel('Jensen-Shannon divergence')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "[1] Oliver Beckstein, Elizabeth J. Denning, Juan R. Perilla, and Thomas B. Woolf.\n", "Zipping and Unzipping of Adenylate Kinase: Atomistic Insights into the Ensemble of OpenClosed Transitions.\n", "Journal of Molecular Biology, 394(1):160–176, November 2009.\n", "00107.\n", "URL: https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164, doi:10.1016/j.jmb.2009.09.009.\n", "\n", "[2] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney, and Oliver Beckstein.\n", "MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations.\n", "Proceedings of the 15th Python in Science Conference, pages 98–105, 2016.\n", "00152.\n", "URL: https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html, doi:10.25080/Majora-629e541a-00e.\n", "\n", "[3] Naveen Michaud-Agrawal, Elizabeth J. Denning, Thomas B. Woolf, and Oliver Beckstein.\n", "MDAnalysis: A toolkit for the analysis of molecular dynamics simulations.\n", "Journal of Computational Chemistry, 32(10):2319–2327, July 2011.\n", "00778.\n", "URL: http://doi.wiley.com/10.1002/jcc.21787, doi:10.1002/jcc.21787.\n", "\n", "[4] Matteo Tiberti, Elena Papaleo, Tone Bengtsen, Wouter Boomsma, and Kresten Lindorff-Larsen.\n", "ENCORE: Software for Quantitative Ensemble Comparison.\n", "PLOS Computational Biology, 11(10):e1004415, October 2015.\n", "00031.\n", "URL: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004415, doi:10.1371/journal.pcbi.1004415." ] } ], "metadata": { "kernelspec": { "display_name": "Python (mda-user-guide)", "language": "python", "name": "mda-user-guide" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }