{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Writing your own trajectory analysis\n",
"\n",
"We create our own analysis methods for calculating the radius of gyration of a selection of atoms.\n",
"\n",
"This can be done three ways, from least to most flexible:\n",
"\n",
" 1. [Running the analysis directly from a function](#Creating-an-analysis-from-a-function)\n",
" \n",
" 2. [Turning a function into a class](#Transforming-a-function-into-a-class)\n",
" \n",
" 3. [Writing your own class](#Creating-your-own-class)\n",
"\n",
"The building blocks and methods shown here are only suitable for analyses that involve iterating over the trajectory once.\n",
"\n",
"If you implement your own analysis method, please consider [contributing it to the MDAnalysis codebase!](https://www.mdanalysis.org/UserGuide/contributing.html)\n",
"\n",
"**Last executed:** Feb 07, 2020 with MDAnalysis 0.20.2-dev0\n",
"\n",
"**Last updated:** February 2020\n",
"\n",
"**Minimum version of MDAnalysis:** 0.19.0\n",
"\n",
"**Packages required:**\n",
" \n",
"* MDAnalysis (Michaud-Agrawal *et al.*, 2011, Gowers *et al.*, 2016)\n",
"* MDAnalysisTests\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import MDAnalysis as mda\n",
"from MDAnalysis.tests.datafiles import PSF, DCD, DCD2\n",
"from MDAnalysis.analysis.base import (AnalysisBase,\n",
" AnalysisFromFunction,\n",
" analysis_class)\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Radius of gyration\n",
"\n",
"Let's start off by defining a standalone analysis function.\n",
"\n",
"The radius of gyration of a structure measures how compact it is. In [GROMACS](http://manual.gromacs.org/documentation/2019-rc1/reference-manual/analysis/radius-of-gyration.html), it is calculated as follows: \n",
"\n",
"$$ R_g = \\sqrt{\\frac{\\sum_i m_i \\mathbf{r}_i^2}{\\sum_i m_i}}$$\n",
"\n",
"where $m_i$ is the mass of atom $i$ and $\\mathbf{r}_i$ is the position of atom $i$, relative to the center-of-mass of the selection.\n",
"\n",
"The radius of gyration around each axis can also be determined separately. For example, the radius of gyration around the x-axis:\n",
"\n",
"$$ R_{i, x} = \\sqrt{\\frac{\\sum_i m_i [r_{i, y}^2 + r_{i, z}^2]}{\\sum_i m_i}}$$\n",
"\n",
"Below, we define a function that takes an AtomGroup and calculates the radii of gyration. We could write this function to only need the AtomGroup. However, we also add in a `masses` argument and a `total_mass` keyword to avoid recomputing the mass and total mass for each frame."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def radgyr(atomgroup, masses, total_mass=None):\n",
" # coordinates change for each frame\n",
" coordinates = atomgroup.positions\n",
" center_of_mass = atomgroup.center_of_mass()\n",
" \n",
" # get squared distance from center\n",
" ri_sq = (coordinates-center_of_mass)**2\n",
" # sum the unweighted positions\n",
" sq = np.sum(ri_sq, axis=1)\n",
" sq_x = np.sum(ri_sq[:,[1,2]], axis=1) # sum over y and z\n",
" sq_y = np.sum(ri_sq[:,[0,2]], axis=1) # sum over x and z\n",
" sq_z = np.sum(ri_sq[:,[0,1]], axis=1) # sum over x and y\n",
" \n",
" # make into array\n",
" sq_rs = np.array([sq, sq_x, sq_y, sq_z])\n",
" \n",
" # weight positions\n",
" rog_sq = np.sum(masses*sq_rs, axis=1)/total_mass\n",
" # square root and return\n",
" return np.sqrt(rog_sq)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Loading files\n",
"\n",
"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": 3,
"metadata": {},
"outputs": [],
"source": [
"u = mda.Universe(PSF, DCD)\n",
"protein = u.select_atoms('protein')\n",
"\n",
"u2 = mda.Universe(PSF, DCD2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating an analysis from a function\n",
"\n",
"`MDAnalysis.analysis.base.AnalysisFromFunction` can create an analysis from a function that works on AtomGroups. It requires the function itself, the trajectory to operate on, and then the arguments / keyword arguments necessary for the function."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"rog = AnalysisFromFunction(radgyr, u.trajectory, \n",
" protein, protein.masses, \n",
" total_mass=np.sum(protein.masses))\n",
"rog.run();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Running the analysis iterates over the trajectory. The output is saved in `rog.results`, which has the same number of rows, as frames in the trajectory."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"(98, 4)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rog.results.shape"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEGCAYAAABlxeIAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3zV1f348de5N/fmZic3OyQhYYYtkLBXUAQVELQDZ4ta6l7VTltFv/ZXW7RV0VIHxT1QqYMComxBIAxZBgIZZJCd3OTm7nvP748bqcjKhSQ34zwfjzzI/az7viG57/s5432ElBJFURRF0fg7AEVRFKVjUAlBURRFAVRCUBRFUZqphKAoiqIAKiEoiqIozQL8HYAvYmJiZFpamr/DUBRF6VR27dpVLaWMPd9xnSohpKWlkZOT4+8wFEVROhUhRFFLjlNNRoqiKAqgEoKiKIrSTCUERVEUBehkfQhn4nQ6KSkpwWaz+TsUvzAYDCQnJ6PT6fwdiqIonVynTwglJSWEhYWRlpaGEMLf4bQrKSU1NTWUlJSQnp7u73AURenkOn2Tkc1mIzo6utslAwAhBNHR0d327khRlNbV6RMC0C2TwXe682tXFKV1dfomI0VRlI5KSsnhika2Hq0hLjyQsb2iiQ4NbPG5JXVW9pWY2Fdaz13ZfQg3tG1foUoIbeS7SXQxMTGEhoZiNpv9HZKiKOdhcbgw210IBEJAVLAereb0u3CX28O63Ere3H6cr4/VEB4UgDFET1SwntDAAIIDA9AI2FFQywnTqU26AxLDSTUGodUINEIQog8gJkxPTGggUkJ+tZmjlWaOVJipbXIAoNdquGJwIpekRLbp61cJQVGUbi+3vIF/bynkP3tLsbs8J7cnhBu4fnQq87JSMIbo2VVUx7rcSj75powTJtvJ/XaXh9omO3VNTsobbFgcbuxON8OSI7n/slgm9I2lssHG1mM1bD1WTWG1BbeUuD2SJruLmiYHbo93sbJwQwB94kKZNiCeIckRDEuOpF9CKIEB2jb/OaiE0ArmzJlDcXExNpuN++67jwULFvg7JEVRfsDp9rDpSBUr9pSy6UgVBp325B3AoRMNGHQarhnRg4FJESAlLo9k/eEqnll7hOe+zCNYr6XB5kKnFYzrHcOjswZx2YA4ArQt64rtERnE8NQo7sruc9o+j0disjpxS0l0iN5vfYNdKiEs/PQgh8oaWvWaA5PCeXTWoHMes3TpUoxGI1arlaysLK699tpWjUFRlFM5XB4qG21UNtqpbLBTWm+luNZCca2F6iYHVocLq9ON0yXRagQBWkG9xYnJ6iQqWMf0QQlohKDe6sBsd/HbKzKYl5VCZLD+lOeZPz6dguom3tlxnHqLg+z+cUzoG0NYK7flazSCqBD9+Q9sY10qIfjLc889x4oVKwAoLi4mLy/PzxEpStfSYHOy4XAVu4vq2FNcz7dlDTjcnlOOCdFrSTEGExduIDHcQLBeS4BW4PaAy+PBEKDl8kHxTOwbiz6g5QMs02NC+P2VA1r7JXVIXSohnO+TfFvYsGEDX3zxBdu2bSM4OJgpU6aoeQGKcpGcbg8F1U3sKzGx+sAJNh2pxuH2EKzXMqRHBPPHp5EeE0JceCBxYQYSIwwY/djU0lV0qYTgDyaTiaioKIKDg8nNzeXrr7/2d0iK0uk02JxsO1bDlrxqdhTUkl9txun2drImRhi4aWxPrhySwLDkyBa32Su+UwnhIs2YMYMlS5YwdOhQ+vfvz5gxY/wdkqJ0CnaXmy8OVfJ+TjFbjlbj9kiC9Vqy0oxMHRBH//gw+ieE0T8+DM0Zhn4qrU8lhIsUGBjIqlWrTtteWFh48ns1B0HpDg6VNfB+TjGf7TtBgEaQYgwiJSqY3nGhZCSEkZEYjtXhJqewlp2FdazLraDO4iQxwsCCSb2Y0i+W4alRPrXvK61LJQRFUS6YlJLPD1WweN1R9pea0Gs1TBsYT5BeS3Gtha/za/hoT+lp5xlD9EzoG8uPRyYzvk/MGSd/Ke1PJQRFUS7I1qPV/HXNYfYW15MeE8JjswYyZ3iP04ZuNtqcHKlo5NsTjei1GjLTokiPCVEdwB2QSgiKovjE6nDz2CcHeS+nmKQIA09dO4RrRySftbM3zKBjZE8jI3sa2zlSxVcqISiK0mJHK83c9dZujlQ2cld2b+6Z2heDru1LKijtQyUERVHOye5ysz2/lnW53hFBBp2W1+aPYlK/WH+HprSyNk8IQoilwEygUko5uHnbMGAJEAoUAjdIKVu35oSiKBeloLqJJRuO8em+MiwON4EBGqZmxPHorEEkRBj8HZ7SBtrjDmEZsBh4/XvbXgEeklJuFELcAjwM/LEdYukUlixZQnBwMDfffLO/Q1G6obyKRp5bd5SV+8rQaTXMHd6DywfFM7ZXDEF61TzUlbV5QpBSbhJCpP1gc39gU/P3a4E1qIRw0u233+7vEJRuqMHm5B9r83htWyGGAA0LJvXm1gnpxIa1bEEXpfPz1wyQA8Ds5u9/DKSc7UAhxAIhRI4QIqeqqqpdgvPFzp07GTp0KDabjaamJgYNGsSBAwdOOebTTz9l9OjRDB8+nMsuu4yKigoA7r33Xh5//HEA1qxZw6RJk/B4PDz22GMsWrQI8BbOGzhwIEOHDmXevHnt++KULstkcbL2UAXv7TzOv78q4O9rjzB10Qb+vbWAn2alsPk3U/ntFRkqGXQz/upUvgV4TgjxJ+ATwHG2A6WULwEvAWRmZspzXnXVb6F8fyuGCSQMgSv+ctbdWVlZzJ49m0ceeQSr1cqNN97I4MGDTzlmwoQJfP311wgheOWVV/jrX//K008/zV/+8heysrKYOHEi9957L//973/RaE7N0X/5y18oKCggMDCQ+vr61n1tSrdR1Whn9/E6dhfVsfVYDQfKTMgf/DWN7BnFsvmjGNwjwj9BKn7nl4QgpcwFLgcQQvQDrvJHHK3lT3/6E1lZWRgMBp577rnT9peUlPDTn/6UEydO4HA4SE9PByA4OJiXX36ZSZMm8fe//53evXufdu7QoUO54YYbmDNnDnPmzGnz16J0HbnlDazYXcqqA+Ucr7UAoNMKLkmJ5L5L+zKudwzJUUEE6bQE6bVq+Kjin4QghIiTUlYKITTAI3hHHF28c3ySb0u1tbWYzWacTic2m40///nPrFy5EoC9e/dyzz338OCDDzJ79mw2bNjAY489dvLc/fv3Ex0dTVlZ2RmvvXLlSjZt2sQnn3zCE088wcGDBwkIUKOFFW+J6MPljeRXN5FfZaaiwYbV4cbm9FBY00RueSMBGsHEvjHcNKYnI3pGMigpQr3xK2fVHsNO3wGmADFCiBLgUSBUCHFX8yEfAf9u6zja0oIFC3jiiScoKCjgN7/5DYsXL+bJJ588ud9kMtGjRw8AXnvttZPbi4qKePrpp9mzZw9XXnklc+bMYfTo0Sf3ezweiouLyc7OZsKECbz99tuYzWYiI9t2oW2l4ztWZeb2N3aRV+ktnCgExIQGEqTTYtBpiArW89isgcwalkR0qOoHUFqmPUYZXXeWXc+29XO3h9dff52AgACuv/563G4348aNY926dUydOvXkMY899hg//vGP6dGjB2PGjKGgoAApJbfeeiuLFi0iKSmJV199lZ///Ofs3Lnz5Hlut5sbb7wRk8mElJIHHnhAJQOFNQfL+dX736AP0LDox8MY0iOCntHB6pO/ctGE/GHPUgeWmZkpc3JyTtn27bffMmBA91je7mzUz6B7sDrc/OOLI/xrUz7DkiN48caR9IgM8ndYSicghNglpcw833GqMVpROoG1hyp47JODlNZbuW5UKo/OGqjuCJRWpxKConRQFoeLL76tZHlOMZvzqukXH8p7C8Ywule0v0NTuiiVEBSlAzhYZuKJzw7hckuCAwMQwI6CWqxON/Hhgfz+ygzmj09Hp9YTVtqQSgiK4md7jtfxs6U7CNRp6RMbisnqxO50c+3IHswamkRWmlGtKay0C5UQFMWPvs6v4dZlO4kJC+St20aTHBXs75CUbkwlBEVpR1JKimos7CysJaewjo+/KSU5Kpi3bhtNfLgqKa34l0oIHZAqf9017S8x8cTKQ+woqAUgIkjHpRnxLLx6EDFq8pjSAaiE0AGp8tddS2m9lWc+P8JHe0owBut55KoBTO4XS+/YUNU3oHQoasjCRfrjH//Is8/+b9L1H/7wh9MK3Kny192PlJKv82u4481dTPrrej79powFk3qx/uEp3DaxF33jw1QyUDqcLnWH8NSOp8itzW3Va2YYM/jNqN+cdf+tt97KNddcw3333YfH4+Hdd99lx44dpxyjyl93fVWNdnYU1HKgzMTh8ka+PdHACZONyGAdt01M5+axaWpWsdLhdamE4A9paWlER0ezZ88eKioqGD58ONHRp04cUuWvu6baJgcvrj/KxiNVJ4vM6bSC3rGhjEo3Mr53DLOGJallJ5VOo0slhHN9km9Lt912G8uWLaO8vJxbbrmFP/zhD6r8dRfm8UjeyynmqdW5mG0uxveJ4ZoRyYztHc3AxHD0AaolVumc1DtLK5g7dy5/+tOfcDqdvP3221x55ZWq/HUXVFpvZX1uJR/sKmFvcT2j0o08OWcwfePD/B2aorQKlRBagV6vJzs7m8jISLTa05sHVPnrzsvicPH+zmLe2VHM4YpGAFKMQSz68TCuHdEDIVTHsNJ1qPLXrcDj8TBixAiWL19O37592/35O8LPoKsxWZy8+lUBr28rpN7iZHhqJFcOTiQ7wztcVCUCpTNR5a/byaFDh5g5cyZz5871SzJQWpfL7eGdncU88/lh6ixOLhsQz+2Te5GZZvR3aIrS5lRCuEgDBw4kPz/f32EoF8jp9lBUY6G41kJRTRPv7iwmt7yRMb2M/GnmIAYmhfs7REVpN+2xpvJSYCZQKaUc3LztEmAJYABcwJ1Syh1nv4qitK5jVWbe3XGcD3eXUtvkOLk9xRjEP28YwYzBCapZSOl22uMOYRmwGHj9e9v+CiyUUq4SQlzZ/HhKO8SidGN2l5vVB8p5e/txthfUEqARTBsYz7SB8fSMDiHVGExMqF4lAqXbavOEIKXcJIRI++Fm4Lt78QjgzIPwFaUVNNicvLj+GO/tPE6dxUmqMZhfz+jPj0emEBumisopynf81YdwP7BGCLEIbz2lcWc7UAixAFgAkJqa2j7RKV2ClJL/7i9n4acHqTLbmTEogetHpzK+d4yqI6QoZ+CvhHAH8ICU8kMhxE+AV4HLznSglPIl4CXwDjttvxD9Iycnh9dff/20AnlKy9Q1OTh0ooFvTzSw4XAVW45WMygpnJdvzmRYiprDoSjn4q+E8DPgvubvlwOv+CmODiczM5PMzPMOF1a+R0rJlqPVLN1SwPrDVSe3x4cH8seZA/nZ2J4EqLWIFeW8WpwQhBBxwHggCbACB4AcKaXnAp63DJgMbACmAnkXcI0OYcmSJSxZsgTwlqhIS0tj/fr1J/cXFhZy00030dTUBMDixYsZN24cK1as4IUXXmDt2rWUl5czefJkNm3aRG5uLosWLeKzzz5j48aN3HefN28KIdi0aRNhYapMwndsTjcf7y1l6ZZCDlc0EhMayL1T+zAqPZoBiWFEq0VnFMUn500IQohs4LeAEdgDVOIdLjoH6C2E+AB4WkrZcJbz38E7gihGCFECPAr8AnhWCBEA2GjuI7hY5X/+M/ZvW7f8deCADBJ+//uz7r/99tu5/fbbcTqdTJ06lQcffPCU/XFxcaxduxaDwUBeXh7XXXcdOTk5zJ07lw8//JAXXniB1atXs3DhQhISEsjN/V/8ixYt4oUXXmD8+PGYzWYMBrXEosPl4ViVmdUHynnz6yJqmhxkJITxtx8NZfYlSQQGqMqiinKhWnKHcCXwCynl8R/uaH5DnwlMAz4808lSyuvOct2RLQ2yM7jvvvuYOnUqs2bNOmW70+nk7rvvZu/evWi1Wo4cOXJy3/PPP8/gwYMZM2YM1113+o9p/PjxPPjgg9xwww1cc801JCcnt/nr6IjqLQ7+ufEY63Mrya9qwuXxdiVdmhHHrRPSGds7Wg0VVZRWcN6EIKV8+By7I6WU/2nFeC7KuT7Jt6Vly5ZRVFTE4sWLWbFiBQsXLgTglVde4bPPPiM+Pp5vvvkGj8dzyqf80tJSNBoNFRUVeDye0xbH+e1vf8tVV13Ff//7X8aMGcMXX3xBRkZGu742f7K73Ly+tYjF64/SYHMyuV8s0wbG0y8+jBGpUaQYg/0doqJ0KT53KgshgvA2F90IDMfbp9Bt7dq1i0WLFrF582Y0Gg1z585l7ty5J/e/9dZbJCcno9FoeO2113C73QC4XC7mz5/P22+/zeuvv84zzzzDQw89dMq1jx07xpAhQxgyZAjbtm0jNze3WySEkjoL7+8s5r2cYioa7EzqF8vvrshgQKIqI6F0T+76ejTh4QhN2w6OOGtCEEJogKF4O489wOV4k8BkvJPKrgY2t2l0ncDixYupra0lOzsb8I4SeuWV/w2auvPOO7n22mtZvnw52dnZhISEAPDnP/+ZiRMnMnHiRC655BKysrK46qqrTrn2P/7xD9avX49Wq2XgwIFcccUV7ffC/GBXUS0vrD/G+sOVAEzuF8uiH6czsW+snyNTlDOTHg/OshMIXQCawEBEUBBC753tLqXEnpuLecsWLDt3EpieTsTVVxM4YMA5mzillNiPHKHpq63YDuzHuv8AzuJieq38jMAzrKrYms5a/loI8T5wHG8fggHYDrwLrAZypZTpbRrZGXTU8tf+1pl/BlJKthfU8vy6PL46WkN0iJ4bRqfyk6wUkqNUk5DSfqTTSdPWrXgsFjTh4WjDIwjslY6m+UPc97kbGzGtWEHdW2/jKCo6dacQ3sQgBJ7m0YX6Xr1wFhcjnU4C+/YlsF8/pNOJdDoRej0BMTEExMbgrquj8ct1OEtKANAlJWEYPBjDkMFEzL4aXXzcBb221ih/3QN4E7gH2AesBzZLKe1CiC4/QUxpW0U1TazYU8p/9pRSWGMhJjSQR64awPWjUwnWqyK8iu9cdXXYDhzEY7EQ2Csdfc+eCL3+1GOqqmj84gssObvQpSRjGDgQfWoqjV98Sf377+OqrDzleE1oKJHXXkPUDTcQEBtL01df0bh2LQ1rv0BaLARdcgnxN9+ECNAh7TY8Vpv3X4sV6XRiGDyYkPHj0MXF4a6vp2H1akyffobtwAGEXgc6HdJmp+nrr/GYTAidjpBx44he8AtCp0xBF3dhCeBCnesv7wbgWmAwoANuArYLIQ4DIUKIYCmlpR1iVLqQA6Umnvsyj88PVSAEjO0VzZ3ZfZg1VC1Gr5xOejze5pNtX2PJycFdW4unqQmPxQJaDZqgYDRBQbgqK3GWlp56slaLLjERbXg4mvBwpNWKdd8+kJKA2Fhcq1dDc58eQhAyYQIJjz2KLjkZT0MDrro6Glevofatt6l9/Q1EYCDSZkMbEUH4lVcQNe86ggYPavFr0UZGEjVvHlHz5p1xv8duBynR+HF4uc8rpgkhJuHtS5iJd2La7LYI7EzO1mSUkZHRbYcdSinJzc3tsE1GdU0O8irNHKloZH1uJV/mVhJuCODn49O5blQKiRFB/g5R6WBc1dU0ffUV5i1f0fTVV7hrawHQ9+xJQFIimpAQNEHB4PHgsVrxWC1oIyIJGjIYw6DBaMJCceTnYz92DGdJKZ7GRtyNjSAlIRPGE3755ej79EE6HNiPHMF+9BjBI0egP0utNGdFJfXvv4+7oYGwqdkEZ2YidLr2/JFctJY2GV3wEppCiEDgSinligu6wAU4U0IoKCggLCyM6OjuNxZdSklNTQ2NjY2kp7d7l8457S2u56+rc9l6rObkNmOInvnj0vjZ+DTCDZ3rD0ppO25zE+aNG7Ds3IklJwfH0WMAaI1GQsaPJ2TsWELGjkGXmOjnSDuvVltCUwhxI/D2D0tUSCntwAohRG8gUUq55YKjvQjJycmUlJRQVVV1/oO7IIPB0GEmrEkp2V9q4sX1x1h9sJzoED2/mtaPoSmR9I0LJTHC0O2SdnckPR5vx+p5/q9d1dXUvvEmde+8g6ehAU1ICEEjRhAx+2pCxo/DMGBAmw+zVE7Vkt67aGCPEGIXsAuowjvqqA/eIajVeEtb+IVOp+twn467m9zyBlbsKWXV/nKO11oIDQzgwWn9uGVCOqGBqoO4O3A3NGDeuJHGz9di3rIF3G60RiNaYxQBkVFoIyPRRkYipQd3bR3u2lqs+/YhHQ7Cpk3D+LObCRo2DBGgfl/8qUVNRkIILd4idOOBRLzF7b4FVp2ppEVbOVOTkeI/e4vrWbwujy++rSRAIxjfJ4arhiQyfVACEcGqSag78FitVL/0ErWvLkU6HATExRE6NRtNcAju2lpcdbW46+ubv0wIQBsdjdYYRWDfvhhvupnAXuoDXVtrtSYjACmlG1jb/KV0Y26PZF1uJa9tLWTL0Woig3U8OK0fN43pSVSI/vwXULoEKSWNaz6n4q9P4So7QfhVV2G86UYMQ4eqZp5OTN2fKedlcbg4UNrA9vwa3t1ZTGm9lYRwA7+ZkcFNY3uqZqFuRjocnFi4ENOHHxGYkUGPv/6VYLWGR5eg/pKVM6posPHh7hI+/eYEh8sbaC4wythe0Txy1QAuGxiPTi060+24amspufderDm7iL7jdmLvuku1+3ch6n9SOam41sKWo9WsOVjOpiNVeCRkpUVx99S+XJISwdDkSGLUojPdkquujqbNm6l69jlc1dX0eOZpwq+80t9hKa3MlxXTAvHOXE77/nlSysdbPyylLTlcHjYeqeKb4nqqzXaqzQ6OVZkpqPbWXekRGcSdU/rwo5HJpMWcXsdF6R48Fgv1K1bQ8OlnWL/5BqREl5REzzffIGjIEH+Hp7QBX+4QPgZMeIee2tsmHKWtWB1u9hTXsWp/OZ/tK6PO4kSrERhD9ESH6OkdG8JNY3oyqV8MvWND1XyBbsxVVUXtm29R9+67eEwmAgcMIObOOwmdMhnDoEGq07gL8yUhJEspZ7RZJEqr8XgkRbUW9pea2F9ST05RHQdKTTjdksAADZcPSuCa4T2Y0DdG9QMoJ7nq6qh55RXq3nobabcTdtmlGOffQvCI4f4OTWknviSErUKIIVLK/W0WjXJRyuqtvPl1Ee/tLKamyQGAXqthSHIEt07oxaj0KLLSjISpshEK4KyowH74MI7jxTjyj2H6+BM8FgvhM2cSe9ed6NPS/B2i0s58SQgTgJ8LIQrwNhkJQEoph57rJCHEUryF8CqllIObt70H9G8+JBKol1Je4mvwitehsgYWr89jzcEKpJRcOiCeSzPiGJIcQb/4MHUXoJzCnpdH9T+X0LB6NXi8FWlEUBChEyYQc8/dGPr183OEir/4khAudLmuZcBi4PXvNkgpf/rd90KIp/H2TSg+yq8y88zaI3y27wRhhgBum5jOjaN7qrWGldNIKbHu3k3tstdoXLsWTXAw0bfMJ3TqVPQpKWhjYlS/kdLyhCClLBJCDAMmNm/aLKX8pgXnbRJCpJ1pn/D+Bv4Eb1kM5SxsTjdHKhrZX2oi90QjhTVNFFQ3UVpvJUin5e7sPvxiYi9VLkI5jcdqpWHlSmrfehv7t9+iCQ8n5s47iLrpJgKiovwdntLB+DLs9D7gF8BHzZveFEK8JKV8/iKefyJQIaXMO8fzLgAWAKSepV55V+PxeKuGbjxSdXJ4qKt5ZlhYYADpsSGMSI1iXlYK80alqrkByimk2409L4/695dj+vRTPI2NBPbtS8LjC4mYNQtNkFqDQjkzX5qMbgVGSymbAIQQTwHbgItJCNcB75zrACnlS8BL4C1udxHP1aF5PJI9xXV8+s0JVu4/QVWjHSFgaI8IbpvYi6HJEQzpEUFyVJC6tVeQHg+m/3xM09atuOvqvF/19bgbG/GYzSAlQq8nbPp0on7yY4IyM9XvjXJeviQEAbi/99jdvO2CCCECgGuAkRd6jc5OSsm+EhOf7Stj5b4TlJlsBAZomJoRx4zBCUzsG4tRFYxTfsCen0/5nx7FkpNDQGIiAXGxaGNjCOzbB014BNqwMALi4wm7fJpqFlJ84ktC+DfeNZW/WyFtDvDqRTz3ZUCulLLkIq7RoXk8kvxqM0cqzFgdbpxuD1anm6IaC8eqzORVmClvsKHTCib1jeXhGf2ZNjBBFYtTzsjT1ETNq0upefllRHAwiU/+HxHXXKM++SutxpdO5WeEEBvwDj8VwHwp5Z7znSeEeAeYAsQIIUqAR6WUrwLzOE9zUWdUVm/ls31lrMut5EBpA2a767RjgvVaesWGMCrdyIS+MUwfqNYPUM5Out2YVqzw1hGqqiJ85kzif/dbAqKj/R2a0sW0ZAnNcCllgxDCCBQ2f323zyilrD3X+VLK686y/ec+RdpB1FscnDDZcLo9ON0Sk9VBWb2NsnorOUV17Cjw/jgGJoYzd3gPhiZHMCAxnDBDADqthsAADcYQvfpUp5yXdLloWLWa6n8twXH0GEHDh9PjuWcJHq5mDittoyV3CG/jnVi2C/h+p65oftyrDeLyGyklJXVWjlaZSTUG0ysmBCEE5SYbSzYe450dx7G7PKedF6AR9IoN4cFp/Zg9LEkVhVN84rHZMP3nY6TDjtDp8Fht1L37Ls7jx9H36U2Pf/yDsOmXqw8SSps6b0KQUs5s/rdLrXNXUN3Es18cYe2hCsIMOowheoL0WvIqGmmw/a+ZJypYR0ZCOLuK6nBLyTXDe5CdEYdeqyFAKwgz6OgRGURsWCBajfpjVXznaWqi+M67sGzffsp2w6BBxD3/HGGXXqoKyintwpd5CF9KKS8937aO6GhlIxUNdqwON1anmy151XywuwSdVjB7WBJSQm2TA7PdxcxhSQxKCqdPbCiFNU3sKqpjf2kD147swR2T+5AarWYBK63H3dBA8YJfYt2/n6Sn/kLo5MlIpxPp8RAQF6fuCJR21ZI+BAMQjLdTOIr/DTUNB5LaMLZW87c1h1lzsOLkY71Ww01jenJndm/iwgxnPW90r2h+mtU9JsMp7UtKie3AQcoffRRbXh49/vF3wqdN87RicpcAACAASURBVHdYSjfXkjuEXwL3433z38X/EkID8EIbxdWq7r+sH7eMTydIr8Wg0xIXFkhksBrfr7Q/e0EB9R98QOPqNThLSxFBQaS8sJjQSZP8HZqitKgP4VngWSHEPRdZpsJvBiSG+zsEpZuz5+VRveRfNKxaBRoNIePGEnPnnYRdOhVtZKS/w1MukNPt5ETTCRJDEtFpO//QcV/mITwvhBgMDAQM39v++tnPUpTuy2O10vjlOkyffEzTps2I5gqjxp//nICYGH+Hp/xAcWMxawrXYHVZCdOFEaYPo7+xP4OiB53Sl1NjrWFTySY2l25ma9lWmpxNaIWW5LBk0iPS6R/VnwHGAaRFpHG49jDby7ezu2I3AEaDkeigaPpG9mVs0lgGxwwmQHP2t+GKpgq2ndjGtrJtPDLmEcL0YW36M/ClU/lRvBPMBgL/xVsOewvfK2utKN2ddDpp2raNhv+uonHtWjxNTQQkJqoKox2MxWnhRNMJypvKKWwo5PPCz9ld6X3TFgjk90bYJ4cmc0X6FYTpw1hfvJ69lXuRSOKC4piRNoMhMUMoayqjwFTAsfpjbCrZhEf+b2h6mD6MkfEjCdQGUmOtIa8ujy+KvuDFb14kVBdKv6h+RAZGEhEYgV6rp9HRiNlppqSxhHxTPgDRhmiONx5nUPSgNv25CClbVi9OCLEfGAbskVIOE0LEA69IKWe1ZYDfl5mZKXNyctrr6RSlRaTTSdOOHTSuXkPj2rW46+vRhIURdvk0ImZfTXBWpho22sqklFRaKiloKKDQVIjZaUYrtGiEhkBtIKH6UEJ1oUQZokgLTyMiMAK3x82W0i0sP7KczaWbT3nTTo9IZ3bv2czsNZO44DgsTgsmh4kdJ3awunA1209sxy3dDDAOIDslmykpU8gwZpxxFJjNZSOvLo+ChgJ6R/Qmw5iBVqM95Zh6Wz07ynew7cQ2Ck2FmBwmTDYTTo/zZOyxwbFkxWcxNmks/aL6XdSIMyHELill5nmP8yEh7JBSjhJC7AKygUbggJSybVPW96iEoHQkrro6qp57jsZVq71JIDiY0Oxswq+6ipAJ49Ho1cCF1iSl5Juqb1iZv5LPiz6n1nbOIgmniAyMRCu01NhqiAmKYVbvWfSP6k9CSAKJIYkkhiSe8w23zlaHw+0gPiS+NV5Ku2tpQvClilqOECISeBnvaCMzsOMC41OUTq1p61bKfvs7XHV1hE+fTviM6YRMmIDGcPZhzIrvnG4nOyt2srF4I+uL13Oi6QSB2kAmJ08mKyGL9Ij0k3cAHunBLd3Y3XbMDjNmp5lqazVFDUUUNRRhdpi5PO1yJqdMRqfxrQM4ytA9mvpalBCaVzb7f1LKemCJEGI1EC6l3Nem0SlKByKlxFFQSN0771D3xhvoe/Ui7Z8vEjSo3W6SuzSzw8zeqr3sr9pPQUMBRQ1FFJgKsLqsGLQGxiSO4Z7h9zA1dSohurOXhgkjjJgg1Wl/IVqUEKSUUgjxH5rXLpBSFrZlUIrSUbhNJiw7d9K0dRvmzZtxFhcDEHX99cQ9/JBafewimOwmdlXsIqcih10Vu8itzcUjPQgESaFJpIWnMbzvcMYkjmF04miCAtTPuq350mT0tRAiS0q5s82iUZQOQDoc1H/0EfXLP8B26JB39bGgIELGjCH6lvmETJyEPrmHv8PscCxOC5tLN7OheAMWpwW9Vo9eqychJIH+Uf3pb+yPyW5ic+lmNpds5lDNISQSvUbP0Nih/GLILxgZP5JhscMI1qkSMf7gS0LIBn4phCgCmmiudiqlHNomkSlKO5JS4mmy0Pj551S/8ALO0lIMgwYRc9ddhIwZTdDQoQjVSXwaKSXby7fz4ZEP2ViyEavLenKsvdPtxO62U2mpxC3/t9iiRmgYGjOUOy65g6z4LIbEDiFQq9YF7wh8SQhXtFkUitLOpJRYduyk9rXXsB04gLu+HulwAN4qowmPPUrIhAmquNxZ2N12Psr7iHdz3yXflE9kYCSze89metp0RsSNOGWYpd1t52j9UQ7XHsagNTAuaRyRBjU7uyPyJSHcAyyVUh5qq2AUpa15HA4a16yh9t/LsB06hNZoJHTyZAKijWijogjs25eQiRNVIjiHfVX7+ONXfyTflM/g6ME8OeFJpqdNP+un/EBtIIOiB7X5pCrl4vmSEHKBl4UQAXjXV35HSmlqm7AUpXU5S0upe+996j/4AHdtLfr0dBIeX0jE7NlqqGgLON1OSswlrDi6gtcOvkZsUCwvXvoiE5Mn+js0pRX5UsvoFeAVIUR/YD6wTwjxFfCylHJ9WwWoKBdDulzUvLqUqsWLwe0mdGo2UfOuI2TcWDV7+DxqbbW8vO/lk+P/v5vZe23fa/lV5q/avK6O0v58uUNACKEFMpq/qoFvgAeFEL+UUs47yzlL8S7BWSmlHPy97fcAdwMuYKWU8tcX9hIU5czsR49S9rvfY9u/n7ArZhD/8MPokjrFEh5+ZXFaePPbN1l6YClWl5XslGxm9Z5FSljKydFCStfkS3G7Z4DZwJfAn6WU381SfkoIcfgcpy4DFvO9InhCiGzgamColNIuhIjzNXBFORNHcTHm9etpXLcey86daMPC6PH3Zwi/Qo2JOB8pJasKVvH0rqeptFSSnZLN/SPup1dkl1o2XTkHX+4QDgCPSCktZ9g36mwnSSk3CSHSfrD5DuAvUkp78zGVPsShKKdwlJTQsGoVDatWYT/0LQCBffsQfeutGG++SZWa/gEpJSWNJRyoOYDD7cAQYEAjNLx56E12V+5mgHEAf5v0N0bEj/B3qEo78yUh7AUyfjD6wgQUXUDncj9gohDiScAGPHS2CW9CiAXAAoDUVLWcpeIlpcSyfTs1r7xK05YtABiGDiXu178m7LJL0avfFdweN5/mf8qKvBUA6LQ6BIIjdUfOWBguKjCKR8c+ytw+c0+rzql0D74khBeBEcA+vJPSBjd/Hy2EuF1K+bmPzxsFjAGygPeFEL3kGUqvSilfAl4Cb7VTH55D6YJcdXWYv/ySunffw3bgANqYGGLvu5fwWbPQJyf7O7x2ZXaYqbBUEBscS5guDCEEHumh0dHIjvIdLN6zmHxTPr0jehMdFI3dZcct3UzoMYFhscMYGjuUEF0Idpcdu9tOz/CehOpD/f2yFD/yJSEUArdKKQ8CCCEGAg8DTwAfAb4khBLgo+YEsEMI4QFigCofrqF0I43r11P3xhs0bd8Bbjf6tDQSFi4kYs7VaAK71yxXl8fF8iPLeWHvC5js3pvzoIAgggOCqbfXn5wVnB6RzjNTnuGy1MvUvAqlRXxJCBnfJQMAKeUhIcRwKWX+Bfyy/QeYCmwQQvQD9HhHLSnKKaTHQ9Xzz1PzzyXoUlKIvu02wqdfTuCAAd3yTe7rE1/z1I6nOFp/lNGJo7m699XU2mopbyo/WTbCaDDSI7QHE5MnnnN5RkX5IV9+Ww4LIf4JvNv8+KfAESFEIOA820lCiHfwLr0ZI4QoAR4FlgJLhRAHAAfwszM1Fyndm9tspuzhX2Nev56Ia68h4dFHu+2iM+VN5SzKWcSawjUkhybzbPazZKdkd8ukqLQdXxLCz4E7gfvx9iFsAR7Cmwyyz3aSlPK6s+y60YfnVroR6XBQ/5//ULPkXzgrKoh/5BGibri+y7/5SSmpslZR3FhMqbkUs8OMw+2gxlbDe4ffwyM93HXJXcwfPF8Vg1PahC8zla3A081fP2RutYiUbktKSf37y6lesgTXiRMYhgwh6am/EJyV5e/Q2txn+Z/x5+1/ptHReMb9U5Kn8OtRvyYlLKWdI1O6E9XAqHQIbnMTJ37/exo//5yg4cNJfPxxQiaM7/J3BS6Pi2d2PcMbh95gRNwIZqTPIDUsleSwZML0YRi0BvRaveoLUNqF+i1T/M6eX0DJPffgKCgg7uGHMd4yv8smgpLGEnaU78DpduKSLr48/iU7y3dyw4Ab+FXmr3xe61dRWtN5E4IQ4g0p5U1CiPuklM+2R1BK92DPL6DurbeoX7ECTWAgqUtfJWTMGH+H1SYKTYW8vP9lVuavPGWxmKCAIJ6c8CSze8/2Y3SK4tWSO4SRQoiewC1CiNfxdiifJKU8fcqjopyD7fARKv/2N+8MY52O8CtmEHf//V2q8JyUkrz6PLaVbWNr2Va+PvE1eo2e6zKu4yf9f0KYPgydRkdQQBB6bfccOaV0PC1JCEuA1UAvYBenJgTZvF1Rzks6HFS//DLVS/6FNjSUmHvvIeonP+mUtYbcHjcbijfw0dGPMDu8YyqEEFicFmpttdTZ6nB4vCuw9YroxS2Db+GGATcQE9T5XqvSfZw3IUgpnwOeE0L8U0p5RzvEpHQxrqoqzJu3UPvaa9gPHyb8qquI/8PvCTAa/R2aT5weJ3l1eew4sYN3D79LqbmUxJBEUsNSkUg80kNccBz9ovphNBhJj0hnbNJYEkIS/B26orSIL8NO7xBCDAO+WyJpk5RyX9uEpXR20u2m/qOPqHvnnZMVSHVJSSS/sJiwSy/1c3QtZ3Fa+PjYx6zMX8m3Nd+e/NQ/Im4ED2U+xJSUKWoEkNJl+LIewr14q45+1LzpLSHES1LK59skMqXTsuzZQ8X/PYnt4EEMgwYR+8ADhE6aSGBGRqcZPVRrq+WV/a+wIm8FZqeZAcYBXJdxHYNiBjE4ZrCaD6B0Sb58tLkNGC2lbAIQQjwFbANUQlCQDgeNGzdi+s/HmL/8koC4OJIWLSL8qis7TRL4TqGpkNu/uJ2Kpgqm9ZzGDQNvYFjsMH+HpShtzpeEIAD39x67+cGII6X78djtVL/wIvXvvYfbZEIbE0P0HbcTc9ttaEJC/B2ez/ZW7uWedfegERpev+J1hsQO8XdIitJufEkI/wa2CyFWND+eA7za+iEpnYUtN5eyh3+NPS+PsOnTibz2GkLGjUMEdJ42dbPDTG5tLscbj5Nfn8+7h98lPjieJZctISVcNQsp3YsvncrPCCE2ABPw3hnMl1LuaavAlI5LOp3ULFtG9XPPo4mMIOXllwidOPH8J3YQTo+Tr0q/4tNjn7KheMPJjuIATQCjEkbx/yb+P4yGzjUCSlFag08f5aSUu4HdbRSL0gk0bd9Bxf89gT3vKGHTLiPh8ccJiIryd1jn5XQ7+frE16wtWsu64nWY7CaiAqP4Ub8fMTF5Ij3De5IYkqhGDCndmvrtV85LSol1zx7q3nyThv+uQtejB8kvvkBodueox7+heAOPbX2MGlsNobpQpqRMYUbaDMb1GKdqBynK96iEoJyV22ym7u13MH30EY7CQkRwMDF33kH0ggVoDAZ/h3deFqeFRTmLWH5kOf2j+rNw3ELGJo1VpSIU5Sx8mYcQAlillJ7mZS8zgFVSyrOulqZ0TtLlon75cqqeX4y7tpbgzEyif/ELwmdM79Ajh5xuJxtLNpJvyqeooYhdFbsoM5cxf9B87h5+t0oEinIevtwhbAImCiGigC+BHLzLaN7QFoEp7c9dX0/D6jXUvvEGjmPHCM7MJO5f/yJoyGB/h3Zeeyv3snDbQo7WHwUgLjiO9Ih0Hh/3OKMSR/k5OkXpHHyahyCltAghbgWel1L+VQihRhl1Us6yMmyHDuGqrcVdW4ft4AHMGzYinU4C+/YlefHzhF56aYfvI7A4LTy7+1neyX2H+JB4ns1+ljGJYwjWBfs7NEXpdHxKCEKIsXjvCG5t6flCiKXATKBSSjm4edtjwC+AqubDfi+l/K8PsSgXwFFSSuPatTSsXoXtm1PLUAXExhJ1/XWEz56NYeDADp8IAPZX7ec3m39DSWMJ1w+4nnuG30OIruM2aSlKR+dLQrgf+B2wQkp5UAjRC1jfgvOWAYuB13+w/e9SykU+PL/iI4/djnX3bsybNmPetAnHsWMAGAYOJPbBBwkZO5aAmGi0RiOawM6zaLvL42LZwWW8sOcFYoJjWDp9KZkJmf4OS1E6PV8mpm0ENn7vcT5wbwvO2ySESLuQ4BTfuerqaPjkE8ybNmPJyUHa7QidjuCsLCJ//CPCsrPR9+zp7zB9UmWpYumBpeyu3E2VpYoaWw0e6WF62nT+OOaPRARG+DtERekSfBlltB7vgjinkFJOvcDnvlsIcTPezulfSSnrzvK8C/BWWSU1NfUCn6rrs+7dS+3bb9O4eg3S4UDfuzeRP/kJIePGEjJqVIceHXQ21dZq/n3g37x3+D1cHhejEkbRP6o/scGxDIoeRHZK55gHoSidhZDytPf4Mx8oxMjvPTQA1wIuKeWvW3BuGvDZ9/oQ4oFqvAnmCSBRSnnL+a6TmZkpc3JyWhRvd2E7coTKRYto2rQZTUgIEVdfTdR18wjs29ffoV2wgzUHefvbt1lVsAq3dDOz10xuH3q7qi2kKBdICLFLSnnedlVfmox2/WDTV0KIjWc8+PzXqvjueyHEy8BnF3Kd7shtMuE4XoyzpBjzli2YVvwHTUgIcQ8/RNS8eZ3yTgCgzlbHmsI1fHrsU/ZV7yMoIIhr+l7DjQNuJC0izd/hKUq34EuT0ferfWmAkcAFrQ0ohEiUUp5ofjgXOHAh1+kOpNuNde9ezBs2YN6wAXve0f/t1Okw3nQj0bff3inqCZ3JsfpjPLf7OTaVbMIlXfSJ7MNvsn7D1X2uJkwf5u/wFKVb8WWU0S68TTwCcAEF/G/46VkJId4BpgAxQogS4FFgihDikubrFQK/9CnqbsBtMlH/wQfUvvUWrrITEBBAcFYmsbNnE5ieji4lBX1KCprgzjne3uK0sGTfEt44+AZBuiBuGngTV/W6iv7G/v4OTVG6LV+ajNIv5AmklNedYbNaR+Es3I2NVC9+gbr330darQRnZRH/0EOETJqENjTU3+FdFI/0cLD6IOuK1/HJsU+otFQyp88cHhj5gCo3rSgdQEsmlk2VUq4TQlxzpv1Syo/OtF3xjZQS08cfU7noadw1NUTMno3x5z/DMGCAv0O7YHsr97Lu+DoqrZVUW6o5ZjpGtbUardCSmZDJosmLGB433N9hKorSrCV3CJOBdcCsM+yTgEoIF0G6XDR+8SU1S5di27cPw7ChpPzzn52iftDZFDcW8/ddf2dt0Vr0Gj2xwbHEBMWQlZDFxB4TmZQ8Sc0dUJQO6LwJQUr5aPO/89s+nK7P43DgLC7GUVSE/fBh6pd/gLOsDF1KColP/h8Rc+ciNBp/h3lBysxlvHbwNZYfWU6AJoA7h93Jzwb9TNUVUpROoiVNRg+ea7+U8pnWC6frkQ4H9sJCmjZvwbx5M5Zdu8D5v4rhQZkjif/977yLzWi1foz0wh2pO8KyA8tYVbAKgNl9ZnPXJXcRFxzn58gURfFFS5qMvhv71x/IAj5pfjwLb0ls5Xs8Fgt1779P46rVOMvKcFVXQ/Pkv8C+fTHedBOGjP7oe/ZE37Mn2shIP0d8YRxuB58Xfc77h99nT+UeggKCmJcxj5sH3kxiaKK/w1MU5QK0pMloIYAQ4nNghJSysfnxY8DyNo2uA5JS4iorw7p/P9Z9+/FYmtCnpKBLScFRWETtsmW4a2sxDB1KyKSJ6BKT0PXoQcjoUeiSkvwd/gWTUlJjq2Fb2TY2lWziq7KvaHQ0khqWykOZD3F176uJNHTO5KYoipcv8xBSAcf3HjuAtFaNpgOTbjcNK1dS/c8lOAoKABB6PSIoCI/JdPK4kAkTiLnzDoJHjPBXqK3CZDexumA164rXUWYuo8JSgdVlBSDaEM1lqZcxI30GYxLHoBGds89DUZRT+ZIQ3gB2CCFW4B1dNJfTS1p3elJKLDt30rh6DWi1aCMj0BgM1H/wIY6CAgIzMoj/4yMEDbsEQ7++CL0ed0MDjuPFCL0OQ79+/n4JF8zpdrKldAuf5n/KhuINOD1OekX0or+xP5OSJxEfHM/I+JEMiB6gkoCidEG+TEx7UgixGpjQvGm+lLLTr5gmXS6cZWU4jhdjO3jQu6B8UREiKAih1eIxmwFv+3+P554l7LLLThsFpA0PJ2jwIH+Ef9GqLFXsr97P1rKtrC5cjcluwmgw8tP+P2V279lkGDNURVFF6SZ8uUNASrlLCFGMt9opQohUKeXxNomsDUmPh6at26h79x3MGzedOupn5EgS77id8OnT0QQFIR0O3GYz2sjITjsc9DtSSooaithRvoOd5TvZU7mHCou3zqBBayA7NZuZvWYyNmksOo3Oz9EqitLefCluNxt4GkgCKvH2KeQCneKjsfR4sB08hHn9ekwrP8NZdBxtVBTG668jsF8/b22gtDR0cacOlRR6PQHGzltWwe62s+PEDjaWbGRTySZONHlrCsYFxTEyfiRDYocwJGYIGcYMDAEGP0erKIo/+XKH8AQwBvhCSjlcCJENnKlOUYdT9fxi6t9/H1dVFQhB8MiRxN59D2HTL0ej1/s7vFZndVnZXLKZtUVr2ViyEavLSlBAEGMSx3Dr4FsZnTianuE9VVOQoiin8CUhOKWUNUIIjRBCI6VcL4R4qs0ia0XS7SJoxAhCs6cQOnlypy0VfT4mu4k3Dr3Bm9++SZOzCaPByMxeM8lOyWZU4igCtZ1n3WRFUdqfLwmhXggRincy2ltCiEq8ZbA7vLj77/d3CG3K6XHy6v5Xef3g6zQ6G5nWcxrz+s9jRPwIAjQ+dRMpitKN+fJucTVgBR4AbgAigMfbIiil5SxOCw9ufJCvSr9iaspU7rzkTrWmgKIoF8SXYadNzd96gNeEEFpgHvBWWwSmnF+drY67vryLgzUHeWzsY1zb71p/h6QoSid23nGUQohwIcTvhBCLhRCXC6+7gXzgJ20fovJDHulhW9k2frb6ZxyuPcwzU55RyUBRlIvWkjuEN4A6YBtwG/AwoAeullLubcPYlGY2l40qSxWV1kr2Vu7lw7wPKW4sxmgw8q9p/yIzIdPfISqK0gW0JCH0klIOARBCvAJUA6nfFblTWp/L4+Kbqm+8cweKN3HMdOyU/SPjR3LXJXdxWc/L1MghRVFaTUsSwslpvFJKtxCiwJdkIIRYCswEKqWUg3+w7yHgb0CslLK6pdfsqo7UHeHjox+zMn8lNbYaAkQAIxNGMiN9BgkhCcQFx5EalkpyWLK/Q1UUpQtqSUIYJoRoaP5eAEHNjwUgpZTh5zl/GbCYHxTCE0KkANOATlf6wld1tjoqLZVUWiqptlZTa6ul3l5Pna2OOnsd9bZ6qq3VlDWVESACmJQ8iSt7Xcn4pPGE6kP9Hb6iKN1ES9ZDuKhlvKSUm4QQaWfY9Xfg18DHF3P9jqjB0cChmkNsKdnC5tLN5JvyTztGr9ETZYjCaDASZYgiOSyZm2Nv5sr0K4kydM2Jc4qidGx+mbXUXBepVEr5zfnKJwghFgALAFJTU9shOt+VN5Xz8dGP2XZiG4WmQmpsNQDoNDoy4zOZ02cOyWHJxAZ5F5s3GowEBQSp0hGKonQo7Z4QhBDBwB+Ay1tyvJTyJeAlgMzMTNmGobVIjbWGo/VHqbHWUGOr4auyr9hWtg2P9DA0ZiiTkieRFpFGn8g+ZMZnqgXmFUXpNPxxh9AbSAe+uztIBnYLIUZJKcv9EE+LeKSHd3Pf5R+7/3Fy5TCA+OB4bhtyG3P6zCElLMWPESqKolycdk8IUsr9wMka00KIQiCzo44ycnqcHG84zpPbn2Rn+U7GJ41n/uD5xAbFYjQYiQiMUE0/iqJ0CW2eEIQQ7wBTgBghRAnwqJTy1bZ+3vPxSA+7Knaxv3o/PcN70j+qP7HBsewq38WGkg1sK9tGlbWKJqe3YkeoLpSF4xYyt89clQAURemS2jwhSCnPuWaClDKtrWNwepzUWGuoslRRZa1iX9U+VhaspLzp1BYqgUAiMWgNjEocxYQeE4gIjCAqMIrJKZNJCElo61AVRVH8plvURn5g/QNsLNl48rFWaBmXNI4HRjzA2KSxFDcWc7juMGXmMi6JvYTRiaPV6mGKonQ73SIh/Kjfj5icMpnYoFhig2JJDksmIjDi5P4oQxRDY4f6MUJFURT/6xYJYUrKFH+HoCiK0uGdt/y1oiiK0j2ohKAoiqIAKiEoyv9v796Do76uA45/j57oCUIP0AMQD4F5GFtYNnYNBIc42MQDuE4z9jS1ncmMp9OmSZN4WnuYSSd1478yfcwkaeLxI07sOI2dTOu4rh2K7dImNkZgDNiIp0AIEJKRhN7SanX6x/kRCSIkEFrtavd8ZnZ297e/n/ZeaXXP73fv3XOdcwEPCM455wAPCM455wIJMcvIOecmXDgELSfg3GE4dxQyp0PRYihYBO1n4Ng7ULsdOpsgLctuqZmQnArJ6ZCUAheyIojArX8BuSURLbIHBOecu5QqdDVDagakXSZjcXcr1P4PnNkLOmDbwn3QfAw+OQwttTDQP/L75JTA9HnQ0Qh9nRDqskAS7rNjVQG1++X3e0Bwzrlx1d8HnY0wZRqkBysSdjTa2Xrtdmg8AJ8cgp5WkGSYsRRm3QI5xbatu9Ver68GDYMk2X5gZ/V55VC4CBbfAwULIb/CGv2uc9B0AJoOQkYezLsD8ucPXgXEAA8Izrn4pGpn66d2w6lqaNgHLceh7TQQLK2Slm2Boa3enqdPheLlsPReKKiA7hY4+T58+G/Q1w4pGZAxDXJLYdXXYcE6KLvZunlGk5UPhQsjVdtx4QHBORdfes5bA179rJ2Rg/XNz7we5q6BaXMgZ6bt13EWOj+xvv15n4LiGyFpmFWDB8LWlZMa3znOPCA45yan/l5orrVB2+ZjNoDbUgt171lffEklbPguzL4VChdD8jU0d0nJwweKOOMBwTkX28IhOL0H6n4H545Yw996AlrrBgdzwbp+8sph+RdgxUNQuiJqRZ6sPCA458ZP+1k4/r+QPQNKb7r8DJ0LQt3QsB8a9tpAbVMNnK+3vvr0bBusPbMH+jps/8wCyJsDJSvg+j+xAduCYNA2Y1rk6xfnPCA458auy8ZBfgAADNdJREFUtx1O7oDjv4Uj/20N+wVJKVB8gw26llTabSAM9Tvtdmq3BQAN2/6pWTY7Z+b11h3U2273NzwA5atgzu2QXRideiYIDwjOuSvX3wcn34Mj2+yLVQ3BHHxJhlkrYd23bDplZxPUvQt1O2D3T2DHDy/+ORl5dgVx3QYbyC1eDlNnxdQUzEQ0EWsqPwvcAzSq6rJg2xPAJmAAaAQeVtXTkS6Lc24MOprg8Jtw8L/g6NsQ6rSz/1krYfWjMOc2uwpIz7n4uIXr7T7cb91Bpz+wgdmym62Lxxv/mCOqGtk3EFkDdAA/GRIQclW1LXj8VWCJqv75aD+rqqpKq6urI1pe5yal7hZorIHMfMgqsDPw4RrcgQFIGiGFWVczHHoDTu2yGTwttXaP2rdqF90FC+6Euav/MAC4mCUiu1S1arT9In6FoKrbRaT8km1tQ55m8ftviTjnrkrbaXj3+7Drx4MDrwBJqZBbDLlllkOnoxHOn7R591lFkL8A8ufZzJzkVOvyqd9pA8ID/ZCeC9Pnwszl1oe/cL099rP6uBa1MQQR+Q7wIHAeuCNa5XBuUurvg23fhh0/sj78ZffZrbfd0jJ0nIW2MxYwzh2xWT/z10F2kb127gjUvG5BJByygd3p8+C2r8CSjTaLxxv/hBO1gKCqW4AtIvI48BXg74bbT0QeAR4BmD179sQV0LlY1XIcXv4SnN4NlX8Gax61+ffXYrSuJJcQYuET8DPgvsu9qKpPqWqVqlYVFvqUM5fAupptxs6P1lg65S/8FDZ979qDAXgwcECUrhBEpEJVDwdPNwI10SiHczEjHLKUC4WLrFvngnNH4cCvbaD35A7rHiqphM8/Z338zo2jiZh2+hKwFigQkXqsa2iDiCzCpp2eAEadYeRcXFKFmv+Erd+C5qO2rWgpzLrZZvo07LNtM5fD6m9CxXpLyZAAeXXcxJuIWUYPDLP5mUi/r3MxK9wPjR/brJ59r1iOnoJFcO9T0HbKFl3Z+wuYsQw++x0b5J3m42cu8vybys5NBFU4/n+w6zk4+IZ9uQssr/49/wSVDw5m41z9jeiV0yU0DwjORVLPedjzM9j5jKVpnjLVsnHOud26habN8emdLmYkRkCo22F9rqU3Df7z9XXaP+mRrZA3F4qWwIwlUFo1eoZG5y7V12XLL/ach5R0ux19y4JBXweU3QKbfwhLN9s6vc7FoMQICO88aYm48sph2ectre7vvgddn1ggaNgPu5+3fZPTLEfL/E/DigctDYBzw+lqhoOv26Dw0begv+fi15PT7PO28hGbGeRcjIt4LqPxNOZcRj3n4cBrsO9lG7DTAcvIuPYxW01J1b7a37DXAsexd+Dsfru8v2MLVH352lZbcvFhIGzpmuvetc9T7Xb7hm9uGVz3OcvcOXWWBYZQj+Xt9xMKFwOuNJdRYgSEodrPWiKwoutG3q+xBt74WwsORUth3lrrbhKxPDAL7oSppX94XF+XBZOmGnuvjrO2nN/C9bDwbkhJG9w33G+phGtet2ySF/LHF99oi35kFVjemZQp0NNq5VaFmcu822GihEM2E+iDFyxb54XB4OnzbfbP4o129u/jAC6GeUAYD6r2paBtfw/tZ+y5hge7BmYssy8S9ffayk/tDRcv+AGWPEySoLsZMqZbIxLqsQHGpkPQ125dC3M/ZeMcZz609xpJUootIjJrpR1Xvgqm5F5c7pEaqIGwLU7SVg85xYO3ocEqEsIhaDxgK2A111p33eyV0cuDr2oLrLfVB3l/TtnfNiPP/latJ6xrsa3epoXOW2vjUKU3Qf58DwJu0vCAECmq0HTQzugPb7VGJCXDBhEz8+0Mv6QSZiyF3BLbPhC2PPJ7XrCc8pn5dgWQX2GN+YJ1F6cSbj9rjVFHoy00Eu6zwJKRZ49P7QpWnaqG/m7LVFm83K44Os7a2EhWYbC84AJbdjB1il1pNB6AQ2/aPpfKzLfAkF1k5UnLgbQsC2gSpDbobrFjO5sssA2ErKFPmWJZNTOmA2pBrb0BuluDK6skC5wDoeDNhN8nuc0ptuNS0iw4ypAvXaVm2NKIGXlWppQM25aZb1dKhYtHD2SqlvTtzId2RXZypyV3O18P4d6Rj519G6z6BlTc6QHATVoeEGLVaGfvV6O/F06+D8fetvu0bGvML6Q7PnfEbt0tg4uRT5lq3V2L7oaChbZf+2k7Q+5osEa84yz0dtjsmL6O4MpIAbXAlFVgt9RMS52clGJXSN0tg91aucWQM9P2v1DvlDS7qiqptC9aNX5sM8BO7bL36e+1BnroZzLUFQShZttnoP/i30FSqnXdhXqs0Q91WcBIzbT7vk7obbv4uIJFNqNsahlMnW2BO7fY8v2nZlj3XFezBfMZS8fnb+VcFMXMegjuEuN5lpmSbguVzF098n6qdhbf323r1sbKAHnxDXa7GuF+q0d7g6V1aNgHrXV2JZOeYw16qNsCQ6jHphCn51ogLFoCZVUWMEeSMW18EsY5N8nESMvgIkrEzs4jPUYwEZJTIDnHGv+CClj2x9EukXNxw3PeOuecAzwgOOecC3hAcM45B3hAcM45F/CA4JxzDvCA4JxzLuABwTnnHOABwTnnXGBSpa4QkSbgxBgPLwCGSeCTELzuiSdR6w1e9+HqPkdVC0c7eFIFhGshItVXkssjHnndE6/uiVpv8LpfS929y8g55xzgAcE551wgkQLCU9EuQBR53RNPotYbvO5jljBjCM4550aWSFcIzjnnRuABwTnnHJAgAUFE7hKRgyJyREQei3Z5IkVEZonI2yJyQEQ+EpGvBduni8hWETkc3OdFu6yRIiLJIvKBiLwWPE+IuovINBF5RURqgr//bYlQdxH5evBZ3y8iL4nIlHitt4g8KyKNIrJ/yLbL1lVEHg/avIMisv5K3iPuA4KIJAPfB+4GlgAPiMiS6JYqYvqBb6rqYuBW4C+Duj4GbFPVCmBb8DxefQ04MOR5otT9X4A3VPU64AbsdxDXdReRUuCrQJWqLgOSgfuJ33r/GLjrkm3D1jX4v78fWBoc84OgLRxR3AcE4BbgiKoeU9U+4OfApiiXKSJU9Yyq7g4et2ONQilW3+eD3Z4HNkenhJElImXA54Cnh2yO+7qLSC6wBngGQFX7VLWVBKg7tgxwhoikAJnAaeK03qq6HWi+ZPPl6roJ+Lmq9qpqLXAEawtHlAgBoRQ4OeR5fbAtrolIOVAJ7ABmqOoZsKABFEWvZBH1z8DfAANDtiVC3ecBTcBzQXfZ0yKSRZzXXVVPAd8F6oAzwHlV/Q1xXu9LXK6uY2r3EiEgyDDb4nqurYhkA78E/lpV26JdnokgIvcAjaq6K9pliYIUYAXwr6paCXQSP90klxX0l28C5gIlQJaIfDG6pYoZY2r3EiEg1AOzhjwvwy4r45KIpGLB4EVV/VWw+ayIFAevFwON0SpfBN0ObBSR41i34KdF5AUSo+71QL2q7giev4IFiHiv+2eAWlVtUtUQ8Cvgj4j/eg91ubqOqd1LhICwE6gQkbkikoYNtLwa5TJFhIgI1o98QFX/cchLrwIPBY8fAv5jossWaar6uKqWqWo59jd+S1W/SGLUvQE4KSKLgk3rgI+J/7rXAbeKSGbw2V+HjZvFe72HulxdXwXuF5F0EZkLVADvj/rTVDXub8AG4BBwFNgS7fJEsJ6rsMvCvcCe4LYByMdmIBwO7qdHu6wR/j2sBV4LHidE3YEbgergb//vQF4i1B34NlAD7Ad+CqTHa72Bl7CxkhB2BfDlkeoKbAnavIPA3VfyHp66wjnnHJAYXUbOOeeugAcE55xzgAcE55xzAQ8IzjnnAA8IzjnnAinRLoBzsUJEwsC+IZs2q+rxKBXHuQnn006dC4hIh6pmj/B6iqr2T2SZnJtI3mXk3AhE5GEReVlEfg38RkSyRWSbiOwWkX0isinYrzxYi+DpIDf/iyLyGRH5bZCr/pZgv6wgr/3OIBFdXGbedZOTXyE4F7iky6hWVe8VkYeBfwCWq2rzhTTLqtomIgXAe1hagDlYiuFK4CMsZcqH2LdJNwJfUtXNIvIk8LGqviAi07B0ApWq2jlxNXVueD6G4NygblW9cZjtW1X1Qh56AZ4UkTVYmu1SYEbwWq2q7gMQkY+whUtURPYB5cE+n8WS8D0aPJ8CzObiRX2ciwoPCM6NbujZ+58ChcBNqhoKsqtOCV7rHbLfwJDnAwz+rwlwn6oejFxxnRsbH0Nw7upMxdZdCInIHVhX0dV4E/irIDsnIlI53gV0bqw8IDh3dV4EqkSkGrtaqLnK458AUoG9wWLpT4xz+ZwbMx9Uds45B/gVgnPOuYAHBOecc4AHBOeccwEPCM455wAPCM455wIeEJxzzgEeEJxzzgX+H4jZJL0dq00ZAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"labels = ['all', 'x-axis', 'y-axis', 'z-axis']\n",
"for col, label in zip(rog.results.T, labels):\n",
" plt.plot(col, label=label)\n",
"plt.legend()\n",
"plt.ylabel('Radius of gyration (Å)')\n",
"plt.xlabel('Frame');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also re-run the analysis with different frame selections. \n",
"\n",
"Below, we start from the 10th frame and take every 8th frame until the 80th. Note that the slice includes the `start` frame, but does not include the `stop` frame index (much like the actual `range()` function)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(10, 4)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rog_10 = AnalysisFromFunction(radgyr, u.trajectory, \n",
" protein, protein.masses, \n",
" total_mass=np.sum(protein.masses))\n",
"\n",
"rog_10.run(start=10, stop=80, step=7)\n",
"rog_10.results.shape"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xc1Zn4/89R712WVS3JliXLsiTbckFu2AFCC2BINiFAQnXIhkDCl2zahhQ2+SW7JLsh5Le8CGEJu5D8vllgNyEBQuIGrsi9SHKTiyzb6l0jjWae3x93NJK7xtZoJM3zfr300sy9VzOPxtZ5zj3n3OcaEUEppZT/CfB1AEoppXxDE4BSSvkpTQBKKeWnNAEopZSf0gSglFJ+KsjXAXgiKSlJsrOzfR2GUkqNK9u2bWsUkeRzt4+rBJCdnU1FRYWvw1BKqXHFGHPsQtt1CEgppfyUJgCllPJTmgCUUspPjas5gAux2+3U1tZis9l8HYpPhIWFkZGRQXBwsK9DUUqNM+M+AdTW1hIdHU12djbGGF+HM6pEhKamJmpra8nJyfF1OEqpcWbcDwHZbDYSExP9rvEHMMaQmJjot2c/SqmrM+4TAOCXjf8Af/7dlVJXZ9wPASml1ETU73BypLGLylPtVJ7q4IvXTiU2fGTn+jQBeMnARWtJSUlERUXR2dnp65CUUmNUS1cflaetht5q8Ns5WN9JX78TgOBAw41FkynNjBvR99UEoJRSo8ThFGoaO9l/qoMqV0NfeaqD0+2D83hJUSHMSI3h/vJsZqRGUzA5hqnJUYQEjfyIvSaAEXDHHXdw4sQJbDYbTzzxBKtWrfJ1SEopH2vrtrt69e1Uneqg8nQ71ac76HX16oMCDNMmRXHN1EQKJkczIzWGGakxJEeHjlqMEyoBfP+P+9hf1z6ir1mYFsN3PzHzkse8/PLLJCQk0NPTw7x587jrrrtGNAal1NjlcArHmrrOGr6pOt3BydYe9zEJkSHMSI3mvoVTmJEaQ0FqNNMmRREaFOjDyCdYAvCV5557jrfeeguAEydOcPDgQR9HpJTyhnabnapTHVS5evb7T3Vw4HQHPXYHAIEBhqnJkZRlx3Pv5CnMSI2m0NWrH4sr9ryeAIwxLwO3AvUiUuTaVgK8AEQBR4F7ROSqu+6X66l7w9q1a/nrX//Kpk2biIiI4Nprr9V1+UqNcyJCbUsP++ra2F/XTuVpq3df2zLYq4+LCGbG5Bjunp/FjFRrCGfapCjCgn3bq/fEaJwBvAI8D7w6ZNtLwFMiss4Y8yDwNeA7oxDLiGtrayM+Pp6IiAiqqqrYvHmzr0NSSnlgYGJ278l29p5sY19dO/vq2mi39QMQYCA3OYrSzDjunp9FoWusPiVmbPbqPeH1BCAi640x2edszgfWux6/D7zHOE0AN954Iy+88ALFxcXk5+ezcOFCX4eklLqI3n4HB890sq+ujb0nrYa+8tTgEE5oUAAFqTHcWpJGUVosM9NiyJ8cPa569Z7w1RzAXuA24H+BTwGZFzvQGLMKWAWQlZU1KsF5IjQ0lHfeeee87UePHnU/1msAlBp93X39VJ5qZ1+d1bPfe7Kdg/Ud2B0CQHRoEDPSrCGcmWkxFKXHMjU5kqDACVEgYVh8lQAeBJ4zxjwN/AHou9iBIvIi8CJAWVmZjE54SqnxpK3bbvXq69rcDf6Rxi7E1WIkRIYwMy2GZfm57p59VkIEAQHjewjnavkkAYhIFXADgDFmOnCLL+JQSo0/9e02dyO/r66dvXVtZ03OpsWGUZgWyycGhnHSY5gcEzbux+u9wScJwBgzSUTqjTEBwD9irQhSSim3gZU4Qxv6fXXtNHT0uo/JSYqkNDOOexZMoSg9hplpsSREhvgw6vFlNJaB/ha4FkgyxtQC3wWijDFfch3yJvAf3o5DKTX2HW3s4m9V9aytrmfXiVb3SpzAAEPepCiW5iW7G/oZqdFEh+mNkK7GaKwCuvsiu37u7fdWSo1tff1OKo42s7qqntVV9Rxp7AJg2qQov1mJ40t6JbBSalQ1dvaytrqB1VVn+OBAIx29/YQEBrBwaiKfu2YKKwpSyEqM8HWYfkETwBj0wgsvEBERwec+9zlfh6LUVRMR9tW1s6aqnr9V1bOrthURmBQdyi3FqawomMSiaUlEhmpzNNr0Ex+DHn30UV+HoNRV6e7rZ8OhJlZXnWFNVYO73HFJZhxfvW46KwomMTMtRlfm+Jj/XPHgJR999BHFxcXYbDa6urqYOXMme/fuPeuYP/7xjyxYsIDZs2dz3XXXcebMGQAef/xxfvCDHwDw3nvvsXTpUpxOJ9/73vd49tlnAavQXGFhIcXFxXzmM58Z3V9OKQ+caO7mNxuP8vmXt1L6g/d55NUK/rjrFLOz4viXTxbz0bev43+/tIjHP5ZHUXqsNv5jwMQ6A3jnG3B6z8i+5uRZcNOPL7p73rx53HbbbfzjP/4jPT093HvvvRQVFZ11zOLFi9m8eTPGGF566SX++Z//mZ/+9Kf8+Mc/Zt68eSxZsoTHH3+cP//5zwQEnJ2Tf/zjH1NTU0NoaCitra0j+7spdRX6HU62HWthdXU9qyvrOVhvXfGekxTJfQunsKJgEvOyE7xyIxM1MiZWAvCRp59+mnnz5hEWFsZzzz133v7a2lo+/elPc+rUKfr6+sjJyQEgIiKCX/3qVyxdupR//dd/ZerUqef9bHFxMffccw933HEHd9xxh9d/F6UupaWrj3UHGvhbVT3rqutpt/UTFGBYkJvAp+dlsqJgErnJUb4OUw3TxEoAl+ipe1NzczOdnZ3Y7XZsNhs/+tGP+NOf/gTAzp07+fKXv8yTTz7Jbbfdxtq1a/ne977n/tk9e/aQmJhIXV3dBV/7T3/6E+vXr+cPf/gDzzzzDPv27SMoaGL9s6mxS0SoPtPB3yrrWVNVz/bjLTjFum3hDTMn87GCSSzOS9L1+OOUtiQjYNWqVTzzzDPU1NTw9a9/neeff54f/vCH7v1tbW2kp6cD8Jvf/Ma9/dixY/z0pz9lx44d3Hzzzdxxxx0sWLDAvd/pdHLixAmWL1/O4sWLef311+ns7CQubmRvDK3UUDa7g42HG621+ZX11LVZE7hF6TE8tiKPFQWTKE6P9fs6OhOBJoCr9OqrrxIUFMRnP/tZHA4H5eXlrF69mhUrVriP+d73vsenPvUp0tPTWbhwITU1NYgIDz30EM8++yxpaWn8+te/5v777+ejjz5y/5zD4eDee++lra0NEeGrX/2qNv5qxIkIB+s7+eBgIx8ebGDTkSZsdicRIYEsmpbE4x/LY3nBJFJiwnwdqhphRmT8FNgsKyuTioqKs7ZVVlYyY8YMH0U0NuhnoDzV0NHLhkONrD/YwIZDjZxpt+rr5CZFsnR6MssLJrEgJ0Gvvp0gjDHbRKTs3O16BqCUH+jpc7D1aDMfHmzgg4ONVJ3uACA+IpjyaUksmZbE4rwkMuL1Clx/oglAqQnI6RT2n2q3hnUONfDR0Rb6+p2EBAYwd0o8/3BjPkumJTMzLUbH8v2YJgClJoi61h4+PNjIB4ca2XCokeYu6z5LBZOj+dzCKSzOS2J+TgIRIfpnryz6P0Gpcaqzt5/Nh5v44GADHxxq5EiDVUkzOTqUa6cnszgvicXTkpikk7fqIjQBKDVO9Duc7Kpt40PXsM6O4630O4Ww4AAW5CTy2flZLM5LIj8lWsssqGHRBKDUGCUiHGvq5oND1vLMjYeb6LD1YwwUpcXyyNJcluQlMXdKPKFBulpHeU4TwBik5aD9V2t3HxsPN/HBwUY+ONjgvtdtelw4t8xKZXFeEuVTk/S2h2pEaAIYg7QctP/o7Xew/VgrHx5q4MODjew+2YYIRIcGsXBqIquW5rJ4WhI5SZE6rKNGnCaAq/Sd73yHpKQknnjiCQC+/e1vk5KSwuOPP+4+5o9//CP/9E//RF9fH4mJibz22mvuY5KSknj66ad57733+OEPf8jatWv5wQ9+QFRUFE899RTPPfccL7zwAkFBQRQWFvK73/3OV7+qGgEiQtXpDtc4fiNba5rpsTsIDDCUZsbx+Io8luQlUZIZR3CgVtFU3jWhEsBPtv6EquaqEX3NgoQCvj7/6xfd/9BDD3HnnXfyxBNP4HQ6+d3vfsfWrVvPOkbLQfu3U209fHDQWpq54VAjjZ3W8sypyZH8XVkGi/OSWZiboAXV1KjzegIwxrwM3ArUi0iRa1sp8AIQBvQDfy8iWy/+KmNXdnY2iYmJ7NixgzNnzjB79mwSExPPOkbLQfuXdpudzYeb2HCo8azlmUlRISyaZi3NXDQtibS4cB9HqvzdaJwBvAI8D7w6ZNs/A98XkXeMMTe7nl97tW90qZ66Nz388MO88sornD59mgcffJBvf/vbWg7aj9gdTnYcb+VD12qdXbVtOJxCeHAgC3ITdHmmGrO83pKIyHpjTPa5m4EY1+NY4MKt3zixcuVKnn76aex2O6+//jo333yzloOewAaqZw6M42850kRXn4MAA8UZcXxx2VQW5yUxOytOl2eqMc1XXcmvAO8ZY57Fui9x+cUONMasAlYBZGVljU50HgoJCWH58uXExcURGHj+H7yWgx7/zrTb2HCo0d3o13dY1TNzkiJZOSedxdOSuSY3kdgIHcdX48eolIN2nQG8PWQO4DlgnYi8YYz5O2CViFx3udcZq+WgnU4nc+bM4fe//z15eXmj/v5j4TOYaDp7+9la0+SevD1wxrrfbUJkCOVTE1mSZ43ja/VMNR6MtXLQnweecD3+PfCSj+K4avv37+fWW29l5cqVPmn81cgYWmZhw6FGth9vod8phAYFMD8ngbvmZLBoWhKFqVo9U00cvkoAdcAyYC2wAjjooziuWmFhIUeOHPF1GOoKnGjuZu2BBtYfaGDz4SY6es8us7B4mlVmQW+Koiaq0VgG+lusFT5Jxpha4LvAI8DPjTFBgA3XGL9S3mSzO9hS08za6nrWHWhwL8/MiA/n1pI0Fk9LonxqIvFaZkH5idFYBXT3RXbN9fZ7K/8mItQ0drHuQANrqxvYfKSJ3n4noUEBLMxN5N4FU7g2P1nLLCi/pQvK1YTS1dvPpsNNVqN/oJ4TzVYxtdykSD67IItl05NZmJuowzpKoQlAjXMDa/IHhnU+qmmhz+EkIiSQ8qmJrFqSy7Lpk8hK1NU6Sp1LE8AYU1FRwauvvspzzz3n61DGrHabnY2HGllb3cC6Aw2carMBMD0livsXZbNsejJl2VojX6nL0QQwxpSVlVFWdt5yXb8mIuyra2fdgQbWVTew7XgLDqcQHRrEomlJPPGxZJZOT9baOkp5aNgJwBgzCVgEpAE9wF6gQkScXoptXHjhhRd44YUXAKvkQ3Z2NmvWrHHvP3r0KPfddx9dXdaKk+eff57y8nLeeustfvnLX/L+++9z+vRpli1bxvr166mqquLZZ5/l7bffZt26de4y08YY1q9fT3R09Oj/kj7Q2t3H+oONrKtuYP3BBhpcV97OTIvhC0tzuTZ/ErOztGSyUlfjsgnAGLMc+AaQAOwA6rGqeN4BTDXG/DfwUxFp92agw3H6Rz+it3Jky0GHzihg8re+ddH9jz76KI8++ih2u50VK1bw5JNPnrV/0qRJvP/++4SFhXHw4EHuvvtuKioqWLlyJW+88Qa//OUveffdd/n+97/P5MmTqaoajP/ZZ5/ll7/8JYsWLaKzs5OwsIl7c2+HU9hzss09lr/rRCtOgbiIYJbkJbNsejJL8/QG50qNpOGcAdwMPCIix8/d4VrHfytwPfDGCMc2rjzxxBOsWLGCT3ziE2dtt9vtPPbYY+zcuZPAwEAOHDjg3veLX/yCoqIiFi5cyN13n79adtGiRTz55JPcc8893HnnnWRkZHj99xhNDR29fHDQGsdff6CBlm47xlVQ7bEVeVybn0xJRhyBeuWtUl5x2QQgIl+7xO44EfmfEYznqlyqp+5Nr7zyCseOHeP555/nrbfe4vvf/z4AL730Em+//TYpKSns2rULp9N5Vi/+5MmTBAQEcObMGZxO53k3g/nGN77BLbfcwp///GcWLlzIX//6VwoKCkb1dxtJ9R02thxpZktNE1uONHOw3qqvkxgZwvL8SSzLT2ZJXrLe71apUeLxJLAxJhxr+OdeYDbWnIDf2rZtG88++ywffPABAQEBrFy5kpUrV7r3v/baa2RkZBAQEMBvfvMbHA4HAP39/TzwwAO8/vrrvPrqq/zsZz/jqaeeOuu1Dx8+zKxZs5g1axabNm2iqqpqXCWAutYed2O/taaZI43WPEhkSCBzsxO4Y3Y6S/OSmZmm9XWU8oWLJgBjTABQjDXZ6wRuwGr0l2HV8r8d+GAUYhzTnn/+eZqbm1m+fDlgreJ56aXB2nZ///d/z1133cXvf/97li9fTmRkJAA/+tGPWLJkCUuWLKG0tJR58+Zxyy23nPXa//Zv/8aaNWsIDAyksLCQm266afR+MQ+JCCeae9hc08TWGquXP3ARVnRYEPOzE/j0vEwW5CZSlBZDkE7eKjUszu5ubNXVhM+ciQkZ2bPji5aDNsb8X+A41hxAGLAF+B3wLlAlIjkjGskwjNVy0L7mi89goMzClppmthxpYktNs3s9flxEMPOzE1iQm8iCnARmpMboOL5Sw9Df0oJt/356Kyux7a/EVlVFX00NiJDz5huEFRZe0eteSTnodOC/gC8Du4E1wAci0muM8f5NBNSY4nQKhxo62XKkic011pDOwNLMpKgQFuQksiA3gQU5ieRNitIhHaUuQUSwnzxpNfZVVVZjX1lJ/5kz7mOC0lIJm1FIzM03EzajgODMzBGP41IJ4B7gLqAICAbuA7YYY6qBSGNMhIh0j3hEakxwOoXK0+3uSduPjrbQ3NUHwOSYMMqnJrob/VwtpqbURYndTu+RGmyVrp59ZRW2qiqc7a6V8wEBhOTmEDF/PmEzZhBWOIPQ/HyC4uO9HttFE4CIHAV+OmTTN4FvGmOWYs0FHDLGVIjIbd4N8fJExG8boJG6o1u/w8m+una2uMbwt9Y0027rB6xyycvzJ7EgN4GFOYlkJoT77eet1KUMjNfbKivdwzi9Bw8ifVbnyYSFEZo/nZibbhps7KdPJ8BH1/h4vApIRNYD640xoVjzAz4VFhZGU1MTiYmJftcoiQhNTU1XdIFYX7+TPSdbXWP4zWw71kJnr9Xg5yRFcvOsVBbkJjA/J5F0LbGg1Hn6m5utBr6q0j2E03f0KLg6ZYGxsYQWziD+nnsIK5xB2IwZhGRnY4LGTgWe4VwJfC/w+rklH0SkF3jLGDMVSBWRD70U4yVlZGRQW1tLQ0ODL97e58LCwi57gZiIcLK1h921ba6vVnYcb6XHbi1JzZsUxR2z05ifY03apujVtkq5DR2vt1VW0uuanL3geP0ttxA2o4CwGTMISk0d853S4aSiRGCHMWYbsA1owFoVNA1rSWgjVqkInwgODiYnZ9QXJI1p9e02dtW2sae2ld0n29hT20aTa/w+ONCQPznaWpKZk8D8nAQSo0J9HLFSY4ezq4vunTvp2baN7m3bse3fj7Ojw9oZEEDo1FyfjNd7w3CuBP65MeZ5rHv3LsK6NqAHqATuu1CJCDV6mrv62F3byp7aNnaftHr3Z9qt1TkBBqanRLOiYBLFmXEUp8dSkBqtZZKVGsLR2kr39u10V2yju6IC27594HBAQABhBQWuVTi+H6/3hmENRomIA3jf9aV8pN1mZ6+rod9T28au2lZqW3rc+3OTI7kmN5HijDiKM2IpTIshImTsjDcqNRbYz9TTs62C7ooKuiu20euqz2WCgwkrLibx4YeJKCsjfHYpgVFRPo7Wu7R1GKO6+/rZV9fO7oGhnNo2dykFgMyEcEoy47hv4RRmZcRSlB5LTFiwDyNWauwREey1tXR/5Grwt1VgP2YNWpiICCJKS4m56UYiysoIKy4mINS/hkO9ngCMMS9jVQytF5Ei17b/D8h3HRIHtIpIqbdjGat6+x1UnupwN/S7a9s4WN+B07XCc3JMGMUZsdw5J53ijDhmpccSrwXTlDqPOJ30HT5sNfYfVdC9bZt7sjYwNpbwuXOJ//RniJhXRlhBASbYvztNo3EG8ArwPPDqwAYR+fTAY2PMT4G2UYhjTLA7nBw40+Eawmljz8lWqk93YHdYrX1iZAjFGbF8vGgyJRmxzEqP1Rr4Sl2E9Pdjq6xyDedU0LNtG47WVgCCkpOJmFdGeFkZEWVlhE6bhgnQGlRDeXJHsFCsK4Ozh/6ciPzgUj8nIuuNMdkXeU0D/B3WBPOE0tLVR01TF0cbuzja1O363kX16Q56+60VtTFhQRRnxPHwklyrsc+IIy02bMwvHVPKV5y9vdj27HH38Ht27MDZbRUkCM7KImr5ciLKyoiYV0ZwZqb+LV2GJ2cA/4vVU98G9I7Q+y8BzojIwYsdYIxZBawCyMrKGqG3HRlt3XZqmro41tRFTaPV2Ne4Gvu2Hrv7uAAD6fHhZCdGct/CKe4VOVMSI/Q/qFKX4OjsomfnTrorPrJW6Oze476qNjQvj5jbb7Ma/LIyglNSfBzt+ONJAsgQkRtH+P3vBn57qQNE5EXgRbCqgY7w+19Wh83O0cbuwd58Y5f7cUv3YCNvDKTFhpOTFMknSlLJToy0vpIiyUwI16WXSg2D2O10b9tO5/r1dG/diq2y0lqSGRhIWGEh8Z/9rDWsM2fOuF17P5Z4kgA2GmNmiciekXhj1+0k7wTmjsTrXY3O3n73EM3Rxi5qGrs56urZN3b2nXVsamwY2YmR3FiUSk5SBNmJkeQkRZKZEEFYsDbySnmqv6GBzvUf0LluHV0bN+Ls7ITgYCJKSkh85GEiyuYRXlpKYFSkr0OdcDxJAIuB+40xNVhDQAYQESm+wve+Duu+ArVX+PMe6e7r56irYR8YrjnWZPXsB8oaD0iJCSU7MZLrZqSQnRTpbuSnJGojr9TVEqcT2969dK5bT+e6ddj27gUgaNIkYm66kcilS4m8plwb/FHgSQK4ottRGWN+C1wLJBljaoHvisivgc9wmeGfkfLY69t5e/eps7YlR4eSkxjJ8vxkspMiyUmMZEpiJNlJEXrxlFIjzNHeTteGDVaj/8EHOJqawBjCS0pIfuJxopYtI3TGDJ0TG2XDbulE5JgxpgRr4hasm8PsGsbP3X2R7fcP972v1oqCSRRMjnb35rOTIokK1UZeKW8REfoOH6Zz3To6166je/t2cDgIiI0latEioq5dRuSSJTqO72OeLAN9AngEeNO16b+MMS+KyC+8EtkIunPOpatlKqWuntNmo3vLFqvRX7ce+8mTAITm55P44INEXbuM8JKSMVUO2d958i/xELBARLoAjDE/ATYBYz4BKKW8w37yJB3r1tG5bh3dm7cgvb2Y8HAiFy4k8ZFHiFq2lODUVF+HqS7CkwRgAMeQ5w7XNqWUnxC7ne4dO6wVO+vX03vwEADBmZnEfepTRC1bRsT8eX5XU2e88iQB/AfWPYHfcj2/A/j1yIeklBpL+puarGWa69fR9eEGqzZ+UBARZWVMuvMuopYtJSQnRydwxyFPJoF/ZoxZi7Uc1AAPiMgObwWmlPINcTqx7a+kc91aOtetx7ZnD4gQmJxE9A3XE7VsGZHl5RO+VLI/GM4tIWNEpN0YkwAcdX0N7EsQkWbvhaeUGg39jY10bdpE14aNdG74EEdDIxhDWPEskr78GFFLlxFWOEOLqU0wwzkDeB2rnPM2YGgpBuN6nuuFuJRSXuTs6aG7YhtdGzfStXEjvdXVgFUyOaL8GqKWLSNqyRKCEhN9HKnypuHcEvJW13e98a5S49TAsM5Ag9+zbRtit2OCgwmfM4fkr36VyPJyq5cfqFe7+wtPrgP4m4h87HLblFJjg/3kSTpdDX73ps3uOvmh06cTf889RC4qJ2LuXAIiInwcqfKV4cwBhAERWKUc4hlc+hkDpHkxNqWUBxzt7XRt2WI1+Bs30XfsGGDV2Im69loiF5UTuXAhQcnJPo5UjRXDOQP4AvAVrMZ+G4MJoB34pZfiUkpdhtjt9OzaZQ3rbNhIz5494HRiIiKInDeP+Hs+S2R5OSFTp+oSTXVBw5kD+Dnwc2PMl8dD2QelJioRoe/IEbo2uIZ1tm617oYVEED4rFkkfmEVUeXlVrmFEL1ntLo8T64D+IUxpggoBMKGbH/14j+llLoa/U1NdG3c5J68HbjBefCULGJuv43I8nIiFywgMCbGx5Gq8ciTSeDvYpV1LgT+jFUe+kOG3OxdKXV1Lrk885priCy/hsjyRYRkpPs4UjUReFIK4pNACbBDRB4wxqQAL3knLKX8w3nLM7dvR/r6rOWZc+fq8kw/IyI02Zqoaauhpq2GI21H3N9fvuFlMmMyR/T9PEkAPSLiNMb0G2NigHr0IjClPHbR5Zn5+dbyzPJyIsrmEhAe7uNIlbc4nA5Odp48r5Gvaauhva/dfVx4UDg5sTmUpZQhjPwt0T1JABXGmDjgV1irgTqBrSMekVITzLCWZ15zDUFJST6OVI20nv4ejrUf40jrEWraazjSeoQjbUc43n6cPufg/cYTwxLJjcvlxuwbyY3LJSc2h9zYXFIiUry6gmtYCcBYEfw/ItIKvGCMeReIEZHdXotMqXHqYsszAyIiiJg/39XLv0aXZ04gLbaW83ryNW011HXWuXvuASaA9Kh0cmNzWZK+hJzYHPdXbGisT+IeVgIQETHG/A8w1/X8qDeDUmo8Gbj9oTWOv+m85ZlJj36ByPJywouLdXnmOOYUJ6e6Trl78QONfE1bDS29Le7jwgLDyI7NpjipmNun3U5urNWjnxIzhdDAsXWfBE+GgDYbY+aJyEdei0apcWJo9cyuTZvcyzNDpkwh9o7brRU7ujxzXBIR6rrqqG6u5kDLAffwzdG2o9gcNvdxcaFx5MbmsiJrhXvIJjcul9TIVALM+Kia6kkCWA58wRhzDOjCVQ1URIov9UPGmJexqonWi0jRkO1fBh4D+oE/icg/eBq8UqPFWp5Z4b4Iq/fAAQAC4+KIuGYhkeXlRJWXE5yuyzPHE1u/jUOth6hurqa6pZrq5moOthykw+LYLGYAACAASURBVN7hPiY9Kp3s2GzmTZ7n7s3nxuYSHzb+b2jvSQK46Qrf4xXgeYZcL2CMWQ7cDhSLSK8xZtIVvrZSXiEOx/nLMweqZ86dS/KTTw4uz9Qa+WOeiNDQ0+Bu6A80H6C6pZqj7UdxihOwVtxMj5/OTTk3kZ+Qz/T46UyPn05E8MQtludJAvgy8LKI7PfkDURkvTEm+5zNXwR+LCK9rmPqPXlNpbyhr7Z2sMzC5s042toACC0oIP6++6zlmXPn6PLMMc7usHOk7Yi7Rz/Q4A8dp0+LTGN6wnSun3I9+Qn55MfnkxGdMW6GbkaKJwmgCviVMSYI6/7AvxWRtit83+nAEmPMDwEb8JTOLajRJiJ0b/2I9nf+TNfGTdiPHwcgKCWFqBUrrDIL1yzU5ZljWIutxd3QH2g5QHVzNYfbDtPv7AcgJCCEafHTuDbz2rN69b5adTPWeFIL6CXgJWNMPvAAsNsYswH4lYisuYL3jQcWAvOA/2uMyRWR8650MMasAlYBZGVlefg2Sp3Pfvo0bf/zP7S++Rb248et5ZkLFpBw771ELionJDdXl2eOMQ6ng2Mdx9xDN9XN1ld9z+DgQVJ4EvkJ+SxKX0R+fD75CflMiZlCUIAn/Vz/4tEnY4wJBApcX43ALuBJY8wXROQzHrxULfCmq8HfaoxxAklAw7kHisiLwIsAZWVlI38pnPIL0tdHx+o1tL75Bl0fbgCnk4gFC0h+7EtEX3+9DuuMId32bqpbqqlsqnT36g+1HnKvwAkyQeTE5TA/dT758flMT5hOfnw+ieF6+0pPeVIM7mfAbcDfgB+JyMBVwD8xxlR7+L7/A6wA1hpjpgMhWAlFqRFlO3CAtjfeoO0Pf8TR0kLQ5MkkfmEVcStXEqJnlD430Njvb9rPvsZ97G/aT017jXtiNjY0lvz4fD6V/yl3rz43NpeQQL2eYiR4cgawF/hHEem+wL75F/shY8xvsaqIJhljaoHvAi8DLxtj9gJ9wOcvNPyj1JVwdHTQ/qc/0/rGG9j27IHgYKJXrCDuk3cRWV6uRdV8pKe/h+rmavY1WQ39/qb9HGk74m7sk8KTmJk4kxuyb6AwsZCChAKvl0Lwd54kgJ1AwTn/GG3AsUtNBovI3RfZda8H763UJQ1M6La9+Qbt7/0FsdkIzcsj5ZvfIOa22wiKH/9rtseTgcZ+f9N+d4M/tLFPDEtkZtJMrptyHYUJhcxMmsmkCF0NPto8SQD/LzAH2I11EViR63GiMeZREfmLF+JT6pLsZ87Q9tZbgxO6UVHE3nE7cXfdRVhRkfYeR4Gt33b2ME7zfo60HsEhDgASwhKYmTiTj2V9jMLEQmYmWo29/tv4nicJ4CjwkIjsAzDGFAJfA54B3gQ0AahRIX19dKxZS+sb/60TuqPM1m/jQMuBs3r2h1sPn9XYFyYWsjxzOTMTZ1KYWKjDOGOYJwmgYKDxBxCR/caY2SJyRP9x1WiwJnTfpO0Pf7AmdFNSdELXi3odvRxoPruxP9R6yN3Yx4fGU5hUyLKMZcxMmsnMxJna2I8zniSAamPMvwO/cz3/NHDAGBMK2Ec8MqUYMqH75pvYdu/WCV0v6bJ3cbDlIFXNVVQ2V1qNfcsh+sW6oCouNI6ZiTNZmrHU3bOfHDlZG/txzpMEcD/w98BXsOYAPgSewmr8l494ZMpv6YSu94gIZ7rPUNVcdVYBtOMdx93HxIbGMjNxJvcX3e9u7FMjU7Wxn4A8uRK4B/ip6+tcnSMWkfJbOqE7suwOO4fbDlPdXE1VcxUHWg5Q1Vx11i0Hs6KzyE/I57apt5GfkK9LL/2MXiOtfOqCE7rz5+uErodaba1Ut5zd0B9pO+KuiRMWGEZefB43ZN9AfrzV0OfF5xEZHOnjyJUvaQJQo05EsO3dS/vbbw9eoZuSQuKqR4i7806d0L0Epzg50XHivCGcM91n3MckhyczPWE6S9KXWJUuE/KZEj2FwACdL1Fnu2wCMMb8p4jcZ4x5QkR+PhpBqYlHRLDt30/Hu+/S/s672GtrByd077qTyEWLdEL3HN32bg62HnQXPqtusSpe9vT3ABBoAsmJzaFscpm7TILWxFGeGM4ZwFxjzBTgQWPMq1gTwG4i0uyVyNS4JyL0VlfT/s67tL/7DvZjxyEoiMhrriHpi18k+rqPERirZXkB2vva2Vm/86xe/bH2Y+4bikcFR5GfkM/KaSspSChgesJ0psVNG3P3mFXjy3ASwAvAu0AusI2zE4C4tivl1nvwIO3vvEP7O+/SV1MDgYFELlhA4sMPE33ddbqKx+V012nWnFjD6uOrqThd4V5ymR6VTn58Pjfn3Mz0hOkUJBSQFpmmE7NqxF02AYjIc8Bzxph/F5EvjkJMahzqPXLE1ei/Q9+hwxAQQMS8eSR8/vNE33A9QQkJvg7R50SEQ62H3I3+vibrusrsmGw+N/NzLE5fTEFCAdEh0T6OVPkLT5aBftEYUwIscW1aLyK7vROWGg/6jh6l3TWm31tdDcYQMXcu8U9/h5jrrycoOdnXIfqcw+lgV8MuVh9fzeoTqznRcQKA4uRinpjzBCuyVpAbqyfRyjc8uR/A41h35nrTtek1Y8yLIvILr0SmxqS+EyfcY/q9+ysBCJ8zh5RvfYvoj3+c4BSt6Gjrt7H51GbWnFjD2hNrabY1ExwQzPzU+dw/836WZy4nOUKTo/I9T5aBPgwsEJEuAGPMT4BNgCaACc5+8qS7p2/buxeA8JISJn3j68TceCPBkyf7OELfa+ttY33telYfX82Gug309PcQFRzFkowlrMhaweK0xUSFRPk6TKXO4kkCMIBjyHMH56wIUhOH/dQp2t97j/Z33sG2yxrpC5s1i0lf+xoxN36c4PR0H0foe3Wdde7x/G1ntuEQB5PCJ3Hb1NtYkbmCeZPnERwY7OswlbooTxLAfwBbjDFvuZ7fAfx65ENSvmI/U0+Hq9Hv2bEDgNDCGST/nyeJufFGQjIzfRyhb4kIB1oOsPrEatYcX0NlszUENjV2Kg8WPciKrBUUJhYSYAJ8HKlSw+PJJPDPjDFrgcVYPf8HRGSHtwJTo6O/oYH2v/yFjnfepXvbNhAhND+f5K88YTX62dm+DtGn+p397Kjfwerjq1lzYg0nO09iMJROKuXJuU+yPHM52bHZvg5TqSviUSkIEdkObPdSLGqU9Dc10fH++7S/8y7dH30ETiehedNIeuxLxNx0E6G5/r0qpae/h011m1h9fDXratfR2ttKSEAIC9MW8sisR1iWuYyk8CRfh6nUVdNaQBOcOBz0HjpMz86d7q++mhoAQnJySHr0UWJuupHQvDwfR+pbLbYW1tWuY/Xx1Wyq24TNYSM6JJplGctYkbWCRWmLiAiO8HWYSo0oTQATjKO9nZ5du+jZ4Wrwd+/G2WlV6w6MiyN89mxi77iDqGuXETp9ut9eXSoi1LTVsL52Petq17G9fjtOcTI5cjJ35t3J8qzlzE2ZS3CATuKqicuT6wAigR4RcRpjpgMFwDsicsm7gRljXgZuBepFpMi17XvAI0CD67BvicifryB+vyZOJ31HjtCzcyfdO3fSs2MnfYcPWzsDAgidPp2YW28hvLSUiNJSgqdM8dsGH6xbHFacrmB97XrW166ntrMWgLz4PB6Z9QgrslYwI2GGX39Gyr94cgawHlhijIkH/gZUYN0W8p7L/NwrwPPAq+ds/1cRedaD9/d7js5Oq3fvaux7du/G2W7d3CMwNpaw0hJib72F8NmzCSuaRWCU1no/03WG9SetBn/LqS309PcQFhjGgtQFPFD0AEvSl5AalerrMJXyCY+uAxCRbmPMQ8AvROSfjTGXXQUkIuuNMdlXGqC/EhH6ao7Ss2OHe+y+99AhEAFjCJ02jZgbbyS8tJTw0lJCcrK154pVemFP4x7W167ng5MfUNVcBUBaZBq3Tb2NpRlLmT95PmFBYT6OVCnf8ygBGGOuwerxP3QFP3+ux4wxn8M6k/g/ItJykTddhVWCgqwJfKMQR2cXtj273cM5tp27cLS1ARAQE0N4SQnRN37cavCLiwmM1oJhA9r72tl4ciPratex4eQGWnpbCDSBlCSX8JU5X2FZxjKmxk3VBKnUOTxpwL8CfBN4S0T2GWNygTVX+L7/DjyDVU76Gaz7DD94oQNF5EXgRYCysjK5wvcbU0QE+7Fj1rj9zp307NxF74ED4HQCEDJtKlHXX0fEQO8+NxcToBcXDRARDrcedg/t7KzfiUMcxIXGsTh9MUszllKeVk5sqN5rQKlL8eRCsHXAuiHPjwCPX8mbioj7/nXGmF8Bb1/J64wXzp4ebHv30r1jp3tIx9FinfAEREURXlxM9KOPEj57NuElxQTGxPg44rHH1m/jo9MfuYd2TnaeBCA/Pp8Hix5kacZSZiXN0tseKuUBT1YBrQHO64GLyApP39QYkyoip1xPVwJ7PX2Nscx++jQ9O3bQvWMHPTt2YqushH7rZh8hOTlELV9OeGkJ4aWlhE6dqrdCvIjTXafdK3a2nNqCzWEjPCicBZMX8NCsh1iSvoTJkVqITqkr5ckQ0FNDHocBdwH9l/shY8xvgWuBJGNMLfBd4FpjTClWQjkKfMGDOMYUsduxVVW7evY76N6xk/5TVm4zYWGEFxeT+NBDhM8uJbykRO+GdQkOp4Pdjbvdjf6BlgOAdYeslXkrWZqxlHmT5+ltEJUaIUbkyofVjTHrRGTZCMZzSWVlZVJRUTFab3dB/S0tgxda7dhBz549SI91k+6g1FQiZs+2hnJmzyYsfzomWC8kupS23jY2nNxgTeDWbaCtt41AE8jsSbNZmrGUpRlLyY3N1Qlcpa6CMWabiJSdu92TIaCh9/QLAOYCE/r8W5xO+mpqzhrO6TtyxNoZFETYjBnEfeqTVqNfWkpwqq4nv5xmWzO7G3azs34n2+u3s6thF05xEh8az7KMZSzJWEJ5WjkxIToPopS3eTIEtA1ryMZgDf3UMLgcdEJwdnfTs2cvPTu2Ww3+zl04XUsxh5ZRiJhdSlhREQHh4T6OeGxzOB0caj3EroZd7GrYxc76nRzvOA5AkAmiIKGAh2c9zNKMpRQlFukErlKjzJNVQDneDGS0iQj9p065e/Y9O3Zgq6oCh3XPm5BpU4m54XrCS63hHL3Q6vLaetvY3bDbauwbdrK3cS9d9i4AEsISKEku4a7pd1GSXEJhYiHhQZpAlfKlyyYAY8wKEVltjLnzQvtF5M0LbR9rpK8PW1WVazjHavD7z1irUU1EhDVZu+oRazinpITAWF1DfilOcVLTVuPu2e9q2MWRNmt4LMAEkB+fz625t1KSXELppFIyojI0gSo1xgznDGAZsBr4xAX2CYM3iR+zTn3/+7S9+RbS2wtAcFoaEfPmET67lIjZs62qmEFaGPVSOvs62d242z2cs7thNx19HQDEhsZSklzCLbm3UJpcSlFSkZZOVmocuGyrJyLfdX1/wPvheEdoTi7xn/mMa3VOKcEpKb4OaUwTEY61Hxscu2/YyaGWQwiCwTA1biofz/44JckllCSXkB2jw2NKjUfDGQJ68lL7ReRnIxeOdyR87j5fhzCmddu72de0zz2Us6thF629rQBEB0dTnFzM9VnXUzKphFlJs4gO0TpESk0Ewxn3GPhrzwfmAX9wPf8EVoloNc7UddZZSzDrrcb+QMsBHGJNfufE5nBt5rWUJpdSklxCblyu3uRcqQlqOENA3wcwxvwFmCMiHa7n3wN+79Xo1IjotndTcaaCD09+yMa6jRxrPwZARFAEs5Jn8dCsh9zDOVpATSn/4cnMZxbQN+R5H5A9otGoESEiHGg5wMa6jWw4uYHt9duxO+2EBYYxb/I87i64m7KUMqbFTdO190r5MU8SwH8CW40xb2Gt/lnJ+Xf5Uj7Samtl06lNbDi5gY11G2nose62mRefxz0z7qE8rZw5KXO0jo5Sys2TC8F+aIx5F1js2vSAiFz2jmDKO/qd/exp3ONu8Pc27kUQYkNjuSb1GsrTyilPKyclUlc8KaUuzKPF7yKyzRhzAqsaKMaYLBE57pXI1HlOdZ5iQ53V4G+u20yHvYMAE0BxUjFfLP0ii9IWMTNxpg7rKKWGxZNicLdh3bkrDajHmhOoAmZ6JzRl67dRcaaCDSc3sKFuAzVtNQCkRKRwQ/YNlKeVsyB1gU7cKqWuiCdnAM8AC4G/ishsY8xy4G7vhOWfBm51ONDLrzhdQZ+zj9DAUOamzOWTeZ9kUfoiLY+slBoRniQAu4g0GWMCjDEBIrLGGPMTr0XmJ9p629h8arN7xc6Zbqs+UW5sLp8u+DSL0hYxN2UuYUFhPo5UKTXReJIAWo0xUVgXf71mjKlnGHcEU2dzOB3sbdrLxpMb2VC3gT2Ne3CKk+jgaBamLWRR2iLK08pJjdJ7CyilvMuTBHA70AN8FbgHiAV+4I2gJoJ+Zz+nOk9xouMExzuOc7zjOCfaT7C9fjvtfe0YDEVJRTwy6xEWpy+mKKmIoAAtSKeUGj2eLAPtcj10Ar8xxgQCnwFe80Zg40Gvo5fajlqrkW8/zomOE+4G/1TnKfpl8AQpLDCMjOgMlmcuZ1H6Iq5JvYa4sDgfRq+U8nfDKQYXA3wJSMeqA/S+6/nXgJ1M8ATQZe9yN/DHO45T21Fr9eY7TnCm6wzC4D2Vo4OjyYzJZGbiTG7MvpHM6EyyYrLIjM4kOTxZJ26VUmPKcM4A/hNoATYBD2M1/CHA7SKy83I/bIx5GbgVqBeRonP2PQX8C5AsIo0exj4iRIS23rbBYZqOE5xoH+zJN9uazzo+ISyBrOgs5k+eT0Z0BlnRVgOfFZ1FbGisNvJKqXFjOAkgV0RmARhjXgIagayBonDD8ArwPOeUjTDGZALXA16/kExEaOhpOG+YZqCx77AP/ioGQ0pkClnRWSzPXH5WLz4zOpPI4Ehvh6uUUqNiOAnAPvBARBzGmBoPGn9EZL0xJvsCu/4V+Afgf4f7Wlfq8dWPs7Z2rft5oAkkPSqdzOhMinOL3Q18VnQW6dHpWi9HKeUXhpMASowx7a7HBgh3PTeAiEiMp2/quqr4pIjsutyQiTFmFbAKICsry9O3AuD2abezKH2RNVwTk0lqZKquuFFK+b3h3A9gRAvLGGMigG8DNwzneBF5EXgRoKysTC5z+AVdN+W6K/kxpZSa0Hxxq6epQA6wyxhzFMgAthtjJvsgFqWU8lujPg4iInuASQPPXUmgzFergJRSyl95/QzAGPNbrCWk+caYWmPMQ95+T6WUUpfn9TMAEblkxVARyfZ2DEoppc7nizkApZRSY4AmAKWU8lOaAJRSyk9pAlBKKT+lCUAppfyUJgCllPJTmgCUUspPaQJQSik/pQlAKaX8lCYApZTyU5oAlFLKT2kCUEopP6UJQCml/JQmAKWU8lOaAJRSyk9pAlBKKT+lCUAppfzUqN8TWCml1AU4+qG7CboaXF+N0FU/+Hz5P0JM6oi+pSYApZTyBhHobXc15A3nNOwN0Fl/9r6e5gu/TkAQRCbDgi9qAlBKKZ/p74PuxvMb76EN+9DHjt4Lv05YnNWoRyZDcj5kL3Y9T4KoSYP7IpOsY43xyq/j9QRgjHkZuBWoF5Ei17ZngNsBJ1AP3C8idd6ORSk1joiAsx8cdnD0DT522l3bBh73WcMn5z22n3P8pV7jIo/7+6ye+UDDbmu7cKyBoa6G29WAp8y0Hrsb8iFfEYkQFDK6n+VFjMYZwCvA88CrQ7b9i4h8B8AY8zjwNPDoKMSilBqL+rrh5DY4scX6qv0Ielq8/74mEAKDISDY+n7u48AQCI+HycVn98ojkwcb/MhkCInyWi/dm7yeAERkvTEm+5xt7UOeRgLi7TiUUmNIex0c3wwntsKJzXB6j9U7B0jKh4JbIDrNaoADg1yN8nAfB1vj5pd9HAwB/r0Q0mdzAMaYHwKfA9qA5Zc4bhWwCiArK2t0glNKjRxHP5zZO9jYn9gKbSesfUHhkD4Xyh+HrIWQMQ8iEnwbrx8xIt7vfLvOAN4emAM4Z983gTAR+e7lXqesrEwqKipGPkCl1MjpaYXaCtdwzmao3Qb2LmtfdCpkLrAa+8z51tBKYLBv4/UDxphtIlJ27vaxsArodeBPwGUTgFJqjBGB5iNn9+7rKwEBEwApRTD7HqvRz5wPsZnjcqx8ovJJAjDG5InIQdfT24AqX8ShlPJQfy/U7RycrD2xxVodAxAaYw3hzFxpNfbpZRAa5dt41SWNxjLQ3wLXAknGmFqsnv7Nxph8rGWgx9AVQEqNTZ0Ng0M5J7ZC3Q5rOSVAfA5M/RhkLbB6+MkFEBDo23iVR0ZjFdDdF9j8a2+/r1LKQ45+aKiylmAO9O6bj1j7AkMgtRQWfMFq7DPmQ3SKb+NVV20szAEopUab0wnNh60e/cntULcdTu2G/h5rf0SS1dDPvd/6nloKwWE+DVmNPE0ASk10ItB63Grs67ZbDf6pXVadGrCWYqYWW419+hxrWWZCrk7W+gFNAEpNNB1nBhv6gUa/u8naFxBslSmY9UlImw1pc6yx+0BtCvyR/qsrNZ51Nw828nU7rUa/w1VWywRYjfv0myCt1OrdpxRBUKhvY1ZjhiYApcaL3g5r6GZoz77l6OD+hKmQvWiwZz95li7DVJfkHwngL9+BYxussc6gUAgOh6Awa1IryPUV7NoXFO7B9iHP/bymiBphdptVH2fouH3jAdxls2IzrV79nM9bPfvUEqtomVIe8I8EEBZj1dTut1kVBjtOW6sd7DZrW78N7D1cVU26wJALJImw85NOUKjrWNd39+Ngq6TsedsHHocM2T9wrGuf+7HreE1GY1d/n1UWoa8b7N3Q1zX4vf3k4Iqc+srB4miRyVaPvuhO63taqVWJUqmr5B8JYOnXLn+MiFUD/EKJob93GNtdz+09Q44Z2O5KPP026xhHn+u73bphxMCFNSMlIMiVDIJdCWNosggZHAMWp/V7i3PwCzln24X2D91+sf0D+y6zfyApus/Kzj3jOves7TJnZO5tlzhju1yCdNitBnlo42zvdjXarsa7r/P8bQPHXvTnugYb9YsJi7WGcMoft76nz4GYdF2Ro7zCPxLAcBhjNZJBIdYf4WgaSD6OXquH6Ogb8nggYQw8truSx9Bj+wYTy3Ae9/dav68JAFzfTcDgtoHHZ+07d//Az17s54e+9rnfB/bjSqK2wcTpTqw2sLVDf/2Fk604r/zzPvdsLSDIeo+Bhtxp9+z1gsIgJBKCIyEkAoIjrOfRqdb3kIgh+1zfL3R8ZJJ1da029mqUaAIYC4YmH12gcXkDd4o6K2mcc9Z13tnaJc7aHHZXIzykMQ6JHHzs3hd5gQY9QssfqHFLE4Aaf4wZvGMTMb6ORqlxS2cLlVLKT2kCUEopP6UJQCml/JQmAKWU8lOaAJRSyk9pAlBKKT+lCUAppfyUJgCllPJTRuQqCqCNMmNMA9ZN5K9EEtA4guGMd/p5DNLP4mz6eZxtInweU0Qk+dyN4yoBXA1jTIWIlPk6jrFCP49B+lmcTT+Ps03kz0OHgJRSyk9pAlBKKT/lTwngRV8HMMbo5zFIP4uz6edxtgn7efjNHIBSSqmz+dMZgFJKqSE0ASillJ/yiwRgjLnRGFNtjDlkjPmGr+PxFWNMpjFmjTGm0hizzxjzhK9jGguMMYHGmB3GmLd9HYuvGWPijDH/bYypcv0/ucbXMfmKMearrr+TvcaY3xpjwnwd00ib8AnAGBMI/BK4CSgE7jbGFPo2Kp/pB/6PiMwAFgJf8uPPYqgngEpfBzFG/Bx4V0QKgBL89HMxxqQDjwNlIlIEBAKf8W1UI2/CJwBgPnBIRI6ISB/wO+B2H8fkEyJySkS2ux53YP1xp/s2Kt8yxmQAtwAv+ToWXzPGxABLgV8DiEifiLT6NiqfCgLCjTFBQARQ5+N4Rpw/JIB04MSQ57X4eaMHYIzJBmYDW3wbic/9G/APgNPXgYwBuUAD8B+uIbGXjDGRvg7KF0TkJPAscBw4BbSJyF98G9XI84cEYC6wza/XvhpjooA3gK+ISLuv4/EVY8ytQL2IbPN1LGNEEDAH+HcRmQ10AX45Z2aMiccaKcgB0oBIY8y9vo1q5PlDAqgFMoc8z2ACnsoNlzEmGKvxf01E3vR1PD62CLjNGHMUa2hwhTHmv3wbkk/VArUiMnBW+N9YCcEfXQfUiEiDiNiBN4FyH8c04vwhAXwE5BljcowxIVgTOX/wcUw+YYwxWOO7lSLyM1/H42si8k0RyRCRbKz/F6tFZML18oZLRE4DJ4wx+a5NHwP2+zAkXzoOLDTGRLj+bj7GBJwQD/J1AN4mIv3GmMeA97Bm8l8WkX0+DstXFgH3AXuMMTtd274lIn/2YUxqbPky8Jqrs3QEeMDH8fiEiGwxxvw3sB1r9dwOJmBJCC0FoZRSfsofhoCUUkpdgCYApZTyU5oAlFLKT2kCUEopP6UJQCml/NSEXwaq1HAZYxzAniGb7hCRoz4KRymv02WgSrkYYzpFJOoS+4NEpH80Y1LKm3QISKlLMMbcb4z5vTHmj8BfjDFRxpi/GWO2G2P2GGNudx2X7aqh/5KrfvxrxpjrjDEbjDEHjTHzXcdFGmNeNsZ85Cq45peVadXYoGcASrmcMwRUIyIrjTH3A/8EFItI80BpYBFpN8YkAZuBPGAKcAirwuo+rBIku4CHgNuAB0TkDmPMj4D9IvJfxpg4YCswW0S6Ru83VcqicwBKDeoRkdILbH9fRJpdjw3wI2PMUqwS0ulAimtfjYjsATDG7AP+JiJijNkDZLuOuQGrAN1TrudhQBYTsM6MGvs0ASh1eUN75/cAycBcEbG7KokO3Cqwd8hxziHPnQz+rRngLhGp9l64Sg2PzgEo5ZlYrHsI2I0xy7GGfjzxHvBlV4VJjDGzRzpApYZLE4BSnnkN2xmrhAAAAFFJREFUKDPGVGCdDVR5+PPPAMHAbmPMXtdzpXxCJ4GVUspP6RmAUkr5KU0ASinlpzQBKKWUn9IEoJRSfkoTgFJK+SlNAEop5ac0ASillJ/6/wGfJpKhMEgq2gAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"for col, label in zip(rog_10.results.T, labels):\n",
" plt.plot(col, label=label)\n",
"plt.legend()\n",
"plt.ylabel('Radius of gyration (Å)')\n",
"plt.xlabel('Frame');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Transforming a function into a class\n",
"\n",
"While the `AnalysisFromFunction` is convenient for quick analyses, you may want to turn your function into a class that can be applied to many different trajectories, much like other MDAnalysis analyses.\n",
"\n",
"You can apply `analysis_class` to any function that you can run with `AnalysisFromFunction` to get a class."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"RadiusOfGyration = analysis_class(radgyr)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To run the analysis, pass exactly the same arguments as you would for `AnalysisFromFunction`."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"rog_u1 = RadiusOfGyration(u.trajectory, protein, \n",
" protein.masses,\n",
" total_mass=np.sum(protein.masses))\n",
"rog_u1.run();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As with `AnalysisFromFunction`, the results are in `results`."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEGCAYAAABlxeIAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3zV1f348de5N/fmZic3OyQhYYYtkLBXUAQVELQDZ4ta6l7VTltFv/ZXW7RV0VIHxT1QqYMComxBIAxZBgIZZJCd3OTm7nvP748bqcjKhSQ34zwfjzzI/az7viG57/s5432ElBJFURRF0fg7AEVRFKVjUAlBURRFAVRCUBRFUZqphKAoiqIAKiEoiqIozQL8HYAvYmJiZFpamr/DUBRF6VR27dpVLaWMPd9xnSohpKWlkZOT4+8wFEVROhUhRFFLjlNNRoqiKAqgEoKiKIrSTCUERVEUBehkfQhn4nQ6KSkpwWaz+TsUvzAYDCQnJ6PT6fwdiqIonVynTwglJSWEhYWRlpaGEMLf4bQrKSU1NTWUlJSQnp7u73AURenkOn2Tkc1mIzo6utslAwAhBNHR0d327khRlNbV6RMC0C2TwXe682tXFKV1dfomI0VRlI5KSsnhika2Hq0hLjyQsb2iiQ4NbPG5JXVW9pWY2Fdaz13ZfQg3tG1foUoIbeS7SXQxMTGEhoZiNpv9HZKiKOdhcbgw210IBEJAVLAereb0u3CX28O63Ere3H6cr4/VEB4UgDFET1SwntDAAIIDA9AI2FFQywnTqU26AxLDSTUGodUINEIQog8gJkxPTGggUkJ+tZmjlWaOVJipbXIAoNdquGJwIpekRLbp61cJQVGUbi+3vIF/bynkP3tLsbs8J7cnhBu4fnQq87JSMIbo2VVUx7rcSj75powTJtvJ/XaXh9omO3VNTsobbFgcbuxON8OSI7n/slgm9I2lssHG1mM1bD1WTWG1BbeUuD2SJruLmiYHbo93sbJwQwB94kKZNiCeIckRDEuOpF9CKIEB2jb/OaiE0ArmzJlDcXExNpuN++67jwULFvg7JEVRfsDp9rDpSBUr9pSy6UgVBp325B3AoRMNGHQarhnRg4FJESAlLo9k/eEqnll7hOe+zCNYr6XB5kKnFYzrHcOjswZx2YA4ArQt64rtERnE8NQo7sruc9o+j0disjpxS0l0iN5vfYNdKiEs/PQgh8oaWvWaA5PCeXTWoHMes3TpUoxGI1arlaysLK699tpWjUFRlFM5XB4qG21UNtqpbLBTWm+luNZCca2F6iYHVocLq9ON0yXRagQBWkG9xYnJ6iQqWMf0QQlohKDe6sBsd/HbKzKYl5VCZLD+lOeZPz6dguom3tlxnHqLg+z+cUzoG0NYK7flazSCqBD9+Q9sY10qIfjLc889x4oVKwAoLi4mLy/PzxEpStfSYHOy4XAVu4vq2FNcz7dlDTjcnlOOCdFrSTEGExduIDHcQLBeS4BW4PaAy+PBEKDl8kHxTOwbiz6g5QMs02NC+P2VA1r7JXVIXSohnO+TfFvYsGEDX3zxBdu2bSM4OJgpU6aoeQGKcpGcbg8F1U3sKzGx+sAJNh2pxuH2EKzXMqRHBPPHp5EeE0JceCBxYQYSIwwY/djU0lV0qYTgDyaTiaioKIKDg8nNzeXrr7/2d0iK0uk02JxsO1bDlrxqdhTUkl9txun2drImRhi4aWxPrhySwLDkyBa32Su+UwnhIs2YMYMlS5YwdOhQ+vfvz5gxY/wdkqJ0CnaXmy8OVfJ+TjFbjlbj9kiC9Vqy0oxMHRBH//gw+ieE0T8+DM0Zhn4qrU8lhIsUGBjIqlWrTtteWFh48ns1B0HpDg6VNfB+TjGf7TtBgEaQYgwiJSqY3nGhZCSEkZEYjtXhJqewlp2FdazLraDO4iQxwsCCSb2Y0i+W4alRPrXvK61LJQRFUS6YlJLPD1WweN1R9pea0Gs1TBsYT5BeS3Gtha/za/hoT+lp5xlD9EzoG8uPRyYzvk/MGSd/Ke1PJQRFUS7I1qPV/HXNYfYW15MeE8JjswYyZ3iP04ZuNtqcHKlo5NsTjei1GjLTokiPCVEdwB2QSgiKovjE6nDz2CcHeS+nmKQIA09dO4RrRySftbM3zKBjZE8jI3sa2zlSxVcqISiK0mJHK83c9dZujlQ2cld2b+6Z2heDru1LKijtQyUERVHOye5ysz2/lnW53hFBBp2W1+aPYlK/WH+HprSyNk8IQoilwEygUko5uHnbMGAJEAoUAjdIKVu35oSiKBeloLqJJRuO8em+MiwON4EBGqZmxPHorEEkRBj8HZ7SBtrjDmEZsBh4/XvbXgEeklJuFELcAjwM/LEdYukUlixZQnBwMDfffLO/Q1G6obyKRp5bd5SV+8rQaTXMHd6DywfFM7ZXDEF61TzUlbV5QpBSbhJCpP1gc39gU/P3a4E1qIRw0u233+7vEJRuqMHm5B9r83htWyGGAA0LJvXm1gnpxIa1bEEXpfPz1wyQA8Ds5u9/DKSc7UAhxAIhRI4QIqeqqqpdgvPFzp07GTp0KDabjaamJgYNGsSBAwdOOebTTz9l9OjRDB8+nMsuu4yKigoA7r33Xh5//HEA1qxZw6RJk/B4PDz22GMsWrQI8BbOGzhwIEOHDmXevHnt++KULstkcbL2UAXv7TzOv78q4O9rjzB10Qb+vbWAn2alsPk3U/ntFRkqGXQz/upUvgV4TgjxJ+ATwHG2A6WULwEvAWRmZspzXnXVb6F8fyuGCSQMgSv+ctbdWVlZzJ49m0ceeQSr1cqNN97I4MGDTzlmwoQJfP311wgheOWVV/jrX//K008/zV/+8heysrKYOHEi9957L//973/RaE7N0X/5y18oKCggMDCQ+vr61n1tSrdR1Whn9/E6dhfVsfVYDQfKTMgf/DWN7BnFsvmjGNwjwj9BKn7nl4QgpcwFLgcQQvQDrvJHHK3lT3/6E1lZWRgMBp577rnT9peUlPDTn/6UEydO4HA4SE9PByA4OJiXX36ZSZMm8fe//53evXufdu7QoUO54YYbmDNnDnPmzGnz16J0HbnlDazYXcqqA+Ucr7UAoNMKLkmJ5L5L+zKudwzJUUEE6bQE6bVq+Kjin4QghIiTUlYKITTAI3hHHF28c3ySb0u1tbWYzWacTic2m40///nPrFy5EoC9e/dyzz338OCDDzJ79mw2bNjAY489dvLc/fv3Ex0dTVlZ2RmvvXLlSjZt2sQnn3zCE088wcGDBwkIUKOFFW+J6MPljeRXN5FfZaaiwYbV4cbm9FBY00RueSMBGsHEvjHcNKYnI3pGMigpQr3xK2fVHsNO3wGmADFCiBLgUSBUCHFX8yEfAf9u6zja0oIFC3jiiScoKCjgN7/5DYsXL+bJJ588ud9kMtGjRw8AXnvttZPbi4qKePrpp9mzZw9XXnklc+bMYfTo0Sf3ezweiouLyc7OZsKECbz99tuYzWYiI9t2oW2l4ztWZeb2N3aRV+ktnCgExIQGEqTTYtBpiArW89isgcwalkR0qOoHUFqmPUYZXXeWXc+29XO3h9dff52AgACuv/563G4348aNY926dUydOvXkMY899hg//vGP6dGjB2PGjKGgoAApJbfeeiuLFi0iKSmJV199lZ///Ofs3Lnz5Hlut5sbb7wRk8mElJIHHnhAJQOFNQfL+dX736AP0LDox8MY0iOCntHB6pO/ctGE/GHPUgeWmZkpc3JyTtn27bffMmBA91je7mzUz6B7sDrc/OOLI/xrUz7DkiN48caR9IgM8ndYSicghNglpcw833GqMVpROoG1hyp47JODlNZbuW5UKo/OGqjuCJRWpxKConRQFoeLL76tZHlOMZvzqukXH8p7C8Ywule0v0NTuiiVEBSlAzhYZuKJzw7hckuCAwMQwI6CWqxON/Hhgfz+ygzmj09Hp9YTVtqQSgiK4md7jtfxs6U7CNRp6RMbisnqxO50c+3IHswamkRWmlGtKay0C5UQFMWPvs6v4dZlO4kJC+St20aTHBXs75CUbkwlBEVpR1JKimos7CysJaewjo+/KSU5Kpi3bhtNfLgqKa34l0oIHZAqf9017S8x8cTKQ+woqAUgIkjHpRnxLLx6EDFq8pjSAaiE0AGp8tddS2m9lWc+P8JHe0owBut55KoBTO4XS+/YUNU3oHQoasjCRfrjH//Is8/+b9L1H/7wh9MK3Kny192PlJKv82u4481dTPrrej79powFk3qx/uEp3DaxF33jw1QyUDqcLnWH8NSOp8itzW3Va2YYM/jNqN+cdf+tt97KNddcw3333YfH4+Hdd99lx44dpxyjyl93fVWNdnYU1HKgzMTh8ka+PdHACZONyGAdt01M5+axaWpWsdLhdamE4A9paWlER0ezZ88eKioqGD58ONHRp04cUuWvu6baJgcvrj/KxiNVJ4vM6bSC3rGhjEo3Mr53DLOGJallJ5VOo0slhHN9km9Lt912G8uWLaO8vJxbbrmFP/zhD6r8dRfm8UjeyynmqdW5mG0uxveJ4ZoRyYztHc3AxHD0AaolVumc1DtLK5g7dy5/+tOfcDqdvP3221x55ZWq/HUXVFpvZX1uJR/sKmFvcT2j0o08OWcwfePD/B2aorQKlRBagV6vJzs7m8jISLTa05sHVPnrzsvicPH+zmLe2VHM4YpGAFKMQSz68TCuHdEDIVTHsNJ1qPLXrcDj8TBixAiWL19O37592/35O8LPoKsxWZy8+lUBr28rpN7iZHhqJFcOTiQ7wztcVCUCpTNR5a/byaFDh5g5cyZz5871SzJQWpfL7eGdncU88/lh6ixOLhsQz+2Te5GZZvR3aIrS5lRCuEgDBw4kPz/f32EoF8jp9lBUY6G41kJRTRPv7iwmt7yRMb2M/GnmIAYmhfs7REVpN+2xpvJSYCZQKaUc3LztEmAJYABcwJ1Syh1nv4qitK5jVWbe3XGcD3eXUtvkOLk9xRjEP28YwYzBCapZSOl22uMOYRmwGHj9e9v+CiyUUq4SQlzZ/HhKO8SidGN2l5vVB8p5e/txthfUEqARTBsYz7SB8fSMDiHVGExMqF4lAqXbavOEIKXcJIRI++Fm4Lt78QjgzIPwFaUVNNicvLj+GO/tPE6dxUmqMZhfz+jPj0emEBumisopynf81YdwP7BGCLEIbz2lcWc7UAixAFgAkJqa2j7RKV2ClJL/7i9n4acHqTLbmTEogetHpzK+d4yqI6QoZ+CvhHAH8ICU8kMhxE+AV4HLznSglPIl4CXwDjttvxD9Iycnh9dff/20AnlKy9Q1OTh0ooFvTzSw4XAVW45WMygpnJdvzmRYiprDoSjn4q+E8DPgvubvlwOv+CmODiczM5PMzPMOF1a+R0rJlqPVLN1SwPrDVSe3x4cH8seZA/nZ2J4EqLWIFeW8WpwQhBBxwHggCbACB4AcKaXnAp63DJgMbACmAnkXcI0OYcmSJSxZsgTwlqhIS0tj/fr1J/cXFhZy00030dTUBMDixYsZN24cK1as4IUXXmDt2rWUl5czefJkNm3aRG5uLosWLeKzzz5j48aN3HefN28KIdi0aRNhYapMwndsTjcf7y1l6ZZCDlc0EhMayL1T+zAqPZoBiWFEq0VnFMUn500IQohs4LeAEdgDVOIdLjoH6C2E+AB4WkrZcJbz38E7gihGCFECPAr8AnhWCBEA2GjuI7hY5X/+M/ZvW7f8deCADBJ+//uz7r/99tu5/fbbcTqdTJ06lQcffPCU/XFxcaxduxaDwUBeXh7XXXcdOTk5zJ07lw8//JAXXniB1atXs3DhQhISEsjN/V/8ixYt4oUXXmD8+PGYzWYMBrXEosPl4ViVmdUHynnz6yJqmhxkJITxtx8NZfYlSQQGqMqiinKhWnKHcCXwCynl8R/uaH5DnwlMAz4808lSyuvOct2RLQ2yM7jvvvuYOnUqs2bNOmW70+nk7rvvZu/evWi1Wo4cOXJy3/PPP8/gwYMZM2YM1113+o9p/PjxPPjgg9xwww1cc801JCcnt/nr6IjqLQ7+ufEY63Mrya9qwuXxdiVdmhHHrRPSGds7Wg0VVZRWcN6EIKV8+By7I6WU/2nFeC7KuT7Jt6Vly5ZRVFTE4sWLWbFiBQsXLgTglVde4bPPPiM+Pp5vvvkGj8dzyqf80tJSNBoNFRUVeDye0xbH+e1vf8tVV13Ff//7X8aMGcMXX3xBRkZGu742f7K73Ly+tYjF64/SYHMyuV8s0wbG0y8+jBGpUaQYg/0doqJ0KT53KgshgvA2F90IDMfbp9Bt7dq1i0WLFrF582Y0Gg1z585l7ty5J/e/9dZbJCcno9FoeO2113C73QC4XC7mz5/P22+/zeuvv84zzzzDQw89dMq1jx07xpAhQxgyZAjbtm0jNze3WySEkjoL7+8s5r2cYioa7EzqF8vvrshgQKIqI6F0T+76ejTh4QhN2w6OOGtCEEJogKF4O489wOV4k8BkvJPKrgY2t2l0ncDixYupra0lOzsb8I4SeuWV/w2auvPOO7n22mtZvnw52dnZhISEAPDnP/+ZiRMnMnHiRC655BKysrK46qqrTrn2P/7xD9avX49Wq2XgwIFcccUV7ffC/GBXUS0vrD/G+sOVAEzuF8uiH6czsW+snyNTlDOTHg/OshMIXQCawEBEUBBC753tLqXEnpuLecsWLDt3EpieTsTVVxM4YMA5mzillNiPHKHpq63YDuzHuv8AzuJieq38jMAzrKrYms5a/loI8T5wHG8fggHYDrwLrAZypZTpbRrZGXTU8tf+1pl/BlJKthfU8vy6PL46WkN0iJ4bRqfyk6wUkqNUk5DSfqTTSdPWrXgsFjTh4WjDIwjslY6m+UPc97kbGzGtWEHdW2/jKCo6dacQ3sQgBJ7m0YX6Xr1wFhcjnU4C+/YlsF8/pNOJdDoRej0BMTEExMbgrquj8ct1OEtKANAlJWEYPBjDkMFEzL4aXXzcBb221ih/3QN4E7gH2AesBzZLKe1CiC4/QUxpW0U1TazYU8p/9pRSWGMhJjSQR64awPWjUwnWqyK8iu9cdXXYDhzEY7EQ2Csdfc+eCL3+1GOqqmj84gssObvQpSRjGDgQfWoqjV98Sf377+OqrDzleE1oKJHXXkPUDTcQEBtL01df0bh2LQ1rv0BaLARdcgnxN9+ECNAh7TY8Vpv3X4sV6XRiGDyYkPHj0MXF4a6vp2H1akyffobtwAGEXgc6HdJmp+nrr/GYTAidjpBx44he8AtCp0xBF3dhCeBCnesv7wbgWmAwoANuArYLIQ4DIUKIYCmlpR1iVLqQA6Umnvsyj88PVSAEjO0VzZ3ZfZg1VC1Gr5xOejze5pNtX2PJycFdW4unqQmPxQJaDZqgYDRBQbgqK3GWlp56slaLLjERbXg4mvBwpNWKdd8+kJKA2Fhcq1dDc58eQhAyYQIJjz2KLjkZT0MDrro6Glevofatt6l9/Q1EYCDSZkMbEUH4lVcQNe86ggYPavFr0UZGEjVvHlHz5p1xv8duBynR+HF4uc8rpgkhJuHtS5iJd2La7LYI7EzO1mSUkZHRbYcdSinJzc3tsE1GdU0O8irNHKloZH1uJV/mVhJuCODn49O5blQKiRFB/g5R6WBc1dU0ffUV5i1f0fTVV7hrawHQ9+xJQFIimpAQNEHB4PHgsVrxWC1oIyIJGjIYw6DBaMJCceTnYz92DGdJKZ7GRtyNjSAlIRPGE3755ej79EE6HNiPHMF+9BjBI0egP0utNGdFJfXvv4+7oYGwqdkEZ2YidLr2/JFctJY2GV3wEppCiEDgSinligu6wAU4U0IoKCggLCyM6OjuNxZdSklNTQ2NjY2kp7d7l8457S2u56+rc9l6rObkNmOInvnj0vjZ+DTCDZ3rD0ppO25zE+aNG7Ds3IklJwfH0WMAaI1GQsaPJ2TsWELGjkGXmOjnSDuvVltCUwhxI/D2D0tUSCntwAohRG8gUUq55YKjvQjJycmUlJRQVVV1/oO7IIPB0GEmrEkp2V9q4sX1x1h9sJzoED2/mtaPoSmR9I0LJTHC0O2SdnckPR5vx+p5/q9d1dXUvvEmde+8g6ehAU1ICEEjRhAx+2pCxo/DMGBAmw+zVE7Vkt67aGCPEGIXsAuowjvqqA/eIajVeEtb+IVOp+twn467m9zyBlbsKWXV/nKO11oIDQzgwWn9uGVCOqGBqoO4O3A3NGDeuJHGz9di3rIF3G60RiNaYxQBkVFoIyPRRkYipQd3bR3u2lqs+/YhHQ7Cpk3D+LObCRo2DBGgfl/8qUVNRkIILd4idOOBRLzF7b4FVp2ppEVbOVOTkeI/e4vrWbwujy++rSRAIxjfJ4arhiQyfVACEcGqSag78FitVL/0ErWvLkU6HATExRE6NRtNcAju2lpcdbW46+ubv0wIQBsdjdYYRWDfvhhvupnAXuoDXVtrtSYjACmlG1jb/KV0Y26PZF1uJa9tLWTL0Woig3U8OK0fN43pSVSI/vwXULoEKSWNaz6n4q9P4So7QfhVV2G86UYMQ4eqZp5OTN2fKedlcbg4UNrA9vwa3t1ZTGm9lYRwA7+ZkcFNY3uqZqFuRjocnFi4ENOHHxGYkUGPv/6VYLWGR5eg/pKVM6posPHh7hI+/eYEh8sbaC4wythe0Txy1QAuGxiPTi060+24amspufderDm7iL7jdmLvuku1+3ch6n9SOam41sKWo9WsOVjOpiNVeCRkpUVx99S+XJISwdDkSGLUojPdkquujqbNm6l69jlc1dX0eOZpwq+80t9hKa3MlxXTAvHOXE77/nlSysdbPyylLTlcHjYeqeKb4nqqzXaqzQ6OVZkpqPbWXekRGcSdU/rwo5HJpMWcXsdF6R48Fgv1K1bQ8OlnWL/5BqREl5REzzffIGjIEH+Hp7QBX+4QPgZMeIee2tsmHKWtWB1u9hTXsWp/OZ/tK6PO4kSrERhD9ESH6OkdG8JNY3oyqV8MvWND1XyBbsxVVUXtm29R9+67eEwmAgcMIObOOwmdMhnDoEGq07gL8yUhJEspZ7RZJEqr8XgkRbUW9pea2F9ST05RHQdKTTjdksAADZcPSuCa4T2Y0DdG9QMoJ7nq6qh55RXq3nobabcTdtmlGOffQvCI4f4OTWknviSErUKIIVLK/W0WjXJRyuqtvPl1Ee/tLKamyQGAXqthSHIEt07oxaj0KLLSjISpshEK4KyowH74MI7jxTjyj2H6+BM8FgvhM2cSe9ed6NPS/B2i0s58SQgTgJ8LIQrwNhkJQEoph57rJCHEUryF8CqllIObt70H9G8+JBKol1Je4mvwitehsgYWr89jzcEKpJRcOiCeSzPiGJIcQb/4MHUXoJzCnpdH9T+X0LB6NXi8FWlEUBChEyYQc8/dGPr183OEir/4khAudLmuZcBi4PXvNkgpf/rd90KIp/H2TSg+yq8y88zaI3y27wRhhgBum5jOjaN7qrWGldNIKbHu3k3tstdoXLsWTXAw0bfMJ3TqVPQpKWhjYlS/kdLyhCClLBJCDAMmNm/aLKX8pgXnbRJCpJ1pn/D+Bv4Eb1kM5SxsTjdHKhrZX2oi90QjhTVNFFQ3UVpvJUin5e7sPvxiYi9VLkI5jcdqpWHlSmrfehv7t9+iCQ8n5s47iLrpJgKiovwdntLB+DLs9D7gF8BHzZveFEK8JKV8/iKefyJQIaXMO8fzLgAWAKSepV55V+PxeKuGbjxSdXJ4qKt5ZlhYYADpsSGMSI1iXlYK80alqrkByimk2409L4/695dj+vRTPI2NBPbtS8LjC4mYNQtNkFqDQjkzX5qMbgVGSymbAIQQTwHbgItJCNcB75zrACnlS8BL4C1udxHP1aF5PJI9xXV8+s0JVu4/QVWjHSFgaI8IbpvYi6HJEQzpEUFyVJC6tVeQHg+m/3xM09atuOvqvF/19bgbG/GYzSAlQq8nbPp0on7yY4IyM9XvjXJeviQEAbi/99jdvO2CCCECgGuAkRd6jc5OSsm+EhOf7Stj5b4TlJlsBAZomJoRx4zBCUzsG4tRFYxTfsCen0/5nx7FkpNDQGIiAXGxaGNjCOzbB014BNqwMALi4wm7fJpqFlJ84ktC+DfeNZW/WyFtDvDqRTz3ZUCulLLkIq7RoXk8kvxqM0cqzFgdbpxuD1anm6IaC8eqzORVmClvsKHTCib1jeXhGf2ZNjBBFYtTzsjT1ETNq0upefllRHAwiU/+HxHXXKM++SutxpdO5WeEEBvwDj8VwHwp5Z7znSeEeAeYAsQIIUqAR6WUrwLzOE9zUWdUVm/ls31lrMut5EBpA2a767RjgvVaesWGMCrdyIS+MUwfqNYPUM5Out2YVqzw1hGqqiJ85kzif/dbAqKj/R2a0sW0ZAnNcCllgxDCCBQ2f323zyilrD3X+VLK686y/ec+RdpB1FscnDDZcLo9ON0Sk9VBWb2NsnorOUV17Cjw/jgGJoYzd3gPhiZHMCAxnDBDADqthsAADcYQvfpUp5yXdLloWLWa6n8twXH0GEHDh9PjuWcJHq5mDittoyV3CG/jnVi2C/h+p65oftyrDeLyGyklJXVWjlaZSTUG0ysmBCEE5SYbSzYe450dx7G7PKedF6AR9IoN4cFp/Zg9LEkVhVN84rHZMP3nY6TDjtDp8Fht1L37Ls7jx9H36U2Pf/yDsOmXqw8SSps6b0KQUs5s/rdLrXNXUN3Es18cYe2hCsIMOowheoL0WvIqGmmw/a+ZJypYR0ZCOLuK6nBLyTXDe5CdEYdeqyFAKwgz6OgRGURsWCBajfpjVXznaWqi+M67sGzffsp2w6BBxD3/HGGXXqoKyintwpd5CF9KKS8937aO6GhlIxUNdqwON1anmy151XywuwSdVjB7WBJSQm2TA7PdxcxhSQxKCqdPbCiFNU3sKqpjf2kD147swR2T+5AarWYBK63H3dBA8YJfYt2/n6Sn/kLo5MlIpxPp8RAQF6fuCJR21ZI+BAMQjLdTOIr/DTUNB5LaMLZW87c1h1lzsOLkY71Ww01jenJndm/iwgxnPW90r2h+mtU9JsMp7UtKie3AQcoffRRbXh49/vF3wqdN87RicpcAACAASURBVHdYSjfXkjuEXwL3433z38X/EkID8EIbxdWq7r+sH7eMTydIr8Wg0xIXFkhksBrfr7Q/e0EB9R98QOPqNThLSxFBQaS8sJjQSZP8HZqitKgP4VngWSHEPRdZpsJvBiSG+zsEpZuz5+VRveRfNKxaBRoNIePGEnPnnYRdOhVtZKS/w1MukNPt5ETTCRJDEtFpO//QcV/mITwvhBgMDAQM39v++tnPUpTuy2O10vjlOkyffEzTps2I5gqjxp//nICYGH+Hp/xAcWMxawrXYHVZCdOFEaYPo7+xP4OiB53Sl1NjrWFTySY2l25ma9lWmpxNaIWW5LBk0iPS6R/VnwHGAaRFpHG49jDby7ezu2I3AEaDkeigaPpG9mVs0lgGxwwmQHP2t+GKpgq2ndjGtrJtPDLmEcL0YW36M/ClU/lRvBPMBgL/xVsOewvfK2utKN2ddDpp2raNhv+uonHtWjxNTQQkJqoKox2MxWnhRNMJypvKKWwo5PPCz9ld6X3TFgjk90bYJ4cmc0X6FYTpw1hfvJ69lXuRSOKC4piRNoMhMUMoayqjwFTAsfpjbCrZhEf+b2h6mD6MkfEjCdQGUmOtIa8ujy+KvuDFb14kVBdKv6h+RAZGEhEYgV6rp9HRiNlppqSxhHxTPgDRhmiONx5nUPSgNv25CClbVi9OCLEfGAbskVIOE0LEA69IKWe1ZYDfl5mZKXNyctrr6RSlRaTTSdOOHTSuXkPj2rW46+vRhIURdvk0ImZfTXBWpho22sqklFRaKiloKKDQVIjZaUYrtGiEhkBtIKH6UEJ1oUQZokgLTyMiMAK3x82W0i0sP7KczaWbT3nTTo9IZ3bv2czsNZO44DgsTgsmh4kdJ3awunA1209sxy3dDDAOIDslmykpU8gwZpxxFJjNZSOvLo+ChgJ6R/Qmw5iBVqM95Zh6Wz07ynew7cQ2Ck2FmBwmTDYTTo/zZOyxwbFkxWcxNmks/aL6XdSIMyHELill5nmP8yEh7JBSjhJC7AKygUbggJSybVPW96iEoHQkrro6qp57jsZVq71JIDiY0Oxswq+6ipAJ49Ho1cCF1iSl5Juqb1iZv5LPiz6n1nbOIgmniAyMRCu01NhqiAmKYVbvWfSP6k9CSAKJIYkkhiSe8w23zlaHw+0gPiS+NV5Ku2tpQvClilqOECISeBnvaCMzsOMC41OUTq1p61bKfvs7XHV1hE+fTviM6YRMmIDGcPZhzIrvnG4nOyt2srF4I+uL13Oi6QSB2kAmJ08mKyGL9Ij0k3cAHunBLd3Y3XbMDjNmp5lqazVFDUUUNRRhdpi5PO1yJqdMRqfxrQM4ytA9mvpalBCaVzb7f1LKemCJEGI1EC6l3Nem0SlKByKlxFFQSN0771D3xhvoe/Ui7Z8vEjSo3W6SuzSzw8zeqr3sr9pPQUMBRQ1FFJgKsLqsGLQGxiSO4Z7h9zA1dSohurOXhgkjjJgg1Wl/IVqUEKSUUgjxH5rXLpBSFrZlUIrSUbhNJiw7d9K0dRvmzZtxFhcDEHX99cQ9/JBafewimOwmdlXsIqcih10Vu8itzcUjPQgESaFJpIWnMbzvcMYkjmF04miCAtTPuq350mT0tRAiS0q5s82iUZQOQDoc1H/0EfXLP8B26JB39bGgIELGjCH6lvmETJyEPrmHv8PscCxOC5tLN7OheAMWpwW9Vo9eqychJIH+Uf3pb+yPyW5ic+lmNpds5lDNISQSvUbP0Nih/GLILxgZP5JhscMI1qkSMf7gS0LIBn4phCgCmmiudiqlHNomkSlKO5JS4mmy0Pj551S/8ALO0lIMgwYRc9ddhIwZTdDQoQjVSXwaKSXby7fz4ZEP2ViyEavLenKsvdPtxO62U2mpxC3/t9iiRmgYGjOUOy65g6z4LIbEDiFQq9YF7wh8SQhXtFkUitLOpJRYduyk9rXXsB04gLu+HulwAN4qowmPPUrIhAmquNxZ2N12Psr7iHdz3yXflE9kYCSze89metp0RsSNOGWYpd1t52j9UQ7XHsagNTAuaRyRBjU7uyPyJSHcAyyVUh5qq2AUpa15HA4a16yh9t/LsB06hNZoJHTyZAKijWijogjs25eQiRNVIjiHfVX7+ONXfyTflM/g6ME8OeFJpqdNP+un/EBtIIOiB7X5pCrl4vmSEHKBl4UQAXjXV35HSmlqm7AUpXU5S0upe+996j/4AHdtLfr0dBIeX0jE7NlqqGgLON1OSswlrDi6gtcOvkZsUCwvXvoiE5Mn+js0pRX5UsvoFeAVIUR/YD6wTwjxFfCylHJ9WwWoKBdDulzUvLqUqsWLwe0mdGo2UfOuI2TcWDV7+DxqbbW8vO/lk+P/v5vZe23fa/lV5q/avK6O0v58uUNACKEFMpq/qoFvgAeFEL+UUs47yzlL8S7BWSmlHPy97fcAdwMuYKWU8tcX9hIU5czsR49S9rvfY9u/n7ArZhD/8MPokjrFEh5+ZXFaePPbN1l6YClWl5XslGxm9Z5FSljKydFCStfkS3G7Z4DZwJfAn6WU381SfkoIcfgcpy4DFvO9InhCiGzgamColNIuhIjzNXBFORNHcTHm9etpXLcey86daMPC6PH3Zwi/Qo2JOB8pJasKVvH0rqeptFSSnZLN/SPup1dkl1o2XTkHX+4QDgCPSCktZ9g36mwnSSk3CSHSfrD5DuAvUkp78zGVPsShKKdwlJTQsGoVDatWYT/0LQCBffsQfeutGG++SZWa/gEpJSWNJRyoOYDD7cAQYEAjNLx56E12V+5mgHEAf5v0N0bEj/B3qEo78yUh7AUyfjD6wgQUXUDncj9gohDiScAGPHS2CW9CiAXAAoDUVLWcpeIlpcSyfTs1r7xK05YtABiGDiXu178m7LJL0avfFdweN5/mf8qKvBUA6LQ6BIIjdUfOWBguKjCKR8c+ytw+c0+rzql0D74khBeBEcA+vJPSBjd/Hy2EuF1K+bmPzxsFjAGygPeFEL3kGUqvSilfAl4Cb7VTH55D6YJcdXWYv/ySunffw3bgANqYGGLvu5fwWbPQJyf7O7x2ZXaYqbBUEBscS5guDCEEHumh0dHIjvIdLN6zmHxTPr0jehMdFI3dZcct3UzoMYFhscMYGjuUEF0Idpcdu9tOz/CehOpD/f2yFD/yJSEUArdKKQ8CCCEGAg8DTwAfAb4khBLgo+YEsEMI4QFigCofrqF0I43r11P3xhs0bd8Bbjf6tDQSFi4kYs7VaAK71yxXl8fF8iPLeWHvC5js3pvzoIAgggOCqbfXn5wVnB6RzjNTnuGy1MvUvAqlRXxJCBnfJQMAKeUhIcRwKWX+Bfyy/QeYCmwQQvQD9HhHLSnKKaTHQ9Xzz1PzzyXoUlKIvu02wqdfTuCAAd3yTe7rE1/z1I6nOFp/lNGJo7m699XU2mopbyo/WTbCaDDSI7QHE5MnnnN5RkX5IV9+Ww4LIf4JvNv8+KfAESFEIOA820lCiHfwLr0ZI4QoAR4FlgJLhRAHAAfwszM1Fyndm9tspuzhX2Nev56Ia68h4dFHu+2iM+VN5SzKWcSawjUkhybzbPazZKdkd8ukqLQdXxLCz4E7gfvx9iFsAR7Cmwyyz3aSlPK6s+y60YfnVroR6XBQ/5//ULPkXzgrKoh/5BGibri+y7/5SSmpslZR3FhMqbkUs8OMw+2gxlbDe4ffwyM93HXJXcwfPF8Vg1PahC8zla3A081fP2RutYiUbktKSf37y6lesgTXiRMYhgwh6am/EJyV5e/Q2txn+Z/x5+1/ptHReMb9U5Kn8OtRvyYlLKWdI1O6E9XAqHQIbnMTJ37/exo//5yg4cNJfPxxQiaM7/J3BS6Pi2d2PcMbh95gRNwIZqTPIDUsleSwZML0YRi0BvRaveoLUNqF+i1T/M6eX0DJPffgKCgg7uGHMd4yv8smgpLGEnaU78DpduKSLr48/iU7y3dyw4Ab+FXmr3xe61dRWtN5E4IQ4g0p5U1CiPuklM+2R1BK92DPL6DurbeoX7ECTWAgqUtfJWTMGH+H1SYKTYW8vP9lVuavPGWxmKCAIJ6c8CSze8/2Y3SK4tWSO4SRQoiewC1CiNfxdiifJKU8fcqjopyD7fARKv/2N+8MY52O8CtmEHf//V2q8JyUkrz6PLaVbWNr2Va+PvE1eo2e6zKu4yf9f0KYPgydRkdQQBB6bfccOaV0PC1JCEuA1UAvYBenJgTZvF1Rzks6HFS//DLVS/6FNjSUmHvvIeonP+mUtYbcHjcbijfw0dGPMDu8YyqEEFicFmpttdTZ6nB4vCuw9YroxS2Db+GGATcQE9T5XqvSfZw3IUgpnwOeE0L8U0p5RzvEpHQxrqoqzJu3UPvaa9gPHyb8qquI/8PvCTAa/R2aT5weJ3l1eew4sYN3D79LqbmUxJBEUsNSkUg80kNccBz9ovphNBhJj0hnbNJYEkIS/B26orSIL8NO7xBCDAO+WyJpk5RyX9uEpXR20u2m/qOPqHvnnZMVSHVJSSS/sJiwSy/1c3QtZ3Fa+PjYx6zMX8m3Nd+e/NQ/Im4ED2U+xJSUKWoEkNJl+LIewr14q45+1LzpLSHES1LK59skMqXTsuzZQ8X/PYnt4EEMgwYR+8ADhE6aSGBGRqcZPVRrq+WV/a+wIm8FZqeZAcYBXJdxHYNiBjE4ZrCaD6B0Sb58tLkNGC2lbAIQQjwFbANUQlCQDgeNGzdi+s/HmL/8koC4OJIWLSL8qis7TRL4TqGpkNu/uJ2Kpgqm9ZzGDQNvYFjsMH+HpShtzpeEIAD39x67+cGII6X78djtVL/wIvXvvYfbZEIbE0P0HbcTc9ttaEJC/B2ez/ZW7uWedfegERpev+J1hsQO8XdIitJufEkI/wa2CyFWND+eA7za+iEpnYUtN5eyh3+NPS+PsOnTibz2GkLGjUMEdJ42dbPDTG5tLscbj5Nfn8+7h98lPjieJZctISVcNQsp3YsvncrPCCE2ABPw3hnMl1LuaavAlI5LOp3ULFtG9XPPo4mMIOXllwidOPH8J3YQTo+Tr0q/4tNjn7KheMPJjuIATQCjEkbx/yb+P4yGzjUCSlFag08f5aSUu4HdbRSL0gk0bd9Bxf89gT3vKGHTLiPh8ccJiIryd1jn5XQ7+frE16wtWsu64nWY7CaiAqP4Ub8fMTF5Ij3De5IYkqhGDCndmvrtV85LSol1zx7q3nyThv+uQtejB8kvvkBodueox7+heAOPbX2MGlsNobpQpqRMYUbaDMb1GKdqBynK96iEoJyV22ym7u13MH30EY7CQkRwMDF33kH0ggVoDAZ/h3deFqeFRTmLWH5kOf2j+rNw3ELGJo1VpSIU5Sx8mYcQAlillJ7mZS8zgFVSyrOulqZ0TtLlon75cqqeX4y7tpbgzEyif/ELwmdM79Ajh5xuJxtLNpJvyqeooYhdFbsoM5cxf9B87h5+t0oEinIevtwhbAImCiGigC+BHLzLaN7QFoEp7c9dX0/D6jXUvvEGjmPHCM7MJO5f/yJoyGB/h3Zeeyv3snDbQo7WHwUgLjiO9Ih0Hh/3OKMSR/k5OkXpHHyahyCltAghbgWel1L+VQihRhl1Us6yMmyHDuGqrcVdW4ft4AHMGzYinU4C+/YlefHzhF56aYfvI7A4LTy7+1neyX2H+JB4ns1+ljGJYwjWBfs7NEXpdHxKCEKIsXjvCG5t6flCiKXATKBSSjm4edtjwC+AqubDfi+l/K8PsSgXwFFSSuPatTSsXoXtm1PLUAXExhJ1/XWEz56NYeDADp8IAPZX7ec3m39DSWMJ1w+4nnuG30OIruM2aSlKR+dLQrgf+B2wQkp5UAjRC1jfgvOWAYuB13+w/e9SykU+PL/iI4/djnX3bsybNmPetAnHsWMAGAYOJPbBBwkZO5aAmGi0RiOawM6zaLvL42LZwWW8sOcFYoJjWDp9KZkJmf4OS1E6PV8mpm0ENn7vcT5wbwvO2ySESLuQ4BTfuerqaPjkE8ybNmPJyUHa7QidjuCsLCJ//CPCsrPR9+zp7zB9UmWpYumBpeyu3E2VpYoaWw0e6WF62nT+OOaPRARG+DtERekSfBlltB7vgjinkFJOvcDnvlsIcTPezulfSSnrzvK8C/BWWSU1NfUCn6rrs+7dS+3bb9O4eg3S4UDfuzeRP/kJIePGEjJqVIceHXQ21dZq/n3g37x3+D1cHhejEkbRP6o/scGxDIoeRHZK55gHoSidhZDytPf4Mx8oxMjvPTQA1wIuKeWvW3BuGvDZ9/oQ4oFqvAnmCSBRSnnL+a6TmZkpc3JyWhRvd2E7coTKRYto2rQZTUgIEVdfTdR18wjs29ffoV2wgzUHefvbt1lVsAq3dDOz10xuH3q7qi2kKBdICLFLSnnedlVfmox2/WDTV0KIjWc8+PzXqvjueyHEy8BnF3Kd7shtMuE4XoyzpBjzli2YVvwHTUgIcQ8/RNS8eZ3yTgCgzlbHmsI1fHrsU/ZV7yMoIIhr+l7DjQNuJC0izd/hKUq34EuT0ferfWmAkcAFrQ0ohEiUUp5ofjgXOHAh1+kOpNuNde9ezBs2YN6wAXve0f/t1Okw3nQj0bff3inqCZ3JsfpjPLf7OTaVbMIlXfSJ7MNvsn7D1X2uJkwf5u/wFKVb8WWU0S68TTwCcAEF/G/46VkJId4BpgAxQogS4FFgihDikubrFQK/9CnqbsBtMlH/wQfUvvUWrrITEBBAcFYmsbNnE5ieji4lBX1KCprgzjne3uK0sGTfEt44+AZBuiBuGngTV/W6iv7G/v4OTVG6LV+ajNIv5AmklNedYbNaR+Es3I2NVC9+gbr330darQRnZRH/0EOETJqENjTU3+FdFI/0cLD6IOuK1/HJsU+otFQyp88cHhj5gCo3rSgdQEsmlk2VUq4TQlxzpv1Syo/OtF3xjZQS08cfU7noadw1NUTMno3x5z/DMGCAv0O7YHsr97Lu+DoqrZVUW6o5ZjpGtbUardCSmZDJosmLGB433N9hKorSrCV3CJOBdcCsM+yTgEoIF0G6XDR+8SU1S5di27cPw7ChpPzzn52iftDZFDcW8/ddf2dt0Vr0Gj2xwbHEBMWQlZDFxB4TmZQ8Sc0dUJQO6LwJQUr5aPO/89s+nK7P43DgLC7GUVSE/fBh6pd/gLOsDF1KColP/h8Rc+ciNBp/h3lBysxlvHbwNZYfWU6AJoA7h93Jzwb9TNUVUpROoiVNRg+ea7+U8pnWC6frkQ4H9sJCmjZvwbx5M5Zdu8D5v4rhQZkjif/977yLzWi1foz0wh2pO8KyA8tYVbAKgNl9ZnPXJXcRFxzn58gURfFFS5qMvhv71x/IAj5pfjwLb0ls5Xs8Fgt1779P46rVOMvKcFVXQ/Pkv8C+fTHedBOGjP7oe/ZE37Mn2shIP0d8YRxuB58Xfc77h99nT+UeggKCmJcxj5sH3kxiaKK/w1MU5QK0pMloIYAQ4nNghJSysfnxY8DyNo2uA5JS4iorw7p/P9Z9+/FYmtCnpKBLScFRWETtsmW4a2sxDB1KyKSJ6BKT0PXoQcjoUeiSkvwd/gWTUlJjq2Fb2TY2lWziq7KvaHQ0khqWykOZD3F176uJNHTO5KYoipcv8xBSAcf3HjuAtFaNpgOTbjcNK1dS/c8lOAoKABB6PSIoCI/JdPK4kAkTiLnzDoJHjPBXqK3CZDexumA164rXUWYuo8JSgdVlBSDaEM1lqZcxI30GYxLHoBGds89DUZRT+ZIQ3gB2CCFW4B1dNJfTS1p3elJKLDt30rh6DWi1aCMj0BgM1H/wIY6CAgIzMoj/4yMEDbsEQ7++CL0ed0MDjuPFCL0OQ79+/n4JF8zpdrKldAuf5n/KhuINOD1OekX0or+xP5OSJxEfHM/I+JEMiB6gkoCidEG+TEx7UgixGpjQvGm+lLLTr5gmXS6cZWU4jhdjO3jQu6B8UREiKAih1eIxmwFv+3+P554l7LLLThsFpA0PJ2jwIH+Ef9GqLFXsr97P1rKtrC5cjcluwmgw8tP+P2V279lkGDNURVFF6SZ8uUNASrlLCFGMt9opQohUKeXxNomsDUmPh6at26h79x3MGzedOupn5EgS77id8OnT0QQFIR0O3GYz2sjITjsc9DtSSooaithRvoOd5TvZU7mHCou3zqBBayA7NZuZvWYyNmksOo3Oz9EqitLefCluNxt4GkgCKvH2KeQCneKjsfR4sB08hHn9ekwrP8NZdBxtVBTG668jsF8/b22gtDR0cacOlRR6PQHGzltWwe62s+PEDjaWbGRTySZONHlrCsYFxTEyfiRDYocwJGYIGcYMDAEGP0erKIo/+XKH8AQwBvhCSjlcCJENnKlOUYdT9fxi6t9/H1dVFQhB8MiRxN59D2HTL0ej1/s7vFZndVnZXLKZtUVr2ViyEavLSlBAEGMSx3Dr4FsZnTianuE9VVOQoiin8CUhOKWUNUIIjRBCI6VcL4R4qs0ia0XS7SJoxAhCs6cQOnlypy0VfT4mu4k3Dr3Bm9++SZOzCaPByMxeM8lOyWZU4igCtZ1n3WRFUdqfLwmhXggRincy2ltCiEq8ZbA7vLj77/d3CG3K6XHy6v5Xef3g6zQ6G5nWcxrz+s9jRPwIAjQ+dRMpitKN+fJucTVgBR4AbgAigMfbIiil5SxOCw9ufJCvSr9iaspU7rzkTrWmgKIoF8SXYadNzd96gNeEEFpgHvBWWwSmnF+drY67vryLgzUHeWzsY1zb71p/h6QoSid23nGUQohwIcTvhBCLhRCXC6+7gXzgJ20fovJDHulhW9k2frb6ZxyuPcwzU55RyUBRlIvWkjuEN4A6YBtwG/AwoAeullLubcPYlGY2l40qSxWV1kr2Vu7lw7wPKW4sxmgw8q9p/yIzIdPfISqK0gW0JCH0klIOARBCvAJUA6nfFblTWp/L4+Kbqm+8cweKN3HMdOyU/SPjR3LXJXdxWc/L1MghRVFaTUsSwslpvFJKtxCiwJdkIIRYCswEKqWUg3+w7yHgb0CslLK6pdfsqo7UHeHjox+zMn8lNbYaAkQAIxNGMiN9BgkhCcQFx5EalkpyWLK/Q1UUpQtqSUIYJoRoaP5eAEHNjwUgpZTh5zl/GbCYHxTCE0KkANOATlf6wld1tjoqLZVUWiqptlZTa6ul3l5Pna2OOnsd9bZ6qq3VlDWVESACmJQ8iSt7Xcn4pPGE6kP9Hb6iKN1ES9ZDuKhlvKSUm4QQaWfY9Xfg18DHF3P9jqjB0cChmkNsKdnC5tLN5JvyTztGr9ETZYjCaDASZYgiOSyZm2Nv5sr0K4kydM2Jc4qidGx+mbXUXBepVEr5zfnKJwghFgALAFJTU9shOt+VN5Xz8dGP2XZiG4WmQmpsNQDoNDoy4zOZ02cOyWHJxAZ5F5s3GowEBQSp0hGKonQo7Z4QhBDBwB+Ay1tyvJTyJeAlgMzMTNmGobVIjbWGo/VHqbHWUGOr4auyr9hWtg2P9DA0ZiiTkieRFpFGn8g+ZMZnqgXmFUXpNPxxh9AbSAe+uztIBnYLIUZJKcv9EE+LeKSHd3Pf5R+7/3Fy5TCA+OB4bhtyG3P6zCElLMWPESqKolycdk8IUsr9wMka00KIQiCzo44ycnqcHG84zpPbn2Rn+U7GJ41n/uD5xAbFYjQYiQiMUE0/iqJ0CW2eEIQQ7wBTgBghRAnwqJTy1bZ+3vPxSA+7Knaxv3o/PcN70j+qP7HBsewq38WGkg1sK9tGlbWKJqe3YkeoLpSF4xYyt89clQAURemS2jwhSCnPuWaClDKtrWNwepzUWGuoslRRZa1iX9U+VhaspLzp1BYqgUAiMWgNjEocxYQeE4gIjCAqMIrJKZNJCElo61AVRVH8plvURn5g/QNsLNl48rFWaBmXNI4HRjzA2KSxFDcWc7juMGXmMi6JvYTRiaPV6mGKonQ73SIh/Kjfj5icMpnYoFhig2JJDksmIjDi5P4oQxRDY4f6MUJFURT/6xYJYUrKFH+HoCiK0uGdt/y1oiiK0j2ohKAoiqIAKiEoyv9v796Do76uA45/j57oCUIP0AMQD4F5GFtYNnYNBIc42MQDuE4z9jS1ncmMp9OmSZN4WnuYSSd1478yfcwkaeLxI07sOI2dTOu4rh2K7dImNkZgDNiIp0AIEJKRhN7SanX6x/kRCSIkEFrtavd8ZnZ297e/n/ZeaXXP73fv3XOdcwEPCM455wAPCM455wIJMcvIOecmXDgELSfg3GE4dxQyp0PRYihYBO1n4Ng7ULsdOpsgLctuqZmQnArJ6ZCUAheyIojArX8BuSURLbIHBOecu5QqdDVDagakXSZjcXcr1P4PnNkLOmDbwn3QfAw+OQwttTDQP/L75JTA9HnQ0Qh9nRDqskAS7rNjVQG1++X3e0Bwzrlx1d8HnY0wZRqkBysSdjTa2Xrtdmg8AJ8cgp5WkGSYsRRm3QI5xbatu9Ver68GDYMk2X5gZ/V55VC4CBbfAwULIb/CGv2uc9B0AJoOQkYezLsD8ucPXgXEAA8Izrn4pGpn66d2w6lqaNgHLceh7TQQLK2Slm2Boa3enqdPheLlsPReKKiA7hY4+T58+G/Q1w4pGZAxDXJLYdXXYcE6KLvZunlGk5UPhQsjVdtx4QHBORdfes5bA179rJ2Rg/XNz7we5q6BaXMgZ6bt13EWOj+xvv15n4LiGyFpmFWDB8LWlZMa3znOPCA45yan/l5orrVB2+ZjNoDbUgt171lffEklbPguzL4VChdD8jU0d0nJwweKOOMBwTkX28IhOL0H6n4H545Yw996AlrrBgdzwbp+8sph+RdgxUNQuiJqRZ6sPCA458ZP+1k4/r+QPQNKb7r8DJ0LQt3QsB8a9tpAbVMNnK+3vvr0bBusPbMH+jps/8wCyJsDJSvg+j+xAduCYNA2Y1rk6xfnPCA458auy8ZBfgAADNdJREFUtx1O7oDjv4Uj/20N+wVJKVB8gw26llTabSAM9Tvtdmq3BQAN2/6pWTY7Z+b11h3U2273NzwA5atgzu2QXRideiYIDwjOuSvX3wcn34Mj2+yLVQ3BHHxJhlkrYd23bDplZxPUvQt1O2D3T2DHDy/+ORl5dgVx3QYbyC1eDlNnxdQUzEQ0EWsqPwvcAzSq6rJg2xPAJmAAaAQeVtXTkS6Lc24MOprg8Jtw8L/g6NsQ6rSz/1krYfWjMOc2uwpIz7n4uIXr7T7cb91Bpz+wgdmym62Lxxv/mCOqGtk3EFkDdAA/GRIQclW1LXj8VWCJqv75aD+rqqpKq6urI1pe5yal7hZorIHMfMgqsDPw4RrcgQFIGiGFWVczHHoDTu2yGTwttXaP2rdqF90FC+6Euav/MAC4mCUiu1S1arT9In6FoKrbRaT8km1tQ55m8ftviTjnrkrbaXj3+7Drx4MDrwBJqZBbDLlllkOnoxHOn7R591lFkL8A8ufZzJzkVOvyqd9pA8ID/ZCeC9Pnwszl1oe/cL099rP6uBa1MQQR+Q7wIHAeuCNa5XBuUurvg23fhh0/sj78ZffZrbfd0jJ0nIW2MxYwzh2xWT/z10F2kb127gjUvG5BJByygd3p8+C2r8CSjTaLxxv/hBO1gKCqW4AtIvI48BXg74bbT0QeAR4BmD179sQV0LlY1XIcXv4SnN4NlX8Gax61+ffXYrSuJJcQYuET8DPgvsu9qKpPqWqVqlYVFvqUM5fAupptxs6P1lg65S/8FDZ979qDAXgwcECUrhBEpEJVDwdPNwI10SiHczEjHLKUC4WLrFvngnNH4cCvbaD35A7rHiqphM8/Z338zo2jiZh2+hKwFigQkXqsa2iDiCzCpp2eAEadYeRcXFKFmv+Erd+C5qO2rWgpzLrZZvo07LNtM5fD6m9CxXpLyZAAeXXcxJuIWUYPDLP5mUi/r3MxK9wPjR/brJ59r1iOnoJFcO9T0HbKFl3Z+wuYsQw++x0b5J3m42cu8vybys5NBFU4/n+w6zk4+IZ9uQssr/49/wSVDw5m41z9jeiV0yU0DwjORVLPedjzM9j5jKVpnjLVsnHOud26habN8emdLmYkRkCo22F9rqU3Df7z9XXaP+mRrZA3F4qWwIwlUFo1eoZG5y7V12XLL/ach5R0ux19y4JBXweU3QKbfwhLN9s6vc7FoMQICO88aYm48sph2ectre7vvgddn1ggaNgPu5+3fZPTLEfL/E/DigctDYBzw+lqhoOv26Dw0begv+fi15PT7PO28hGbGeRcjIt4LqPxNOZcRj3n4cBrsO9lG7DTAcvIuPYxW01J1b7a37DXAsexd+Dsfru8v2MLVH352lZbcvFhIGzpmuvetc9T7Xb7hm9uGVz3OcvcOXWWBYZQj+Xt9xMKFwOuNJdRYgSEodrPWiKwoutG3q+xBt74WwsORUth3lrrbhKxPDAL7oSppX94XF+XBZOmGnuvjrO2nN/C9bDwbkhJG9w33G+phGtet2ySF/LHF99oi35kFVjemZQp0NNq5VaFmcu822GihEM2E+iDFyxb54XB4OnzbfbP4o129u/jAC6GeUAYD6r2paBtfw/tZ+y5hge7BmYssy8S9ffayk/tDRcv+AGWPEySoLsZMqZbIxLqsQHGpkPQ125dC3M/ZeMcZz609xpJUootIjJrpR1Xvgqm5F5c7pEaqIGwLU7SVg85xYO3ocEqEsIhaDxgK2A111p33eyV0cuDr2oLrLfVB3l/TtnfNiPP/latJ6xrsa3epoXOW2vjUKU3Qf58DwJu0vCAECmq0HTQzugPb7VGJCXDBhEz8+0Mv6QSZiyF3BLbPhC2PPJ7XrCc8pn5dgWQX2GN+YJ1F6cSbj9rjVFHoy00Eu6zwJKRZ49P7QpWnaqG/m7LVFm83K44Os7a2EhWYbC84AJbdjB1il1pNB6AQ2/aPpfKzLfAkF1k5UnLgbQsC2gSpDbobrFjO5sssA2ErKFPmWJZNTOmA2pBrb0BuluDK6skC5wDoeDNhN8nuc0ptuNS0iw4ypAvXaVm2NKIGXlWppQM25aZb1dKhYtHD2SqlvTtzId2RXZypyV3O18P4d6Rj519G6z6BlTc6QHATVoeEGLVaGfvV6O/F06+D8fetvu0bGvML6Q7PnfEbt0tg4uRT5lq3V2L7oaChbZf+2k7Q+5osEa84yz0dtjsmL6O4MpIAbXAlFVgt9RMS52clGJXSN0tg91aucWQM9P2v1DvlDS7qiqptC9aNX5sM8BO7bL36e+1BnroZzLUFQShZttnoP/i30FSqnXdhXqs0Q91WcBIzbT7vk7obbv4uIJFNqNsahlMnW2BO7fY8v2nZlj3XFezBfMZS8fnb+VcFMXMegjuEuN5lpmSbguVzF098n6qdhbf323r1sbKAHnxDXa7GuF+q0d7g6V1aNgHrXV2JZOeYw16qNsCQ6jHphCn51ogLFoCZVUWMEeSMW18EsY5N8nESMvgIkrEzs4jPUYwEZJTIDnHGv+CClj2x9EukXNxw3PeOuecAzwgOOecC3hAcM45B3hAcM45F/CA4JxzDvCA4JxzLuABwTnnHOABwTnnXGBSpa4QkSbgxBgPLwCGSeCTELzuiSdR6w1e9+HqPkdVC0c7eFIFhGshItVXkssjHnndE6/uiVpv8LpfS929y8g55xzgAcE551wgkQLCU9EuQBR53RNPotYbvO5jljBjCM4550aWSFcIzjnnRuABwTnnHJAgAUFE7hKRgyJyREQei3Z5IkVEZonI2yJyQEQ+EpGvBduni8hWETkc3OdFu6yRIiLJIvKBiLwWPE+IuovINBF5RURqgr//bYlQdxH5evBZ3y8iL4nIlHitt4g8KyKNIrJ/yLbL1lVEHg/avIMisv5K3iPuA4KIJAPfB+4GlgAPiMiS6JYqYvqBb6rqYuBW4C+Duj4GbFPVCmBb8DxefQ04MOR5otT9X4A3VPU64AbsdxDXdReRUuCrQJWqLgOSgfuJ33r/GLjrkm3D1jX4v78fWBoc84OgLRxR3AcE4BbgiKoeU9U+4OfApiiXKSJU9Yyq7g4et2ONQilW3+eD3Z4HNkenhJElImXA54Cnh2yO+7qLSC6wBngGQFX7VLWVBKg7tgxwhoikAJnAaeK03qq6HWi+ZPPl6roJ+Lmq9qpqLXAEawtHlAgBoRQ4OeR5fbAtrolIOVAJ7ABmqOoZsKABFEWvZBH1z8DfAANDtiVC3ecBTcBzQXfZ0yKSRZzXXVVPAd8F6oAzwHlV/Q1xXu9LXK6uY2r3EiEgyDDb4nqurYhkA78E/lpV26JdnokgIvcAjaq6K9pliYIUYAXwr6paCXQSP90klxX0l28C5gIlQJaIfDG6pYoZY2r3EiEg1AOzhjwvwy4r45KIpGLB4EVV/VWw+ayIFAevFwON0SpfBN0ObBSR41i34KdF5AUSo+71QL2q7giev4IFiHiv+2eAWlVtUtUQ8Cvgj4j/eg91ubqOqd1LhICwE6gQkbkikoYNtLwa5TJFhIgI1o98QFX/cchLrwIPBY8fAv5jossWaar6uKqWqWo59jd+S1W/SGLUvQE4KSKLgk3rgI+J/7rXAbeKSGbw2V+HjZvFe72HulxdXwXuF5F0EZkLVADvj/rTVDXub8AG4BBwFNgS7fJEsJ6rsMvCvcCe4LYByMdmIBwO7qdHu6wR/j2sBV4LHidE3YEbgergb//vQF4i1B34NlAD7Ad+CqTHa72Bl7CxkhB2BfDlkeoKbAnavIPA3VfyHp66wjnnHJAYXUbOOeeugAcE55xzgAcE55xzAQ8IzjnnAA8IzjnnAinRLoBzsUJEwsC+IZs2q+rxKBXHuQnn006dC4hIh6pmj/B6iqr2T2SZnJtI3mXk3AhE5GEReVlEfg38RkSyRWSbiOwWkX0isinYrzxYi+DpIDf/iyLyGRH5bZCr/pZgv6wgr/3OIBFdXGbedZOTXyE4F7iky6hWVe8VkYeBfwCWq2rzhTTLqtomIgXAe1hagDlYiuFK4CMsZcqH2LdJNwJfUtXNIvIk8LGqviAi07B0ApWq2jlxNXVueD6G4NygblW9cZjtW1X1Qh56AZ4UkTVYmu1SYEbwWq2q7gMQkY+whUtURPYB5cE+n8WS8D0aPJ8CzObiRX2ciwoPCM6NbujZ+58ChcBNqhoKsqtOCV7rHbLfwJDnAwz+rwlwn6oejFxxnRsbH0Nw7upMxdZdCInIHVhX0dV4E/irIDsnIlI53gV0bqw8IDh3dV4EqkSkGrtaqLnK458AUoG9wWLpT4xz+ZwbMx9Uds45B/gVgnPOuYAHBOecc4AHBOeccwEPCM455wAPCM455wIeEJxzzgEeEJxzzgX+H4jZJL0dq00ZAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"for col, label in zip(rog_u2.results.T, labels):\n",
" plt.plot(col, label=label)\n",
"plt.legend()\n",
"plt.ylabel('Radius of gyration (Å)')\n",
"plt.xlabel('Frame');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating your own class\n",
"\n",
"Although `AnalysisFromFunction` and `analysis_class` are convenient, they can be too limited for complex algorithms. You may need to write your own class.\n",
"\n",
"MDAnalysis provides the `MDAnalysis.analysis.base.AnalysisBase` class as a template for creating multiframe analyses. This class automatically sets up your trajectory reader for iterating, and includes an optional progress meter. \n",
"\n",
"The analysis is always run by calling `run()`. `AnalysisFromFunction` actually subclasses `AnalysisBase`, and `analysis_class` returns a subclass of `AnalysisFromFunction`, so the behaviour of `run()` remains identical.\n",
"\n",
"### 1. Define `__init__`\n",
"You can define a new analysis by subclassing AnalysisBase. Initialise the analysis with the `__init__` method, where you *must* pass the trajectory that you are working with to `AnalysisBase.__init__()`. You can also pass in the `verbose` keyword. If `verbose=True`, the class will set up a progress meter for you.\n",
"\n",
"### 2. Define your analysis in `_single_frame()` and other methods\n",
"Implement your functionality as a function over each frame of the trajectory by defining `_single_frame()`. This function gets called for each frame of your trajectory. \n",
"\n",
"You can also define `_prepare()` and `_conclude()` to set your analysis up before looping over the trajectory, and to finalise the results that you have prepared. In order, `run()` calls:\n",
"\n",
" - `_prepare()`\n",
" - `_single_frame()` (for each frame of the trajectory that you are iterating over)\n",
" - `_conclude()`\n",
"\n",
"Class subclassed from AnalysisBase can make use of several properties when defining the methods above:\n",
"\n",
" - `self.start`: frame index to start analysing from. Defined in `run()`\n",
" - `self.stop`: frame index to stop analysis. Defined in `run()`\n",
" - `self.step`: number of frames to skip in between. Defined in `run()`\n",
" - `self.n_frames`: number of frames to analyse over. This can be helpful in initialising result arrays.\n",
" - `self._verbose`: whether to be verbose.\n",
" - `self._trajectory`: the actual trajectory\n",
" - `self._ts`: the current timestep object\n",
" - `self._frame_index`: the index of the currently analysed frame. This is *not* the absolute index of the frame in the trajectory overall, but rather the relative index of the frame within the list of frames to be analysed. You can think of it as the number of times that `self._single_frame()` has already been called.\n",
" \n",
"Below, we create the class `RadiusOfGyration2` to run the analysis function that we have defined above, and add extra information such as the time of the corresponding frame."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"class RadiusOfGyration2(AnalysisBase): # subclass AnalysisBase\n",
" \n",
" def __init__(self, atomgroup, verbose=True):\n",
" \"\"\"\n",
" Set up the initial analysis parameters.\n",
" \"\"\"\n",
" # must first run AnalysisBase.__init__ and pass the trajectory\n",
" trajectory = atomgroup.universe.trajectory\n",
" super(RadiusOfGyration2, self).__init__(trajectory,\n",
" verbose=verbose)\n",
" # set atomgroup as a property for access in other methods\n",
" self.atomgroup = atomgroup\n",
" # we can calculate masses now because they do not depend\n",
" # on the trajectory frame.\n",
" self.masses = self.atomgroup.masses\n",
" self.total_mass = np.sum(self.masses)\n",
" \n",
" def _prepare(self):\n",
" \"\"\"\n",
" Create array of zeroes as a placeholder for results.\n",
" This is run before we begin looping over the trajectory.\n",
" \"\"\"\n",
" # This must go here, instead of __init__, because\n",
" # it depends on the number of frames specified in run().\n",
" self.results = np.zeros((self.n_frames, 6))\n",
" # We put in 6 columns: 1 for the frame index, \n",
" # 1 for the time, 4 for the radii of gyration\n",
" \n",
" def _single_frame(self):\n",
" \"\"\"\n",
" This function is called for every frame that we choose\n",
" in run().\n",
" \"\"\"\n",
" # call our earlier function\n",
" rogs = radgyr(self.atomgroup, self.masses,\n",
" total_mass=self.total_mass)\n",
" # save it into self.results\n",
" self.results[self._frame_index, 2:] = rogs\n",
" # the current timestep of the trajectory is self._ts\n",
" self.results[self._frame_index, 0] = self._ts.frame\n",
" # the actual trajectory is at self._trajectory\n",
" self.results[self._frame_index, 1] = self._trajectory.time\n",
" \n",
" def _conclude(self):\n",
" \"\"\"\n",
" Finish up by calculating an average and transforming our\n",
" results into a DataFrame.\n",
" \"\"\"\n",
" # by now self.result is fully populated\n",
" self.average = np.mean(self.results[:, 2:], axis=0)\n",
" columns = ['Frame', 'Time (ps)', 'Radius of Gyration',\n",
" 'Radius of Gyration (x-axis)',\n",
" 'Radius of Gyration (y-axis)',\n",
" 'Radius of Gyration (z-axis)',]\n",
" self.df = pd.DataFrame(self.results, columns=columns)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because `RadiusOfGyration2` calculates the masses of the selected AtomGroup itself, we do not need to pass it in ourselves."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Step 98/98 [100.0%]\n"
]
}
],
"source": [
"rog_base = RadiusOfGyration2(protein, verbose=True).run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As calculated in `_conclude()`, the average radii of gyrations are at `rog.average`."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([18.26549552, 12.85342131, 15.37359575, 16.29185734])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rog_base.average"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results are available at `rog.results` as an array or `rog.df` as a DataFrame."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"