{ "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:** May 18, 2021 with MDAnalysis 1.1.1\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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:55.164483Z", "iopub.status.busy": "2021-05-19T05:59:55.163573Z", "iopub.status.idle": "2021-05-19T05:59:56.540827Z", "shell.execute_reply": "2021-05-19T05:59:56.541182Z" } }, "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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:56.546291Z", "iopub.status.busy": "2021-05-19T05:59:56.545647Z", "iopub.status.idle": "2021-05-19T05:59:56.785608Z", "shell.execute_reply": "2021-05-19T05:59:56.786089Z" } }, "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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:56.791388Z", "iopub.status.busy": "2021-05-19T05:59:56.790617Z", "iopub.status.idle": "2021-05-19T05:59:58.964292Z", "shell.execute_reply": "2021-05-19T05:59:58.964732Z" } }, "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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:58.970797Z", "iopub.status.busy": "2021-05-19T05:59:58.969750Z", "iopub.status.idle": "2021-05-19T05:59:58.974402Z", "shell.execute_reply": "2021-05-19T05:59:58.975378Z" } }, "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": 1, "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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:58.991618Z", "iopub.status.busy": "2021-05-19T05:59:58.990939Z", "iopub.status.idle": "2021-05-19T05:59:59.088137Z", "shell.execute_reply": "2021-05-19T05:59:59.088508Z" }, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Jensen-Shannon divergence')" ] }, "execution_count": 1, "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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:59.094046Z", "iopub.status.busy": "2021-05-19T05:59:59.093346Z", "iopub.status.idle": "2021-05-19T05:59:59.095367Z", "shell.execute_reply": "2021-05-19T05:59:59.094973Z" } }, "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": { "execution": { "iopub.execute_input": "2021-05-19T05:59:59.099109Z", "iopub.status.busy": "2021-05-19T05:59:59.098472Z", "iopub.status.idle": "2021-05-19T06:00:01.094005Z", "shell.execute_reply": "2021-05-19T06:00:01.094590Z" } }, "outputs": [ { "data": { "text/plain": [ "(9, 3)" ] }, "execution_count": 1, "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": { "execution": { "iopub.execute_input": "2021-05-19T06:00:01.112454Z", "iopub.status.busy": "2021-05-19T06:00:01.108370Z", "iopub.status.idle": "2021-05-19T06:00:01.241556Z", "shell.execute_reply": "2021-05-19T06:00:01.242120Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1RU19rH8e+mgwVEsKKCCnbAir131DTNtaQZu0k05ibRtJvclJvEmNcWS+xpxhSj0VhiL7EiRuwNbNgbVen7/eOgQQUcmBmKPJ+1ZunMnLNnY1bmxzl772crrTVCCCGKLpv87oAQQoj8JUEghBBFnASBEEIUcRIEQghRxEkQCCFEEWeX3x3IKQ8PD+3t7Z3f3RBCiEIlNDT0mtbaM7P3Cl0QeHt7s2fPnvzuhhBCFCpKqTNZvSe3hoQQooiTIBBCiCLOqkGglOqqlDqmlDqplBqXyfttlVLRSql96Y//WLM/QgghHmS1MQKllC0wDegERAIhSqllWuvD9x26VWvdw1r9EEIUDMnJyURGRpKQkJDfXXmkOTk54eXlhb29vcnnWHOwuAlwUmsdAaCUWgQ8BtwfBEKIIiAyMpISJUrg7e2NUiq/u/NI0lpz/fp1IiMj8fHxMfk8a94aqgicy/A8Mv21+zVTSoUppVYppepk1pBSaqhSao9Sas/Vq1et0VchhJUlJCRQunRpCQErUkpRunTpHF91WTMIMvuvfX+p071AFa11ADAVWJpZQ1rrWVrrRlrrRp6emU6DFUIUAhIC1pebf2NrBkEkUCnDcy/gQsYDtNYxWuu49L+vBOyVUh7W6MyN+CT+u/wQCcmp1mheCCEKLWsGQQjgq5TyUUo5AH2BZRkPUEqVU+nxpZRqkt6f69bozLaT11iw/TR9Z+3kWlyiNT5CCFHAvfjii5QpU4a6deve8/obb7xBzZo18ff354knniAqKsrkNl944QV+/fXXHPfl9OnTLFy4MMfnWYPVgkBrnQK8DPwJHAF+1lofUkoNV0oNTz+sN3BQKRUGTAH6aivtlNMzoAIzBjTk6KUYnpi+jZNXYq3xMUKIAuyFF15g9erVD7zeqVMnDh48yP79+/Hz8+PTTz+1el9yEwSpqda5o2HVdQRa65Vaaz+tdTWt9Sfpr83UWs9M//tXWus6WusArXVTrfV2a/ana91yLBrajNtJqTw5fTvbw69Z8+OEEAVM69atcXd3f+D1zp07Y2dnTKJs2rQpkZGRmZ4/fvx46tWrR0BAAOPGPbA0Cm9vb65dM75X9uzZQ9u2bQHYvHkzgYGBBAYGUr9+fWJjYxk3bhxbt24lMDCQiRMnkpqayhtvvEHjxo3x9/fn66+/BmDTpk20a9eO/v37U69ePeLj4wkODiYgIIC6devy008/mf3vUuhqDZkrsJIbS0a24MUFITw3dzefPeVP74Ze+d0tIYqU/y4/xOELMRZts3aFkrzfM9OJhzkyb948/vWvfz3w+qpVq1i6dCm7du3CxcWFGzdumNzmhAkTmDZtGi1atCAuLg4nJyc+++wzJkyYwB9//AHArFmzcHV1JSQkhMTERFq0aEHnzp0B2L17NwcPHsTHx4fFixdToUIFVqxYAUB0dLTZP3PRKjGRmgxAJXcXfh3RnKCq7rz+Sxj/t+YYsnezEOKTTz7Bzs6OAQMGPPDeunXrGDhwIC4uLgCZXllkpUWLFrz22mtMmTKFqKiou1cfGa1Zs4Zvv/2WwMBAgoKCuH79OidOnACgSZMmd9cF1KtXj3Xr1jF27Fi2bt2Kq6trbn7UexSdK4IT6+CPMTBwBbhVxtXZnvkvNOHdpQeYsuEkZ27cYnxvfxztbPO7p0I88izxm7ulffPNN/zxxx+sX78+0ymYWuuHTs20s7MjLS0N4J65/OPGjSM4OJiVK1fStGlT1q1bl2n7U6dOpUuXLve8vmnTJooVK3b3uZ+fH6GhoaxcuZK33nqLzp0785//mFedp+hcEbj7wO2b8OuLd68MHOxs+Pwpf97oUoPf913g2Tm7uRmflM8dFULktdWrV/P555+zbNmyu7/x369z587MmzePW7duAWR6a8jb25vQ0FAAFi9efPf18PBw6tWrx9ixY2nUqBFHjx6lRIkSxMb+M2mlS5cuzJgxg+Rk4/vp+PHjxMfHP/AZFy5cwMXFhWeeeYbXX3+dvXv35v4HT1d0gqB0Neg1GSJDYP2Hd19WSvFSu+pM6VeffZFRPDljO6evPfiPL4Qo/Pr160ezZs04duwYXl5ezJ07F4CXX36Z2NhYOnXqRGBgIMOHD3/g3K5du9KrVy8aNWpEYGAgEyZMeOCY999/n9GjR9OqVStsbf+5uzBp0iTq1q1LQEAAzs7OdOvWDX9/f+zs7AgICGDixIkMHjyY2rVr06BBA+rWrcuwYcNISUl54DMOHDhAkyZNCAwM5JNPPuHdd981+99FFbZ7440aNdJmbUzzxxjYMw/6/wx+916C7Tl9gyHfGm3Peq4Rjb1NvwcohMjekSNHqFWrVn53o0jI7N9aKRWqtW6U2fFF54rgji6fQtl6sGQ4RJ+/561G3u4sGdkCNxcHBszexbKwC1k0IoQQj46iFwT2TtBnAaQmpY8X3Hvp5e1RjN9GNCewkhujfvybrzackBlFQohHWtELAgCP6tBjEpzbCRs/fuDtUsUc+G5wEx4PrMCENcd589f9JKWk5UNHhRDC+orO9NH7+feB01vgr4lQpSX4drznbUc7Wyb+K5DKpYsxZf0JzkfdZsYzDXF1Nn2zByGEKAyK5hXBHV0/hzK1YclQiHlwPEApxWud/PiyTwAhp2/w1IztnLtxKx86KoQQ1lO0g8DBBfp8A8kJsHjwA+MFdzzV0ItvXwziSkwCT0zfxr5zplcmFEKIgq5oBwGApx/0+D84sw02f5blYc2qlea3kS1wdrCl76wdrD54MQ87KYSwhKioKHr37k3NmjWpVasWO3bsMPnc4sWL5+ozly5dyuHDBXuHXgkCgIC+EPgMbJkA4RuyPKx6meIsGdmCWuVLMuKHvczeEiEzioQoREaPHk3Xrl05evQoYWFhebKuITdBkNlCMmuSILij+3jwrAGLh0DspSwP8yjuyI9DmtK9bnk+WXmEd5ceJCVVZhQJUdDFxMSwZcsWBg0aBICDgwNubm4PHHf58mWeeOIJAgICCAgIYPv2e6vjb9q0iR49etx9/vLLL7NgwQLAqClUu3Zt/P39ef3119m+fTvLli3jjTfeIDAwkPDwcMLDw+natSsNGzakVatWHD16FDD2Snjttddo164dY8eOzbR0tbUU3VlD93MoZqwvmNXOGC947newybwAnZO9LVP71aeSuwszN4dzPuo2X/VvQHFH+ecUwiSrxsGlA5Zts1w96Jb17d2IiAg8PT0ZOHAgYWFhNGzYkMmTJ99T0A1g1KhRtGnThiVLlpCamkpcXJxJH3/jxg2WLFnC0aNHUUoRFRWFm5sbvXr1okePHvTu3RuADh06MHPmTHx9fdm1axcjR45kwwbjTsTx48dZt24dtra29OzZ84HS1dYiVwQZlakFwRPg9FbYPD7bQ21sFOO61eTTJ+ux9cQ1es/YzsXo23nUUSFETqWkpLB3715GjBjB33//TbFixfjssweDY8OGDYwYMQIAW1tbk8s8lyxZEicnJwYPHsxvv/2WafG6uLg4tm/fTp8+fQgMDGTYsGFcvPjPeGOfPn3u1igypXS1pcivsPcLHACn/4LNn0OV5lC1TbaH92tSmYpuzoz8YS+PT9vG3OcbU7ei+fXBhXikZfObu7V4eXnh5eVFUFAQAL179840CB4mY6lp+KfctJ2dHbt372b9+vUsWrSIr7766u5v+nekpaXh5ubGvn37Mm0749VJZqWra9asmeP+mkKuCO6nFHSfAB6+8NsQiLvy0FNa+3ny64hm2CrF01/vYP2Ry3nQUSFETpQrV45KlSpx7NgxANavX0/t2rUfOK5Dhw7MmDEDMPYIjom5dye1KlWqcPjwYRITE4mOjmb9+vWA8dt+dHQ03bt3Z9KkSXe/7DOWmy5ZsiQ+Pj788ssvgLEHQVhYWKb9zax0tbVIEGTGsbgxXpAQbYRB2sM3jK5ZriRLX2pBNc/iDPl2D99sP231bgohcmbq1KkMGDAAf39/9u3bx9tvv/3AMZMnT2bjxo3Uq1ePhg0bcujQoXver1SpEk8//TT+/v4MGDCA+vXrAxAbG0uPHj3w9/enTZs2TJw4EYC+ffvyxRdfUL9+fcLDw/nhhx+YO3cuAQEB1KlTh99//z3TvmZWutpail4Z6pwI/QaWj4J270CbN0065VZSCqN+3Me6I5d5sYUP7wTXwtYm+12NhCgKpAx13pEy1JbU4Dmo1wc2fWqMG5jAxcGOr59tyMAW3szbdoph34VyKylv5wQLIUROSBBkRynoMRHcq8KvgyDuqkmn2doo3u9Zhw961mbD0cv86+udXIlJePiJQgiRDyQIHsaxhDFecPumUZwuzfTFYy+08GH2c40IvxrHE9O3c+yS9RaECCFEbkkQmOLOQpXwDbBtYo5O7VCrLD8Pa0Zyahq9Z2xny3HTriqEECKvSBCYquFAqPMkbPgEzmx/+PEZ1K3oytKXWlCxlDMDF4Tw4+6zVuqkEELknASBqZSCnpOhVBVjvCD+eo5Or+DmzC/Dm9Gyugdv/XaAz1YdJS2tcM3YEkI8miQIcsKppDFecOsaLBmWo/ECgBJO9sx9vhEDgiozc3M4r/z4NwnJD1+jIIQwX0JCAk2aNLk7f//999/P0flShlr8o3wAdPkfnFwL26fk+HQ7Wxs+frwu73SvxcqDF+k/eyfX4xKt0FEhREaOjo5s2LCBsLAw9u3bx+rVq9m5c6fVP/eRKEOtlPJTSq1XSh1Mf+6vlHrX+l0rwBoPhtqPwfoP4eyuHJ+ulGJI66pM79+AQxdieGL6ds5ely0whbAmpdTd3+qTk5NJTk5GqQcXe0oZ6szNBt4AvgbQWu9XSi0EPrZarwo6paDXVLgYBr++CMO3got7jpvpVq885VydGLgghH6zd7JoaFMquT9YsVCIR83nuz/n6A3L1s6p6V6TsU3GZntMamoqDRs25OTJk7z00kt3C9BlJGWoM+eitd5932smXbcopboqpY4ppU4qpcZlc1xjpVSqUqq3Ke0WCE6u0Hs+xF2GpSMgl6U66lcuxfeDgohLTKHvrJ1E3pQrAyGsxdbWln379hEZGcnu3bs5ePDgA8dIGerMXVNKVQM0QPqX9UM37FVK2QLTgE5AJBCilFqmtT6cyXGfA3/msO/5r2ID6PwxrB4LO76C5q/kqpm6FV35flAQA+bsTL8yaEZFN2cLd1aIguNhv7lbm5ubG23btmX16tXUrVs3R+cW1TLUL2HcFqqplDoPvAqMMOG8JsBJrXWE1joJWAQ8lslxrwCLgYfXey6IgoZBzR6w7gM4F5LrZup5ufLdoCCibiXTb9ZOLkTJJjdCWNLVq1eJiooC4Pbt21l+sUoZ6kykf5F3BDyBmlrrllrr0ya0XRE4l+F5ZPprdymlKgJPADOza0gpNVQptUcptefq1QK2MlcpeOwrKFnBGC+4fTPXTQVUcuO7QUHcjE+i3+ydXIqW+kRCWMrFixdp164d/v7+NG7cmE6dOt0z6HuHlKHO7ACl/geM11pHpT8vBfxba53tzCGlVB+gi9Z6cPrzZ4EmWutXMhzzC/Cl1nqnUmoB8IfW+tfs2s3TMtQ5ERkK87qAb2fo+4MRELm09+xNnpu7G88Sjiwa2pSyJa03SCREXpEy1HnHGmWou90JAQCt9U2guwnnRQKVMjz3Ai7cd0wjYJFS6jTQG5iulHrchLYLHq+G0Om/cGwF7Mr2AuehGlQuxTcvNuZKTAL9ZknlUiGEdZkSBLZKKcc7T5RSzoBjNsffEQL4KqV8lFIOQF9gWcYDtNY+WmtvrbU38CswUmu91OTeFzRNR0KN7rDmPTgfalZTDau4882LTbgUk0Df2Tu5EithIISwDlOC4HtgvVJqkFLqRWAt8M3DTtJapwAvY8wGOgL8rLU+pJQarpQabk6nCyyl4LFpUKIc/PIC3I566CnZaeTtzoKBTbgUnUD/2bu4GisrkEXhVth2RCyMcvNvbNJWlUqpbkAHQAFrtNb5NtWzwI4RZHRuN8zvBjW6wdPfmTVeALAz4joD54fgVcqZH4c2xaO4KRdkQhQsp06dokSJEpQuXTrTFb3CfFprrl+/TmxsLD4+Pve8l90YgexZbC3bpsDa96DbFxA01OzmdoRfZ+CC3VRxL8bCIUGUljAQhUxycjKRkZF3590L63BycsLLywt7e/t7XjcrCJRST2Is+CqDcUWgAK21LmmRXudQoQmCtDT4sS9EbIRBa6FCoNlNbj95jYELQvDxKMbCIU1xL+ZggY4KIYoCc2cNjQd6aa1dtdYltdYl8isEChUbG3hiJhTzNMYLEmIeesrDNK/uwdznG3PqWjwD5uziZnyS+f0UQhR5pgTBZa31Eav35FHk4g6950HUWVg+Ktf1iDJq6etxdx/kAXN2EXVLwkAIYR5TgmCPUuonpVQ/pdSTdx5W79mjonJTaP8uHFoCe+ZapMnWfp7MerYhJ6/E8czcXUTfSrZIu0KIosmUICgJ3AI6Az3THw+uyxZZa/EqVO8Iq9+Gi/st0mTbGmX4+tmGHL+UHga3JQyEELkjs4bySvw1mNkS7F1g2GZwLGGRZtcfuczw70OpXb4k3w0OoqST/cNPEkIUOWYNFssOZRZSzAOemgs3T8HyVy0yXgDQoVZZZgxoyOGLMTw3dzexCXJlIITIGVNuDc0G3gKSwdihDKNchMgp7xbQ7m04+CvsfejibJN1rF2Waf0bcPB8NM/PkzAQQuSMVXcoE5lo+W+o2g5WjYVLD+6OlFud65Tjq/4N2B8ZzQvzQ4hLlP9EQgjTmBIEudqhTGTBxgaenGVsdfnLC5Bo2n6opuhatxxT+9Vn37koBs7fTbyEgRDCBNbcoUxkpXgZeGoO3AiHFa9ZbLwAoFu98kzpW5+9Z6MYuCCEW0kSBkKI7D10z2KtdQTQUSlVDLDRWsdav1uWd/zmcVafWp3f3biXf1c4uxqXP1+mb/vPKe5Q3CLNBvuXJ01rRi/6mxcXhDDvhca4OFhv42shROH20G8HpdRr9z0HiAZCtdaZ78BcAJ2OPs38g/PzuxsPcnMl5fIWVixszVclG1CxTB1wrwru1Yw/HVxy1WzPgAqkac2Yn/YxaMEe5r3QGGcHWwt3XgjxKDCl6NxCjJ3Elqe/FIyx6UxN4Bet9Xir9vA+hXYdQVbirrJjzWv8OyYMe53G5EuXCUzMUDaiRPn0UPCB0tUyhIQPOBR7aPNL/o7ktZ/DaF6tNHOfb4yTvYSBEEWRudVH/wSe0lrHpT8vjrGb2BMYVwW1LdzfbD1yQZAuIjqCl9e/zOX4y3xUezDdHcvAjQi4HmH8eSMc4q/ee1KJ8unBkP64GxRV7wmJxaGRvP5rGC2rG3WKJAyEKHrMDYIjQIDWOin9uSOwT2tdSyn1t9a6vsV7nI1HNQgAohKieHXTq4ReDmVEwAhGBIy4dwOPhJj0UEgPhhun4Hq48Tz+yr2NFS+XHgw+4F6N7Tdd+XhnIhWq1uGr51tKGAhRxJgbBO9h/Pb/e/pLPTH2Hv4SmKW1HmDBvj7UoxwEAMmpyXy480OWnlxKN+9ufNjiQ5zsnB5+YkKMsWr5TjDceVwPfyAkbtq641qxJjalM4xFlK4GpXzA0TID1kKIgiXXQaCMX0e9MDalaYmxKc1fWut8+yZ+1IMAjO3m5h+az6TQSdTzqMfk9pPxcPbIfYOJsXeDYX/YXo4eDiOw2HV87a+i4i7fe2zxcunBUBXK1IZGL4K9s3k/kBAi35l7RRCqtW5olZ7lQlEIgjvWn1nPW3+9hZujG1PbT6WGew2LtLtw11neXnKAjrXKML23Hw4xZ9JvNd03JhF3GWr1hD7fGgvhhBCFlrlBMA1YoLUOsUbncqooBQHA4euHeWX9K8Qlx/FFmy9o7dXaIu1+t/MM7y09SKf0OkUOdpl80e+YBn++Dc1fgc4fW+RzhRD5w9ytKtsBO5VS4Uqp/UqpA0opyxTVFw9Vu3RtFgYvpErJKryy4RW+O/wdligd/mzTKnz4WB3WHr7MKz/uJTk17cGDmo6ExoNh+1QIscymOkKIgseU5abdrN4Lka2yxcqyoOsC3vnrHcaHjOdU9CneCnoLexvz9h54rpk3aWmaD5YfZtSPfzOlX33sbTP8bqAUdP3c2Gpz5RvgVgV8O5r50wghCpqHXhForc8AlYD26X+/Zcp5wrJc7F34su2XDK43mF+O/8KIdSOITow2u90XWvjwXo/arDp4iVcX7SPl/isDWztj3+UytY0ieRasmCqEKBhM2ZjmfWAsxp4EAPbA99bslMicjbJhdIPRfNziY0Ivh/LMymc4G3PW7HYHtfTh3eBarDhwkVd/yiQMHEtA/5+MqaULn4YYKT4rxKPElN/snwB6AfEAWusLgGX2WRS58lj1x5jTeQ5RiVH0X9mfkEvmj+MPblWVt7vX5I/9F3nt57AHw8C1ohEGt6OMMLBg+WwhRP4yJQiStDE6eWc/gocXuBFW17BsQxZ2X4i7kztD1w5lyYklZrc5tHU1xnatybKwC7z+SxipafcNSpcPgD7z4fJBWDwY0lLN/kwhRP4zJQh+Vkp9DbgppYYA6zC2rxT5rFLJSnzf/Xsal23Mf7b/h4mhE0nTmcz+yYERbavxRpcaLN13gTd+yeTKwK8LdBsPx1fBn++Y9VlCiILBlP0IJiilOgExQA3gP1rrtVbvmTBJSYeSTOs4jc92fca8g/M4HX2aT1t9iot97spXA7zUrjppaZov1x7nzI1bTPpXIJXcM7TXZIix6GzndKOWUdAwC/wkQoj8YsqCsjEY5aYj86ZL2StqC8pMpbVm4dGFjA8ZT41SNZjafipli5U1q83f953n3SXGLKFPnqxHr4AK/7yZlgo/PWtcGfRdCDVklrEQBZm5C8pKAn8qpbYqpV5SSpn37SKsQinFgFoDmNp+Kmdjz9J/RX8OXT9kVpuPBVZk5ehW+JUrwagf/+a1n/cRd2cfZBtbeGo2lPOHX1+EC4VmjyIhxH1MWUfwX611HYy9iysAm5VS60xpXCnVVSl1TCl1Uik1LpP3H0tfrbxPKbVHKdUyxz+BuEdrr9Z82+1b7GzseGHVC6w7Y9J/qixVcnfhp6FNGd3Bl6V/nyd4ylb2nYsy3nQoZswkcnaHhf+C6AJx0SiEyKGcLAy7AlwCrmNUI82WUsoWmIaxMrk20E8pdf8mNusx9joIBF4E5uSgPyILfqX8+CH4B/zc/RizaQxzDswxqyyFna0NYzr58dOwZqSkanrP2M60jSeNWUUlysGAnyEp3giDxEK5pbUQRZopC8pGKKU2YXxpewBDtNb+JrTdBDiptY5I39RmEfBYxgO01nH6n2+oYqRPURXm83D2YF6XeXTz6cbkvZN5d9u7JKUmPfzEbDT2dmfl6FZ0rVuOL/48Rv/ZO7kQdRvK1oGnF8CVI8bq49QUi/wMQoi8YcoVQRXgVa11Ha31+1rrwya2XRE4l+F5ZPpr91BKPaGUOgqswLgqeIBSamj6raM9V69ezewQkQlHW0c+b/U5IwNHsix8GUPWDOFmwk2z2nR1tmdqv/pM6BPAgfPRdJu8lVUHLkL1jhD8JZxcB6veAAsUxhNC5I0sg0ApVTL9r+OBs0op94wPE9pWmbz2wLeD1nqJ1rom8DjwUWYNaa1naa0baa0beXp6mvDR4g6lFCMCRjC+9XgOXjtI/xX9iYiKMLvN3g29WDmqFd6lXRjxw17GLd7PLf9nocVo2DPPKGEthCgUsrsiWJj+ZyiwJ/3P0AzPHyYSo1jdHV7AhawO1lpvAaoppczYiktkpZtPN+Z1ncetlFs8s/IZtl/Ybnab3h7F+HVEc0a2rcZPe87RY8pfHKw1Bmr1gjXvwpHlFui5EMLasgwCrXWP9D99tNZV0/+886hqQtshgK9Sykcp5QD0xdjr+C6lVPX07TBRSjUAHDAGo4UVBHgG8GPwj5QrXo6R60by87GfzW7T3taGN7vWZOHgptxKSuWJGTuYW2YcumJDWDwEIkMt0HMhhDVluaAs/Ys5S1rrvQ9tXKnuwCTAFpintf5EKTU8/fyZSqmxwHNAMnAbeENr/Vd2bcqCMvPFJ8fz5pY32RK5hWdqPcPrjV7H1sbW7HZvxifx1m8HWH3oEsFVbZkS9zq2qQkweD2UqmKBngshcitXW1UqpTam/9UJaASEYdz39wd2aa3zZc6/BIFlpKal8mXol3x3+DtaVWzF+NbjKe5Q3Ox2tdYsCjnHh8sPU9v+Iots38PerSK8+Cc4u1mg50KI3MjVymKtdTutdTvgDNAgfbC2IVAfOGmdroq8Ymtjy5uN3+S9pu+x/cJ2nl31LOfjzpvdrlKKfk0qs/yVliS4Vee5+FGkXj1B6s/PQWqyBXouhLA0U6aP1tRaH7jzRGt9EAi0XpdEXnq6xtPM6DiDy/GX6b+iP/uuWKZURPUyxfltZHPqtujBuORB2J7aTNQvL8u0UiEKIFOC4IhSao5Sqq1Sqo1SajZwxNodE3mnWYVmfB/8PcXsizHoz0GsjFhpkXYd7Wx5J7g2vV54k7k2T+F2dBF7fviPWauchRCWZ0oQDAQOAaOBV4HD6a+JR0hV16os7L6Qep71GLt1LNP3TbfYF3YrX08eHzOdncXa0ejkFGZMm8C1uESLtC2EMN9Dy1AXNDJYbF3Jqcn8d8d/+T38d7p5d+Ojlh/haOtokbZ18m2uTuuK681DDLf7gOeffpq2NR5atkoIYQHmlqEWRYi9rT0ftfiIMQ3HsOr0Kl5a/xK3km9ZpG1l70yZIb+hXCswMe1z3lvwBx8uP0xiimx5KUR+kiAQD1BK8WLdF/mk5SeEXAph2NphxCTFWKbxYqVxeO43XJ1sWVJyIou3HeDxads5cVmqlgqRXyQIRJZ6VevFhDYTOHj9IIP/HMyNhBuWadijOqrvQjxSLrGp0nc8MQMAACAASURBVBxuRMfSY+pffL/zjAwkC5EPTClD7aeUmq2UWqOU2nDnkRedE/mvU5VOTG0/lYjoCAauHsjl+MuWabhKc3hsGqWu7mZTzSUE+bjz7tKDDP0ulBvx5pXLFkLkjClXBL8Ae4F3gTcyPEQR0bJiS2Z2nMnlW5d5fvXzRMZaaCcy/6eh7ds4H/6ZBdU28V6P2mw+dpWuk7aw7eQ1y3yGEOKhTAmCFK31DK31bq116J2H1XsmCpRG5Roxp/McYpNieX7V82aXsr6rzZvg3xebTf9jUMkQlrzUnBJOdjwzdxefrjpCUkqaZT5HCJElU4JguVJqpFKqfA73IxCPmLoedZnfdT6pOpUXVr/AkesWWFeoFPSaAlVawu8vUSfpIH+80op+TSrz9eYInpqxnYirceZ/jhAiSw9dR6CUOpXJy9rEUtQWJ+sI8t+ZmDMMXjOY+KR4pnecTmAZC1QcuXUD5naGW9dg0DrwqM6fhy4xdvF+EpPT+G+vOvRp5EV61XIhRA6ZtY7gvn0IcrIfgXhEVSlZhW+7fkspp1IMXTuUnRd3mt+oizsM+BmUDSzsA/HX6VKnHKtHt6Z+ZTfeXLyflxbuJfqWFK4TwtJMmTVkr5QapZT6Nf3xslLKPi86Jwqu8sXL8023b6hYvCIvrXuJTec2md+oe1Xo+yNEn4dF/SE5gXKuTnw/KIhx3Wqy5tBluk3ewq4I2btICEsyZYxgBtAQmJ7+aJj+mijiPJw9WNB1AX6l/BizcQyrTq0yv9HKQfDETDi3E34fCWlp2NgohrepxuIRzXGws6Hv7J28s+SATDMVwkJMCYLGWuvntdYb0h8DgcbW7pgoHFwdXZndeTYBZQIYu2Usi48vNr/Ruk9Ch/fh4GLY9L+7LwdUcmPFqFa80NybRSHnaPvFRuZvO0VyqswsEsIcpgRBqlKq2p0nSqmqgBSHEXcVdyjOjI4zaF6xOR/s+IDvDn9nfqMtx0D9Z2HLF/D3D3dfLuZox/s967B6dCsCKrnx3+WH6TZ5K1uOXzX/M4UookwJgjeAjUqpTUqpzcAG4N/W7ZYobJztnJnSbgqdqnRifMh4ZobNNK9chFLQYyL4tIHloyBi8z1v+5YtwbcvNmHu841ISU3juXm7GfxNCKeuxZv5kwhR9JhUhlop5QjUwNiz+KjWOt+Kycv00YItJS2F97e/z7LwZQysM5AxDceYN+XzdhTM6wIxF2HwWvCs8cAhiSmpLNh2mqkbTpKYksqLLXx4uX11SjjJnAYh7sjV5vX3NdAc8Abs7rymtf7WUh3MCQmCgi9Np/G/Xf/jp2M/8bTf07zT9B1slBn1DW+egTkdwd4JBq+H4pnvYXAlNoEJfx7jl9BIShdz4M0uNend0AsbG1l7IIRZQaCU+g6oBuzjn7EBrbUeZdFemkiCoHDQWjNp7yTmHZxHj6o9+KjFR9jZ2D38xKxEhsKCYChbB174A+ydszx0f2QU/11+mNAzN6lX0ZX3e9amkbcshhdFm7lBcASorQtIfWAJgsJl9v7ZTPl7Ch0qd2B86/E42DrkvrEjy+GnZ6FWT+izAGxsszxUa82ysAt8tuooF6MT6BVQgXHdalLBLesAEeJRZu4OZQeBcpbtkigqhvgPYVyTcaw/u55XNrzC7ZTbuW+sVk/o/BEcWQYzmsPhZZDF7ydKKR4LrMj6f7dhVAdf/jx0ifZfbmLyuhMkJMukNyEyMuWKYCMQCOwG7g4Sa617WbdrmZMrgsJpyYklfLDjAwI9A/mqw1eUcCiRu4a0NoJgw8dw7TiUD4QO70G1DsZMoyxE3rzFp6uOsmL/RSq6OfN291p0r1dOaheJIsPcW0NtMntda705s9etTYKg8Fp9ejVvbXkLP3c/vu74NW5ObrlvLDUF9v8Emz6D6LNQpQW0fw+qNMv2tF0R1/nv8sMcvhhDE293/tOzNnUruua+H0IUEmbPGipIJAgKty2RWxizcQyVS1ZmVqdZeLp4mtdgSiLs/RY2j4f4K+DbGdq/C+UDsjwlNU3z855zfPHnMW7eSqJv40q83rkGpYs7mtcXIQowc68IngQ+B8pgrCNQGLOGSlq6o6aQICj8dl/czcsbXsbD2YM5nedQoXgF8xtNiofds+CvSZAQBbUfh3bvgKdflqdE305myvoTfLP9NM4Otozu4MtzzbxxsJOtvMWjx9wgOAn01FpbYBcS80kQPBrCroYxYt0IXOxcmNN5Dt6u3pZp+HYU7JhmPFJuQ0B/aDsW3CpnecrJK3F8vOIwm45dpapnMd7rUZt2NTJfqyBEYWVuEGzTWrewSs9yQYLg0XH0xlGGrR0GwKxOs6jh/uCq4VyLuwp/TYSQOaDToNFAaPU6lCib5Skbj17hoz8OE3EtnnY1PHm3R22qeRa3XJ+EyEfmBsFkjOmjS7l31tBvluykqSQIHi0R0REMWTOE2ym3mdlxJv6e/pb9gOjzsGU87P0ObB2g6XBoPsrYCCcTSSlpfLvjNJPXneB2ciovNPfmlQ6+uDpLuQpRuJkbBPMzeVlrrV804YO7ApMBW2CO1vqz+94fAIxNfxoHjNBah2XXpgTBoycyNpIha4ZwI+EGX3X4isblrFDl/Hq4McPowC/gWMIIg6bDjb9n4lpcIl+uOcaikHO4uzjwepcaPN2oErZSrkIUUvkya0gpZQscBzoBkUAI0E9rfTjDMc2BI1rrm0qpbsAHWuug7NqVIHg0Xbl1haFrhhIZF8n/tf0/Wnu1ts4HXT4EGz6BYyvAxQNa/RsavWjUMcrEwfPRfLj8MLtP36B2+ZK837M2QVVLW6dvQliRuVcETsAgoA5w9/+Wh10RKKWaYXyxd0l//lb6eZ9mcXwp4KDWumJ27UoQPLpuJtxk2NphnIg6wWetPqOLdxfrfVjkHtjwEURsgpIVoc2bEDgAbB+8BaS1ZsWBi/xvxREuRCcQ7F+et7rVxKuUi/X6J4SFmVti4juMMYIuwGbAC4g14byKwLkMzyPTX8vKIMACex2KwqqUUynmdplLPY96vLnlTZaeXGq9D/NqBM/9Ds8tg5IVYPlomNYEDvwKaffueKaUood/Bdb/uy1jOvqx/shlOny5mf9be5zbSVKuQhR+pgRBda31e0C81vobIBioZ8J5md1MzfTyQynVDiMIxmbx/lCl1B6l1J6rV2UnqkdZCYcSzOw4kyblmvDetvdYeGShdT+wahsYtBb6/QT2LrB4EMxsCUdXPlDHyNnBltEdfdnw77Z0qVOOKetP0P7LTSwLu2DeJjxC5DNTgiA5/c8opVRdwBVjb4KHiQQqZXjuBVy4/yCllD8wB3hMa309s4a01rO01o201o08Pc1ciSoKPBd7F77q8BXtKrXj092fMufAHOt+oFJQoysM2wpPzYWUBFjUz9gDIeLBSioV3JyZ0q8+vwxvRuniDoz68W/6zNzBgcho6/ZTCCsxZYxgMLAY8AfmA8WB/2itZz7kPDuMweIOwHmMweL+WutDGY6pjLH15XNa6+2mdFjGCIqO5LRk3v3rXVaeWsmguoMY3WB03hSJS02BsIWw6XOIiQSf1tD+P1DpwdlMqWmaxaGRjP/zKNfjkwiuV55RHXzxK5vLonpCWEm+1RpSSnUHJmFMH52ntf5EKTUcQGs9Uyk1B3gKOJN+SkpWHb1DgqBoSU1L5eNdH/Pr8V/pV7Mf45qMM2+3s5xIToDQ+bBlAty6BjW6G2UrytV94NDYhGRmbg5nwbbT3EpOpXtdIxBqlJNAEAWDubOGHDG+rL25d6vKDy3YR5NJEBQ9Wmu+3PMl3xz+hseqPcYHzT8wb7eznEqMg10zYdsUSIyBuk9Bu7ehdLUHDr0Zn8Tcv06xYPtp4hJT6F6vHKM6+FKzXL6U5hLiLnODYDUQDYTyz1aVaK2/tGQnTSVBUDRprZkZNpPpYdPpVKUTw/yHUdW1KvaZTPe0mts3jTDYNdOoelr/GWPaqavXA4dG3Upi3l+nmL/tNLGJKXStYwRC7QoSCCJ/mBsEB7XWD14L5xMJgqLtm0PfMGHPBADslB0+bj7UKFWDGqVq4Ofuh18pPzycPazbidjL8Nf/wZ55gILGg6Dla1D8wYkM0beSmbftFPO2nSI2IYXOtcsyqoOv7IEg8py5QTALmKq1PmCNzuWUBIE4G3OWg9cOcvzmcY7dPMbxG8e5cvvK3fdLO5WmhrsRDr6lfKnhXgMfVx/sbSx89RB1FjZ/DvsWgp0zNB0BzV7KtI5R9O1k5m87xby/ThGTkELHWmV5taMEgsg7uQoCpdQBjHn/doAvEIFRdO7OfgQWrg5mGgkCkZmbCTeNYLhxzAiHm8cJjwonOc2Y/WxvY081t2r4lTKuGu4ERSmnUuZ/+LUTsPF/cOg3IxAC+kLTkZnuhRCTkMyCbaeZszUiPRDKMLqDH/W8JBCEdeU2CKpk16jW+kx271uLBIEwVXJaMqejT9+9arhzBXHt9rW7x5RxLoOvu+/d20s13GtQpWSV3A1GXzkCO6dD2E+Qmgi+XaDZSPBp88B+yjEJyXyz7TRz/jpF9O1k2tcsw+gOvgRUMmP7TiGykdsgcAGStdbJ6c9rAN2BM/lVghokCIT5rt++zrGbxzhx88TdK4iI6AhS0lIAcLBxoJpbtbtXDTXca+BXyg9XRxN/a4+7CnvmGnshxF+FsnWNK4R6vcHu3u0wYxOS+XbHGWZvjSDqVjJta3gyuoMv9Stb4EpFiAxyGwRbgEFa6xNKqerAbuAHoDawW2v9lrU6nB0JAmENyanJRERHPHB76UbCjbvHlHUpezcc/Er54efuR5USVbC1sc2i0QSj7PXO6XDlMBQrA02GGtVOi91bwTQuMYVvd5xm9pYIbt5Kpo2fJ6M7+tJAAkFYSK7HCLTW9dL//hHgrrV+SSnlAITeeS+vSRCIvHTt9rV7guHYjWOcjj5NijauHpxsnajuVh0/dz9qlKpBm0ptqFj8vtqKWkPERmP7zJPrwM4pwzjCvbuyxSem3L1CuBGfRCtfD17t6EvDKplvpCOEqXIbBPvvDAgrpbYBX2itl6Y/D9NaB1irw9mRIBD5LSk1ifCo8HtmLR27eYyoxCgAGpRpQHDVYDpX6Yyb0333/K8cNa4Q9v9k1DSq3skYR6ja7p5xhPjEFL7feYZZWyK4Hp9Ey+oejO7oS2NvCQSRO7kNgu+BSxh1gsYBPlrrW0opN2CzBIEQ/9BaExkbyerTq1kRsYLw6HDsbOxoWaElwdWCaevVFie7DJvfxF8z1iHsng3xV6BMHSMQ6va+Z5OcW0kp/LDzLF9vCedaXBItqpdmdAc/mvhIIIicyW0QOAOjgfIYdYLC0l9vDlTTWn9npf5mS4JAFHRaa47dPMaKiBWsjFjJldtXKGZfjA6VOxBcNZigckH/jCukJBp7IOyYBlcOQTFPaDzEWKRW7J+FcbeTUvlh1xlmbo7gWlwizaqWZnRHX5rKbmnCRBYrOqeUaqC13muxnuWCBIEoTFLTUtlzeQ8rIlaw9sxa4pLj8HT2pKtPV4KrBlPbvbZRUVVrOLXZCIQTa8DWEQL+ZYwjlKl1t73bSaks3H2WmZvDuRqbSJCPO6929KNZNQkEkT1LBsFerXUDi/UsFyQIRGGVmJrI5nObWRGxgi3nt5CSloJ3SW+CqwYTXDWYSiXSt++4ejx9PcKPxjhCtQ7GiuVq7e+OIyQkp/Lj7rPM2BTOldhEmvi482oHX5pVK503pbpFoWPJIPhba13fYj3LBQkC8SiIToxm7Zm1/BHxB6GXQwEI8AwguGowXby74O7kDvHXITR9HCHuMnjWMsYR6j19dxwhITmVRbvPMmNzOJdjEmnsXYrRHfxoUV0CQdzLkkHw+J2ZQ/lFgkA8ai7GXWTlqZWsOLWCEzdPYKfsaF6xOcE+wbSt1BYXZQsHfzNuG10+AC4e0HiwMY5QvAxgBMLPe84xfWM4l2ISaFilFK929KVldQ8JBAFYIAiUUhWBKty7H8EWi/UwByQIxKPs2I1jrDhlDDJfvnUZZztnOlTuQI+qPQgq1wS7szuMQDi+2hhH8O8DTV+CsrUBSExJ5ec9kUzfeJKL0Qk0qOzGqA6+tPHzlEAo4sytPvo58C/gMP/sR6C11r0s2ksTSRCIoiBNpxF6OZQVEStYc2YNsUmxuDu5082nG8E+wdTFEbVrplH5NOW2sQ6h2ctQvQMoRWJKKr+kB8KF6AR8yxRnQFBlnmjghatzHu7hIAoMc4PgGOCvtU60RudySoJAFDVJqUlsjdzKilMr2HxuM0lpSVQpWYVgn2CCK7Sg8rG1sGsWxF0CjxrGOIL/v8DemaSUNJbuO88Pu84Sdi4KJ3sbevpXYEDTKgR4ucpVQhFibhCsAvporeOs0bmckiAQRVlMUgzrzqxjRcQKQi6FoNHU86hHcJUudElIwWPPAri0H1xKQ6NBxlhCibIAHDwfzQ+7zvL7vvPcSkqlToWS9A+qzGOBFSnumIdbf4p8YW4QLAYCgPUY+xEAoLUeZclOmkqCQAjDpfhLrD61mj8i/uDYzWPYKlualm9KcInqdIjYjcvxNWBrD/X6QJuxUMqoLB+bkMzv+y7ww66zHLkYQzEHWx6rX5EBQZWpU0H2RXhUmRsEz2f2utb6Gwv0LcckCIR40MmbJ+8OMl+Iv4CznTNtyzSiR1w8zQ79ib1S0P49CBoG6auatdb8fS6KhbvOsjzsAokpaQRWcqN/UGV6+lfA2SGLqqqiULLErCFnoLLW+pilO5dTEgRCZC1Np7Hvyj7+iPiDNWfWEJ0YjbuDKx+kFKdd+A6o2BB6TYWyde45L/pWMr/9HckPu85y8kocJZzseKqBFwOCKuNbtkQ+/TTCksy9IugJTAActNY+SqlA4EOZNSREwZacmsxf5/9i5v6ZHL9xnC99nqL9jnmQEA0tx0Cr1+8pcAfGVcLuUzdYuPssqw5cIik1jSbe7gxoWpmudcvhaCdXCYWVuUEQCrQHNt1ZVZxxr4K8JkEgRM7EJsUyfO1wDl8/zJfNPqD9wVVG+YrSvtBrClRpnul51+MS+TU0koW7z3Lm+i1KudjTp1El+jWpjI9HsTz+KYS5zA2CXVrroIzlJTLuVZDXJAiEyLl7wqDtl7RP0vDHqxB11tgxreMH4JT5QHFammZ7+HV+2HWGtYcvk5KmaVG9NAOCqtCpdlnsbW3y9GcRuWNuEMzFmDE0DngKGAXYa62HW7qjppAgECJ37obBjcN82eZL2pcLgo3/MwrcFS8HwV9Cze7ZtnElJoGf95zjx93nOB91G88SjjzdyIu+jStTyd0lj34SkRvmBoEL8A7QGVDAn8BHWusES3fUFBIEQuRebFIsw9YO48iNI/xfm/+jXeV2EBkKy14x9kOo/Th0G3937UFWUtM0W45f5YddZ9hw9AoaaOvnSf+gKrSr4YmdXCUUOJYsOmcLFNNax1iqczklQSCEeTINg9Rk2DYZNo83BpA7fwL1n7ln+8ysXIi6zaKQcyzafZYrsYmUd3Wib+PK/KtxJcq5Oj30fJE3zL0iWAgMx6gzFAq4Av+ntf7C0h01hQSBEObLNAwArp2A5aPhzDbwaQ09JkHpaia1mZyaxvojV1i4+yxbjl/F1kbRoWYZ+gdVprWvJzY2Us4iP5kbBPu01oFKqQFAQ2AsECqDxUIUbjFJMQxfO5wjN44wse1E2lZqa7yRlgZ7v4G1/4HUJGj3tlHh1Nb0MhRnr9/ix5Cz/BxyjuvxSVRyd6Zfk8r0aVgJzxKO1vmBRLbMDYJDQCCwEPhKa71ZKRUmm9cLUfhlGQYAMRdgxetwbAWU8zcWolUIzFH7SSlp/HnoEgt3nWVHxHXsbRWd65RjQFBlmlWVzXPykrlBMArjKiAMCAYqA99rrVtZuqOmkCAQwrJikmIYtmYYR28efTAMtIYjy2DlGxB/DZq/DG3GgUPOZwidvBLHj7vP8mtoJNG3k6nqUYyR7arzZP2KctsoD1hssDhDg3Za6xQTjusKTAZsgTla68/ue78mMB9oALyjtZ7wsDYlCISwvGzDAOD2TeNW0d5voZQP9JwMVdvk6rMSklNZeeAi87ed5sD5aGqXL8m7wbVoXt3D/B9EZMncKwJHjPUD3ty7Q9mHDznPFjgOdAIigRCgn9b6cIZjymDsfPY4cFOCQIj8kzEMJrWdRJtKmXzRn9piDCbfiDBmFXX+GJxL5erz0tI0fxy4yOerjnI+6jbta5bh7e41qV5GahtZQ3ZBYMpk39+Bx4AUID7D42GaACe11hFa6yRgUXo7d2mtr2itQ4BkE9oTQlhRSYeSfN35a2qWqsmrm15l87nNDx7k0xpGbIcWr8K+H+GrJnBoqXELKYdsbBS9Aiqw/t9teKtbTUJO36DLpK28s+QAV2MLxD5YRYYpVwQHtdZ1c9ywUr2BrlrrwenPnwWCtNYvZ3LsB0BcVlcESqmhwFCAypUrNzxz5kxOuyOEMFFMUgxD1wzl+M3jTGw7MfMrA4CLYcZCtIthUCMYgidAyQq5/twb8UlMWX+C73eewdHOhhFtqzGoZVUph20h5l4RbFdK5abAXGajPzn/tQHQWs/SWjfSWjfy9PTMTRNCCBOVdCjJrM6z8Cvlx5hNYzK/MgAoHwCDN0CnjyB8A0wLgpC5xvTTXHAv5sAHveqwZkxrWlT3YMKa47T/chOLQyNJS8vVV4cwkSlB0BIIVUodU0rtV0odUErtN+G8SKBShudewIXcdFIIkbfuD4MtkVsyP9DWDlqMgpHboUJ9WPEaLAg2FqblUlXP4sx6rhE/DW2KZwlH/v1LGD2/+ovt4ddy3abInim3hqpk9rrWOtv7M0opO4zB4g7AeYzB4v5a60OZHPsB2dwaykgGi4XIO9GJ0QxbO4zjN48zqd0kWnu1zvpgrWHfD/DnO5B8C9q8Cc1Hg51Drj8/LU2zfP8Fxq8+xvmo23SoWYa3ZEA5VyyxQ1lLwFdrPV8p5QkU11qfMuG87sAkjOmj87TWnyilhgNorWcqpcoBe4CSQBoQB9TOrpaRBIEQeStHYQAQexlWj4VDS6BMHWMhmldDs/qQkJzK/G2nmb7xJLeSU+nXpBKvdvTDo7isUjaVudNH3wcaATW01n5KqQrAL1rrFpbv6sNJEAiR96IToxm6dignbp4wLQwAjq6EFf+G2IvQdAS0ewcci5vVj+txicaA8q6zONvbpg8o++BkLwPKD2N2rSGgPrBXNqYRoujKVRgkxMD6/0LIHHCtDD0mgm9Hs/sSfjWOz1YdZe3hy1RwdeL1LjV4PFBWKGfH3FlDSdpIC53emOxRJ0QR5OroyqxOs/At5curG1/NegA5I6eSxoY3A1cb5a1/eAp+Gwrx183qSzXP4sx+rhE/DmlK6eKOvPZzGL2m/cWOcPPaLapMCYKflVJfA25KqSEYu5XNsW63hBAFUa7CAKBKMxi2FVq/CQcXw7TGsP/nXC1Ey6hZtdL8/lILJv4rgBtxSfSbvZPB3+zh5JU4s9otakwdLO6EsUMZwJ9a63VW7VU25NaQEPkvOjGaIWuGcDLqJJPbTaaVVw5qUF4+BMtGwfk94FIaKgUZj8rNjOqmdrkbAE5ITmXetlNM3xjO7eRU+jepzKsdfSktA8pALscIlFKx/LMA7P4bbwlAOEahuPWW6qgpJAiEKBjMCoO0VOPKIGITnN1h1C4CsHWEig3+CYZKTcDFPUf9uhaXyOR1J1i42xhQHtmuGi+2kAFla1QftQXqAj/kpvyEOSQIhCg47oRBeFQ4k9tPpmXFlrlrKO4KnNsFZ3caj4v7IC29wLFHDajc9J9HKR+TttA8eSWOz1YdYd2RK1R0c+aNLjXoFVChyA4oWzwIMjQ8TGv9da4byAUJAiEKFouFQUZJt+DC3n+C4dxuSIw23itW5t5gKOcPtvZZNrU9/Br/W3mEg+djqFfRlXeCa9G0amnz+1jIWC0I8oMEgRAFj1XCIKO0NLh6JD0Udhm3k6LOGu/Zu0DFhkYoVGoKlRqDk+t9p2uW7jvPF38e42J0Ap1ql2Vct5pU8zRvXUNhIkEghLA6q4fB/WIuZAiGnXBpP+g0QEHZOv+MM1QOAtdKoBQJyanM/esU0zeeJCEljQFBlRndoWgMKEsQCCHyRJ6HQUaJccZMpDu3kyJDICl9GmnJivcEw1UXXyZtCGdRyDlc7G0Z2a46A1t4P9IDyhIEQog8kzEMprSfQouK+VKNBlJT4MohOJt+K+nsTohNL4DsUBy8GnPdvT7zz5Vj3hkPSrmV4s2uNejp/2gOKEsQCCHyVIEJg4y0huhzRjCcS79quHwI0Ghly0kbb/5KrM6VUvXp2OtZGvp65XePLUqCQAiR5wpkGNwvIRrOhcC5neizO0k9G4JdWgKXtRtryg2j7dOjqFT60RhQliAQQuSL6MRoBq8ZTERUBFPbT6V5xeb53aXspSaTEL6Vm8vfo3zsQQ5pbw7UeYMej/WluKNdfvfOLOYWnRNCiFxxdXRldqfZVHWryisbXmH7+e353aXs2drj5Nee8q/9xc1uMyhvf5u+h19iz2ddWblp6yO7ZaYEgRDCqtyc3ApXGAAoRamg/riPDeN8wzcI4iCdNj7GsvHPE3o0Ir97Z3ESBEIIq8sYBqM2jiocYQBg70zFnu/iOCaM895P0jNhOdV+bMkvX73NuatR+d07i5EgEELkiTth4F3Su3CFAWBTsizeA+eQPHgz0aXq0OfaNFK+CuK3hV8Tl5Cc390zmwwWCyHyVFRCFIPXDOZ0zGn61+xPy4otqV+mPvbZ1AsqULTmRtgfJK98m7JJZ9mj6nCjxQd0bN+pQK8/kFlDQogCJSohine3vcu2C9tISUvB2c6ZoPJBtKzQkhYVW+BVohDM4U9N5ty6Gbju/ILiabGsd+xA6V4f0aBu7fzuWaYkCIQQBVJ8cjy7L+5mOx219QAAC7NJREFU24Vt/9/evQdHVd5hHP8+SbgkKzcTCCRcrRYUqYAIkqC2Uq230dZq1Rlt6x9VO+poW8eprZepbR3baa12nKl10KpT74oO0zpFq1UkAYpiKnKRxlAkIQkhKEhCgGV//eOchCUGBNzNWTy/zwzDnr0+YcI+e95z9n1Z2LCQhm0NAIwdOJbK8koqyyqZNnwahQWFESfdN9v+EbXP38mY2sdIWj6vFl/GlEtuY2RpSdTR9uJF4JzLeWbGuq3rukphadNSduzeQb/8fkwrnRYUQ3kl4waOQwewHkFv62iuZf2zN3PMpldpsiN56+jrOe2iaxlQmBsT2nkROOcOOx3JDpY1L2PhhoVUNVRRtyU4bbMsUdZVCjOGz+CIvrn1zd9NK1+nfd7NjO54n1UcRePJt3LamReSH/HxAy8C59xhr2FbA1UNVVQ1VLG4cTHtyXYKVMDkYZOpLK9kVvksxg8Znxt7C6kU6954lMSbv6IktYnqPjMpPPfXTJl8YmSRvAicc18ou3bvoqalJiiGDVWs3rwagJLCEirKKphVPouZI2YyuP/gSHPazjZWv/gbxqz8MwW2iwVDLmT8xXcyqrys17N4ETjnvtBa2luo3lBNVUMV1Y3VbNmxBSEmlUzq2luYWDyR/Lxo1hvYvrmBuqdv4dimeWwhwVtjr+bk79zEgERRr2XwInDOxcbu1G5WtK6gqiE46Lx803IMY1C/QVSMqOg6vlBS2Ptn9WyqfZvNc2/iy+3LWEcZ6068hcpzLic/P/vf7fUicM7F1scdH7OocRELG4KDzq0drQBMOHIClWVBKUweNpk+eb30hTYz6qrn0u+12ynfXc87BSeQ9427OOGk7K7m5kXgnHNAylKs+WhNVynUbKwhaUkSfRLMGD6jaxip7Ijsj+FbcifvzbuX0e/exwBro3rg2Yy9+C5Gjh6XldfzInDOuR5s27mNJU1LuoaRGtsaAUj0SVBaVEppUSnDE8MpTZR2bZcmgusG9BmQkTOUOra2surp2zi+/il2UsBbo77H1EtuZcCAQZ/7udN5ETjn3GcwM9ZuXcuiDYuo/6SeprYmmtubaW5rpmV7C8be75WFBYVdxZBeGsMTw7u2B/UbdMBlsWndKjY8dzNf+WQBTRSz9oSfMP38a8jPz8wB7siKQNJZwH1APjDHzO7udrvC288B2oHvm9my/T2nF4FzrrftSu2idXsrTW1NNLU30dzW3FUSze3NNLU10bK9hZSl9npc//z+n9qb6NrLCLeH9BuyV1l8sHQ+zP85X0r+lzX5x7Bj9i+ZVHH25/4ZIikCSfnAGuAMoB5YClxmZivT7nMOcD1BEcwA7jOzGft7Xi8C51wuSqaStG5vDQoiLIeuwghLY2P7RpKW3OtxffP6Mqxo2N5DUIXDaF9dw4TVLzAh2crawgrKvv1byo+aeMj59lcE2VyEczpQa2Z1YYingAuAlWn3uQB4zII2WixpsKQRZtaYxVzOOZdxBXkFwRt5onSf90lZak9ZtDUHexedl9uaqNlYQ3N7M8lUWBblhcBICmwdw167iMqqqdx+xeOZz57xZ9yjHFiftl1P8Kn/s+5TDuxVBJKuAq4CGD16dMaDOudcb8hTHkOLhjK0aCjHlxzf431SluKjjo/2Kou1zWtY+/4rFA8ck5Vc2SyCno6QdB+HOpD7YGYPAg9CMDT0+aM551xuylMexYXFFBcWc1xxuLbBBOC0O7L3mll75uDT/ai07ZHAhkO4j3POuSzKZhEsBY6RNE5SX+BSYF63+8wDvqvAycAWPz7gnHO9K2tDQ2aWlHQdMJ/g9NGHzWyFpGvC2x8AXiI4Y6iW4PTRK7OVxznnXM+yeYwAM3uJ4M0+/boH0i4bcG02MzjnnNu/7E9555xzLqd5ETjnXMx5ETjnXMx5ETjnXMwddrOPSmoB1h3iw0uATRmMkym5mgtyN5vnOjie6+B8EXONMbOhPd1w2BXB5yHprX1NuhSlXM0FuZvNcx0cz3Vw4pbLh4accy7mvAiccy7m4lYED0YdYB9yNRfkbjbPdXA818GJVa5YHSNwzjn3aXHbI3DOOdeNF4FzzsVcbIpA0lmS3pdUK+mnUecBkPSwpI2S3os6SzpJoyT9S9IqSSsk3RB1JgBJ/SX9W9J/wly/iDpTOkn5kt6R9Leos3SS9D9JyyXVSMqZxb7DZWmfk7Q6/D2bmQOZxof/Tp1/tkq6MepcAJJ+FP7OvyfpSUn9M/r8cThGICkfWAOcQbAYzlLgMjNbud8HZj/XqcA2gnWbe163LgKSRgAjzGyZpAHA28A3c+DfS0DCzLZJ6gMsBG4ws8VR5uok6cfANGCgmZ0XdR4IigCYZmY59eUoSY8Cb5rZnHC9kiIz+zjqXJ3C94wGYIaZHeoXWDOVpZzgd/04M9su6RngJTN7JFOvEZc9gulArZnVmdlO4CnggogzYWYLgM1R5+jOzBrNbFl4+RNgFcFa0pGywLZws0/4Jyc+yUgaCZwLzIk6S66TNBA4FXgIwMx25lIJhGYDH0RdAmkKgEJJBUARGV7JMS5FUA6sT9uuJwfe2A4HksYCU4Al0SYJhMMvNcBG4BUzy4lcwL3AzUAq6iDdGPCypLclXRV1mNBRQAvwl3AobY6kRNShurkUeDLqEABm1gD8DvgQaCRYyfHlTL5GXIpAPVyXE58kc5mkI4DngRvNbGvUeQDMbLeZTSZY33q6pMiH1CSdB2w0s7ejztKDSjObCpwNXBsOR0atAJgK/MnMpgBtQE4ctwMIh6rOB56NOguApCEEIxjjgDIgIenyTL5GXIqgHhiVtj2SDO9afdGEY/DPA4+b2dyo83QXDiW8DpwVcRSASuD8cDz+KeB0SX+NNlLAzDaEf28EXiAYJo1aPVCftjf3HEEx5IqzgWVm1hx1kNDXgbVm1mJmu4C5QEUmXyAuRbAUOEbSuLDtLwXmRZwpZ4UHZR8CVpnZPVHn6SRpqKTB4eVCgv8gq6NNBWZ2i5mNNLOxBL9br5lZRj+xHQpJifBgP+HQy5lA5GeomVkTsF7S+PCq2UCkJyJ0cxk5MiwU+hA4WVJR+H9zNsFxu4zJ6prFucLMkpKuA+YD+cDDZrYi4lhIehL4KlAiqR64w8weijYVEHzCvQJYHo7HA/wsXIM6SiOAR8MzOvKAZ8wsZ07VzEGlwAvBewcFwBNm9o9oI3W5Hng8/GBWB1wZcR4AJBURnF14ddRZOpnZEknPAcuAJPAOGZ5qIhanjzrnnNu3uAwNOeec2wcvAuecizkvAuecizkvAuecizkvAuecizkvAhdrkv6QPsOkpPmS5qRt/17S7Qc7Y62kRyRdlMmszmWLF4GLu2rCb2lKygNKgIlpt1cA883s7giyOdcrvAhc3FWx5+v6Ewm+efuJpCGS+gHHAidIuh+6Pun/UVK1pLrOT/0K3C9ppaS/A8M6X0DS7HByteXhGhT9JE2XNDe8/QJJ2yX1DddcqOvFn985LwIXb+FcPElJowkKYRHBTKszCdYWeBfY2e1hI4BZwHlA557Ct4DxwCTgB+zZy+gPPAJcYmaTCL7h+0OCb4lOCR97CkEBnQTMIEdmenXx4UXg3J69gs4iWJS2Xd3D/V80s1S4UE9peN2pwJPh7KgbgNfC68cTTBi2Jtx+FDjVzJJAraRjCSaCuyd8jlOANzP9Azq3P14Ezu05TjCJ4JP5YoI9ggqCkuhuR9rl9CnOe5qvpacp0Du9STDT5S7gnwR7GbOABQca3LlM8CJwLnizPw/YHH6i3wwMJiiDRQf4HAuAS8OFc0YAXwuvXw2MlXR0uH0F8EbaY24EFplZC1AMTAAinxDRxUssZh917jMsJzhb6Ilu1x1hZpvC2Ts/ywvA6eHj1hC+2ZtZh6QrgWfDZQaXAg+Ej1lCMLTUuQfwLsECNz4TpOtVPvuoc87FnA8NOedczHkROOdczHkROOdczHkROOdczHkROOdczHkROOdczHkROOdczP0fE1TU1qcQfc8AAAAASUVORK5CYII=\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": { "execution": { "iopub.execute_input": "2021-05-19T06:00:01.247319Z", "iopub.status.busy": "2021-05-19T06:00:01.246438Z", "iopub.status.idle": "2021-05-19T06:00:03.689365Z", "shell.execute_reply": "2021-05-19T06:00:03.689783Z" } }, "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": { "execution": { "iopub.execute_input": "2021-05-19T06:00:03.696587Z", "iopub.status.busy": "2021-05-19T06:00:03.695214Z", "iopub.status.idle": "2021-05-19T06:00:03.699634Z", "shell.execute_reply": "2021-05-19T06:00:03.699083Z" } }, "outputs": [ { "data": { "text/plain": [ "array([[0.5480449 ],\n", " [0.40686249],\n", " [0.32580422],\n", " [0.25577172],\n", " [0.1872541 ],\n", " [0.11802003],\n", " [0.06266304],\n", " [0.03212954],\n", " [0. ]])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dres_conv" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T06:00:03.710016Z", "iopub.status.busy": "2021-05-19T06:00:03.709181Z", "iopub.status.idle": "2021-05-19T06:00:03.815587Z", "shell.execute_reply": "2021-05-19T06:00:03.816008Z" }, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Jensen-Shannon divergence')" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3jV9fn/8edNEsLeYQiEvRFkiwxBBUFRRLCK1rbWqlgRtVVr+21rrb+2ropbVFq1DlBxfkVFUGSL7L3DFEvC3iPk/v1xDjRfGuBknHxOcl6P6zoX+Zz50gtyn/c2d0dEROJXiaADiIhIsFQIRETinAqBiEicUyEQEYlzKgQiInEuMegAuVWtWjWvX79+0DFERIqUefPmbXf3lJweK3KFoH79+sydOzfoGCIiRYqZbTzdY+oaEhGJcyoEIiJxToVARCTOqRCIiMQ5FQIRkTinQiAiEudUCERE4lzcFIJtew/z0P8u49jxrKCjiIjElLgpBAs27eLVGRt4fMKqoKOIiMSUuCkE/VrX4sbz6/Hy1DS+WrEt6DgiIjEjbgoBwP9c3oKWtSrw6/cWsXX3oaDjiIjEhLgqBKWSEnj+hvYcy8zizjELNF4gIkKcFQKABtXK8rfBbZi3cRdPfKnxAhGRuCsEAFe2PYfru6Ty0pQ0Jq9MDzqOiEig4rIQAPxxQEua1yzPr95dyA97NF4gIvErbgvBifGCI5lZjBizgEyNF4hInIrbQgDQKKUcfx10LnM27OLJiauDjiMiEoi4LgQAV7WrzXWd6vLCN+uYsjoj6DgiIoUu7gsBwINXtKJZjfLc885C/r3ncNBxREQKlQoBULpkaLzg8LHjjBir8QIRiS8qBGGNq5fj/13Vmu/W7+Tpr9YEHUdEpNCoEGRzdfs6XNOhDs9NXsv0NduDjiMiUiiiWgjMrJ+ZrTKztWb2QA6P9zKzPWa2MHz7YzTzROKhga1onFKOu99ZQPpejReISPEXtUJgZgnA80B/oCUw1Mxa5vDUae5+Xvj252jliVSZkom8cEN7Dhw5zl1jF3I8y4OOJCISVdFsEXQG1rp7mrsfBcYCA6P4eQWmSY3y/HlgK2al7eAZjReISDEXzUJQG9ic7XpL+L5TdTWzRWb2uZm1imKeXLmmY12ubl+bZ75ew8y1Gi8QkeIrmoXAcrjv1H6W+UA9d28LPAt8lOMbmd1qZnPNbG5GRuEt+np4YGsaVivLiLELSd+n8QIRKZ6iWQi2AHWzXdcBtmZ/grvvdff94Z8/A5LMrNqpb+TuL7t7R3fvmJKSEsXI/1fZ5EReuKED+48c4553NF4gIsVTNAvBHKCJmTUws5LAdcAn2Z9gZjXNzMI/dw7n2RHFTLnWrGZ5HrqyFTPW7uD5yWuDjiMiUuASo/XG7p5pZsOBCUAC8E93X2Zmw8KPjwKGALebWSZwCLjO3WPua/ePOtZl1rodPDVpNZ3qV6Fro6pBRxIRKTAWg793z6hjx44+d+7cQv/cA0cyueK56ew/nMlnd/WgWrnkQs8gIpJXZjbP3Tvm9JhWFkeobHIiz1/fnj2HQuMFWRovEJFiQoUgF1rUqsCDV7Ri2prtvDhlXdBxREQKhApBLg3tXJcr2p7D379cxXfrdwYdR0Qk31QIcsnM+Oug1tSrWpY7x8xnx/4jQUcSEckXFYI8KF8qieeub8eug8f41buLNF4gIkWaCkEetTqnIn8Y0JIpqzN4aWpa0HFERPJMhSAfftwllcvPrcUTX65izgaNF4hI0aRCkA9mxt8Gn0udyqUZMWYBuw4cDTqSiEiuqRDkU4VSSTx/fXt27D/Kr9/TeIGIFD1nLQRm1tTMvjKzpeHrNmb2++hHKzpa167I/1zegq9XpjN6usYLRKRoiaRF8ArwW+AYgLsvJrSBnGTzk6716N+6Jo9+sYp5G3cFHUdEJGKRFIIy7v7dKfdlRiNMUWZmPDqkDedUKsWIMQvYfVDjBSJSNERSCLabWSPCh8qY2RDgh6imKqJOjBek7zvMve8toqht6Cci8SmSQnAH8BLQ3My+B+4Gbo9qqiKsTZ1K/LZ/CyatSOcf09cHHUdE5KzOeh6Bu6cBl5hZWaCEu++Lfqyi7aZu9fk2bQePfL6SDvUq0y61ctCRREROK5JZQ381s0rufsDd95lZZTP7f4URrqgyMx4f0paaFUsx/O0F7Dl4LOhIIiKnFUnXUH93333iwt13AZdFL1LxULFMEs8Obce2vYe5b5zGC0QkdkVSCBLM7ORxXGZWGtDxXBFol1qZB/o358vl23h1xoag44iI5CiSQvAm8JWZ3WxmPwcmAq9HN1bxcXP3BlzSojp/+3wFizbvPvsLREQK2VkLgbs/BvwFaAG0Ah4O3ycRMDOeuKYt1cuXYviY+ew5pPECEYktEe015O6fu/u97v5rd58Q7VDFTaUyJXlmaDt+2H2YB95frPECEYkpkcwautrM1pjZHjPba2b7zGxvYYQrTjrUq8z9/Zrx+dJ/869ZG4OOIyJyUiQtgseAK929ortXcPfy7l4h2sGKo190b8hFzavzl/ErWPr9nqDjiIgAkRWCbe6+IupJ4kCJEsbfr2lL1XIluePt+ew9rPECEQleJIVgrpm9Y2ZDw91EV5vZ1VFPVkxVLluSZ4e2Y8uuQ/z2/SUaLxCRwEVSCCoAB4G+wBXh24BohiruOtavwr19mzF+yQ+8OXtT0HFEJM5FstfQTYURJN7c1rMhs9fv4OFPl9M+tRKtzqkYdCQRiVNRPaHMzPqZ2SozW2tmD5zheZ3M7Hh4i+u4cGK8oHKZJO54az4Z+44EHUlE4lTUTigzswTgeaA/0BIYamYtT/O8R4G4W59QtVwyzw5tzw97DnP5M9P4bv3OoCOJSByK5gllnYG17p7m7keBscDAHJ53J/A+kB7BexY7nRtU4aM7ulGmZAJDX/mWV6amaQBZRApVNE8oqw1szna9JXzfSWZWGxgEjIoobTHVolYFPrmzO31a1OAvn61g2JvzNLVURApNNE8osxzuO/Wr7lPAb9z9+BnfyOxWM5trZnMzMjIi+Oiip0KpJF78cXt+f3nodLMrn53O8q1awC0i0RfJpnNp7n4JkAI0d/fu7r4hgvfeAtTNdl0H2HrKczoCY81sAzAEeMHMrsohw8vu3tHdO6akpETw0UWTmfGLHg0Ze+v5HDx6nEEvzOC9uZvP/kIRkXw46/RRM/vVKdcAe4B57r7wDC+dAzQxswbA94QGmK/P/gR3b5DtfV8DPnX3jyINX1x1ql+F8SN6cNfYBdw3bjHzNu7iT1e2olRSQtDRRKQYiqRrqCMwjFD/fm3gVqAX8IqZ3X+6F7l7JjCc0GygFcC77r7MzIaZ2bD8Bi/uUson88bNXbijdyPGztnM4BdnsmnHwaBjiUgxZGeboWJmE4DB7r4/fF0OGEdokHeeu//XlNBo6tixo8+dO7cwPzJwX63Yxj3vLMSBJ390Hn1a1gg6kogUMWY2z9075vRYJC2CVOBotutjQD13PwRoFVQhuLhFDcaP6EG9qmW45V9z+dvnK8g8nhV0LBEpJiIpBG8D35rZg2b2IDADGGNmZYHlUU0nJ9WtUoZxwy7g+i6pvDQljRtGzyZ93+GgY4lIMXDGriELjQzXAaoD3QlNCZ3u7oH1zcRj19CpPpi/hd99uITypZJ4dmg7zm9YNehIIhLj8tw15KEq8ZG7z3P3p939qSCLgIRc3b4OH9/RnfLJidwwejajpqzTamQRybNIuoa+NbNOUU8iudKsZnk+Ht6NS1vV4JHPV3LrG/PYc0irkUUk9yIpBL0JFYN1ZrbYzJaY2eJoB5OzK18qieevb88fB7Rk8sp0rnxuOsu26ghMEcmdsy4oI7R7qMQoM+Pn3RvQtm5F7nhrAYNemMnDA1txbafUoKOJSBERyRYTGwltFXFR+OeDkbxOCleHelX4dER3OtWvzG/eX8J97y3i8LEzbuEkIgJEdjDNg8BvCJ1JAJAEvBnNUJI31col86+fd+HOixrz3rwtDHphJhu2Hwg6lojEuEi+2Q8CrgQOALj7VqB8NENJ3iWUMH7dtxmv3tSJH/Yc4opnpzNh2b+DjiUiMSySQnA0PI30xHkEZaMbSQpC72bV+fTO7jRMKcttb8zjr5+t4JhWI4tIDiIpBO+a2UtAJTO7BZhE6PhKiXF1Kpfh3WFdufH8erw8NY0bXpnNtr1ajSwi/9dZN50DMLM+QF9CK4snuPvEaAc7Ha0szpuPF37PA+8voWxyIs8ObUfXRlqNLBJP8rXpnJndA6xw9/vc/d4gi4Dk3cDzavPx8G5UKJ3IDaO/5YVv1pKVpdXIIhJZ11AFYIKZTTOzO8xMeyAXUU1rlOeT4d257NxaPPbFKm59Yy57Dmo1ski8i2QdwUPu3orQ2cXnAFPMbFLUk0lUlAt3DT10ZSumrM5gwHPTWPq9ViOLxLPcLAxLB/4N7CC0G6kUUWbGTy+ozzu3deX4cefqF2cy5rtN2rhOJE5FMkZwu5l9A3wFVANucfc20Q4m0dc+tTKfjuhBlwZV+O0HS7j3vcUcOqrVyCLxJpK9huoBd5/loHopoqqULclrN3Xmma/W8MzXa1i2dQ8v/rgDDappuYhIvDhti8DMKoR/fAzYZGZVst8KJ54UhoQSxj19mvLaTZ3ZtvcwVzw7nS+W/hB0LBEpJGfqGno7/Oc8YG74z3nZrqWYubBpCp+O6EHj6uUY9uZ8/t+nyzmaqdXIIsVdRAvKYokWlEXf0cws/vrZCl6buYHmNcvzxDVtaV27YtCxRCQfzrSg7LSFwMzan+lN3X1+AWTLNRWCwvPVim387sMlbN9/lF/2asTwixqTnJgQdCwRyYMzFYIzDRb/PfxnKaAjsIjQFhNtgNmEDrOXYuziFjX4sl4VHh6/nGe/XsuXy7bx+DVtaFOnUtDRRKQAnXaMwN17u3tvYCPQ3t07unsHoB2wtrACSrAqlkniiWva8urPOrHn0DEGvTCTxyes5EimppmKFBeRLChr7u5LTly4+1LgvOhFkljUu3l1JtzTk8Hta/P85HUMeGY6izbvDjqWiBSASArBCjMbbWa9zOxCM3sFWBHtYBJ7KpZO4rEhbXntpk7sP5LJoBdm8OgXK3UkpkgRF0khuAlYBtwF3A0sD98ncapXs1Dr4Ecd6/LiN+sY8Ox0FmzaFXQsEcmjSDadO+zuI919UPg20t0jOt3EzPqZ2SozW2tmD+Tw+EAzW2xmC81srplpALqIqFAqiUcGt+H1n3fm4JFMBr84k799vkKtA5EiKDebzuWKmSUAzwP9gZbAUDNrecrTvgLauvt5wM+B0dHKI9FxYdMUJtzTk2s7pfLSlDQuf2Ya8zaqdSBSlEStEACdgbXunubuR4GxwMDsT3D3/f6fhQxlCZ+LLEVL+VJJ/O3qc3nj5s4cPpbFkFEz+cv45WodiBQR0SwEtYHN2a63hO/7P8xskJmtBMYTahX8FzO7Ndx1NDcjIyMqYSX/ejQJtQ6u75zKK9PWc9nT05i3cWfQsUTkLCLZhrqpmb1iZl+a2dcnbhG8t+Vw339943f3D929OXAV8HBOb+TuL4fXMXRMSUmJ4KMlKOWSE/nLoHN56xddOJKZxZBRs3j40+Xa3lokhkWyDfV7wCjgFSA3/5q3AHWzXdcBtp7uye4+1cwamVk1d9+ei8+RGNStcTUm3NOTRz9fyT+mr+frlek8NqQNnepr41qRWBNJ11Cmu7/o7t+5+7wTtwheNwdoYmYNzKwkcB3wSfYnmFljM7Pwz+2BkoROQJNioFxyIg9f1Zq3b+lCZlYWP3ppFg/97zIOHs0MOpqIZBNJIfhfM/ulmdXKzXkE7p4JDAcmEFqA9q67LzOzYWY2LPy0wcBSM1tIaIbRtdkGj6WYuKBRNb64qyc/Ob8er87YQP+npzE7TfVeJFacdRtqM1ufw93u7g2jE+nMtPto0fZt2g7uH7eYTTsP8rML6nN/v2aUKRlJD6WI5EeetqGOVSoERd/Bo5k89sUqXpu5gdQqZXh0cBu6NqoadCyRYu1MhSCSWUNJZjbCzMaFb8PNLKngY0q8KFMykT9d2Yp3b+tKCYOhr3zLHz9eyoEjGjsQCUIkYwQvAh2AF8K3DuH7RPKlc4MqfH5XT37erQFvfLuRS5+aysx1mjAmUtgiGSNY5O5tz3ZfYVHXUPE0d8NO7hu3mPXbD/Dj81N5oH8LyiVr7ECkoOSrawg4bmaNsr1ZQ3K3nkDkrDrWr8JnI3rwi+4NeGv2Ji4dOZUZa9U6ECkMkRSC+4DJZvaNmU0BvgZ+Hd1YEo9Kl0zg9wNaMm5YV5ITS3DD6Nn87sMl7Dt8LOhoIsVaRLOGzCwZaEZo24iV7n4k2sFOR11D8eHwseM8OXE1o6elUatiaR4ZfC49mmh7EZG8ym/XEIQGiFsDbYFrzewnBRVOJCelkhL43WUtGHf7BZRKKsGN//iO336wWK0DkSiIZProG8ATQHegU/iWY1URKWjtUyszfkQPbruwIe/M2cylI6cyZbV2oBUpSJHMGloBtIyVrR/UNRS/FmzaxX3jFrM2fT9DOtThD5e3pGIZLWkRiUR+u4aWAjULNpJI7rVLrcynd3bnjt6N+HDB91wycgpfLP0h6FgiRV4khaAasNzMJpjZJydu0Q4mkpNSSQncd2lzPhnejerlkxn25nxuf3Me6fsiOkZbRHIQyYqdP0U7hEhutTqnIh/d0Y1XpqXx1KQ1zFy3gz8MaMng9rUJ72wuIhHSpnNS5K3L2M8D7y9mzoZd9Gyawl8HtaZO5TJBxxKJKfnddO5qM1tjZnvMbK+Z7TOzvQUfUyRvGqWU451bu/Lnga2Yt2EnfUdO5fWZG8jKKlpfckSCEskYwWPAle5e0d0ruHt5d68Q7WAiuVGihPGTrvWZcE9POtavwoOfLONHL81iXcb+oKOJxLxICsE2d18R9SQiBaBO5TK8flMn/n5NW9ak76f/09N4fvJajh3PCjqaSMyKZLB4rpm9A3wEnNxawt0/iFoqkXwwMwZ3qEPPpik8+MlSHp+wivGLf+CxIW1oXbti0PFEYk4kLYIKwEGgL3BF+DYgmqFECkJK+WReuKEDo37cnoz9Rxj4/Awe/WIlh49p81yR7DRrSOLCnoPH+Mtny3l37hYaVivLo0Pa0Kl+laBjiRSafJ1ZbGalgJuBVkCpE/e7+88LMmSkVAgkP6av2c4DHyxmy65D/KRrPe7v11wH4EhcyO8WE28Q2mLiUmAKUAfYV3DxRApP9ybVmHB3T27qVj90PKY2sROJqBA0dvc/AAfc/XXgcuDc6MYSiZ6yyYk8eEUrxg27gNIlE/jpP7/jV+8uZPfBo0FHEwlEJIXgxAbwu82sNVARqB+1RCKFpEO9yowf0Z07L2rMJwu3csmTU/hsiTaxk/gTSSF42cwqA38APgGWE1pkJlLkJScm8Ou+zfhkeHdqVSzNL9+az7A35pG+V5vYSfzQrCGRsMzjWYyevp6RE1eTnFiC3w9oyTUd6mgTOykW8jtrKBkYTKg76OT0Cnf/cwFmjJgKgURbWsZ+Hnh/Cd9t2EmPJtX466BzqVtFm9hJ0ZbfWUMfAwOBTOBAtlskH9zPzFaZ2VozeyCHx28ws8Xh20wzaxvJ+4pEU8OUcoy99Xwevqo18zfu4tKnpvLqjPUc1yZ2UkxF0iJY6u6tc/3GZgnAaqAPsAWYAwx19+XZnnMBsMLdd5lZf+BP7t7lTO+rFoEUpu93H+J/PlzCN6syaJ9aiceGtKFx9fJBxxLJtfy2CGaaWV6mi3YG1rp7mrsfBcYSalmc5O4z3X1X+PJbQmsURGJG7UqlefVnnRh5bVvSth/gsqen89zXa7SJnRQrpy0EZrbEzBYD3YH54S6exdnuP5vawOZs11vC953OzcDnp8lyq5nNNbO5GRla/COFy8wY1K4Ok351IX1b1eCJL1dzxbPTWbJlT9DRRArEmdbW53djuZymWuTYD2VmvQkVgu45Pe7uLwMvQ6hrKJ+5RPKkWrlknru+PVe2/Te//2gpV70wg1t6NOTuS5pQKikh6HgieXamrqEMYKu7b3T3jYT2Gboa6BC+PpstQN1s13WArac+yczaAKOBge6+I+LkIgHp26omE391Idd0qMOoKevo//Q0Zqfpr64UXWcqBF8QXkFsZo2BWUBD4A4z+1sE7z0HaGJmDcysJHAdoQVpJ5lZKvABcKO7r859fJFgVCydxCOD2/DWL7qQmZXFtS9/yx8+Wsr+I5lBRxPJtTMVgsruvib880+BMe5+J9CfCLqN3D0TGA5MAFYA77r7MjMbZmbDwk/7I1AVeMHMFpqZpgNJkdKtcWgTu5u7N+DN2drEToqm004fNbPF7t4m/PMM4HF3/yh8vcjdA5nzr+mjEqvmb9rF/eMWszZ9P9d0qMPvL29JxTJJQccSAfI+fXSxmT1hZvcAjYEvw29WKQoZRYq89qmV+fTO7tzRuxEfLPiePiOnMHH5tqBjiZzVmQrBLcB2QuMEfd39YPj+lsATUc4lUiSVSkrgvkub8/Ed3ahStiS3/GsuI8YsYOcBbXEtsStXm86ZWXt3nx/FPGelriEpKo5mZjFqyjqe/XoNFUol8dDAVlx+bi1tYieByO/K4uxGF0AekbhQMrEEIy5uwqd39qBO5dIMf3sBw96cR/o+bXEtsSW3hUBfZURyqVnN8rx/+wX8tn9zJq/KoM+TU3l/3haK2hbwUnzlthA8FJUUIsVcYkIJbruwEZ/f1YMm1cvx6/cWcdNrc9i6+1DQ0UQiKwRmVju8U+hOM+tpZj2jnEukWGqUUo53b+vKn65oyey0nfQdOZW3Z29S60ACdaa9hgAws0eBawkdUXk8fLcDU6OYS6TYKlHC+Fm3BlzUvAYPfLCY3324hE8Xb+WRq9uQWlUH4Ejhi+Q8glVAG3c/UjiRzkyzhqQ4cXfGztnMX8av4HiWc3+/Zvy0a31KlNBwnBSs/M4aSgO0PFIkCsyMoZ1T+fKennRpWIWH/nc5P3ppFusy9gcdTeJIJIXgILDQzF4ys2dO3KIdTCSenBM+AOfJH7VlTfp++j89jVFT1pGpA3CkEJx1jIDQjqGfnPVZIpIvZsbV7evQvUk1/vDRUh75fCWfLfmBx4a0oXnNCkHHk2IsopXFZlYaSHX3VdGPdGYaI5B44O6MX/IDD368jL2HjzG8dxNu79WIkom5nfEtEpKvMQIzuwJYSOh8AszsPDNTC0EkisyMAW3O4ct7etK/dS1GTlrNlc/peEyJjki+XvyJ0EH0uwHcfSHQIIqZRCSsarlknhnajpdv7MDOA0e56oUZPPbFSg4fO372F4tEKJJCkOnup34N0eoXkULUt1VNJt5zIVe3q80L36zj8memMW/jrqBjSTERSSFYambXAwlm1sTMngVmRjmXiJyiYpkkHr+mLa//vDOHjh5nyKiZPPzpcg4dVetA8ieSQnAn0Ao4AowB9gJ3RzOUiJzehU1TmHBPT27okso/pq+n39NTmbVuR9CxpAjL7XkECUBZd98bvUhnpllDIv8xa90OfvP+YjbtPMiPz0/lgf4tKJccyaxwiTf5nTX0tplVMLOywDJglZndV9AhRST3ujaqyhd39+Dm7g14a/YmLh05lSmrM4KOJUVMJF1DLcMtgKuAz4BU4MaophKRiJUpmcgfBrRk3LALKJVUgp/+8zvue28Rew4eCzqaFBGRFIIkM0siVAg+dvdjaNaQSMzpUK8y40f04Je9GvHBgu/pM3IKE5dvCzqWFAGRFIKXgA1AWWCqmdUjNGAsIjGmVFIC9/drzke/7EaVsiW55V9zGTFmAel7dTymnF6uBotPvsgs0d0zo5DnrDRYLBKZo5lZvPjNOp6bvIaEEsZN3RpwW8+GVCpTMuhoEoAzDRZHch5BMjAYqE+2Terc/c8FmDFiKgQiubNh+wGemrSajxdtpVxyIrf1bMhN3RpQVrOL4kp+zyP4GBgIZAIHst1EpAioX60sT13Xjs9G9KBLg6o88eVqLnx8Mq/OWM+RTC1Gk8haBEvdvXUh5TkrtQhE8mfexl08PmEl36btpHal0tx1SROublebxATtbFqc5bdFMNPMzs3jB/czs1VmttbMHsjh8eZmNsvMjpjZvXn5DBHJnQ71KjPmlvN58+YuVCtXkvvHLabvU1MZv/gHsrI0ITAeRdIiWA40BtYT2mbCAHf3Nmd5XQKwGugDbAHmAEPdfXm251QH6hGamrrL3Z84W2C1CEQKjrszYdk2/v7lKtak76d17Qrc27cZFzZNwUznJhcnZ2oRRDJa1D+Pn9sZWOvuaeEQYwmNNZwsBO6eDqSb2eV5/AwRyQczo1/rmvRpWYOPFnzPyEmr+dmrc+jcoAr3X9qMjvWrBB1RCsFZu4bcfSNQF7go/PPBSF4H1AY2Z7veEr4v18zsVjOba2ZzMzK0fF6koCWUMAZ3qMPXv+7Fnwe2Ii3jAENGzeKmV79j2VYdhlPcRbLX0IPAb4Dfhu9KAt6M4L1zalfmqQPS3V92947u3jElJSUvbyEiESiZWIKfdK3P1Pt78Zt+zZm/aTeXPzOd4W/PJy1jf9DxJEoi+WY/CLiS8JRRd98KlI/gdVsItSROqANszW1AESl8ZUomcnuvRky9vzfDezfm65Xp9Bk5lQfeX8zW3YeCjicFLJJCcNRDI8oOEN6FNBJzgCZm1sDMSgLXATrrWKQIqVg6iXsvbcaU+3pz4/n1+GD+9/R64hse/nQ5O/YfCTqeFJBIZg3dCzQhNPvnb8DNwNvu/sxZ39zsMuApIAH4p7v/xcyGAbj7KDOrCcwFKgBZwH7+s9tpjjRrSCQ4W3Yd5OlJa3h//hZKJyVwc/cG/KJnQyqUSgo6mpxFvraYCL9BH6Bv+HKCu08qwHy5okIgEry16fsZOXE145f8QKUySdx+YSN+ekF9SiUlBB1NTiNPhcDM9vGfwd1TB34PA+uA/3H3rwoqaCRUCERix9Lv9/D4hFVMWZ1BjQrJ3HlRE67tVJckrVKOOfluEeTwhglAa+Ctwt5+QoVAJPbMTtvB4xNWMXfjLlKrlOFXfZpyRdtzSJ1+OgEAAAwPSURBVCihRWmxIr9bTPwXdz/u7ouAZ/OVTESKhS4Nq/LesK68+rNOlE1O5O53FnLZ09P4ctm/ycuXTSlceWoRBEktApHYlpXljF/yA09OXM367Qc4r24l7r+0GRc0rhZ0tLhW4C0CEZHTKVHCuKLtOUy8pyePDj6XbXsPc/3o2fx49GwWbt4ddDzJgVoEIhJVh48d563Zm3h+8lp2HjhK35Y1uPfSZjStEcm6VCkoBT5YHCQVApGiaf+RTP45fT2vTE1j/9FMBp1Xm+u7pHJe3Uo6C6EQqBCISMzYdeAoo6au47UZGziSmUXF0kn0bJpC72YpXNg0harlkoOOWCypEIhIzNl7+BjT12xn8sp0Jq/KYPv+I5hBmzqV6N0shd7NqnNu7YqU0BTUAqFCICIxLSvLWbZ1L5NXpTN5VToLN+/GHaqVK8mFTavTu3kKPRqnULGMtrLIKxUCESlSdh44ytTVGUxelc6U1RnsPniMhBJGh9TK9Goeai00r1lep6jlggqBiBRZx7OchZt38024tbD0+9CelDUrlKJ38xR6NatO98bVKJscyYGL8UuFQESKjW17DzNlVai1MG3NdvYfySQpwejSoCq9mqXQu3l1GlYrq9bCKVQIRKRYOnY8i7kbdp1sLazeFjpFLbVKGXo3S6FX8+p0bVhVu6KiQiAicWLLroN8syqDySvTmbFuO4ePZZGcWIILGlWld/Pq9G5WnbpVygQdMxAqBCISdw4fO87s9TvD01PT2bjjIACNUspyUbgodKxfhZKJ8bGYTYVAROLe+u0HThaF2Wk7OXo8i7IlE+jepBq9m1WnV7Pq1KxYKuiYUaNCICKSzYEjmcxatyO0bmFlOlv3HAagRa0K9G6WwsUtqnNe3crF6jwFFQIRkdNwd1Zv23+yKMzduIvjWU61ciW5qHl1LmlRg+5NqlGmZNGenqpCICISoT2HjjFldQaTlm9j8qp09h3OJDmxBN0aV+OSFjW4uEV1alQoel1IKgQiInlw7HgWc9bvZOKKbUxasY3NOw8B0LZORS5pUYNLWtYoMiucVQhERPLpRBfSpBXbmLh828lDdmpXKk2fljW4pEUNOjeI3VlIKgQiIgUsfd9hJq9MZ+LydKavzeDwsSzKJydyYbMU+rSsQa+m1WNqkzwVAhGRKDp09Dgz1m5n0optTFqRzvb9R0goYXSuX4VLWtagT4sapFYNdiGbCoGISCHJynIWbdkdKgrL01m1bR8ATWuUOzmucF6dSoV+zoIKgYhIQDbtOBhuKWxj9vqd4ampyVzcvDqXtKxB98bVKF0y+nshqRCIiMSAPQeP8c3qdCatSOeblensOxKamtqjSWhq6kUtqlO9fHSmpp6pEER1hYSZ9QOeBhKA0e7+yCmPW/jxy4CDwM/cfX40M4mIBKVimSQGnlebgefV5mhmFnM27GTi8m0nxxYAzqtb6eQspKY1yhXK1NSotQjMLAFYDfQBtgBzgKHuvjzbcy4D7iRUCLoAT7t7lzO9r1oEIlLcuDurtu1j0vJtTFyRzqLw1NS6VUpzSYvQYHOnBlVISsj71NSgWgSdgbXunhYOMRYYCCzP9pyBwL88VI2+NbNKZlbL3X+IYi4RkZhiZjSvWYHmNSsw/KImpO89zFcr05m0fBtvz97EqzM2UL5UIndd3IRf9GhY4J8fzUJQG9ic7XoLoW/9Z3tObeD/FAIzuxW4FSA1NbXAg4qIxJLqFUoxtHMqQzuncvBoJtPXhKamRmtri2gWgpw6tk7th4rkObj7y8DLEOoayn80EZGioUzJRPq2qknfVjWj9hnRXAu9Baib7boOsDUPzxERkSiKZiGYAzQxswZmVhK4DvjklOd8AvzEQs4H9mh8QESkcEWta8jdM81sODCB0PTRf7r7MjMbFn58FPAZoRlDawlNH70pWnlERCRnUV1H4O6fEfpln/2+Udl+duCOaGYQEZEzi839UkVEpNCoEIiIxDkVAhGROKdCICIS54rc7qNmlgFszOPLqwHbCzBOQYnVXBC72ZQrd5Qrd4pjrnrunpLTA0WuEOSHmc093aZLQYrVXBC72ZQrd5Qrd+Itl7qGRETinAqBiEici7dC8HLQAU4jVnNB7GZTrtxRrtyJq1xxNUYgIiL/Ld5aBCIicgoVAhGROBc3hcDM+pnZKjNba2YPBJ0HwMz+aWbpZrY06CzZmVldM5tsZivMbJmZ3RV0JgAzK2Vm35nZonCuh4LOlJ2ZJZjZAjP7NOgsJ5jZBjNbYmYLzSxmDvsOH0s7zsxWhv+edY2BTM3C/59O3Paa2d1B5wIws3vCf+eXmtkYMyvQo8riYozAzBKA1UAfQofhzAGGuvvyM74w+rl6AvsJndvcOsgs2ZlZLaCWu883s/LAPOCqGPj/ZUBZd99vZknAdOAud/82yFwnmNmvgI5ABXcfEHQeCBUCoKO7x9TiKDN7HZjm7qPD55WUcffdQec6Ifw743ugi7vndQFrQWWpTejvekt3P2Rm7wKfuftrBfUZ8dIi6Aysdfc0dz8KjAUGBpwJd58K7Aw6x6nc/Qd3nx/+eR+wgtBZ0oHykP3hy6TwLSa+yZhZHeByYHTQWWKdmVUAegL/AHD3o7FUBMIuBtYFXQSySQRKm1kiUIYCPskxXgpBbWBztustxMAvtqLAzOoD7YDZwSYJCXe/LATSgYnuHhO5gKeA+4GsoIOcwoEvzWyemd0adJiwhkAG8Gq4K220mZUNOtQprgPGBB0CwN2/B54ANgE/EDrJ8cuC/Ix4KQSWw30x8U0ylplZOeB94G533xt0HgB3P+7u5xE637qzmQXepWZmA4B0d58XdJYcdHP39kB/4I5wd2TQEoH2wIvu3g44AMTEuB1AuKvqSuC9oLMAmFllQj0YDYBzgLJm9uOC/Ix4KQRbgLrZrutQwE2r4ibcB/8+8Ja7fxB0nlOFuxK+AfoFHAWgG3BluD9+LHCRmb0ZbKQQd98a/jMd+JBQN2nQtgBbsrXmxhEqDLGiPzDf3bcFHSTsEmC9u2e4+zHgA+CCgvyAeCkEc4AmZtYgXO2vAz4JOFPMCg/K/gNY4e5PBp3nBDNLMbNK4Z9LE/oHsjLYVODuv3X3Ou5en9Dfra/dvUC/seWFmZUND/YT7nrpCwQ+Q83d/w1sNrNm4bsuBgKdiHCKocRIt1DYJuB8MysT/rd5MaFxuwIT1TOLY4W7Z5rZcGACkAD8092XBRwLMxsD9AKqmdkW4EF3/0ewqYDQN9wbgSXh/niA34XPoA5SLeD18IyOEsC77h4zUzVjUA3gw9DvDhKBt939i2AjnXQn8Fb4i1kacFPAeQAwszKEZhfeFnSWE9x9tpmNA+YDmcACCniribiYPioiIqcXL11DIiJyGioEIiJxToVARCTOqRCIiMQ5FQIRkTinQiBxzcxGZt9h0swmmNnobNd/N7M/5nbHWjN7zcyGFGRWkWhRIZB4N5PwKk0zKwFUA1ple/wCYIK7PxJANpFCoUIg8W4G/1mu34rQytt9ZlbZzJKBFkBbM3sOTn7Tf8bMZppZ2olv/RbynJktN7PxQPUTH2BmF4c3V1sSPoMi2cw6m9kH4ccHmtkhMysZPnMhrRD/+0VUCCS+hffiyTSzVEIFYRahnVa7EjpbYDFw9JSX1QK6AwOAEy2FQUAz4FzgFv7TyigFvAZc6+7nElrhezuhVaLtwq/tQagAdQK6ECM7vUr8UCEQ+U+r4EQhmJXtemYOz//I3bPCB/XUCN/XExgT3h11K/B1+P5mhDYMWx2+fh3o6e6ZwFoza0FoI7gnw+/RA5hW0P+BImeiQiDyn3GCcwl9M/+WUIvgAkJF4lRHsv2cfYvznPZryWkL9BOmEdrp8hgwiVArozswNdLgIgVBhUAk9Mt+ALAz/I1+J1CJUDGYFeF7TAWuCx+cUwvoHb5/JVDfzBqHr28EpmR7zd3ALHfPAKoCzYHAN0SU+BIXu4+KnMUSQrOF3j7lvnLuvj28e+fZfAhcFH7dasK/7N39sJndBLwXPmZwDjAq/JrZhLqWTrQAFhM64EY7QUqh0u6jIiJxTl1DIiJxToVARCTOqRCIiMQ5FQIRkTinQiAiEudUCERE4pwKgYhInPv/ZIt/dPUSpoEAAAAASUVORK5CYII=\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": { "execution": { "iopub.execute_input": "2021-05-19T06:00:03.820777Z", "iopub.status.busy": "2021-05-19T06:00:03.820060Z", "iopub.status.idle": "2021-05-19T06:00:03.821891Z", "shell.execute_reply": "2021-05-19T06:00:03.822264Z" } }, "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": { "execution": { "iopub.execute_input": "2021-05-19T06:00:03.826293Z", "iopub.status.busy": "2021-05-19T06:00:03.825752Z", "iopub.status.idle": "2021-05-19T06:00:05.200317Z", "shell.execute_reply": "2021-05-19T06:00:05.200683Z" }, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(9, 3)" ] }, "execution_count": 1, "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": { "execution": { "iopub.execute_input": "2021-05-19T06:00:05.216927Z", "iopub.status.busy": "2021-05-19T06:00:05.216371Z", "iopub.status.idle": "2021-05-19T06:00:05.325075Z", "shell.execute_reply": "2021-05-19T06:00:05.325585Z" }, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1iV5f/A8ffNRpYKuEBlOBAHuLeIiuLKtHKWOcowzW197Zt9KxuWVo7KUT9XDrRpmXvvLbhxKzhBVARF1v3740EjQznCec5h3K/rOpeec57z3B+8is+51+cWUkoURVGUosvC3AEoiqIo5qUSgaIoShGnEoGiKEoRpxKBoihKEacSgaIoShFnZe4AnpWbm5v08vIydxiKoigFyoEDB+KklO7ZvVfgEoGXlxf79+83dxiKoigFihDi4pPeU0NDiqIoRZxKBIqiKEWcSgSKoihFnK5zBEKIUGAqYAn8IKWc+Nj7Y4E+WWKpBrhLKeP1jEtRFCU7qampxMTEkJycbO5Qcs3Ozg5PT0+sra0N/oxuiUAIYQl8C4QAMcA+IcQfUsrjD6+RUk4CJmVe3xkYqZKAoijmEhMTg5OTE15eXgghzB3OM5NScvPmTWJiYvD29jb4c3oODTUAzkgpz0kpU4BwoMtTru8FLNExHkVRlKdKTk7G1dW1QCYBACEErq6uz9yj0TMReADRWZ7HZL72L0KIYkAo8MsT3h8khNgvhNgfGxtr9EAVRVEeKqhJ4KHcxK/nHEF20Typ5nVnYMeThoWklLOB2QD16tXLVd3sM7fOsObiGhytHXGwdsDRxhFH6yyPzOfFrIthIdQcuqIoRYeeiSAGKJ/luSdw5QnX9kTnYaGz8SeYFTkL+cRc9DcHa4dHCcLBxgEna6d/JQ8HawecbP79+sO/21vZF/hvFoqimN6AAQNYsWIFpUqV4ujRowD069ePLVu24OzszP3792nUqBGfffYZHh7ZDrI8Mz0TwT6gshDCG7iM9su+9+MXCSFcgCDgZR1joV2KIOT8Re6VqUFixUYkedTmrltlkmQad1PvkpSSRGJqovZISSQpNenR3xNSEriSdIWklCTupt7lftr9HNuzEBaPEso/kkaWZPHw7y62Lvi7+uPt7K2Sh6IUcf369WPo0KH07dv3H69PmjSJF198ESklU6ZMITg4mKNHj2JjY5PnNnVLBFLKNCHEUGAN2vLROVLKY0KIsMz3Z2Ze2hVYK6VM0isWAEr7Y9FqPI7nNuO4fwHs+QEsbaB8Q/AJAp9gKBsIljn/k6RlpHEv7R6JKX8njsRULXncTbn7jz+zJpZbybeIvhv96PoH6Q/+cV8XWxcC3AMIdA8ksFQg1V2rU8y6mF7/Ioqi5EMtWrTgwoULT3xfCMHIkSP57bffWLVqFV26PG0NjmF03UcgpVwJrHzstZmPPZ8HzNMzDgBK+kCLMdoj5R5c2gXnNmuPjR9rD1tn8GoOPi21h1tlyOYbupWFFc42zjjbOOcppNT0VJJSk4i7H8eRuCNExEYQcSOCrTFbAbAUllQtWfVRYgh0D6SMQxnVa1AUE/jwz2Mcv5Jg1Hv6l3Pmf52rG+VederU4eTJk/k/EeRbNsWgUmvtAZAUB+e3ZiaGTRD1l/a6U7m/k4JPEDiVMWoY1pbWFLcsTnG74lQqUYmulbsCcOfBHSJjI4m4EUFkbCS/nfmNxScXA1CqWKl/JAa/kn5YWxq+cURRlMLBmOfNF81E8DgHN6jRTXsAxJ//u7dwahVEar+Eca/2d2Lwagq2TrqE42LrQgvPFrTwbAFoQ1Gnbp0i4kYEEbERRN6IZO3FtQDYWtpS3bX6o8QQUCqAknYldYlLUYoSY31z18uhQ4do3bq1Ue6lEkF2Snprj3r9ISMDrh/5OzEcmAt7ZoCFFXjU+zsxeNYDnb6ZW1lY4e/qj7+rP72rafPt15Oua72GzMSw4PgC5mTMAaCic0VtriEzOfgW91VLYhWlkJBSMn36dK5evUpoaKhR7imM2b0whXr16kmznkeQmgwxe/9ODFcOgcwAawetl+DTUnuU8s92fkEvyWnJHL95/NE8Q2RsJPHJ2rYMJ2snarnXIqCUNhFd060mjjaOJotNUQqKEydOUK1aNbPG0KtXLzZv3kxcXBylS5fmww8/ZNu2bY+Wj967d+/R8lFPT89s75HdzyGEOCClrJfd9SoR5NX9W3Bh+9+J4eYZ7XWHUpmrkVqCdxAUL//ke+hASkn03ehHiSEiNoIzt84gkVgICyoXr0xgqcBHPQdPR081Ca0UefkhERiDSgTmdicGzm35OzEk3dBed62UZX6hGdiXMHlod1PuciRWW50UGRtJZGwkSanaql1XO9dHQ0mBpQKp5loNW0tbk8eoKOZUVBOBmiMwNhdPqN1He0gJN078nRQilsC+H0BYaHsWfFpqj/INwdpO99CcbJxo4tGEJh5NAEjPSOfsnbOPhpIibkSw4dIGAKwtrKnuWp3Ovp3p5NNJ7WdQlEJM9QhMKS0FLh/4OzHE7AOZDpa24FEHPOtrSaF8A3AsZZYQb96/+WgSesflHZy6dQpHa0e6VOpCj6o98HYxvLStohQ0RbVHoBKBOT24Cxd3ansYovfA1UhIT9HeK+EFng20pFC+AZSqbtCuZ2OSUhIZG8mSk0tYe3EtaRlpNCrbiJ5+PQnyDMLKQnUolcJFJYIColAlgselJmvJIGavlhii90Lide09awet1/Cwx+BZH4qZbr9A3P04fj39K8uilnH93nXKOpSle9XudKvcTe1bUAoNlQgKiEKdCB4nJdy+pA0hPUwM145ow0kArpUzE0PmkJJbVbDQd79AWkYam6M3E34ynD3X9mBtYU07r3b08utFTbeaauWRUqAV1USg+vb5mRBQoqL2qPmi9lpKkrZ3IXoPRO+DqJUQsVB7z9ZF29j2cDjJox7Y5a0e0uOsLKxoU7ENbSq24ezts4SfDOePs3+w4twK/F396Vm1J+2922Nnpf/kt6IURtHR0fTt25dr165hYWHBoEGDGD58uK6lqFWPoKCTEuLP/d1jiNkH14+hnQEktI1tD3sMng3A1dfoG92SUpP48+yfhJ8M5+yds7jYutCtUje6V+2Op1P2G14UJT/KDz2Cq1evcvXqVerUqcPdu3epW7cuv//+O1988QWdOnX6RynqGTNmZFuKWvUIihohtF/urr4QmHncQ3KCtjopeq8233D0NzgwT3uvmGvm6qQGWmLwqAM2DnkKwcHagZ5+PelRtQf7r+9nycklLDi+gHnH5tHcszk9q/akqUdTVeZCUQxQtmxZypYtC4CTkxPVqlXj8uXL/7jG2KWoVSIojOycwTdYe4BWLynulNZriNmrJYhTq7X3hCWUqfF3j6F8AyheIVe9BiEE9cvUp36Z+lxLusbPp37m51M/82bMm1RwqkD3qt15vtLzuNi6GPGHVRSdrPqPNidnTGVqQvuJBl9+4cIFDh06RMOGDVm0aNG/3jdWKWqVCIoCCwso5ac96r6qvXYvHmL2/71C6dAi2Dtbe8+xtJYQKjTWehm52AVdxqEMQ2sP5Y1ab7Du4jrCo8KZvH8y3xz6ho4+Henp1xO/kn5G/CEVpXBJTEzkhRdeYMqUKTg7Zz/XZ6yhfZUIiqpiJaFKW+0BkJ4GN47/3WOI3gsn/oQtX2iH+dR/PVe7n60treng04EOPh04GX+S8JPh/HXuL345/QuB7oH09OtJ24pt1ZkKSv7zDN/cjS01NZUXXniBPn360K1btydeZ6xS1GrQVtFYWkHZWlD/Neg2G4ZHQNgObRXS2vfgm/pweJk2zJRLfiX9+KDJB6x/aT1j640lPjme/2z7DyE/hzD90HSuJV0z4g+kKAWTlJKBAwdSrVo1Ro0a9cRrpk2bZrRS1GrVkJKzs5tg3XhtvLRsAIRM0Cqr5lGGzGDXlV2EnwxnS8wWLIQFweWD6eXXi/pl6qs9CYrJ5YdVQ9u3b6d58+bUrFkTi8x9QZ9++inLli0zuBR1vtpQJoQIBaaiHV7/g5TyX30tIURLYApgDcRJKZ/6G0YlAjPJyIAjP8HGCXAnGiqFQMhHUNrfKLePuRvDslPL+O30b9x+cBtfF196+PXgOd/ncLDO26omRTFUfkgExpBvEoEQwhI4BYQAMcA+oJeU8niWa4oDO4FQKeUlIUQpKeWNp91XJQIzS03WJpW3TdZqJQX2huD/gnM5o9w+OS2ZNRfWsOTkEo7dPEYxq2I85/scPf164lvc1yhtKMqTFNVEoOccQQPgjJTynJQyBQgHHl/j1Bv4VUp5CSCnJKDkA9Z20HQYDIuARm9q8wbT6sCGCdr+hTyys7KjS6UuhHcKZ3GHxbSp2IZfT//K88ufZ+Cagay7uI60jDQj/CCKojykZyLwAKKzPI/JfC2rKkAJIcRmIcQBIURfHeNRjKlYSWj3CQzdB9U6aT2EaYGwZ7ZWbtsIarrX5JNmn7DupXUMrzOc6LvRjNo8itBfQtkWs80obSiKom8iyG6m7/FxKCugLtARaAeMF0JU+deNhBgkhNgvhNgfGxtr/EiV3CvhBS/8AK9v0spZrBoL3zWE48u18hdGUNKuJK/VfI1V3VYxLXgaLrYuDNkwhBkRM8iQuV/FpCiKRs9EEANkPajXE7iSzTWrpZRJUso4YCsQ8PiNpJSzpZT1pJT13N3ddQtYyQOPOvDqn9D7J+2gnWV94f9C4NJuozVhaWFJcIVgFnZYSGffznwX+R1DNwzlzoM7RmtDUYoiPRPBPqCyEMJbCGED9AT+eOya5UBzIYSVEKIY0BA4oWNMip6E0DaoDd4Bz32jnd88px2E94G400Zrxt7Kno+bfsx7Dd9j19Vd9FzRk5PxJ412f0UpanRLBFLKNGAosAbtl/syKeUxIUSYECIs85oTwGrgMLAXbYnpUb1iUkzEwhLqvAJvHYBW78G5LfBtQ1gxChKNsx5ACEEPvx7MC51HSkYKL698mT/P/mmUeyuKOSUnJ9OgQQMCAgKoXr06//vf/wDo168f3t7eBAQEUKVKFfr27fuvYnS5pevOYinlSillFSmlr5Tyk8zXZkopZ2a5ZpKU0l9KWUNKOUXPeBQTs3GAFmNh2CGoPxAOzodptWHz59q5CkYQ4B7Ask7LqOVei3e3v8vHuz8mNT3VKPdWFHOwtbVl48aNREZGEhERwerVq9m9WxtinTRpEpGRkURFRVG7dm2Cg4NJScn74owcE4EQoooQYoMQ4mjm81pCiPfy3LJSdDi6Q4dJMGQv+LaCzZ9qCeHAPK3GUR652rsyO2Q2/ar3Y2nUUvqt6cf1pOt5j1tRzEAIgaOjI6DVHEpNTf3XLvuHZajLlCnDqlWr8tymIUXnvgfGArMApJSHhRCLgY/z3LpStLj6Qo8ftYJ2a9+DP4fDru8g5EOoEpqnA3OsLKwYXW80Nd1qMn7HeLqv6M7koMnUL1PfiD+AUpR8vvdzo889+ZX0450G7+R4XXp6OnXr1uXMmTMMGTKEhg0bMmPGjH9dZ6wy1IYMDRWTUu597DW1o0fJvfINYMAa6LFIO395SU+Y10k7TCeP2nq1ZUnHJTjbOPP62teZf2y+0Ur1KoqpWFpaEhERQUxMDHv37uXo0eynTk1ZhjpOCOFL5h4AIcSLwFWjtK4UXUJoG9GqtNPmDjZPhO9bQfVu0Pp9KOmd61v7FPdhScclvL/zfSbvn8zh2MN81PQjVbNIeSaGfHPXW/HixWnZsiWrV6/O9n1TlqEegjYs5CeEuAyMAAbnuWVFAbC01kpfDzsEQe9oJ6d9Ux9Wj9MOz8klRxtHvgz6klF1R7H+0np6/9Wbc3fOGTFwRdFHbGwst2/fBuD+/fusX78eP79/HuJk7DLUOSaCzFpBbQB3wE9K2UxKeSHPLStKVrZOEPyulhACe8OemTA1ELZ/Dan3c3VLIQT9a/Rndshsbj+4Te+/erP+4nojB64oxnX16lWCg4OpVasW9evXJyQkhE6dOgEwduzYR8tH9+3bx6ZNm/51cH1u5Fh9VAjxKfCFlPJ25vMSwGgppVlWDqnqo0XEjZOw/gM4tQqcPbX9CLW6a3sUcuFa0jVGbR7Fkbgj9K/Rn2G1h2FloQ7oU/5JVR99svYPkwCAlPIW0CFPUSpKTkr5Qe9w6PcXOJaC38NgVhCc2ZCr25VxKMO80Hl0r9KduUfnErYujPjk3A89KUphYkgisBRC2D58IoSwB2yfcr2iGI9XM3h9I7w4Bx4kwMJu8GNXuPvs+wRsLG0Y33g8E5pOICI2gu5/dudI7BEdglaUgsWQRLAQ2CCEGCiEGACsA+brG5aiZCEE1HhBK3kdOlErZLekJ6Tcy9Xtnq/0PD+2/xErCyteXf0qP536SS0xVR4p6P8t5CZ+QyaLvwA+AaoB1YEJma8pimlZ2UKjwVrv4MohbbgoI3dlqKu5VmNpp6U0KNuAj3Z9xPs73yc5LdnIASsFjZ2dHTdv3iywyUBKyc2bN7Gzs3umz6nD65WCaec3sPa/0Hy0tu8gl9Iz0pl5eCYzI2dSrWQ1vg7+Gg/Hx89PUoqK1NRUYmJiSE4uuF8K7Ozs8PT0xNra+h+v5+nMYiFEN+BzoBTaYTMCkFJKZ6NE/YxymwjS0jPYFBVLiH9pHaJSTE5KWDFCq1f0/AxtyWkebInewrht47CwsODz5p/T1KOpceJUlHwir6uGvgCek1K6SCmdpZRO5koCefHTgRheX7Cfz1aeICOjYPWClGwIAR0mg09L+GMYXNiRp9sFlQ8ivFM4pYuVZvD6wcyKnKVOP1OKDEMSwfXMcwMKtO71yvNKo4rM2nqOEUsjeJCWbu6QlLyytIaX5mvlKJb2gZtn83S7Cs4VWNhhIR18OvBNxDcM2ziMhJQEIwWrKPmXIYlgvxBiqRCilxCi28OH7pEZmaWF4KMu1Xk7tCp/RF6h/9x9JCSruvUFnn1x6L0UELC4O9y/lbfbWdnzWbPPGNdgHDsu76Dnip5ExUcZJ1ZFyacMSQTOwD2gLdA589FJz6D0IoTgzZaV+Kp7AHvPx9N95i6u3Sm4k0JKppI+0HMR3LqonZWcx4NphBD0rtabuaFzSU5L5uWVL7Pi3AojBaso+U+RXTW07XQsYT8ewMXemvkDGlC5tJMRolPMKmKJtqS0Tl/oPC1P5xs8FHc/jjFbxnDg+gF6+/VmTL0xWFta5/xBRcln8jRZXFhPKGte2Z2lbzQmNUPywoyd7D2vyg0UeIG9oPkYOLgAdk43yi3d7N34vu339PXvy+KTixmwZgA37hnn3GVFyS8MGRr6HhgHpIJ2QhnQU8+gTKWGhwu/Dm6Cm5MtL//fHlYeUccsFHjB/wX/52Hd+3DyL6Pc0trCmrH1xzKpxSSibkXR/c/u7L+m9rIohYeuJ5QJIUKFEFFCiDNCiP9k835LIcQdIURE5iP3O4NyqXzJYvwS1oSaHi4MWXyQuTvOmzoExZgsLKDrTPCoA7+8BlcjjXbrUO9QFndYjJONE6+tfY0FxxYU2B2oipKVIYkgVyeUCSEsgW+B9oA/0EsI4Z/NpduklIGZj48MD914SjjYsOi1hoRUK82Hfx5Xew0KOmt76LkEirnC4p6QcMVot65UohKLOy4myDOISfsn8fbWt7mXmruaR4qSX+h5QlkD4EzmwTYpQDiQtxOWdWRnbcmMl+uqvQaFhVNp6BWuVSxd0hNSkox3axsnpgRPYUSdEay9uJY+K/tw4c4Fo91fUUxNzxPKPIDoLM9jMl97XGMhRKQQYpUQonp2NxJCDBJC7BdC7I+NjTWg6dxRew0KmTI14MW5cO0I/Doo1wXqsiOEYGDNgcxsM5Ob92/S86+ebLiUu7MSFMXcDFk1NEoIMQp4A3g98/lAIURgTh/N5rXHx1sOAhWllAHAdOD37G4kpZwtpawnpazn7u6eU8h5ovYaFDJV2kK7z+DkCtjwgdFv37hcY5Z2WoqXsxcjNo3ghyM/qHkDpcAxZGioHhCG9m3eAxgEtAS+F0K8/ZTPxQDlszz3BP4xWCulTJBSJmb+fSVgLYRwMzh6HXWr48nc/vWJjr9Ht+92cPr6XXOHpORWwzeg/muwY6q2tNTIyjqWZX77+bT3bs/Ug1P56sBXKhkoBYohicAVqCOlHC2lHI2WGNyBFkC/p3xuH1BZCOEthLBBW3L6R9YLhBBlhNB2/QghGmTGc/OZfwqdqL0GhYQQEPo5+LaGFSPh/FajN2FracvE5hPpWbUn847N4/2d75OWYdDiOkUxO0MSQQUgJcvzVLThnPvAgyd9SEqZBgwF1gAngGVSymNCiDAhRFjmZS8CR4UQkcA0oKfMZ1+l1F6DQsLSCl6aC66VYOkrEHfG6E1YCAvebfgugwMG8/uZ3xm1eRQP0p/4v4ii5BuGnEcwHugKLM98qTPaN/svgdlSyj66RvgYcx1McysphdcW7OfgpVu838mf/k29TR6DYgS3LsD3rcHOGV7bAMVK6tLMohOLmLh3IvXL1Gda8DQcbRx1aUdRDJXrEhOZwzbzgNeB28AdIExK+ZGUMsnUScCc1F6DQqKEF/RaAncuaz2DtJQcP5Ibfar14bPmn3Ho+iEGrBnAzfv5ZsRTUf7lqYkgc5jmdynlASnlVCnlFCllkd1br/YaFBLlG8Dz38HF7dopZzqNRnby6cTUVlM5f+c8/Vb340qi8Ta2KYoxGTJHsFsIUV/3SAoItdegkKj5IrQcBxGLYPvXujXTwrMFs9vO5mbyTV5Z9Qpnb+ft8BxF0YMhiSAYLRmcFUIcFkIcEUIc1juw/EztNSgkgt6Bmi/Bhg/h+PKcr8+l2qVqM7fdXDJkBq+ufpUjsUd0a0tRcsOQRNAe8AFa8fehNJ31DKqgUHsNCjgh4LlvwLMB/PoGXD6oW1NVS1ZlQegCnKydGLh2ILuu7NKtLUV5VoaUmLiItjGsVebf7xnyuaJC7TUo4KztoOdicHTXahLdidGtqfLO5VnQfgGeTp68ueFN1l5Yq1tbivIsDCkx8T/gHbQzCQCsgYV6BlXQqL0GBZyjO/ReBqn3tWqlDxJ1a8q9mDtz282lpltNxmwZw0+nftKtLUUxlCHf7LsCzwFJAFLKK4A61/Ex6lyDAq5UNW3D2Y3j8MtAyNBvNZiLrQuzQmbR1KMpH+36SNUnUszOkESQkrmM9OF5BA76hlRwqb0GBVylNtD+czi1GtaO17Upeyt7prWaRgfvDkw9OJXJ+yeTIY1XHVVRnoWVAdcsE0LMAooLIV4HBqAdX6lk4+Fegw/+OMasree4eieZSS/VwtbK0tyhKYZo8DrcPAO7vwW3SlBvgG5NWVtY81nzz3CxdWHB8QXcfnCbD5t8iJWFIf9bKorx5PhfnJRyshAiBEgAqgLvSynX6R5ZAfZwr0HZ4nZ8sTqKuMQHzHylLs521uYOTTFEu08h/hz8NUbbiezbSremLIQF4xqMo4RdCb6L+I6ElAQmtZiEnZWdbm0qyuMMmSweCZyQUo6VUo5RScAwaq9BAWZhCS/OAXc/WNYPYqN0bU4IweCAwbzb8F22RG9h8PrB3E1RS5EV0zFkjsAZWCOE2CaEGCKEKK13UIWJ2mtQQNk6Qe+lYGULi16CpDjdm+zl14uJzScScSOCgWsGqvpEiskYso/gQylldbSzi8sBW4QQ63WPrBBRew0KqOLltXOPE69DeB9I07+kdAefDkxrNY3zd87z6upXuZx4Wfc2FeVZNobdAK6hHRxTSp9wCi+116CA8qwLXWdC9G744y3dCtRl1dyzOd+3/Z745Hj6rurLmVvGPztBUbIyZI5gsBBiM7ABcANel1LW0juwwkjtNSigqneFVu/B4aWwdZJJmgwsFci80HlIKem3ph+HY4t0eS9FZ4b0CCoCI6SU1aWU/5NSHtc7qMJM7TUooJqPgYBesOkTOPqLSZqsUqIK89vPx9nGmdfWvsbOyztN0q5S9DwxEQghnDP/+gVwSQhRMuvDNOEVTtmda3A/RZ1rkK8JAZ2nQoUm8NtgiN5nkmbLO2n1iSo4VWDIxiGsvrDaJO0qRcvTegSLM/88AOzP/PNAludKHjx+rkHHads4dOmWucNSnsbKFnosBOdyEN4Lbl00SbNu9m7MCZ1DLbdavL3lbZZFLTNJu0rR8cREIKXslPmnt5TSJ/PPhw8fQ24uhAgVQkQJIc4IIf7zlOvqCyHShRAvPvuPUHA93Guw6LWGJKem88KMnUxeE0VKmio1kG85uGoF6tJTtGqlyQkmadbZxpmZITNp7tmcCbsn8P3h71V9IsVonnh4vRCiztM+KKV8avF2IYQlcAoIAWKAfUCvx+cYMq9bByQDc6SUPz/tvuY6vF5vCcmpfPTncX4+EIN/WWe+7hFI1TKqtl++dW4zLHwBfIK1JaaWpikLkZqRyvs73mfFuRW84v8KY+qNwUKoqvBKznJ7eP2XmY9vgT3AbLQaQ3uAaQa02wA4I6U8J6VMAcKBLtlc9xbwC9ry1CLL2c6ayS8FMPuVuty4m0zn6duZueUs6WoiOX/yaQkdv4Qz62DNuyZr1trCmk+afUKfan348fiPjN8xntQMdVSqkjdPGxoKllIGAxeBOlLKelLKukBtwJCFzR5AdJbnMZmvPSKE8EArcz3zWQMvrNpWL8OaES1o5VeKiatO0mPWLi7EJZk7LCU7dftB46GwdxbsNV0dRgthwTv132FI4BD+OPsHozaNIjlNlS9Rcs+QPqWflPLRIatSyqNAoAGfE9m89vjX2ynAO1LKpy6ZEUIMEkLsF0Lsj42NNaDpgs3V0ZYZL9fh6x4BRF2/S/up2/hx90U1JpwfhXwEVTvAqrfhtOk23AshCAsI472G77ElZgth68NUfSIl1wxJBCeEED8IIVoKIYKEEN8DJwz4XAzaEZcPeQJXHrumHhAuhLgAvAh8J4R4/vEbSSlnZ/ZI6rm7uxvQdMEnhKBrbU/WjmxBPa8SjP/9KH3n7OXqnfvmDk3JysISun0PpavDT/3gSoRJm+/h14PPW3xO5I1IBqwZQNx9/WsiKYXPEyeLH10ghB0wGGiR+dJWYIaU8ql9USGEFdpkcWvgMtpkcW8p5bEnXB+7cYQAACAASURBVD8PWFFUJ4ufRkrJwj2X+PSvE1hbCj7qUoMugeUQIrtOl2IWdy7DnHbaKqKXf4byDUza/PbL2xm1eRTu9u7MbjsbD0ePnD+kFCm5nSwGQEqZLKX8WkrZNfPxdU5JIPNzacBQYA1aD2KZlPKYECJMCBH2rD9EUSaE4JVGFVk1vDmVSzsxYmkEby46yM1E/YugKQZy8YD+q7TlpQueh/NbTdp8M49mzA6Zze0Ht+m7si+nb502aftKwZZjjyC/KYo9gqzSMySzt57j63WncLa34rNutQjxV5XB842717REcOs8dP8RqrQ1afOnb53mjXVv8CD9Ad+2/pbAUoZM5ylFQZ56BEr+YmkhGNzSlz/eaoq7kx2vL9jP2J8iuZuslhDmC05loN9f4F4VwnvD8eUmbb5yicosaL+A4rbFGbRuEDsu7zBp+0rBpBJBAeVXxpnlQ5oyNLgSvxyMIXTKNnaeUROF+YKDK7z6J3jU0SaQI5aYtHlPJ0/mt59PBacKDN04lNXnVX0i5ekMKUNdRQjxvRBirRBi48OHKYJTns7GyoIx7aryy+Am2FpZ0PuHPXzwxzFVwC4/sHOBV34Dr+bwexjs+z+TNu9m78bc0LlafaKtb7Pi3AqTtq8ULIasGopE2/B1AHj0G0ZKeUDf0LJX1OcInuR+Sjqfrz7JvJ0X8HFz4MvuAdSuUMLcYSmpybCsL5xeA20/gSZDTdp8cloyQzYM4eCNg8xqM4sGZU27mknJP542R2BIIjiQuaM4X1CJ4Ol2nIlj7E+RXEtI5s2WlRjWujI2VmoE0KzSUuDX1+H479DyXQh6WytrbSIJKQn0XdmXG/du8GOHH/Et7muytpX8I6+TxX8KId4UQpRV5xHkf00rubF6ZAu61fHkm01neP7bHURdUztOzcrKBl6cA4F9YPOnsO59kxx5+ZCzjTPftfkOWytb3lz/ptp0pvyLIT2C7M5TlIaWojY21SMw3Npj13j3tyMk3E9jVNsqvN7cB0sLtQnNbDIyYNVY2PcD1H8N2k8CC9P11o7FHaP/mv54u3gzt91cilkXM1nbivnldUOZdzYPsyQB5dmoAnb5jIUFdJgMTYdryWD5EEhPM1nz1d2qM6nFJE7Gn+Sdre+QnqEWFSgaQ1YNWQshhgkhfs58DBVCWJsiOCXvVAG7fEYIaPMhBP8XIhfDLwO1OQQTCSofxLgG49gcs5mJeyeq/w4UAAw5TWMGYA18l/n8lczXXtMrKMW4Hhawa+Tjyts/H2b870dZe+waX7xYi7Iu9uYOr+gRQpswti4Ga/8Lqfeh+wKwtjNJ8z39ehJzN4b5x+fj6eTJq9VfNUm7Sv5l0PJRKWVATq+ZipojyBtVwC6f2T8HVowC7+bQcwnYOpqk2QyZwZgtY1h/cT1ftvySkIohJmlXMZ+8rhpKF0I8Wm8mhPAhy34CpWBRBezymXoDoOssuLAdFnaD+7dN0qyFsODTZp9Sy70W47aNI+KGactnK/mLIYlgLLBJCLFZCLEF2AiM1jcsRW9ebg4se6Mx74T6seHEDdpN2cq649fNHVbRFNADXpoHlw/Cgucg6aZJmrWzsmN6q+mULlaaYRuHcSnhkknaVfIfg6qPCiFsgapop46dlFKa7eujGhoyvpPXEhi5NJITVxN4qa4n73f2x8lOrQcwudPrYOnLUMIL+i7XCtiZwMWEi7y88mVcbF34sf2PlLBTO9ILI2NUH60L1AACgB5CiL7GCk4xP1XALp+oHAJ9fobb0TAnFG6b5ht6ReeKTG81nauJVxm+aTgP0tUwYVFjyPLRH4HJQDOgfuYj26yiFFxPKmCX9MB069wVtEnjvsvhXjzMaQ83z5qk2cBSgXza/FMO3TjEu9veJUNmmKRdJX8wZNXQCcBf5pMFx2poSH9ZC9iVc7FjfCd/QmuUUSuLTOlqJPzYFSys4JXfobS/SZqdd3QeXx74kv41+jOq7iiTtKmYRl6Hho4CphmsVPIFextLPniuOj+HNcbZ3prBiw7Sb+4+tSvZlMoGaEdfImBeR7hyyCTNvlr9VXpU7cHco3NZFrXMJG0q5mdIj2ATEAjsBR4NHkopn9M3tOypHoFppaVnsGDXRb5ad4qU9AzCgnx5s6UvdtaW5g6taIg/B/O7QPJt6PMTVGike5NpGWkM3zSc7Ze3M73VdFp4ttC9TUV/eS1DHZTd61LKLUaI7ZmpRGAeNxKS+WTlCZZHXKF8SXs+fK46rfzUWckmcScG5j8Hd69CryXg01L3Ju+l3qPf6n5cSLjAvNB5+LuaZmhK0U+eEkEeGw4FpgKWwA9SyomPvd8FmABkAGnACCnl9qfdUyUC89p5Jo7xy49yNjaJEP/S/K+zP54lVBVL3d29Dj8+r00ed18AVUN1bzL2Xix9VvYhNSOVxR0WU9axrO5tKvrJa4+gG/A5UAptH4FAK0PtnMPnLIFTQAgQA+wDekkpj2e5xhFIklJKIUQtYJmU0u9p91WJwPxS0jL4v+3nmbbhNBLJW60q81pzb2yt1HCRru7Fa7uPrx2BF36A6l11b/LMrTP0XdWX0g6lmd9+Ps42T/3fXsnH8jpZ/AXwnJTSRUrpLKV0yikJZGoAnJFSnpNSpgDhQJesF0gpE7OsRnIA8sXKJOXpbKwsGNzSl/Wjg2hZpRST1kTRfuo2tp9Wew90VayktrTUsz78PAAOLdK9yUolKvF18NdcSLjAqE2jSE1P1b1NxfQMSQTXpZQncnFvDyA6y/OYzNf+QQjRVQhxEvgLGJDdjYQQg4QQ+4UQ+2NjY3MRiqIHj+L2zHylLvP61yc9Q/Ly/+1h6OKDXLuTbO7QCi87F3j5F/AOguVvwt7vdW+yYdmGfNjkQ/Zc28MHuz5QpasLIUMSwX4hxFIhRC8hRLeHDwM+l92i83/9FySl/C1zOOh5tPmCf39IytlSynpSynru7u4GNK2YUsuqpVgzogUj2lRm7fHrtP5yMz9sO0dqutqUpAsbB+gVDlU7wMoxsGOq7k0+5/scbwa8yR9n/2Bm5Ezd21NMy5BE4AzcA9oCnTMfnQz4XAxQPstzT+DKky6WUm4FfIUQbgbcW8ln7KwtGdGmCutGtqCBd0k+/usEnaZtZ+/5eHOHVjhZ22mTxtW7aWcgb/pU93OQwwLC6OLbhe8iv2P5meW6tqWYlm6rhoQQVmiTxa2By2iTxb2llMeyXFMJOJs5WVwH+BPwfNouZjVZnP9JKVl7/Dof/Xmcy7fv062OB+PaV8PdydbcoRU+Genw5zA4tBAaD4W2H2sH3+gkNT2VwRsGc+DaAWaEzKBRWf33NSjGkddVQ3bAQKA68OgIJSlltuP5j322AzAFbfnoHCnlJ0KIsMzPzxRCvAP0BVKB+8BYtXy08LiXksY3G8/w/bZz2Flb8na7qvRuWBFLC1WqwqgyMmD1f2DvLKjbHzp+pZ2PrJO7KXfpu6ov15Ous6D9AiqVqKRbW4rx5DUR/AScBHoDHwF9gBNSyuHGDtQQKhEUPGduJPL+8qPsPHuTmh4uTHi+BoHli5s7rMJFStjwIWz/Gmr1hC7fgqUhJ9HmztXEq/RZ2QcrCysWdViEezE1d5ff5XX5aCUp5Xi09f7zgY5ATWMGqBRulUo5sui1hkzrVZvrCcl0/W4H4349wq0k0x3aXugJAW0+gFbj4XA4/NwP0vT79y3rWJZvW3/L7Qe3GbJhCPdS7+nWlqI/QxLBw4XDt4UQNQAXwEu3iJRCSQjBcwHl2DA6iAFNvVm2P5pWX25m6b5LZGSo5YhG02IMtPsMTvwJ4b0h9b5uTVVzrcbkoMlE3Ypi7NaxpGWokuUFlSGJYLYQogQwHvgDOI62yUxRnpmTnTXjO/mz4q1mVCrlyDu/HOHFmTs5duWOuUMrPBq/CZ2nwpn1sOgleHBXt6ZaeLbgvw3/y9aYrUzcO1HtMSigdK01pAc1R1B4SCn55eBlPlt5glv3Uujb2ItRbavgrI7JNI7DP8Fvb4BHHa1yqb1+R1B+deAr5h6dy6i6o+hfo79u7Si597Q5ghxnkzLPK34BbTjo0fVSyo+MFaBSNAkheLGuJyHVSjN5bRTzd11gxeGrvNexGl0Cy6mDcPKq1ktgbQ8/9dMOuXnlN92SwYg6I7iaeJWvDnxFOcdytPNqp0s7ij4MGRpajlYjKA1IyvJQFKNwKWbNhOdrsHxIUzyK2zFiaQS9vt/N6ev6DWkUGdU6QY+FcP2Ylgzu39KlGQthwcfNPqZ2qdq8u+1dDt0wzUE6inEYsnz0qJSyhoniyZEaGirc0jMk4fsu8cXqKJIepDGwmTfDWlfGwVa/pZBFQtRqWPYKlK6ua8/gdvJtXl71Mnce3GFhh4VUdK6oSzvKs8vr8tGdQgi1XFQxCUsLQZ+GFdk4OohudTyYtfUcbb7awsojV9VEZF5UDYXuP+reMyhuV5wZrWcgEAxeP5j4ZFVipCB4YiIQQhwRQhwGmgEHhRBRQojDWV5XFN24OtryxYsB/DK4MS721ry56CB95+zlvDo3OfdMlAzKO5dneuvp3Lh3g2Ebh5GcpqrR5ndPHBoSQjy1TyelvKhLRDlQQ0NFzz/OTU7LICzIh7CWvhSzUcNFuWKiYaJ1F9cxevNo2lRsw+SgyVgI/cpeKDnL7dBQLHBFSnkx85e+HdANqGuuJKAUTVaWFgxo5s3G0UG0r1mGaRvPEDxZ24yWrjajPbusPYMFz+vWMwipGMLoeqNZd3EdX+3/Spc2FON4WiJYTeYO4swqobsAH2CIEOIz/UNTlH8q5WzH1J61+SmsMeWK2/POL0doP3Urm07eUPMHz6pqqLaa6MZxXZNBX/++9PLrxfzj81lycokubSh597REUEJKeTrz768CS6SUbwHtMew8AkXRRX2vkvw6uAnf9anDg7QM+s/bR+/v93AkRu1OfiZV2umeDIQQvFP/HVp6tmTi3olsjt5s9DaUvHtaIsj6FasVsA4g8/xhdfSUYlZCCDrULMu6kUF80NmfqOt36fzNdoaHHyI6XhVAM5gJkoGlhSWft/icaiWr8fbWtzkWdyznDykm9bTJ4oXANbRDZf4DeEsp7wkhigNbpJQBpgvzb2qyWMlOQnIqs7ac5Ydt55ESXm1SkSHBlShezMbcoRUMp9bA0pehlD/0/V2XCeS4+3H0+asPD9IfsKjjIjwc/3WEuaKj3E4Wvw7Eoc0TtJVSPvya5Q9MNmqEipJHznbWjG3nx+axLekSWI4ftp8naNJmvt96juTUdHOHl/+ZoGfgZu/GjDYzSMlI4c31b5KQkmD0NpTceaaic0KIOlLKgzrGkyPVI1AMceJqAhNXnWTLqVg8itvzdmhVOtcqh4U6He3pTNAz2HdtH4PWDaJOqTrMbDMTa0tVZNAU8rqzOKsfjBCPouiuWlln5g9owMKBDXGxt2Z4eARdvt3BzrNx5g4tf6vSDnosyuwZdNGlZ1C/TH0+avIRe6/tZfD6wdxOvm30NpRn86yJQH2dUgqUZpXdWPFWM77uEUB8Ugq9v99D/7l7ibqmCto9UZW2mcnghG7JoLNvZz5u+jEHbxyk11+9OH3rdM4fUnTzrIngQ12iUBQdWVgIutb2ZMPoIMa192P/xVu0n7qVd34+zLU7qvxBtkyQDLpU6sLc0Lkkpyfz8sqX2XBpg9HbUAxjUCIQQngIIZoA8UKIFkKIFgZ+LjSzRtEZIcR/snm/T2b9osNCiJ1CCLOsRFKKBjtrS94I8mXr2GD6N/Xm10MxtJy8iclroribnJrzDYoaEySDAPcAwjuG4+Piw4hNI5gZOVNtDjQDQ8pQfw70QDui8uHyCymlfC6Hz1kCp4AQIAbYB/SSUh7Pck0T4ISU8pYQoj3wgZSy4dPuqyaLFWOJjr/HF2ui+DPyCq4ONoxoU5meDSpgbalq4vzDqbWwtA+UqgZ9l+sygZyclsyHuz5kxbkVhFQM4eOmH1PMupjR2ynKnjZZbEgiiAJqSSkfPGOjjdF+sbfLfD4OQEqZbXmKzHORj0opn7q4WCUCxdgio2/z6coT7Dkfj7ebA++EVqVd9TLqhLSssiaDV36HYiWN3oSUkvnH5vP1wa+pVLwS01pNU3sNjCivq4bOAblZ3+UBRGd5HpP52pMMBFZl94YQYpAQYr8QYn9sbGwuQlGUJwsoX5zwQY34v1frYWkhCFt4kBdn7uLARVVL/5Gsw0Q/Pg/3jP9vI4SgX41+fNv6W64mXqXXil7su7bP6O0o/2ZIIrgHRAghZgkhpj18GPC57L5OZdv9EEIEoyWCd7J7X0o5W0pZT0pZz93d3YCmFeXZCCFoXa00q4c357NuNbkUf48XZuxi8MID6gyEh6q0hZ6LdU0GAM08mrG442JcbF0YtHYQy6KW6dKO8jdDEsEfwARgJ3AgyyMnMUD5LM89gSuPXySEqIW2P6GLlPKmAfdVFN1YWVrQq0EFtoxtycg2VdhyKpaQr7bw/vKjxCU+0+ho4VQ5xCTJwMvFi8UdF9O4XGMm7J7AhF0TSE1XE/p6MWhnsRDCHqggpYwy+MZCWKFNFrdGq1e0D+gtpTyW5ZoKwEagr5RypyH3VXMEiinF3n3A1A2nWLI3GntrS8KCfBjYzAd7G0tzh2Zep9dBeG9d5wwA0jPSmXZoGnOOzqFu6bp81fIrStrp01Zhl6c5AiFEZyAC7XwChBCBQog/cvqclDINGAqsAU4Ay6SUx4QQYUKIsMzL3gdcge+EEBFCCPUbXslX3J1s+fj5mqwZ0YImvq5MXnuKlpM3sWxfdNE+FMdEPQNLC0tG1h3JxOYTORp3lJ4renIy/qQubRVlhqwaOoBWhnqzlLJ25mtHpJRmOdBe9QgUc9p3IZ5PV57g0KXbVC3txH86+NGyinvRXWFkop4BwLG4YwzbNIy7KXeZ0HQC7bza6dZWYZTXVUNpUsrHT/wowl+FlKIs66E4yWnp9J+7jz4/7OHo5SJ6KE7WnsGCLrr1DACqu1VnaaelVClRhTFbxvDNoW/IkOpoFGMwJBEcFUL0BiyFEJWFENPRJo4VpUh6/FCcE1cT6DR9OyPCD3EuNtHc4Znew2QQG6V7MnCzd2NOuzl0rdSVWYdnMWLTCJJS1aquvDJkaKgY8F+gLdqS0DXABCmlWYq0qKEhJb9JSE5l5uaz/N/286SkZxBavQxhQb4ElC9u7tBM6/R6bZjIvaq2A1nHYSIpJYtPLmbSvkl4u3gzLXga5Z3L5/zBIixPO4sfu5El4CClNNuJEioRKPlV7N0HzNt5nh93XSQhOY0mvq6EBfnSvLJb0ZlDMGEyANh9dTejN49GCMHkoMk0KttI1/YKsryuGloshHAWQjgAx4AoIcRYYwepKAWdu5MtY9v5seM/rXi3gx9nYxPpO2cvnaZv58/IK6SlF4Hx7MptTDZMBNCobCPCO4bjbu9O2LowFp1YpIrW5YIhQ0MRUspAIUQfoC7a7t8DUspapgjwcapHoBQUD9LS+f3QZWZtPce52CQqlCzGoBY+vFjXEzvrQr4PwcQ9g6TUJMZtG8em6E10rdSV9xq9h42lOq86q7yuGrIWQlgDzwPLpZSpqFVDipIjWytLetSvwPqRQcx8uS4lHGx47/ejNPt8I99uOsOd+4V4p6yJewYO1g5MCZ7CG7Xe4LczvzFwzUDi7qvT6AxlSI9gGFovIBLoCFQAFkopm+sf3r+pHoFSUEkp2X0unhlbzrL1VCyOtlb0bliBgc28Ke1sZ+7w9GHingHAmgtrGL9jPM42zkwNnkp1t+q6t1kQGG2yOMsNrTJ3DpucSgRKYXDsyh1mbTnHisNXsLKwoGttDwYF+eDr7mju0IzPDMngZPxJhm0cRnxyPB82+ZCOPh11bzO/y+t5BLbAC4AXYPXwdSnlR0aM0WAqESiFyaWb9/h+2zmW7Y8mJT2Ddv5lCGvpS2BhW3pqhmQQnxzPqM2jOHD9AANqDGBY7WFYWhTyuZmnyGsiWA3cQas4+vCEMqSUXxozSEOpRKAURnGJD5i34wILdl0gITmNxj6uhLX0pUVhWnpqhmSQmp7KxL0TWXZqGc09mvN5i89xsnHSvd38KK+J4KiUsoYukeWCSgRKYZb4II0ley7xf9vPcy0hGf+yzoS19KVDjTJYFYYjNM2QDACWRS3jsz2f4enkyfRW0/Fy8TJJu/lJXhPBbGC6lPKIHsE9K5UIlKIgJS2D3yMuM3PL2UdLT19v4cNLhWHp6aNkUAV6/wTOZU3S7L5r+xi9eTRpGWl8EfQFzTyamaTd/CKvieA4UAk4DzxAKzMh1T4CRdFfRoZk3YnrzNh8lojo27g52tCviRevNPLCpVhuTpDNJ86shyW9QQio/xo0HQGO+p8+eDnxMsM3Duf07dOMrDOSV6u/WniG3nKQ10RQMbvXpZQXjRDbM1OJQCmKpJTsOR/PjM1n2XIqFgcby8ylpz6UcSmgS0/jz8OWL+BwOFjZQ8NB0GSY7sNF91Lv8d6O91h3cR2dfDrxQZMPsLW01bXN/CDPy0eFEM2AylLKuUIId8BRSnneyHEaRCUCpag7fiWBWVvP8mfkFSwthLb0tIUvlUoV0KWncWdgy0Q48jPYOEKjwdB4CNjrt3JKSsmsw7P4NuJbarjWYErwFEo7lNatvfwgrz2C/wH1gKpSyipCiHLAT1LKpsYPNWcqESiKJjpeW3q6dJ+29LStf2nCgnypXaGEuUPLnRsnYPNEOP472LpAk6HQMAzsnHVrcuOljYzbNg4Hawe+Dv6aAPcA3doyt7wmggigNnAwywllh9UcgaLkD3GJD5i/8wILdl3kzv1UGvmUJCzIl6CCenLatSNaQji5AuxLaMNFDQaBrT49ntO3TjNs4zCu37vO/xr/jy6VuujSjrnlNRHslVI2EEIclFLWyaxCukslAkXJXxIfpBG+9xI/bNOWnlYr60xYkA8da5YtmEtPrxyCTZ/C6bVQzA2ajYB6A8GmmNGbup18mzFbxrDn2h5e8X+FUXVHYWVhlfMHC5C8JoIxQGUgBPgMGAgsllJOM6DhUGAqYAn8IKWc+Nj7fsBcoA7wXynl5JzuqRKBojxdSloGyzOXnp6NTaJ8SXv6N/HmpXqeONkVwJVG0ftg0ydwbhM4loZmo6BuP7A27iR5WkYak/dPZtGJRTQp14QvWnyBi62LUdswJ2NMFoegnVAGsEZKud6Az1gCp9ASSAywD+glpTye5ZpSQEW0yqa3VCJQFOPJyJCsP3GdWVvPceDiLRxtrXixrif9mnjh5eZg7vCe3cWdWg/hwjZwKgctRkPtvmBl3HLTv57+lQm7J+Dp6Mm0VtPwdvE26v3NJVeJQAhxl7/LTT8+0JgMnEX7Fr/hCZ9vDHwgpWyX+XwcgJTys2yu/QBIVIlAUfQRGX2buTvO89eRq6RlSFpVLUX/pt40reRa8OYRzm3RegjRe8ClPLQYC4G9wdJ4vZ2D1w8ycvNIUtNTmRQ0iaYeZlkbY1R6VB+1BGoAi55UfkII8SIQKqV8LfP5K0BDKeXQbK79gKckAiHEIGAQQIUKFepevGiWLQyKUuDdSEhm4Z5LLN5zkbjEFKqUdqRfE2+61vbA3qYA7ViWEs5u0HoIlw9ACS8IegdqdgdL44ztX0m8wrCNwzh9+zSj6o6ir3/fgpc0s8jrwTT/IqVMl1JGAtOf1m52H81le7OllPWklPXc3fXffagohVUpZztGhVRh+zutmPxSANaWFrz72xEaT9zAxFUnuXL7vrlDNIwQUKkNvLYBei0FW2f4fTB81xAO/wQZ6TnfIwflHMuxoP0CWpVvxeT9kxm/Yzwp6SlGCD7/yVWPwKAbq6EhRcn3pJTsPR/P3B0XWHv8GkIIQquXoX9TL+pWLFFwvgFLqS033fQZ3DgG7n7QchxUew4s8rZiKkNmMDNyJjMiZxDgHsCU4Cm42bsZKXDTMfrQkIGNWqFNFrcGLqNNFveWUh7L5toPUIlAUcwqOv4eP+6+SPjeSyQkp1HTw4X+Tb3oWKsstlYFZNgoI0PbkLZ5IsRFQemaEDwOqnbQehF5sPbCWt7b8R7ONs5MazUNf1d/IwVtGmZJBJkNdwCmoC0fnSOl/EQIEQYgpZwphCgD7AecgQwgEfCXUiY86Z4qESiKvu6lpPHLwcvM23Ges7FJuDna8nKjCvRpWBF3pwJSkycjHY7+oiWE+LNQNhCC/wuVQ/KUEE7Gn+StjW9xO/k2E5pNINQr1IhB68tsiUAPKhEoimlkZEi2nYlj7o7zbI6KxcbSgk4BZRnQ1JsaHgVkfX16GhxeCls+h9sXwbM+BL8LPsG5Tghx9+MYuWkkEbERvFHrDd4MfBMLkf837KlEoChKnpyNTWT+zgv8fCCGeynpNPAqSf+mXoT4ly4Yu5bTU+HQQtg6GRJioEITLSF4N8/V7VLSU/h498f8duY3WldozafNPqWYtfF3PBuTSgSKohjFnfup/LQ/mnk7LxBz6z4exe3p27giPeqXp3gx427s0kXaAzi4QEsIidfAu4U2ZFSh0TPfSkrJwhMLmbx/MpWKV2Jaq2l4OHroELRxqESgKIpRpWfuWp674zy7z8Vjb21J1zoe9G/iReXSBeBM4NT7sH8ubP8KkmLBt7WWEDzrPvOtdlzewdgtY7GysOLr4K+pW/rZ72EKKhEoiqKb41cSmLfzPL9HXCElLYPmld0Y0NSboCruWFjk8+WnKUmw7wfYPgXux0OVUG3ZabnAZ7rNhTsXeGvjW8QkxvBew/d4ocoLOgWceyoRKIqiu5uJD1iy9xI/7r7I9YQH+Lg58GoTL16o64mjbT6v5PngLuyZBTunQ/JtqNIemo+G8vUNvkVCSgJvb3mbHVd20NuvN2Prj81XFUxVIlAUxWRS0jJYdfQqc3dcICL6Nk62VnSvX55XG3tRwTV/T6iSfAd2z4Q9M+D+LW0OClWBVQAADXdJREFUoflo8A4yaJVRWkYaXx/4mgXHF9CobCMmB03ONxVMVSJQFMUsDl26xdwdF1h55CrpUtKmWmn6N/WisU8+L3b3IBEOzNN6CInXwKOulhCqtDdop/Jvp39jwu4JlHUoy/RW0/Ep7qN/zDlQiUBRFLO6dieZhbsvsnjvJeKTUvAr40T/pl50qFk2f5+RkJoMkUtgxxS4dQHcq0HzUVC9W47F7SJuRDB803BS0lP4vMXntPBsYZqYn0AlAkVR8oXk1HT+iLjCnB3nOXntLtaWggbeJQmuWorW1UrjnV/PSUhPg2O/wravIPaEVu206XAI7ANWT95tfTXxKsM2DSMqPoqRdUfSr3o/s/WEVCJQFCVfkVJy4OIt1p24zqaTNzh1PRH+v737D7KqvO84/v6wP2GXX7vLAvJDUAhrCLYYZQUUs5CkmpDYzrQJzsRO7RTbVKnGzDi1TmozyUzS0Rh/JaUWjWYSdKLBxklt1QIqIBCEoIgsFCHICuWyu4Vld2F/fvvHeRaumwX5ce+eC+f7mtnhnnvPvfd7md37Oc9zzvM8wMSKEuZWVTK3qpKrJpRRmJ9jg9W6u2HHf8KqH0TTX5eOglmLohXTTrKmcmtHK99a8y1e2fMKX7rkS9w36z6K8vp/qg4PAudcTtvb2MrK7SmWb0uxdlcD7Z3dlBblM+cTFdRMqaSmqpKK0hya58gMdr8eBcLuN2DgcKj+OsxYCIPK+tjdWPzOYn68+cdcXnE5D9U8xIhB/TulvgeBc+680dreyZqdDayoPcCK2hQHmtqQ4PKxw5gXWgtTLxqSOyeb926IBqZtfwkKS+HKv4SZt8HgUb+366t7XuXe1fcyuHAwj9Q8wtSKqf1WpgeBc+68ZGZs3dfEytoUy2tTvF13CDMYOaSImilRKMyeVEFJLoxTOLAVVv8wmvV0QAFM/xrM/rvofEKa7Y3bWbRiEY3HGvnO7O9ww8Qb+qU8DwLn3AWhvrmN17YfZEXtAVbtqOdIWyeFeQO4+tLy462FcWUxj1VoeB/WPAybl4J1w7Q/g2u+AZVVJ3Y52sBdr93FptQmFk5byO3Tb8/6DKYeBM65C057Zzdv/a6RFbUpVtSm2FXfAsDkytLjJ5w/ffHw+GZHbdoHbz4GG38CHa1QNT8aizDmCgA6ujr47vrvsux/llEzrobvXfs9Sgqyd9WUB4Fz7oK3u76FFbUpVtamWL+7gY4uY0hxPtdNqWRu1Qiu+0QlZSUxzJDa0gDrF8Nv/jUauXzp3CgQLp6NAUtrl3L/hvuZOHQij859lLGDx2alDA8C51yiHDnWwZqd9SzflmLl9hT1ze0MEEwfP/x4a6Fq1OD+PeF8rAneehLW/ghaUjCuOgqEyZ9n7f51fPP1b5KnPB78zINcNer05zg6XR4EzrnE6u42tnx4mOWhtbDlw8MAXDS0mJqqSuZdVsmsSysoLuindZk7jkaL5Kx5BA5/ACM/BdfexZ6x01n02p3sbdrLPdX38JUpX8no23oQOOdckGo6dnzMwuqd9bS2d1GUP4DZkyqoCa2FMcMGZr+Qrg7Y8nx06Wn9Dii7hCMz/5a7D21k9b41LJiygLtn3E3BgMxMweFB4JxzfWjr7GL9rhMnnD9obAWgpDCPstJCykqKqCgppKykkLLSQipKij56u7SQ8pLCc2tNdHdD7a9h1QOw/226hozhoUmf5qn/20z1qGoeuO4BhhUPO+fPGlsQSLoeeBjIA5aY2fd7Pa7w+BeAVuAvzGzTqV7Tg8A5lw1mxvsHW3h9x0H2HTpKQ3MbDS3tNDS309gS/bR3dff53J7gKC8pojwER3lp+u3wWGm03WdwmMH7y6P5jPas4VdlI/n20IGMLBnFo/N+xKThk87p88USBJLygB3A54A6YANwk5m9l7bPF4BFREFQDTxsZtWnel0PAudcHMyMI22dNDa3h4Boo7GlPS0sPhocDS1tdHT1/f1aUphHeWnUuigPQVEWQqS8tJCJre8wafu/sSu1mjtGVnIsv5B/nvlPXDf5xrOu/1RBkM3heDOAnWa2KxTxLHAj8F7aPjcCP7UojdZJGiZptJntz2Jdzjl3xiQxpLiAIcUFTDiNWVI/GhxtNIQAaQxh0dASBcn+w8d4d99hGlvaewXHX/FJfZaFbS/wq4v2smjNvXx13XPce/PPMv7ZshkEY4C9adt1REf9H7fPGOAjQSDpVuBWgPHjx2e8UOecy7SzCY6mY52hG6qN+uZ2Glum0djyR8w8uJXBh+5n5IhLs1JrNoOgrwt0e7eTTmcfzOxx4HGIuobOvTTnnMstkhg6sIChAwv6WJdhElEHSnZkc+x1HTAubXsssO8s9nHOOZdF2QyCDcBkSRMlFQILgBd77fMi8OeKXA0c9vMDzjnXv7LWNWRmnZJuB14munz0STPbKulvwuOLgZeIrhjaSXT56C3Zqsc551zfsjqJt5m9RPRln37f4rTbBtyWzRqcc86dWo4tCOqcc66/eRA451zCeRA451zCeRA451zCnXezj0o6COw5y6dXAPUZLCdTcrUuyN3avK4z43WdmQuxrovNbERfD5x3QXAuJL11skmX4pSrdUHu1uZ1nRmv68wkrS7vGnLOuYTzIHDOuYRLWhA8HncBJ5GrdUHu1uZ1nRmv68wkqq5EnSNwzjn3+5LWInDOOdeLB4FzziVcYoJA0vWStkvaKenv464HQNKTklKS3o27lnSSxklaKWmbpK2S7oi7JgBJxZJ+I+ntUNe3464pnaQ8Sb+V9Ou4a+kh6XeStkjaLClnFvsOy9I+L6k2/J7NzIGapoT/p56fJkl3xl0XgKRvhN/5dyU9I6k4o6+fhHMEkvKAHcDniBbD2QDcZGbvnfKJ2a9rDtBMtG7zp+KsJZ2k0cBoM9skaTCwEfjjHPj/ElBiZs2SCoDVwB1mti7OunpIugu4EhhiZvPjrgeiIACuNLOcGhwl6WlglZktCeuVDDKzQ3HX1SN8Z3wIVJvZ2Q5gzVQtY4h+1z9pZkcl/QJ4ycyeytR7JKVFMAPYaWa7zKwdeJZsrvt2mszsDaAx7jp6M7P9ZrYp3D4CbCNaSzpWFmkOmwXhJyeOZCSNBb4ILIm7llwnaQgwB3gCwMzacykEgnnA+3GHQJp8YKCkfGAQGV7JMSlBMAbYm7ZdRw58sZ0PJE0ApgPr460kErpfNgMp4FUzy4m6gIeAu4HuuAvpxYBXJG2UdGvcxQSXAAeBn4SutCWSPn519/61AHgm7iIAzOxD4AHgA2A/0UqOr2TyPZISBOrjvpw4ksxlkkqBXwJ3mllT3PUAmFmXmf0h0frWMyTF3qUmaT6QMrONcdfSh9lmdgVwA3Bb6I6MWz5wBfAvZjYdaAFy4rwdQOiq+jLwXNy1AEgaTtSDMRG4CCiR9LVMvkdSgqAOGJe2PZYMN60uNKEP/pfAz81sWdz19Ba6El4Dro+5FIDZwJdDf/yzwFxJP4u3pIiZ7Qv/poAXiLpJ41YH1KW15p4nCoZccQOwycwOxF1I8Flgt5kdNLMOYBkwK5NvkJQg2ABMljQxpP0C4MWYa8pZ4aTsE8A2M3sw7np6SBohaVi4PZDoD6Q23qrAzO4xs7FmNoHod2uFmWX0iO1sSCoJJ/sJXS+fB2K/Qs3M/hfYK2lKuGseEOuFCL3cRI50CwUfAFdLGhT+NucRnbfLmKyuWZwrzKxT0u3Ay0Ae8KSZbY25LCQ9A3wGqJBUB9xnZk/EWxUQHeHeDGwJ/fEA/xDWoI7TaODpcEXHAOAXZpYzl2rmoJHAC9F3B/nAUjP7r3hLOm4R8PNwYLYLuCXmegCQNIjo6sK/jruWHma2XtLzwCagE/gtGZ5qIhGXjzrnnDu5pHQNOeecOwkPAuecSzgPAuecSzgPAuecSzgPAuecSzgPApdokn6YPsOkpJclLUnb/oGkfzzTGWslPSXpTzNZq3PZ4kHgku5NwihNSQOACmBq2uOzgJfN7Psx1OZcv/AgcEm3hhPD9acSjbw9Imm4pCLgMuAPJD0Gx4/0H5H0pqRdPUf9ijwm6T1J/wFU9ryBpHlhcrUtYQ2KIkkzJC0Lj98o6aikwrDmwq5+/PzOeRC4ZAtz8XRKGk8UCGuJZlqdSbS2wDtAe6+njQauAeYDPS2FPwGmANOAhZxoZRQDTwFfNbNpRCN8v040SnR6eO61RAF0FVBNjsz06pLDg8C5E62CniBYm7b9Zh/7/7uZdYeFekaG++YAz4TZUfcBK8L9U4gmDNsRtp8G5phZJ7BT0mVEE8E9GF7jWmBVpj+gc6fiQeDcifME04iOzNcRtQhmEYVEb21pt9OnOO9rvpa+pkDvsYpopssO4L+JWhnXAG+cbuHOZYIHgXPRl/18oDEc0TcCw4jCYO1pvsYbwIKwcM5ooCbcXwtMkDQpbN8MvJ72nDuBtWZ2ECgHqoDYJ0R0yZKI2Ued+xhbiK4WWtrrvlIzqw+zd36cF4C54Xk7CF/2ZnZM0i3Ac2GZwQ3A4vCc9URdSz0tgHeIFrjxmSBdv/LZR51zLuG8a8g55xLOg8A55xLOg8A55xLOg8A55xLOg8A55xLOg8A55xLOg8A55xLu/wGrnnrcuPvoQAAAAABJRU5ErkJggg==\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 }