Tica_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import pyemma\n",
    "import glob\n",
    "import numpy as np\n",
    "import pyemma.plots as mplt\n",
    "%pylab inline\n",
    "import mdtraj as md\n",
    "import pyemma.coordinates as coor\n",
    " \n",
    "def average_by_state(dtraj, x, nstates):\n",
    "    assert(len(dtraj) == len(x))\n",
    "    N = len(dtraj)\n",
    "    res = np.zeros((nstates))\n",
    "    for i in range(nstates):\n",
    "        I = np.argwhere(dtraj == i)[:,0]\n",
    "        res[i] = np.mean(x[I])\n",
    "    return res\n",
    " \n",
    "def avg_by_set(x, sets):\n",
    "    # compute mean positions of sets. This is important because of some technical points the set order\n",
    "    # in the coarse-grained TPT object can be different from the input order.\n",
    "    avg = np.zeros(len(sets))\n",
    "    for i in range(len(sets)):\n",
    "        I = list(sets[i])\n",
    "        avg[i] = np.mean(x[I])\n",
    "    return avg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['/net/jam-amaro-shared/bccgc4/CatSMD/md1/1.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md1/2.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md1/3.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md1/4.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md1/5.nc'], ['/net/jam-amaro-shared/bccgc4/CatSMD/md2/1.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md2/2.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md2/3.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md2/4.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md2/5.nc'], ['/net/jam-amaro-shared/bccgc4/CatSMD/md3/1.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md3/2.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md3/3.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md3/4.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md3/5.nc'], ['/net/jam-amaro-shared/bccgc4/CatSMD/md4/1.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md4/2.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md4/3.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md4/4.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md4/5.nc'], ['/net/jam-amaro-shared/bccgc4/CatSMD/md5/1.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md5/2.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md5/3.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md5/4.nc', '/net/jam-amaro-shared/bccgc4/CatSMD/md5/5.nc']]\n",
      "['RES_DIST (closest-heavy) PHE146 - TRP187', 'RES_DIST (closest-heavy) GLU116 - SER136', 'RES_DIST (closest-heavy) ARG142 - HIS143', 'RES_DIST (closest-heavy) PHE71 - PHE212', 'RES_DIST (closest-heavy) TRP27 - PHE29', 'RES_DIST (closest-heavy) LEU117 - TRP187', 'RES_DIST (closest-heavy) ASN64 - THR160', 'RES_DIST (closest-heavy) TYR119 - GLN20', 'RES_DIST (closest-heavy) PRO216 - PHE191', 'RES_DIST (closest-heavy) GLY24 - PHE146', 'RES_DIST (closest-heavy) TYR114 - TYR119']\n",
      "<pyemma.coordinates.data.fragmented_trajectory_reader.FragmentedTrajectoryReader object at 0x7f1d2887e400>\n",
      "trajectory length =  90000\n",
      "number of dimension =  11\n"
     ]
    }
   ],
   "source": [
    "#Makes a list of nested lists for 5 mds, and 5 .nc files in each\n",
    "trajfile = []\n",
    "for i in range(1,6):\n",
    "    nc_list = []\n",
    "    for k in range(1,6):\n",
    "        path = \"/net/jam-amaro-shared/bccgc4/CatSMD/md\"+str(i)+\"/\"+str(k)+\".nc\"\n",
    "        nc_list.append(path)\n",
    "    trajfile.append(nc_list)\n",
    "    \n",
    "#print(trajfile)\n",
    "\n",
    "#assigns the features\n",
    "feat = coor.featurizer(\"/net/jam-amaro-shared/bccgc4/CatSMD/md1/system.prmtop\")\n",
    "feat.add_residue_mindist([[146,187],[116,136],[142,143],[71,212],[27,29],[117,187],[64,160],[119,20],[216,191],[24,146],[114,119]], scheme=\"closest-heavy\")\n",
    "print(feat.describe()[:13])\n",
    "#selecting distances at minimum & closest heavy atoms\n",
    "    \n",
    "#feat.add_minrmsd_to_ref(traj[0])\n",
    "inp = coor.source(trajfile, feat)\n",
    "print(inp)\n",
    "print('trajectory length = ',inp.trajectory_length(0))\n",
    "print('number of dimension = ',inp.dimension())\n",
    "  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c16276aed4b4cf18078ea2b492b9fcc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='calculate covariances'),), layout=Layout(max_width='35%', min_width…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "807907bafc744624b9973d8a0fc128a3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='getting output of TICA'),), layout=Layout(max_width='35%', min_widt…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "ename": "AttributeError",
     "evalue": "'list' object has no attribute 'shape'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-66-02c617ceadfa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m# have a list of trajectories, so we just get the first element.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mY\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtica_obj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_output\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Projected data shape = '\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Mean values: '\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'list' object has no attribute 'shape'"
     ]
    }
   ],
   "source": [
    "lag=1000\n",
    "tica_obj = coor.tica(inp, lag=lag, var_cutoff=0.95, kinetic_map=False)\n",
    "# here we get the data that has been projected onto the first 2 IC's. It's a list, because we could generally\n",
    "# have a list of trajectories, so we just get the first element.\n",
    "Y = tica_obj.get_output()\n",
    "print('Projected data shape = ', Y.shape)\n",
    "\n",
    "print('Mean values: ', np.mean(Y, axis=0))\n",
    "print('Variances:   ', np.var(Y, axis=0))\n",
    "print(-lag/np.log(tica_obj.eigenvalues[:5]))\n",
    " \n",
    "# histogram data\n",
    "z,x,y = np.histogram2d(Y[:,0],Y[:,1], bins=50)\n",
    "# compute free energies\n",
    "F = -np.log(z)\n",
    "# contour plot\n",
    "extent = [x[0], x[-1], y[0], y[-1]]\n",
    "contourf(F.T, 50, cmap=plt.cm.hot, extent=extent)\n",
    "\n",
    "#making clusters from kmeans\n",
    "cl = coor.cluster_kmeans(data=Y,k=10)\n",
    "\n",
    "#cl=coor.cluster_uniform_time(data=Y, )\n",
    "# for later use we save the discretetrajectories (dtrajs) and cluster center coordinates:\n",
    "dtrajs = cl.dtrajs\n",
    "cc_x = cl.clustercenters[:,0]\n",
    "cc_y = cl.clustercenters[:,1]\n",
    " \n",
    "#plotting out the results and clustercenters \n",
    "contourf(F.T, 50, cmap=plt.cm.hot, extent=extent)\n",
    "plot(cc_x,cc_y, linewidth=0, marker='o')\n",
    " \n",
    "print(dtrajs)\n",
    "print(np.size(dtrajs))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ac6e2e1316fe49d9ac076a59cb6634ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='initialize kmeans++ centers'),), layout=Layout(max_width='35%', min…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "52a9071ef6cb4b8680213558506b57b8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='kmeans iterations'),), layout=Layout(max_width='35%', min_width='35…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10-11-18 21:39:37 pyemma.coordinates.clustering.kmeans.KmeansClustering[76] WARNING  Algorithm did not reach convergence criterion of 1e-05 in 10 iterations. Consider increasing max_iter.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a6b4f0ce73445f3baaf4367a82e32d2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HBox(children=(Label(value='getting output of KmeansClustering'),), layout=Layout(max_width='35…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([0, 0, 0, ..., 8, 8, 5], dtype=int32), array([0, 0, 0, ..., 7, 7, 1], dtype=int32), array([0, 0, 0, ..., 6, 6, 6], dtype=int32), array([0, 0, 0, ..., 8, 2, 2], dtype=int32), array([0, 0, 0, ..., 6, 5, 9], dtype=int32)]\n",
      "450000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnX+sXdl11787Mx478/rsV/nN2MRjzwNR0akH0mKrIEB6QKdMJmlLWorUVsCDCFmV6NBooCkhQJEACSYiApUiYbUjXDW0IEpJxQSlCT9eiVBK4yrTjuuEtjCjjFO7eab+UWf8MtNs/jhn3bPOumvtvc+P+959966PZN17z697nu+5n7Pu2mvvHWKMcBzHcRaHt+z3CTiO4zjj4mJ3HMdZMFzsjuM4C4aL3XEcZ8FwsTuO4ywYLnbHcZwFw8XuOI6zYLjYHcdxFgwXu+M4zoLx4H686fr6etzY2NiPt3YcxzmwXL58eSfG+Ehuu30R+8bGBj796U/vx1s7juMcWEIIr5Zs56kYx3GcBcPF7jiOs2C42B3HcRYMF7vjOM6C4WJ3HMdZMFzsjuM4C4aL3XEcZ8FwsTuO4ywYLnbHcZwFw8XuOI6zYLjYHcdxFgwXu+M4zoLhYnccx1kwXOyO4zgLhovdcRxnwXCxO47jLBgudsdxnAXDxe44jrNguNgdx3EWDBe74zjOguFidxzHWTBc7I7jOAuGi91xHGfBcLE7juMsGC52x3GcBePBoQcIIRwB8PMADtfH+/cxxh8aelzHWUxeSKx7z56dhbPYDBY7gF0AfzrG+DshhEMAPhlC+M8xxk+NcGzHWTBO7PcJOEvAYLHHGCOA36lfHqr/xaHHdZzF5OR+n4CzBIwRsSOE8ACAywB+P4AfiTH+grLNBQAXAODMmTNjvK3jHEDOJtbdB3Bkr07EWWBGaTyNMf5ujPHrATwG4BtDCE8q21yMMZ6PMZ5/5JFHxnhbx+nA/cS/vTwHx5k9o1bFxBhvAfjvAN4x5nEdZ3bsRYRcegPxaN0ZhzGqYh4B8EaM8VYI4a0AngLwTwafmeOMynVl2cYevK8U+n6dh7NMjJFj/z0ALtV59rcA+Hcxxv80wnEdZ0Rusufn9uD9LKHfFMv34lycZWOMqphfBvANI5yL4xTQp4HxBTRlhu8a8TxKuY5G6CT4sc7DcaYZpSrGcWbPfeV5ieBfRCP1GyOfRwp+bhvwdIuzl7jYnQOAJdNU9H65fpR148+jipo/NOJ5cLwB1Nl/XOzOHFMiUi16fwXAcWN7itqfQ7nc+0TpjrN/uNidOaVrzTffPtW78+0AXqqfb6GK3j828Dxc6M584WJ35pBZd+Q5AWAHTVS/VT8/iapn6En2OsdYUn8R1a8JHwjMGY6L3Vkw+tSJvwxgveDYUvRU7XId/atcLqM5552ex3CcNi52Zw65gvKImeAlhRxeJ75RH1tjp15Hon4STTUNnQcd/7h4DVQllV2j7VfY8zEqdhynwsXuzCEvoZLqTaQHzQLaUTOXoyVZiqw1we/U/9bZ6/X6uNpwuxRp36i3fa4+l0uZc77P9qXz3oH+a8NxuuNid+YQmZKwoncpdZJjSbXL+9A0nmrvT9H7cbaMhE+iv8HW8eM8jSZfv15vS7l7jvYLw3GG42J35ogt9vw4KmGSREmyJ2ELnZY9jXSlC/FyYh1F70AlZ56GOcnW8aEC6PlJtm6drefLrfN3nOG42J05gaJc4iYamXPBUwqFhEivr6Mt4nNoOilJHutxfprk6X0tbqK6eWiNoi/Bhe7MChe7s888nVlPggf0iJ2kSaKnnDhQCXwH7fLJPuWJMg3DJZ+DS5uncKTUPb/ujIeL3dlHclKXkORJ6Eyy93aBWwDWrlWvVw43u90LzXO+fLJ+t+C9+XG52E+gLXitLYDO+yz0KN2l7oyLi30hyXXw2Yueko9BryQpiXI5ssabNViSkG8BuFsvpsdVJus1driUxG8VnM4tftxrwApF81LwGrISRxP6FvJVNY6TxsW+kJC492sqNpJ6icR3xPN1Zd2N6WVS6lLKt9AI/S6AVbQFLymRusa93VruJfDGU4Lvu44q3fQ8qqodx+mHi32hOYJ2JxhiY0bvR0Ln3fIBu3HQKjUE2qWEAkvqd8V2JHROSu6SSfSvrJv6FSCjdw7drCyha7xYP/q47U53XOwLz0b9eFlZNibn0ETpcqwVavAE2pKn9VxyNzAlPS19wtMvgB6103ZraCRfEr3T8fhzTe6STtF7CSdR3Zg3Rjymswy42OeaLhNK5DiHJgrs0/2dekveRFOqdwVVOR9Ve0ip8x6cJZ1zaqnnGjO51HnULiN2IiVmTfDacbrIHdf0RtoJWtsDYA+h4CNMOt1wsc8tfWcMsngebdF2zePKtIl8nRrrhG9rVYMwqZfmu6XUb3fYR0btd8U2HJmr17aZgt2cVmQenf6vTojlgJ22yo2b02fKQGdRcbHPHanorO+X9/n6UaYJnqsfS7rgb9SP/5UdR1asUI6ZHq0BtwgqWayfk9StyJvTVeq363/H2DIrCi95/05IyecalXlZJw1GxoUvJxEh6Y/5C885yLjY5wZrVnuiz5f3ufwmuIluJXaUwqHu/zJSl13xOTIFIVIvMsWSQpO6tQ+PzEnwudQKHYvy86i358vlthpyn0keXg4uZtW/b6NqjOZR/nVje36NuNyXGRf7XKBJnfdI5I2QfB/ry0sTR6SQP/mfRiVkqxu+5BImXfPNnLjsLCSFL6ROpEQpc+sl+2ipFS54QhN9TvDa+aXWA8CKXFBa/07DK6RGnPTo3RlB7CGE0wB+HNUV9RUAF2OM/3zocRePkgYwbUApHhHzgbDomNqAUsdRNv6IVm6YGmNFcK0W9xraIptKZew20lwDJo2LUuqp0kVCrr+d2SclawlF/sfE8pyoJdp5rInHJPJz0dI3Kbk7y84YEfubAP5GjPGXQgirAC6HED4eY/zVEY69QKQ6DclSQC5poJ2LlvlWa4RDwmp0S9VUa2OsCK6xbvqaWK0OQ5OIdre9T0rO8lhDpK4du0vpYwl0PH4zo8balcNoS7pkgg2t2oincI6z5xyP1peVwWKPMf4mgN+sn98NIVwFcAqAi13lCOyIWE7cQKI+jiYlwyUP6ELnKY/cOCRS6uzY9wKwEturPxvyDY5a9M6jXi7TVGNpSuh8/dDGTil32cgqSXVckiSlTtD/ubZOKx89gUbi2lj1LvRlZ9QcewhhA8A3APgFZd0FABcA4MyZM2O+7QGEpmt7USyXM+mQrEnqcmo2vo3W/Z4jpaF1pFF6e94LlfhOxUrqQFqkqV6gqdclefUuQi+RLh2DbjqlvVJvs8djxnuRzKekzifo4Ghy7xqlu9CditHEHkL4KgA/DeC9McY7cn2M8SKAiwBw/vz5KNcvJ+9C1VmIkNG3HEcFmI7CjcGxJqMQlozXoohG5r8/G6a3kZSkZLoghW69RwmyHYAfQ5O7jNpzlTdAI/ik1MEelTFwAPRLvbjUnYZRxB5COIRK6h+OMf6HMY65HLwA/UvO4ctPoGgs8FavRxmZ595PieT7pjo0CXY5libT0v1luqTrTYcqZ1Lb8ghfypx3gGoNNSAj9h1lmy54iaMzzRhVMQHAjwG4GmMs6eniTHqB8lwpkMx3t1gXz3fQCINeS3pIQ5YgaqkKEmeqO38Oq3NRaXTe5VeBtW3fXxaUiuEyP61s1xL3Tnv55P1L5M5nk+JtMCfhJY4OMUbE/scB/EUAvxJC+Ey97G/HGD86wrEXjOfQ/mnNf1YDjch51YuVTpHzapbIHSiK+Dm8BDGXh9bSHRazkHnJTSSXl5fHTY0xUyp1guTO6/6n/t7SyJ3Xs7+9XuY17E7FGFUxnwRQkIBddiyp84idp1ospJBp5MQCud/Txi9RBK+N16Llobn0aOREQhOvJnMr720dQyMl65KKGz4WDJAXPP/7SeqnIfLpSlROz7NpIU3uvKOarJR6CdMTfgOLnZ5Z5L9tON7zdOZQt34+S30Jddf8qS7oQFvGiRmGCK1352T8Ej4SoZA877ZP+3apILHGcLFk3jUdMlYjrRwmoCSqp5p8vo/2f6L93/cuz5RyT9HlWjtIvILluHENw8U+U+RYLbxruIzeX0I1BC7NbH+lie5u1UPBrsGItI3osLSb/a3dpkcoF7nWE5S/5g2Ecr9cw2Ou7r0U62/MHcfqkVpa8gg0f2Prl8tuM+9q6teIZKqKRoPLmlfGaKm9RRTei2j/H3i7goWLfSbQaIp8ggkOjd5HY39cqbe7Uv8zhq+dSB55eZTWcfNjp4RuRd5y7BarLDBVt56raZeyHVptI7fv0kNVpqTo7/08WzfK6JB8HlVN6IAeld9E009iUeBDTgP6+EkevXNc7KPzQn6TCSR4oJmwQjSupapNUtEqj6Zz4ioVunx+TNnOEnYX0Uv6RuWl8JEeZdSujQJpyf0u2o2nJdG/zN1P1bxztPp1vu2iCR1oBrTjPa+BpjLIq4I0XOyd4eOnyFEYS8b94PsC7QvVGO2QkBFyCtqWC4hLREqcHzuXRkmtGyrzIbLO1Zp3OY62H1+uyf0Yqsh9VWyjnYc2F+tUCoZHqccxDa2/ge4zYh0Enhavec9r7ZewR++Ei70zR1A2UiOgX3yFQpeCtAbLssoGKZrWJEOv+bHkcfrmvftG5n2EXtJgKhtGS0ht2yUPf0vZXt4YpuCy5q8pMpWcwOJOeG01AJPgt6G3Vx1n2y3ir5g8LvZedI0GeFQPFM0cZEXsQF7qAPAq2pGkJvaUzPn7lDAk1dKFsY7D6doe0ZVcm0ErBcPhvwCl5EhkZwef3vzCBW3B05lc8kD1//MCqu9bl2kgDz4u9kFsoGzcc0CV+r1rejQu6SpYbQwUoGzWnyFC7rvv0J6lfH3X1IuWEhmCdZPoWnEzgeQuxU7R6aKWNQLtv01+z6xfw1R5JvuJPIeyKSAXAxf7YOinnhypUY6xDjRSF1G6jNY1ctG6dYPgUi+9QYwZufelT2MqrdMiYo2+kbr2Hqlhfvl+xBp7NIfzJW6w9WfRHiRsLylNQWr0yXlTRy9Zu2/NJcDnMOCNyvT/JXP2QNmN0Ur3lOy7P7l+F/socKnzC41fiIbUP4/y8r5U46ZW0cLpUt6n7d9H2NYNQaaExngvuX9JdGzmuWu0dVbFTG4Md42W1M+iEsXLbANt9MfjaEb43EuxD5E67X8FZTnvLbRTTNeV5+L7BLSDm9VrVUlwqyewdvO05o8lSOD0+ZQKHdjPBlwX+2Auo6lX58PuAm2ho3p+bVcvL+RYEbNVrZKTumRo6mPIsbVtur7fWo99UudSOmGGhVYtYx2DjymzchiVMJ5EI7LjmPRlmCp5pD4PN9n2s24g7FookJPe86j+Bmvy9KfrY1xhy+R8A5kJ0DunvCy5H0f12XSN0oH9rspxsffmFfaccns8x5doIE1JnVNSfqhJfQzpjVEj3vU8lFTFfzy1iQ+e3cIX3rqOt72+gx+4cgnvvrbd7vE6xjlKGVtlirnjaA3VdJyk1EnUlDag60mOGbOD6fTDrMReKvUraI9XU9Kgew6tmcTuBWDlVP1C/n31MinzVDCj3lz5wHkSKXcu9bM4KEInXOy9uI8mSpdSLxA6UNZ4WdKLc6jUx+rk0/W9CxpE/+PZTbz/G57F6w9WX5ZrDz+K9//hZ4EV4N3/e7varm/0rrVDyHFwuqZzuNw5XOqrUKT+TWhkSNfROpohJjTBv1wvn0VVTJco/Qqq6/3l+t+T7LwyMrwXmrr/NVTFBACA+pEmPJfpxpJrnbZZKfxTWrn8gy11wMXeA7ro+UUtpW5EF6kLsiTdIklJvW9D6SzoWe3zwSe2JlInXn/gCD74+7YasWvk0ibyeUkaxULbno8hI6W+BlTSXkeTL5ecRHUdlQjlCqr0xcc6nPTQfDmlG19Ce0gMYDpNVHDj+Xz9b+pzM9ItslRXa9tYrf/d263lLqcc1PqT0K8lOVxBjvkROuFi7wxdEDxSkaVY62ywrt32Kqths4/Q+frSHHtJTr+UkvJJjcKbzhdW9EqRLxxmy3Pvm5J8LtrPNcL2KZWcTFl4xdiAejDzRtSSRtKt+tHKXQPtkRH7wsc0UgIanEV7nl7+SPv9lyo6p8IBWZ4r4YGR1ndDSp6KEmjf09dquWspHs4VtBupgWase4uNzPr9wcXeGfpCcqnznJ0YD33lMFpyv4tpoXSZAs6K0nOdm7RjDGXGaZ+33dvBta96dHr56zv2++fkKv/vte2t1AwnJXVtAo6p99pB1XNyB9U1RFUbXOic0goYq177codjaGhRujK+0Qqlkejv4bnr/1LtR/03NKyGdS34kZ+LjNxpYDagqpI5VTIH8I36H/3SsPaZ7yEcXOydociEV79YKHI/jSaaALpVtKRKG+UXZa8bUPuUUGb4gU9ewvuferaVjnnrm/fxA1culUndknxpXp538MpF/jnM/x8euVvi7Rpl0+ii78P0ULd94NE5j9Lrc2+NRFoPMT2ZQ4CLcbuqCgPy/Tfk9V06fhHBP7tbqAdn2xUlkCnob3wS05/L/Hd0crF3hiIqrbMEIScsRiN3fsHxKKRE6jLPSM/ldpIxG0gtZvAe7/61beCtwAfPbeELK+t4270d/MDVS3j3lUR+fSxI5prctdmUiK617ACatB4fppdkwrvId4VPlt6HwihdE/QtGpeeTeQipa6REjpfXwJP7fBhldd26zp3ZCRP0Tv9H142tps/XOydsYReME+ldtGmohZtH0vqBWmcL39petlDD9tvrW1vkTrOEN79y9t49/8tFHkusrbos4+2H0l9VVmXRdRpT43D3lfuUk5y+jwLLUpnOWitWgXssXVTZELn26TSiSVC7/p5qzcFPkOZxglUUXuq7WL+cLF3RramF6D1NO1az576EvHtFFKC7iLvFHScEsHz9xz9hsC/7NoXf1ZVQFzq2jlN2E3MkMTR5jzVIEGX5o6JXCRPbUiycRS21Pn1N8ltK8v6iNyqVeefsdb7V2uQnTqWJXfqX3DwxuNxsXemROhGHTuXempkxlQDkib1RMSvibukqs86ZG6/L3/JlrV2LqntWyfTpQxRyt1CSr/kPbSOS1LqqQZZAP3lnqM0qt+B3UhL0BSNyuQv2q9GK0ApkTjfTm5rbUOUyF0eg35N0Ou7AE7vipp3kvqTyFfGzB8u9s6QtHMRD5M65fe41GW03qVOvaCx1RK69rwLWiBmvTcJO3dzGSx3/sdYMyFJuqZKrO1LpM6ZnCtVkWgTifPUSVe5y2NxqNGfl/VZGEIH8mnBkk522rbyfVLb81FLc3JPlVHSvqtgUTuX+p/G0kbsIYQXAHwLgN+KMT45xjHnE94Rgf90U8b0KJE6YZUqassTUrfSKl0lnvoxQd+XLoK3zoUfo0sqZ+pAFprcNfGmShdXledraEoaU8dNnQ89v1cSvackzK9DPhkHl7vswXqlLWyCn0dO6KnXtxPb5Ci9WLmUge6Ru3a8ezxqP4n29/pgzco0VsT+rwH8CwA/PtLx5ozH6kcq3yJ22CP7IvWJ1FM/RTNSzwmUk5J2CbS/FDxRElT3JneQlKC7Cp2vXxPPHzO273r8ltxLaqxT8HpxaywUPtyFgpS9dr3JKF1Wr5Q0jGqM2fahyT2XzqN1964xufPetLzkcWP4Oc6YUcQeY/z5EMLGGMeaP85BLV+cQjQs3UVa6iVCl+sVclIfKvI77PlRcVwtIBoq8KK0jIYVPQ8ROq2j14+j3fGoyx9L7yHf25Q7T8d0gSTEe3wq/S1oHBYN7Zq0ApISoXdJw0hSMpb5cska0nK30mf3rrGe40B7hEfgIEy5t2c59hDCBQAXAODMmTN79bYD2EL1pZBRuoGs6ZUNpSmpWxG6toxtK6VeKvQ7iXUpaD8SvHyPXiXc6N422qI0Ss+NrW4di0Su9SZdQ1mjXy6yH03uJHMu947kBM0b/q0bAN9eW9cH6wLpKndtvVaeem8XWOE9g6lElCbboDkY5nO+2T0Te4zxIoCLAHD+/Pk47tGHDGpEEQ2N0UENTAXICaiB5oLuKvXURV/whZABFNBf4Dmk4InSXwe5G0Bx1N5X6iVC52hSp4g9J20rUpeMmpYB2lG6vJ6VYCUXpWtCt351anIfK9WSKie1BC/lrrW9aNzbrYdI4J3HzqKden0B8zi8wIJUxRzB8BHrJImJp4nUl2EMqRvRoBWp0/duiNBzcuZStgTflc7Rei76LvnSyuNYKRvrWJZAUu9RQlEZpIRPBycnvOhwHafSLlqjf0n6pTTH3pVcf4FVsVxeHyWfy1SvVPo1RAy9+c6OBRE7UMn9lR77yWidT3AgOmVwchd1H6lrgk8Ina8ujZbHbjwFKsEPlXsRJemU0mjcyq+WFvinbh4l5zkTNMHTcmPCCiIVpefahCyhZ67dHA89jHb1yxp7HJS/q7EE3ypBpVQYFftRSma+Gavc8ScB/EkA6yGE1wD8UIzxx8Y4djc26scuYzpoUq/rfK3onNCqA2RUw7Gknjh2Ko9OlKRfhgrdOh4JvjR6l2kY7XulpmH6DJEr0aRuyYG27dJwUFrLbqHWtKfQGvRPYLpNKCH1rkLPpV8KhV5STTUpgdUPkd6Zw9tIrPVEK0rnY/cAez95eH/Gqor57jGOMx7n0J5gOoUcfreW+jWjYsD6KWr9bJWiz/0kFULPbZ6Tuib0vr+Kte+FJnhL7n0bV4vrxUuEmpJ6XyFLofdKp3SFyRpAu8FVIyN1Ks0tFboWnCjBSOm1Zm1HfRwm0XuKkgi+U5SuSX3+o3VgoVIxEmqtfsFYL3v2Ufql/gKU1N1qFzvvnCG3116LL4r1pchF3FLqYwpd219+N3jb1Fi5906kUir8y2zJO9Vgqh1Pvs+epFqAtri1PhXEenu5JnUrSueTVBCa0I1rd8xyW2KVd2DT0jMauSBgSvCa1HnqReuQmCuD3h8WWOwpeE0qfTDUKLIO4FrzoXcJOW6L16m0zIjkpL4Hp1BESWfAbDWMVtGgiZYvtxpBrTwtSf20cnyOWQdt/NoDOo4Rsw5dHEa7T4t63lDrVyWtSwlda/hMROiW0K30YMmN/6ef2MQLm1u4cXQdv+fODp7bvoTvfHV7Wu59WQNLv2hSfxJtkfNG6fksdQSWTuxc6MfZMnp9vXq9sgOg7mikdUSRDZ7ywr8LO4oYkYMidE7Xnt5ZUlIviaal3LnUtWg99f4lFA0hQHDBs4i8pGMRkG+4LxW61sB/d1ro8vrLVWfl1v+vJzbx/DvfizcfPAQA+MKxR/H+d74X+Ch0ueckT9fEVLSuSf0sW87HxweqIHB+pQ4sjdjpQ5MTGABNpC6nKLvWXCya3OXr0uECRiIl9ZTQU6dS6qihxQicMYobAExLXP0CszfVOIZK6FLsgP4f1zdSVOUuUictMpVZHJlioWVaDXoqmufbyoZSVFLXIvQ7yrIStJv9Dz91YSJ14o0HD+EfPnUB3/bD2+20TA6rUqn1OUipvx3tVMx1VAUaGwVvuL8sgdh5lE7w1xSp03Pj568VuWuNpRYjhdCpL4/2FqX3Fdquj6+6ROGjRe1alQu9XlO2A5obNWCXZmiNoDR+SDL9MSNy1VlAOmdO64FBMgfsKN26JlOXPHet2ib0Vj1Zc6te3urIlorWiy42nn7hETx3xUbJgeaCJRC7lLocN+MmKpFTD9S6AVWilYX1xWg41d4OsPOVuS9Q31McmrYcFS1nXpIiKT02/ceV3mmsERA5qf9AMw2TGLaC9rHeLyX11C/IDtG5WDQl9NS1aAUMqcb4ElrtMWtoXx/yJm2mYYgdNK7YwfRQDgdnZEdgKcSuweX+EppyR5pYoMZqbAJbblXBjBSd77XUx2D7iU38xOYWdo6u45E7O/ir25fw1NVmeruUQ7PDCYwldSLVuzSZD1+v2mK6RPC9yyApTXNtelUXqWtRuJE/B6aF3kfmqWX8v13LqKy+fgd3H56+WtZev9NcI/Lz08pZrRLXW0AzexJQ/T/zXPrB5S37fQKz5yz0D+s6mkl6Re9TbYRG2Qglq2AsOuTbU73zxpL6rIW//cQm/uUzz+KLxx5FDG/Bbx17FP/0mWfxiSc2hx2YR2RDpD74p8g6GtGuT8t69MZz1llp5XD+/LU+FFLw8lq81d72y19qp1zoUr9T/7td/7uL6a/ILUz/6dZ/hdxWXsvf84mLeODNL7eWPfjml/F3PnGxeqE1lstoPTVlIZ3EvV00v9Qp6JvPMsZSliRilw2l1hACaEsdyEftEmpwzZHYxqoysLYrOawMzHJ+6JuO+YnNLeweav9s3T10BD+6uYWnrm4XZTwmUTv/gpZKfdV4zin9jAARtRsTq6yIYIAfeypfbzSSUuRf8l7YnZa3vC6tNAzQFjra67UcOo/QtcNYr1PrZLMHveaR+x+7uo2HAXx4cwtfPLqOE3d28Dep3JEfiEsdmI7WOwUE19Hk2en1weiUxFlgsctBwaTcDVYOYzKrekrq1kh3XVF+9kpKKw1KpE78zyc28cnNLdw5uo5jd3bwzduX8PUsXdKXnaN6d/gvHl0vGk4A6CF1eQdaS6zTlsvRG/uwchiT6e7AjkPWat0gEqWLk+2oOiszdK+VgqH3hvJaWZdqFL3NNtUOwU+FKBnMEmg+Bkvum1e38a0sIFgFAH59SKiqaQ3VhChS+hJafm8XWOkyFMl8s6Bit0Z6PIl2bzJteN4dvcyNX8UpqQ/4CW41mGrLutxHZID6q09s4ueeeRZv1pH17WOP4iPPPAsALbnLL18J63d28MVjj04tP3GnHaEmx4jp0ptTy6/K57mfCakvPpCQLRskDpiWO9GSe8EYRLd261I8GdVfsyN0QF+npV3Y8lyVC4/Sta8C3ze1rMv9krbN/rqT0foq+kl9ZeRRxOeABcux30fZ8L1c7jRxLQ2iv16WyxwaqbN9rZpgYDZjqn9yc2sideKNQ0fw8c0tdfvcF5l/Af/C9iUcfqP9GRx54z6+b/vS5LX8kj/0cELqVrTOI3m5XHteQq4hFUAjciF1gq4deRz6T+Rd+2XeW+a/r+1W29+7pkud58YNkzW4AAAgAElEQVRlJKBF7uz9KI/Oc+gyj/6aeBvrlEsYo5bAvDmQwKnvAZe6tjP/fBZQ6sBCRexdx2On6F1CX1ZlWAHt524fmCllCmZWUudR+x0jXXLbWN6Fb726jbcC+FGWF/2+7Ut4pv4loEl9aoWUukQunyqnYNukRvXjxysJKVuRe4JJykWk9JSc9gTrlx7PS0ip07LUoHP8vZVAwqp0sdIupafN0VIsfP/U/Tc57ICM0rnU5ZvL1wsqdGJBxD5kko3jKOqgVEKXNMzd6aep/PkYgymR3I/e2cEdJV1y7E5GWIU8dXVbbSjtJHVrx5TUpaBzv+X5r4HS6H6qkZPDb4w7zbyi/K5qpe1ygQLfX0qdo0X09fJcwyhfl7sXdbnUZf4c6Jaa4R/j5JrRftHxDZM9TceelGf+WBCxd0Hm1WnqKxqytx62F0iXrWk9UUckFa0Peds1AH9m+xI+8syzeIOlYw69cR/fzNIlKSx3HjWWy32AnhNWl0p9DXaZm8zJp4SeCyenyExppzViliSotcZQ3s6jVcckpG41jPK30E5R+ypop685NfdfKTNuyetH1qfTidxG+9eNevdYB/B0/fxjiTM62CyB2EnkshqGlpPMxQQbPH/Jr2g+b+KM5U50KZcvgRpIP765hdsjV8XsGdowAvw1X6Ytl/tpMl2F3YLcZTAvLach3zP3AWrVWAVCB6ZTL7MSeh9yaZipaJ2nXSS32A7a571yGE2F0VkAzwN4X/eTPgAssNhpiAB6zuGTFFxB0+NU6ZxkQT+vu8i9i31rxhwLhvP1V7f3TeStaL30N7kUsbaM59W1fWnbNYgvO5rPHZiuKCHBD5H7XeN5STrGGsCLBx2fhyp0flgu9T5C79JQOqQP2TG0f/2tosMvPPq/Sv7S4uOrv1i/PtfhDOefBRQ7j9DpuZxthkvfkLp2xUuk3NegR2cdwxuty3aPw8yEkjRMap/iL6gmcmsd/xJbKRjZrXwi9Tp1Qh1/CCliK3pPyV0bjVGTuvxgbye2kRE6LRPHTqVetMs79eNhr647+ZFNtXvzaF27LuhEj6H5XgLis6ZU2Un2eBxV3v1gjQeTYkHEzj+QDfH4PFsno3g2YTWP2ErRvnT89chYPw7k/WRMrMAn1S5prUtG6tobaTKWy603lmGjGUaKBmPtxryG6eidv7/V+YiOcVdsLxtACTnukHZddRA6HVLWpCvZmtbbyOXaev7nlGBlyOQ6Ga0DIgXDD8CvC3p8Eu1OSiun2JH43Av89eKxIGJPQTm059gyrWMSdEN2taaMrAZwDNNRe0ru9LZjIb+41s/rLo2mUwvlm2iv5ZdXYt1JeKTN5TwhEVFLwXMpr2rb7Oq166XXgNVrFJi+poTQc6MvyhLGrlK3rim+PCV5S+pyOf8YzcBhjT3y6+JxAF8H4A9CmRGJ4HOY8kbugzdkQI4lEDvxIQC8A06itE9Geql8uvalg7Id0P0mkTg9yxdjCD4VPAMDUzBWtK3dRTSpdx3Inf9HqXIX22qfn4zS5TE00ctjpM6hdK5RxcqlPUdLUi99yxlp+z4RvNxHza3LFAx/PI0mSv/aw2gmybA4jkUUuWSJxA4AVM5HDSVG70GCR2kpKZdIvQNHke+clGuz7Sv4LlKXWGXESalbb2hJXf4E6NKiJ5E20lJrWhpGnnPXz5uOQwlwem6l9LiVOwid1ltSt6L0kvtT1zZvvo+UelG0Lq+H02ii9NMATp1CJfVvYjvxX+ayKo44biw/2Iwi9hDCOwD8cwAPAPjRGOM/HuO4s+MKqouAfqopnUm4OeVzwK5BHCB1LfWiLSOsSj1OqeBzQqdz4RxVlhd1QurypjnkryuJlsu6pTy39tcY427Jj3UXTVWLVa7SQ+iJw2SrYHL/FanKF215idR5eeNUtC4f6Tk/2GRYEKA9cQbQTsHweUwXM3ofLPYQwgMAfgTAN6MaXuIXQwg/G2P81aHHng2PQf+5xuTO4ZG6jNqtL6z1reggBB6108WeEnzuiygFn/vZnJK6lVPvJPXUm/IvbZ+cUqrwnxtJPpdw+3TBStcRJPBXMS31PRL6q09s4jObW7h/dB1H7uzg9PYlPFKXv5bc30r+S+R/X2+p5y7WW6irmq6jqnKzZL0cUgfGidi/EcCvxxj/DwCEEH4KwJ8FMGdiP4Lmg+XzG15h29RpmVTUzi8yLhDt4ispLUCZmIG04EuPUZIHLZF6kdBzlCZmR2iXMN879Z/WR+gl5/oampw6r3TpmUPn60ty6Z99YhO/8syz+Erd8/j+sUfxG888i9cBfHVB34ZZSZ22zebVtRO4twusXIE9BpQ2OfXiMobYT6Fp9gGqy/aPyI1CCBcAXACAM2fOjPC2XZBSp9EcCZZ/06Y7s0TPkcvl9oXw1IuVa7cE3+PtpvbX3oeQUi+qTx8SrXNKhFnaeierZaxzKUWrhMoNnG9J3RA6P2TJkLqpxtHPbW5NpE585dAR3NjcSoo99d+i/bjRpE7LNKlPTXeXkznQtE2sAVVgRnn1s2yj5ZI6MI7Yg7Jsaui0GONFABcB4Pz583s0tBpdvJR64R8qT8XcbL9e2cFksg0NmZIBplMH8mbQM+pMNaRqgi/Ju3Os78ueSN0K6yy0/8PVxPpS0WvHkvtan59sNJe9RK3tDamXCl2+dYnQadl9YxTPN5TlQ6Jz/lxKnV9PrdQLbaxJnS+TF+ItACtA8wucJH6WPS4PY4j9NVTt0sRjAL4wwnEH8jQAmmczNTwvbzlnnZVSaF/aVHljx3xxqsG0dPuuQac8HpGqUe+F9pucm6AkLLS2s0Su5Y2sY2mhp1UVpUldsynB1xlSn6XQiUN3dvCGMrrnoR6je8r/cu1j5dvI6pciqcuDHkO7E9Jp1LXrZ1HVPh6vH7UofXF6l6YYQ+y/COBrQgi/F9U06t8F4HtGOG5PnkP1QVLZE+9tStwUz9kAYBraN0V+kQFd3FLuYK8Lw+qS8seuNwPrGPJ9tXVFN42UeGUUlkOTrVxHWHefrkLn22mfkxWpa9vL64ek/kqzrHSMdNq9VOh8Oa07sX0J1555FpGlY8Ib93GicHRPIF3GmLtPTqVfclLnQpclj2uoe5eeQCVzPnEOj9KXQ+jEYLHHGN8MIXwfqjEwHwDwQozxSma3GfEi2nWsQPXh8omrCZI7GysG0AeD4s/pS9llnjoe7dGVXr9+6OH0fKdEqdyBtOBLo+5Rpa5F6fx16uAl6Rrtjyo5Vu74VCWVugmn0i/SqJnUy9hC104BaBpIb2xu4Y2j6zh0Zwcnti8VNZwCZVK3PiZV6vKG31nqJPPljtI5o9Sxxxg/CuCjYxyrHzRwvpQ6QRH5OiqRJ6QukRF5qidHCrp4+f4JYQyJwIekTGRP0tyxshUwpVJP7Sv3KZ0VO3XjSEl+Mu7L+vTgYDI/zpcDbQsD05G9kXrJTXrRN0KHsc1XX90uFjnR9ccOX56Ueheh0+Mp3suUS315o3TOEvU8PYlK4ifRpF8KZkqSuVVLxl0krxy3JDND0p3FPKj8+EQuUp+Set9IPbXvXtH3/bRfbqlcer28dIz0oULXtivFakueqdRPswNYUm9hVbksr9SBpRA75dhfQpOOoQZTmmm+AN4AuorqWyilD5R3TmKvKR0j5W5F7X0En5w7UiE7PACMDQC9LEK2MUiR3oV+99C4LU7Q2le+D9+O1t1l64hVitCvNetuoSlRpHOg4/D9tYg+kXrJlS5aPxDltpIhZa9AukBILpft39y/ptQ3oEu7ROinFnu+0jFYcLFTCuYG2imYjpDAuRSkdTXJA9PfRrmNsHmHNlVV1ncS60rRpN5b6Np2+4kmdwnvlcElnRsnXW7PhA40UXpq5EXt0LllnDGFDtjZLE3qMkoHMlKnPLkUOuBSH8iCip2idEq3cKlrQ/auIzvzvJS7pLTlyoA3onaRuyQldJma1n4NFNWplwpd25awpJqC7yOj9pJ9CPrPXYVtR25SWf0kjyOj9I5pF0voqR+D2qX1209s9m4QBcaVulqnTiKnUkV6fUwcTD4HXOgdWQCxyxnHKdVCjaKpKfJoRpVa6iuon+/qJYlS7qVpl1SNO5r30eSecleucTXnvVxlINBhcoxSoVuUbm/JPZXKsW4ich8ZnctciHZcuV+PtEtfofN1v/3EJl575/cDDz4EAHjj2KN47Z3fjy8BWLm6nb2HDpW6lnqh7UypPy4OrEXrALDiUu/KARe7lLoWjX8oc4ynUbWos1HhVm5Ukr/HcqyEFDIw/W2TaRl+g+DrDLkDes5dMkqHIYGZesk1jGrbaNtbb1iyrbZOyt06Lv88ZN6d1kuhWwlsvu8IaZcSocv1ct0XnrowkfqEBx/Cbz91AStXt817Wx+h89epfDqQidTlgeXzyYQZ5wBcVs7esTigYpdC13hP4bE+hmoCDiqTuolmHkw0Y8fwb+AapmUvv2l0kQpxT9ZZcgda0fuQ9HRpOke+Ryep96lmsf6oEvNoaGkZ649PpV9SUbpmWaNx1CphtKL0lLBLfvzdAvCVt+pJuGgsB9IfXa7yRaZegILKF/6a2z8pdep4RDOg5QI1BziQYi+R+rs6HpN63D2HSvDH2bo6goeI3mU6Rv6kJ5SofGq5EeFzwZZ0YpJoTkxlLHqVL44t8hSp7Ut7ZVnCt3oRW6bNpF2A9uTRfFctSu/aKJr7MSHJZcqkvLVlVpQOFJYz0nOr6mWSeuFSX0fT4WgdVSfErt/v5eMAib1E6BtoJrHuw4dQTX79JFtWR/A8eudVE6uwB2rRJJ6SO6BG70CHoXAV+E0hG51bG+akXhJZdxX5WNU0JT29NMNaZSmFQueba4eQb1ManWu85fU7+MrD08m5t7yeLorVPsKU0Pny5GBeUuo8YoeybhKlA9PDa59ky86iccFy16qnOEBi32t29OfaN0t+n3ISt5bTa3kMwI40C0jeFLpE01pYlzqGdZw+pI7TNXy10AxrCB0YL+2SSsXI9RZrn7iI//fO9wIPHmoWvvkG1j5xcWpbKWxtWUrq2uQYxy2h0068AoYP4KWmXoBmMK/lHaFxCAdI7HR31hpR6AIYcienHN5JNOWR1zHp1MTz7PwRSFdjlMhdQ0vPILNPCX1y4nw7ywDW9n3RcgKzQpY3AtO2NaTep3E0V58u15WwUpc13t7cwu8eXccDd3ZwbPvSVEVM6uNbU9ZrUTowYHgAdVRGQJ8vwRpH3SP1HAdI7MQ5AC+IZSdQibiv4LfqfW+iLfQbVWWM9u2U9JE72DLt2NryPoJPSbKLNEtTI0NEnLrxlIxy1oWU0KV9O6ZetCh9LKFr8+ICldxXRN16l6i8l9RPwBY6YEt9FUounWY1exLNFHY+9ksfDqDYgabihaLsWsKTn3F0d7+P9MXwNNuehu2lIQZ2gGu76QhOUip3vr1saJXvZS2X76UZwhK69m230BqHLXF3OW4J1vuMMUax9svLasE0pJ5KvVj3iT5CH+M+1kXqKaED1YQLU1H6abaTlDvdBfi6U1LqHqWPyQEVO0GlT7xckeByB5qOS/R4s16+g0bobJIN+aUHyqJk+jbcFq+53NeUZZyc4KGs10RqyVxun4qw+S+MlLjHlrqkpGg/VfLDt5HmtX6NMQtrUk+lXkqFPkuZl0bmsk2Tb1M0NICM1ulAtIMcOmAyfrpH6bPigIud4IPra3f6V4z9eAPpierfCvVYFUP4koSBtOD5b3QOCTIlHkvwQ9COYUldOze6EXWRdYlgS+A3FPo/5dG6dtPVfhXxbTTzWkKv13fNp/eVeh+Za/e7LkLnr/nHLKP0otQLT7sAldDpn5l28Sh9Fiyg2Al+UWygPclGCiZ43jkJmG4szSVDJTK/bsmvNH8+ND8uv+Gl+w05jxJk+kcTfCqnoR2L9klJXRE6bWaNl95H6tpl01Xq1g+YsaReXMbI8+Z8BxJ6a9x0KXOgKWF0oY/JAok9d2FQp4ac4NfZo+icRHDRyKhRiyZl46k8lhUpS7RIW7NELmK2onRL0GNV5GjnIM/fqiNfhR6hl6CJ3DJtYdUL341nc7T1KanPKkrXluXy6cm0C+0g0yo8SueyPw0xu5GUOX1fPUqfBQsidrowSi6K96HKyQPNhcaPoSCHFiAoTSG/2V2QMr2rLLO21Zbl0j0lUs/l64k+f3Mqus+luLqg3TCscJptkxqJUds1l663pK4JvSSDJaWufWxDUy+9ZzmaitL5xNK8sOE4e0640MdkQcTe9aKgIQS2xPLjckOGGNaXvh08bZCqcCmJeukLRHJO5YsJLS8vhSVl2if/rf0ySP1NlqVuKedTcrxSclF9wswlUbqxq5rRKZW6Fu13/YiktCFel6RekkLnO0mp88G9WqkXGY1buNTHZkHE3peUyCVi/kstD8y/jfTt7Sv4UrTjpfYvMUbXxk/tJsTPSR4r1zYxi9SPZuERpZ4rbwTSQpeUfASpJpJVZb2UeucoHWLZY2jn2afy6Tlc6LNiycVecvEJtMmNtZJG2ZlmVoLnx0sdK2cJGUn3kTvtJyn55cFJ/WrpgxZO95S6la4fS+gpeBomlZXj/70l+XS1Ll3eGVL5dFXqqaDpJqqOhs6sWHKxr+c3SU2np1VtSCFpguff+q5lhBJ+vJL0DTdCLnLuU7ZYWhLKtyXWjOd8n1vGeus9NPtmpC4PQbvIQ+Wi9LGELinNq6ekbkbppzEtdP5c1qZPShlJ6hwtcHKh7wWDxB5C+PMA/j6AJwB8Y4zx02Oc1HyhRB4rO/VEHHXkTkLlUtcidOubnRMsYeXYpRAt+fH9U420qffj7yPRblh98g0lPUutXL3WxqCVPMKWuja1adfUS1epp5oeAL0SZlSpb8BOu9BzNfXCa9OB6Wid1t9A+RwJzlCGRuwvA/gOAP9qhHOZU1LVMkzw9CWQwpQyHZpikNKyEqyS0gi3q8y19VLwXf/W0hpA64ZoGZhtn5O6EeB3lnrpn566B2spmNEj9VQunRpIk0IH2t8VPtTuSfgY6nvLILHHGK8CQAhhnLOZa7ScYR2J0FypPIoHphtUKULlZZKSnNRK8u4cqyeL7MHZ52Zg0ac9gbbXUkVWeaJGLuGtnIM15gs9DonU+2DdW2cmdfk5c6nzKH1Smy6FzqFKGC71LkUKzhjsWY49hHABwAUAOHPmzF697UicRDW+jHaB0mTYXPLARPSce7vttE2KoTNVp0RdmvoZSh/Bk9ylOTm5/fk2yt/Ko3Wr45GWmh8j/ZJLuciPTX7MM5G6FrWT1NUepID+S9aSeo8iBWcQWbGHED4B/ZP5QIzxI6VvFGO8COAiAJw/f34Oph2XQ/92gf93XIcevVCHjBvNIqqo4WkboC37PpUoklJzkED5uWjzh5Zi3Yzkr5OSCh1rmWVLKOtvGc8VhkrdOp52mjm5E1YVzMyl/jjEOC8UpfOZxSQu9XkiK/YY41N7cSJ7y4toGnRKuN7zfRLHl99AnqbhaA2MfVMkdDztWJrkx6LvLwRZmy/PTx43J/W7erQ+htRL8ur0J1gfX0le3ZL6mnhNx8hWvxRF6Tmk1OUQH85es2TljpcLt+O9TEnqibLH0mPe29VXnTrcpGk+r6znJkhFu0OFnLth9Dm+JfVcWsgq7eTrrXYKKXWlCmYsqXetfhkidbltL6lvoGOUzrmOaWEfR3vIXRmleyek/WBoueO3A/hhAI8AeDGE8JkY49OZ3faJV9DMkpSCpM6j9D5S59G6IvXJMKZAazTJtd1pucv0CV/GTSOlXypibodUVFzaCEpo0qX3ocZbS84p6Bj8PLTzraVujaWek/pnn9jE5za3cP/oOg7d2cGJ7Uv46qvbg6pfJKmcei7lokm98/jpp9EzSgeanLsVpbvU94uhVTE/A+BnRjqXGXE/vwmA2UTp9XG51NVxqYHWaJKnr7Xlqn3DwZ5bMpYVORzNFrScR8NaY29pI6iES50/yu2t4+be95Z4VFIvQLnUf+WZZ/GVQ5Wc3jj2KK498yy+hGoKunmVujnc7gbaUj8NI0rngubBDU1I4+OoHwSWLBWjIX9u0sXbJ69OEfiOvroVpSdYOdWMJsm/wZOJfxlUQ88lrEXBXfL0Wp6db2/dOFKNrqUNsjmpS+TfmojSgfL0y+c2tyZSJ+KhI7i9uYU363lFc/eX3H8xp4/UrUZSADiuNZLKvDqvTV85jPZwukD7e0DXrVUdBrjQ54clELuMvq+j3QvuRVQXMF3EJPp1VP2vtKg9J30hb6pzby2XDas7Ytl6M1zw2EjjyAZaLUKXDZkcLQ1UmraRSeqSctBbmDYxt7Tx1l1y6veP6jfg3zWWS/pIXWtKSUm9V+qFN5KustdFw2t8rGAbZx5YArFzoQPTXZupR9xlNLnFm2gilx2Uz77EKfmi3EAzCbdVQSOPI0S/cthulJWkbJMTcKpaJnVciqi75iu0SJwgExd0QNLGVNcOxXe9BeDwnR3sHnt06rTCneb/f4wCIin1XJ26lDqP0lchpq7jsxvRoxzrZQ3sl2DqmnWpHySWQOyW0CXnUDWw8uj97QBeQiV8krsVrcuGI/4TlnVgUtGWy959VuReL081QqbEq2Edp7TGXr6fjMJ5zXxXLKmzc6Y0jLYrIe818nCPbF/CtWeeReTpmDe/jHjoMO6872cR7uzg8PYloE7LcEqj9RKpy+ea1Iuj9GPKgSbIlCThQj+ILIHYSZolAxBt1I/30Yhayp0LNtULjzdAyUiIzkk7Fk8F5WqB6+PwqJ3n2TmalHOdfbqg5RHkeh7idg11UzWIyjnL3DrYpkpKvnWYuwC+uhb2b25u4XePriO8fhfxoYeBhysdx2OP4v4zzwIAgiJ3DV6KSOTy6bQuK3WeP5e5dC50+TjV7kP58w8V/U3OfLIEYu8zotwR6HKnml2tntdC206L4q2ICWiifxm1n2iWrxwGsNtuRNXKJFOkasNz8BZD2aOW1vNoO5Ub0c7LkrrYjk9tR1Bu3UrBWISr23hbLe3XvveFidQnHDoCbG6pUTsnFaXT+pzge0n9cUyLnDgVUV3n2kTTzkFnCcTeF9nCTxJ9ucexjrPH0hJKKXH5WsmH0pAFJHf+RdYqarQJQ1Jy58K20PLifYWeWq68hzVfqdUhlR9e3jck0Wo0zTSmlkbp8nWqnHEidZ5P7yL1r6URPe6jSkGeRLsm3TnouNizbKDdmAq05Z4SNQm9JArScuhWXn1deX4DU9P3AemGMaq4oZJJQpO7lpDOwW15WyzLHavLr4aE1LXDlZw+3/42qkbTqDSm4o5dtSSlXlL1QutkGpxLPVvK+LhYRnytNkTT5frxefhojIuDi72Ic/Xji/Xjk2iky2t9b6L9peCdjywoiufHyQndGoSMy13uk4jEWtuzOV37VLQMkbm2be697rYbS7URG7VUvPVcQsc5vH2pyqnzxtQ37gPbl6b2GRKl03MepQNNgYs6iBefzUhWvRCq1DnvQ7syzDnILIjYqXdpaQcJ3hu1S6eKd6EZFZJ3RpLD+nLx8p+3tXgn+4DtJ38V0LH5cejYMheqyL0TdE5ocvWSkpSIrHZJ3RRkrWAq9SPXUQvo3fam1jC8/BCpHx6pP/GhOo++u7mFeHR9UhXzkJFf75t64euKRmbkkTqXOi1XU3AW5/KbOAeCBRE7USL40iEGLKgxlgueixpsOdDuzSd7tGo3BR69a7PT8O7cQBVhydy9Vj5J9fLWNgVyB9L20wq7U/l4vn1Jo63y3nImJImWgunyo4DPR/LQ1W1T5HIfoKyBlK+zUi9qPl0+8vTL1Jgv6wCehpcuLg8LJnbiPqblPlTokvdAj95l6iOXs8xF71SBw4V+QqzXSMk91RmKyV3m3jkpwacG+NKOZ1XCANOhthKtE6mJqIfQZ3j60RpIS8dPV6X+JNq/Iul69blHF50FFTvQjt7HljqhRe9AW7ylJWSp6J2OcRZ6KobeW6vYkRIvHYOekSqVLMm/y/3JlLfZupJqHLHO6oikkbpPjM1MUi/W+Oky/TKR+ibav+z4NfkKmj4bziKywGInZiV1Dhc8T3dYUuepE4qg6QvIomYA+VlpZGOXlLvWezUldyVqBxrpypJI6zWP2jnHxHPKc2hpHK2FU8mtc2QlTGkjqTyFksyQdb9LCZ1eaw2kR1EN4VI8dZ0p9SfRBAJA+1cjXTvar1pnUVgQsecu0L2QO1AJ/sXsVvn9eO/T1AQG9HfxVE9q0DA6zk1UN4BM9C4n5gZmG+pqGDkVWQlj7aK95kiBS7l3oUTo9HpQ6iUr9eNo91zWAgyX+iKzIGLPMct0jIQPKla6LT2nfWSeXZM6fyROQpe1LLvcqY95BUUDnOWknuvYxBkyp6qAC1sbOqAUTe7y+Ba8bbhE6PQ4SOqPiddrQPXZrqM9o5FWpeVCXwaWROxAc0HvleDPocplptZby15BFX3JoQusL+VG/e8y2qkgrSGXeIkt49G7EfF3idRlOkZyO7Euh7DtbfHcajQtSato96sSpMBpmSZ02q5Y6nyERm3Ml9Y0djSAHX/OcakvC0skdmJo9N4lN7nB9pHLcvvxsWpQ+J7n0KR0UvlV4gSaiH29fl7n2CkNw3PaKTtqOfIUfaReIxtNU71Mu0DiLd3fEjpfJhtIgQKpp4bblesnE2RQtM7HfSFc6MvGEood6H+h970h9Hm/vuco69WB6R6xfChjgg9Ktt5MDrKCSvLSdpb9rEZUC5nmKSRVIVlC6j5UmlsvScMUS30D6VLGqSiddzwqSfs5y8SSir0rr9SP8zzy3fOYrnF/CW15c+FTfv06skMY0ExOJHjKpxMkZ6tCphRL8lTekkjD5LBSLR3vJ1PpmS5CBzIjM2qpFy16n0h9r9KKzkHDxZ7lMuZ/UCQpdTpf2RmJnr+M6bHlCWtMmRPVRNtS8EC+0dSK2jWR87QPH+ilRkvDaNFtjDEAAAxnSURBVIfU6Cr3VI69j9DptSp1a2TGqbQL4VJ3bAaJPYTwQQDfCuDLAH4DwF+JMe5VIdyMoVz1PEfpzyFd534T0+mWKygbeliON0Pvs9MI/tquPoRAqq6d0CpttOd8W7ELMUbkXrKftk8XodN6c/o6GsRLLmtNhvFawdk6y87QiP3jAN4fY3wzhPBPALwfwA8OPy2nO1aPVGogvY7p8kYtOtfGp+GvmeBP3QDuXUvOZFSEVb7CnxfmTEqi8pJGUi3PruXU18T6pNC1nqRW6qU1iNdZVDdol7pTxiCxxxh/jr38FIDvHHY688S78pvsOx9ClYaxOI7psWR4esbKvxPWTFFM8Ct1iubetWoVj9KtqD2VgpHPC24W8scAvS2XvFabnmsktfLp8rXWMErLH3qYvZBROsnbrHihMsZ1eAOp04Uxc+zvAfBvRzyeU8T70IxVY7GORtJWQylV0/CBzHJwwQM4rchdI5eTN8iNDaO9pXwrLmvt9EplzreVlS4Ai85pJy2Xbkn9FKVeSOrznA505pGs2EMIn4B+ZX0gxviRepsPAHgTwIcTx7kA4AIAnDlzptfJOhY0Vs0riW1Oosq583Hbh3K8/bhyAvjaK1Xu/fMDDz2wpYYL3bqPWLn1NeO53M+sdFlFO1fDq1q0KJ0/TkmdBn5znHKyYo8xPpVaH0LYAvAtAL4pxmhO0xJjvAjgIgCcP38+N52L04sNTFdLkBRoGAE5zjtxA/moXevoxLerOzidqqfn+zy6z4jUkaOoKmMoHcOjdp5Ll4KWb5WSubwBZIUuZW1F6GoZI5e5HHvfccoYWhXzDlSNpZsxxg6DqDqzg/espVQJCZl+1l9HM7rkdbZOg+fZqaMTX6akdvjgYT3TLi2MlA7JXJM7301rLC1JsdB7EEfF8qIInT/KMV7UXLo2sfQYv66cZWJojv1fADgM4OMhBAD4VIzxeweflTMycgYnPjQwH/sd0Cfu0JbJ3qxUF18f1yqDHBlN7kBa8BqazKXI+badhG5Wu2gpl3UAb2fv5tPVOd0ZWhXz+8c6EWfW8OidBC/nYKVJPrQRH3fYI0X3MnK/og8/sEeQ3AFd8Cm0yJxXtwCiwgXQ0y5gz/nMRlM16SfQDLErJ1EBqhunS93ph/c8XTqk4IFpyQOVWChS38F0D9YTbB2QTBfwdEzfUR0T8EOSlKXgc3SWOS3ThM6j9KlKFy5wGl6XT4oBVOm0jcIzd5xpXOxLC8+b03OSOwmez/RE67XnQGtEyFK0nkOylIXlTx5CVfKo1aRzgVNqhnNHWcbpLXS5jE+AMYnSKYfO4Q3QPhKjMy4u9oVkjLHnSfBAJXnquSobSwsa9nhqxorWueRlByfW0emhehNL8IAdpaei96TQ6VHLodO6bD064J2MnL3Cxb7QjDVzFEme8vBc8DwSvTadX9d6mXbBaPHUBE9ob8U7D6nHywkdmO5cxA9s9holXOrO3uFiX3j6Ru981iXeoPcSdMGjahi8tduMykiGpUcrWtdSMFaahq9HW/CEJvCJuC1So3rJKP1x2LKfjO1CuNCdvcfFvjSURu98PHdtULG3o13+eBPJ0SJzwwqkUjApuaO9TVbcdHzrXPij1hia6mA0Ve1CbRQudWd/cLEvFSVyp5w5DUFAz4nrbBteOQO98ZQnwmULp4i+J0i50zIg3Z30rnIsiZba6Sp0iPUApksYP5Q5EceZHS52R4EPLyB7qd5AFaGT1KkUUjSqcoFKeUvJ8+VW338teqf1KVnnlssUjBQ2n/RC7kP7TeYdfRJNr1HH2T9c7E4CObwvVcdIqQssCQPTUTgVoWtCX1Vea8M30vrUOAHy/OR6rYRRi9L5sVfBcuokdR+0y9l/XOxOARSNS6mz9VoaRkuVpPr1A3pqRltnCV4ij6VF3jmhy+PQ9qc0qbvYnf3Hxe6g3RhaOvAU76W6Xk2VBwBQBG+lUzRkbl1rXOXb5Y5FaFG6XG5F6NpxVKlTr17H2V9c7EsH79m4Ab3r+mW0hxuQk2InWAGqeVFRRfE8VdKXVUxXzADpG4WWO+dwWfMid9nxSDuWmX6Z90nPnWXBxe4onEMl9+NohhWQI0NqIz4KSPKnd9sStmTM8+3ajUCL4i209IqF1pMU4vVUOSON78Jr1h1nPnCxOwY0suB9TNe1y6F/JWLdSj3xBoeLOdUoamHVutPxeGoldzyrMgYwhC4nwPBo3ZkvXOxOBjlxh+Qs2tG7MXGHlLuVVikZ0IX2typqtE5EOdR0S07oqYG9HGf/cLE7BXQZloAi+gw8Z05oKRitsZPvd1rso3UuWjlVr6x/RWgVPCuH6ycpmXORv8v6yxxn33GxOx0gwVtph5vG8g6kSgy150Ajdy77qflDieusgoc3CJO0NaHzGY0AnwDDmXdc7E5HNupHKXE5vnsmai/JoQO23CWn68dJxQqXOjVwUg9Zqsfn6SItOpe5843MSTvOfOBid3pyDsCLaM+klJH7ymFM8uxWZyVN5K00SQo+fyjNTiQj9itohkbg1T08OtcGP/MJMJyDg4vdGQDlmV/IbHcCkxvAyo0qFZLMcxOykTIndorOSdBaKeJZTFf30HvIKQMBF7pzEHGxOyOiVcdIasmvAE2OW1aXAO1SQqBJiaQqTyiFIqN02fjL530FXOjOouFid2ZIqs79JKp0iYUUuaxK0dAmhYbyWivfdKk7i8MgsYcQ/gGAPwvgKwB+C8BfjjF+YYwTcw4qVmWM7MyTq/nmIpfRdQkpOVu1+S50ZzEYGrF/MMb4dwEghPDXAfw9AN87+KycGUJCG1tiPP99U1meqgm36JKGIUr/Lhm9u9SdxWGQ2GOMd9jLFQBx2Ok4e8eYQnsPqgZUKXESvNZjs0svzdL6+I0OxyRc6M7iMTjHHkL4RwD+EqrO3X9q8Bk5e8x9zEbuHBmll4ytosncErx3GHIcTogxHWSHED4BPbz6QIzxI2y79wM4EmP8IeM4FwBcAIAzZ86ce/XVV3uftDOE1LAAY0WvWvljbpJsCZ+Sj/OSeP2ebqfmOAeYEMLlGOP57HY5sXd4w8cBvBhjTJU6AADOnz8fP/3pT4/yvk5XcuO9jC33IQ2g2nyr1KnIJ4t2lo89EXsI4WtijL9WP38WwGaM8TsL9vsiACtkT40Hu6gs29/sf+9i43/v7Hg8xvhIbqOhYv9pAH8AVbnjqwC+N8Z4rfcBq2N+uuSOtEgs29/sf+9i43/v/jO0KubPjXUijuM4zji8Zb9PwHEcxxmXeRT7xf0+gX1g2f5m/3sXG/9795nRqmIcx3Gc+WAeI3bHcRxnAHMp9hDCPwgh/HII4TMhhJ8LIbxtv89ploQQPhhC+Gz9N/9MCKFk+uUDTQjhz4cQroQQvhJCmKuKgrEIIbwjhPC5EMKvhxD+1n6fz6wJIbwQQvitEMLL+30ue0EI4XQI4b+FEK7W1/L37/c5EXMpdlSDi/2hGOPXA/hPqAYXW2Q+DuDJGOMfAvC/Abx/n89nL3gZwHcA+Pn9PpFZEEJ4AMCPAHgGwNcB+O4Qwtft71nNnH8N4B37fRJ7yJsA/kaM8QkAfxTAX5uXz3guxb5sg4vFGH8uxvhm/fJTAB7bz/PZC2KMV2OMn9vv85gh3wjg12OM/yfG+GUAP4VqiOuFJcb48wD+336fx14RY/zNGOMv1c/vArgK4NT+nlXF3E60scSDi70HwL/d75NwBnMKwOfZ69cA/JF9OhdnxoQQNgB8A4Bf2N8zqdg3secGF4sxfgDAB+rBxb4PgDq42EGhZDC1EMIHUP28+/BentusKB1AbkEJyrKF/uW5rIQQvgrATwN4r8g27Bv7JvYY41OFm/4bAC/igIs99/eGELYAfAuAb4oLUoPa4TNeRF4DcJq9fgyAzy62YIQQDqGS+odjjP9hv8+HmMscewjha9jLbwPw2f06l70ghPAOAD8I4NtijF/a7/NxRuEXAXxNCOH3hhAeAvBdAH52n8/JGZEQQgDwYwCuxhjnarjRueygNIvBxeaZEMKvAziMZozaT8UYF3qKwRDCtwP4YQCPALgF4DMxxqf396zGJYTwTgD/DMADAF6IMf6jfT6lmRJC+EkAfxLVaIc3APxQjPHH9vWkZkgI4U8A+B8AfgWVqwDgb8cYP7p/Z1Uxl2J3HMdx+jOXqRjHcRynPy52x3GcBcPF7jiOs2C42B3HcRYMF7vjOM6C4WJ3HMdZMFzsjuM4C4aL3XEcZ8H4/xcEBozt6OyVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cl = coor.cluster_kmeans(data=Y,k=10)\n",
    "\n",
    "#cl=coor.cluster_uniform_time(data=Y, )\n",
    "# for later use we save the discretetrajectories (dtrajs) and cluster center coordinates:\n",
    "dtrajs = cl.dtrajs\n",
    "cc_x = cl.clustercenters[:,0]\n",
    "cc_y = cl.clustercenters[:,1]\n",
    " \n",
    "#plotting out the results and clustercenters \n",
    "contourf(F.T, 50, cmap=plt.cm.hot, extent=extent)\n",
    "plot(cc_x,cc_y, linewidth=0, marker='o')\n",
    " \n",
    "print(dtrajs)\n",
    "print(np.size(dtrajs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.477317    0.46967688  1.3487787  -1.0906587  -0.78296345  0.654191\n",
      " -0.26333734 -0.20036615  1.33917     0.4688478 ]\n",
      "[ 0.4292571  -0.3879486   0.6676682   0.5185451   1.3539321  -0.03853653\n",
      " -1.2127491  -1.386519    0.68939865 -0.40003198]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvX18HFd97//+ZtfWONo4GyQnIjZBYIU4RGkEMRRDL3GJCgHSXugDuG0MhuZS2hvAt4j7C/2Fx6a9oT+nr9CGQnMDV8QCzFPpbQM1kLQOBcyDVZSyIQ6xiUrsRMRyotjreG2vcn5/nHNmzpw5s1rHcmzJ89FrX6udmZ05Mzvzme98vk+ilKJAgQIFCswfnHaiB1CgQIECBWYXBbEXKFCgwDxDQewFChQoMM9QEHuBAgUKzDMUxF6gQIEC8wwFsRcoUKDAPENB7AWOG0SkV0SUiJRP9FjmGkTk90XkGy3mbxGRq5/OMRWYOyiIfRYgIr8nIttEpC4iD4vIP4vIr7T5XSUifd60PxWRB8z6donI53O+W3deT4rIQefz74vIB0VkxFleROSdIlITkQNm3V8UkYu99X7QjOvFbYz/eWYdkyLyuIj8h4j8iYiU2tn/dnC8SExEXiUi3xKR/SKyR0TuEpHfmIX1rhORbx/LOpRSn1FKvfJYx1Lg1ERB7McIEfkT4CbgL4BzgPOAvwX+61Nc35uBtcCgUqoCrATuDC2rlKrYF/Bz4NedaZ8JfOWjwLuAdwLPAJ4H/APwWmf7Yrb/KPDmGca6HPg+8CBwsVLqTOB3zJjPaH+vjy9CNxkR+W3gi8BtwDL0b/d+4NdP1JgKFJg1KKWK11N8AWcCdeB3WizzYmArMAU8DNwMLDTzvgUo4IBZzxvN/JuewljG0TcDd9oHgRHz//nANPDiGdbzcuAgcBWw1441Z9kR4Kst5vea/SuHxuiNLzLr22uO1Q/RZPvnZtwNc4xuNsuvAL6JvgHdB7zBWe8w8HHga+bY+sdF0DfC98xwLN4K3As8BnwdeLYzTwFvB+438z9m1nuhGeu0Ge9U3pjM+XMbsAf4T+A64DSz/Drg2872fg3YDjxuzpG7gKvNvD7z+XFgEvj8ib42iteJfRUW+7FhFZqQvtJimWngfwDdZvnLgT8GUEq93CxzidJW9ueB7wFvEpH3iMjKWbTsLgd2KaV+MMNybwb+CbDyz5Utlh0EvjQLY7PbPRN4FtCFJs2DSqn/F/g34BpzjK4RkU40qX8WOBv4XeBvReQiZ32/h74pnAH4ssgFZju5YxeR1wF/CvwmsMSM4XPeYlcCLwIuAd4AvEopda8Z+1Yz3mqLMf2N2efnApcBbwLeEhhLN/BlNPF3AzuBlzmL/BnwDeAs9NPH3+TtV4FTAwWxHxu6gEmlVDNvAaXUqFLqe0qpplJqHPg79EWct/wI8A7gVWgr7BERuXaWxvpwqwVE5HS0lPJZpdQRNPG1kmNmXOdR4IhZX59Satoct305y14JjCul/o85rv+OJr7fdpb5v0qp7yilnlRKNQLjZoax/yHwv5RS95rf9y+AARF5trPMDUqpKaXUz4F/BQZm2Md4TGZ/3wi8Vym135wbN6JlMB+vAX6ilPqS+V1uAiac+UeAZwPnKqUaSqlj0vcLzH0UxH5s2At0t4r6MM7F20VkQkT2oQmiu9VKlXacDQJVtPX3YRF51SyM9ZkzLPN6oImWCwA+A7xaRJYcwzrbxUa03LFJRB4Skb8UkQU5yz4b+GURmbIv4PeBHmeZB1tsa695bzX2ZwMfddb/KFpqWeos45LrE0Clxfr8MXUDC9ESjMV/euu3ONf9rlJKeev6n2ZsPxCRe0TkrTOMo8A8R0Hsx4ataD31dS2W+ThaGz1fKbUY/Xgv7axcKXVEKfVF4D+A/mMc653AMhFZ2WKZN6PJ6eciMoF2Li5ASx0h3AH81lGM4QBwuvM5JmKzrx9SSj0feCnaKn+Tne2t50HgLqVU1XlVlFJ/5CzTqmzpfWYdrcb+IPCH3jYWKaW+23oXW27bnT5JYmlbnAfsDnzvYbR0BMQO7vizUmpCKfXflFLnop80/taPtCpwaqEg9mOAUupxdCTFx0TkdSJyuogsEJFXi8hfmsXOAPYBdRFZAfyRt5pfoDVWIA6Ve62InCEip4nIq4GL0NEnxzLW+9HROp8TkdUislBEIhFZIyLXishStA5/JVpSGEBrxx8hX475APBSEfn/RKTHjL9PREZEpBpYfgxYY47RShzpRER+VUQuNj6FfWjSmzazU8cIuB14noisNetaICIvEpEL2zwWCvgT4H0i8hYRWWyO9a+IyC1msU8A77W6vYicKSK/0876zXiXicjCFmOYBr4A/Ln5rZ9txjQSWPyrwEUi8pvm6fCdODdFEfkdEVlmPj6GvoFMZ1dT4FRBQezHCKXUX6EvyOvQ0Q0PAtegwwgBhtBOs/3A/yZxSlp8EPi0eeR/A5rU/hQdtTEF/CXwR7Okm74THVHxMbPunWj55Z/Q2u6YUuobxgKcUEpNAH8N/JKIZJ4YlFI70Q7hXuAeEXkcrXVvM/vr433AcjT5fAjt/LToQWv6+9CRKHeRkNxHgd8WkcdE5K+VUvuBVwJrgIfQkshHgI52D4RS6ktojfutZh2/AK4H/q+Z/xWzzk1GQqsBr25z9f8C3ANMiMhki+XegX6K+RnamfpZ4FOBsU6ifR83oGWk84HvOIu8CPi+iNSBfwTepZR6oM2xFpiHEG28FChQoECB+YLCYi9QoECBeYaC2AsUKFBgnqEg9gIFChSYZyiIvUCBAgXmGU5IOdXu7m7V29t7IjZdoECBOYbR0dFJpVReklxbeOmLXqSm9uUlMie496c//bpS6opj2dbJgBNC7L29vWzbtu1EbLpAgQJzDCLynzMv1RpT+/ax8W//dsblVg4OtswKnysopJgCBQoUmGcoOtsUKFBg/uOJJ2Bs7ESP4mlDYbEXKFCgwDxDQewFChQoMM9QSDEFChSY/yikmAIFChQoMJdRWOwFChSY/ygs9gIFChQoMJdREHuBAgUKHAVEpCoiXxKR7SJyr4isylnuRSIyLSJuQ5m/NO0L7xWRvzbdsGYdhRRToECB+Y9GA2q12VrbR4HNSqnfNl2yTvcXMJ3APoLu42unvRR4GfBLZtK30Y3tt8zWwCwKi71AgQIF2oSILAZeDnwSQCl1WCk1FVj0HehuYo840xQQoZuYd6D7Cf/ieIyzIPYCBQrMewhanpjpBXSLyDbn9TZvVc9Ft8D8PyLyIxG5VUQ6U9vS/YNfj+6bG0MptRX4V3Rz8oeBryul7p3N/bQoiL1AgQIFEkwqpVY6r1u8+WXghcDHlVIvQPesvdZb5ibg/zENy2OISB9wIbAMWAq8QkRefjx2otDYCxQoUKB97AJ2KaW+bz5/iSyxr0Q3QQfoBl4jIk10E/LvKaXqACLyz8BLgG/N9iALi71AgQKnBNqUYlpCKTUBPCgiF5hJlwM/8ZZ5jlKqVynViyb+P1ZK/QPwc+AyESmLyAK04/S4SDGFxV5gzmL0zjvj/y+9/PITOJICpxjeAXzGRMT8DHiLiLwdQCn1iRbf+xLwCuDHaEfqZqXUPx2PARbEXqBAgQJHAaXUGFpucREkdKXUOuf/aeAPj9/IEhTEXmDOo7DWC8wEGxVzqqDQ2AvMaRSkXqBAFqfSTazAPENB6gXaRWGxFyjwNGP0zjtTjtACBQocGwpiL3DSoCD3AgVmBwWxn6I4WUn0qYyrsPgLzITT0EVaZnrNFxTEfgpjdO1a7jk+VUPb274l5FlqgFCQe4ECGsdM7CISicgPRORuU2f4Q7MxsAJPH+4RYfTiizXR3nhjat5uEXYHyH+nCHeLMHrjjWFCPbv1DSP+zvAwlaEhWL/+KY8/uN4CBRwIhcV+tDgEvEIpdQkwAFwhIi+ZhfUWOA5oR7YYNUTuEvoB5/+9AaJ313mfCAf2tB4DAIOD9IyMUAHKtVpsuR8NORdEXqBAFsccAaSUUkDdfFxgXupY11vg6DF6550tQwCPhgR3iqQsmAYQGUKvOtOaY2MwMJCsf2ys5UllxxANDtJjpk0BzcByoX0piLxAgZkxK6GdplvIKNAHfMypfOYu8zbgbQDnnXfebGy2QACtyN2d7ksuLkInRQMYN/9bYvfJGIhJfgp9IygpfY8fFYENGwCoDA2xzKyz4Wyz6Wvtof1oR48v4tsLeLDO01MFotTsGdciUgW+ArxDKZXbh2rlypVq27Zts7bdAhpHWxTrbpGYxBv9/XDTTTA8TDQyQi/JhTCBJurGhg3aIjfySR24xBK3a0kPDlIBKnbdZtky+qZQNevz4VvulwTOzbtncPaGvlNgbkNERpVSfm2Wo8ILRdS/tbFcBY55WycDZjUqxrSI2gJcMZvrLXD0sFr6aI7zEzQJVgPTy2gyt68Ih3AHBmiaz8uc77S6kTTMshetgqXBtr8aVXTxavvyx71bJJ4XKrlakHqBPJxq4Y7HbLGLyBLgiFJqSkQWAd8APqKUuj3vO4XFfvwQ1KCHh+N/L924EciSZt28LIk3+/uhVuNSc34cEGEcTaDdJBdBhCbuHUDzjjsAKA8OxlLLFNC46ioqIyNcsAq4AqY/oG8YOOuw6+tUCl6qx7Z3ayLVWFjrv3nVVcH9t/tXYP5gNiz2lSJZfTiA8jyx2GdDY38m8Gmjs58GfKEVqRc4vrj08suz5L5uXfzv6Nq1lEdGAE3QFq5s0rzqKhgYIBoaiud3KsUKEUpLAht1I2DGxugGuj4P02+ESdBWvtkmy9DrMN+pAp1L0N6Z7xoj47sKPiV07YDdezSR2/HZpwUGBmI930W87+vXE5mu9GXggsKaLzBLEJFxYD8wDTT9G4GIvAf4ffOxjG6Ht0Qp9aiRq28F+tFBJm81vVBnd4yzqbG3i8JiP35wiS2Ecq2mibG/n2qtlvv4aTXx5R+CQ78LHT9Cd3K0uALYnP7ObnN6Lv0qHLgYOn8Me1+rrfMVQOnzwPbEYo/QxF5alV3v9FbtrK2jnx7KtVp8I4qtdvOE4CJ1YxsehrExLv3xj3P2ssBcwKxY7AtEbTurjW3tmdliN8S+Uik1OeP6RH4d+B9KqVeYz58G/k0pdatp1HG6kbBnFadSwbN5iaD0sn49lVrWdx1bu8YibwwN0QN0rkrLHhGwdC2wGm1J+7gaDv0X6FhBiuyXfkh/58BzzIRubbl33QTY+8yO9KpKS0g8MpuBm+HAHiPhYCShdetoDg8zVavRA/SYebsGB6lfdVX8RGJ1/pjc163LaP/2eBWVIdMYvfFGGBiYlePi3lijkREuOrWfln4X+ByAiCwGXg6sA1BKHQYOH4+NFhb7PMDo2rXJBxOxsixn2Vj3xljRX0VrMjX0A2IfsEZPO3QWNKrQWAjn3GeWqcCBVbDjbOjdB2f+EG1+9wERTC+CpnkMKDeg9BiJoG5DaW6H6Y0OqW82733E1vqUGadrsZfRN6bIDHkCE80DbVvlRTu9MEbXrs3cIJ/yulxjY3j4mP0es2KxV0Rt+6U2trWV/8QoiAa3KKVu8cbzAPAYWkr5O3++s9zp6ObXfUaGGQBuQfdIvQQdIv4updSBo9+jGfajIPaTB22TzoeN4/P9yW9na75EQC9QsoUdPAsZ0AS6DKZfpIm7GWkSLnueyoYJmWks1O/VR/SyE8+AehkqTeh+As58ILuJ6UX6vXSQ9KMA6Mtmmx7bgY2aoHuATvuUcKuZtych+CaJ9s/wMJEh+qPVzv0nnILcNe4RSUJeaX1cjjYR7liP8dNM7O1IMecqpR4SkbOBb6LDu78VWO6NwFVKqV83n1cC3wNeppT6voh8FNinlHrf0e9RaxRSzByDzQjtAUo3Czyiic2GLUaQODhXmJeLHjj0zMQSt6ifniZri8ZCPR2gfq7+P/5s/q9cBD2PQmSUwmaU3CwAWJRMAyifBR0VYEdC3BHQacbH1cAW6NwM7NHzJoH6yIiOo6/V4v09IKIjadpASnsHWLduRpI6VXCRUprch4fj4wJZUh69+OKY/E9VKKUeMu+PiMhXgBcDGWJHP/t+zvm8C9jlJHB+Cbj2eIyxIPaTEO7FtFMkCe8bGICrroqjWtgDBGLUp/dAaTPaWl9NElLSCweeqQm2sRCigLpX8dJJ62Woder/Qxln/SWgAyYi6HlG8p16GXoaen2VfXq6JfpmBDwTOjCJShs2MDk0xNIdJKE6fcAVmtyjPfrG1QCmAr4DzhZ27oEpo8df+u53p2Zb/ZixMaKhoVjSaYyNwU03FeRucJFSaVkPczN0Ioww0lc7OKmOaQdhf5GPGeJTRKQTOE0ptd/8/0rgw4HlzgQuA+K4XKXUhIg8KCIXKKXuAy5HyzKzjoLYTzasXw9GL95trPMGToKQSamfRPN1W0kVFZh+FtR7YPJ0M6kZJvd6OZk2ebom7BrpV795gUP2JYzJ7aBTmyy9kSF5NLk3I6gvhnNWQ7RRE2wEWmcfNzu2gthRUDJyUrQ1KUNgj0cDqBqrHkM+d4ukkpVior/8cvjIEAf26O9N1mrUBwd1qYOTiYROJJzQ2BhHIXudAs7pc4CviDaoysBnlVKbReTtAEqpT5jlXg98I6CfvwP4jImI+RnwluMxyEJjP4ngWkd+Rqibbt80NVcYGkql7oMXQngFsAKmL9SkPr44kVEsrIVeL8OWEqyeTiztiUhPnyxpAt9E2mpfk7Mflvwt7I2g95Bet4tzvgNcTxIZY58ybIUwx9k6BnG9mRjeMWiYY9ADcZ2aEPaKMGmWLzJWNTIRVuYJpwfoOgpSh9kl9lnR2JeK2vaHbWzrA0WCUoFZhBuV4MLqz1WS7NAYJjs0F8uA3oTUb+/Iyin9Jf2+Cag1YE0E10RZSQY0Ocffb8CmKCHtvFHUnHn9HbCmI03wv3gZVD6vY97ZhCb2PjRT95j31TC10exvICnJIpXsBNrJ/P4wIe3KXcupiQypm1BFG1paYG6hIPaTAPeZYlxNE6rYTSI5NK66ioa5wLrRXD1lMkJTRGZfy8CasIfOh6mztdU93pG1uGPYGEhgUxVqEawxhO9KMPEV7rzX0Mv7qAH9znRXuql1wOoy9Jm7VDPSY+1Yg75z9cD0WSZUEmBXgFyGh6mOjISlKBNRw1bofH929ujatTqW3/lclCJIcOnllzM6OEiTJOR04uKL2w4pnccyzJxBQewnADY0sYzRnUlqtFSArlVwYKuxKo3Tj1qNrlXA1UbK9jUYE3fuRp9Mnp44MmNS93Pc3Nq5Zn6tCtdFzjyfVf1C7VMExf5aVU8Pudu2lPT4reVejoDztQZvQzCjyGS8ehmurF9PT63G0iWwd08iv1jYMElIy/5xREdIRz6FMHrjjTojN3Azi5O8lOI+kSSpzYmUcZfL+368rZNFc19AfnLHPEShsT+NmBZhB+hsSUjXFjcSgy2J65bMjXBS8ntIKm9NANvRJ2yfzvisL9bfs4RutfPrIMlMsgiRNiRVuXwLPQQ3Rt0nfFuj15m+xnnvnjZRM+blOnKjKei43+zjrXra3q3aadxnjwWkk53sU8vmNLnHu+5FdJyqpQZGL74Y0OUlIoy858Wwj158cVzGoe4uY9B2QpgTHvlUyX1WNPZeUdsCT2+Zbf1BobEXQEdgNJ209hh+xt1LhQnMBdLCYmyOjVE39VwqzqtkotCmTb2Lki2YshlNZj0mnHBxQujjHYnTM9f6ttOq3rRQ9YrQdy1cAm8Q3p7R5dcErHs75u4nDKk/bPbPEcO7VkHXFcBqI9UcBK6EktXUl5nlr0hi4O0QypAqs5AX4TGThWm7RPkhlcFlb7yxreXawVMJybTNVPwx+OUmqkC1VmNicJDGhg2Mgi4CZ0o4TGEaoTjfu0+krSiZHrNe7rgjaVw+NBRXDX3aYJ4QTxUUxH6MaG7YkDj0nPK4jI2xW4SlSsGbhANbPcens1z8vYEBXaSrv596rRaTelzD/CYoXaf16MZzIDoLSsbBeOA5OiN0ItJRLFZ6cZ2dMVlH3jQ3S6gVQpa71UGmAtP85ayOb546+ktJHHz3tJ5XaTpf7XG+b55Kpk2ZA5sMRS9JqQK7XMBid6NmptCWZCurc8bsSpFccoqzOAGOkdifaivAURHtRxgby4zBrerpPpT1ov039lSwpO7+nO6TT8tj6OYWuE8CAIHibQVmFwWxHysMqV96+eVxLPQ9xhlaBXipJvUJzEXkygHDw1RGRpKEGVt50Vj0jaEhekHLLBv1RbV0E/A+uONZ0H0uVJ6XyC02soX9Zv0uuboyjEvs9t0l7RDBtzN/wpufdyOYSMgdTGRNSb9XIv3qsMv3EMtPNrmq3PBKFVjWthb+jmTWUod8dxvN2MoO7aTGt5qfZ7Xap4NjNRDbIfWdIlSBriX6ZjZBUk8Hc2756LKZyY+kj419cHMjYTz3S6rsQDzG4eFYQrTq2zhZUrfvRd/a449CYz9OuMdcKNYpaqMHy5gs0rExqrUay5foTFFrXdooBNBEV7oLfbXepItjlT4PO18Lmzodh2hIN3clFmumurVyrZVtpx+JoNyr1yHb0ztzxGHyBc7GXF3dbi9yXj2ke+HZnXRlG2f5NWafr34cun+uZ1tncDPSIZugi49VJpyoGVu8zJZPuEk7n4FMuYFQSrwlnVDJ45Ndh3fPMys55cblvzTQSeu7igNmHaDlLkCXdbhW18MH7dto4py/LhzjxEZ1TaAjuiDd/OSpOFNnRWM/X9S2j7axrdcWGnuBFrDp2Q1PT2+uX0/ZhOn1AFyjMyu7NmuCr2IcgTfoIl2HTPp9+e+0lfqLC7TckkHk/W8t6Cnn3X2mdondknq5F5rj+vOCRkLolvBB17MDsxzQMJRgl7XT7BgssbtPC+44ArLN986E/ufq/632PmnCL7undbRPuQqdj6HJZ6teRdcSfTxZD5156eOBOicxoY+NxY22ZxvHq2zBRUoxevHFNEy/2kvywjYtqfehE8C2oHMGzpb43ppqorIFuAF6/sA8BZjEsPLQEGWHyOOXmV8fHo5zK/yom8JSf/pQEPvxxNhYOvLFhC3a3JvOtSSZln0mdf5KUjVdLNzQRZsJmoJvtVuMk5D6kSghYrtcg4TUXQI/Mh6e3hzX7/YmYInePZOOjCeDCElAM5B7jXR9GpvJeu0hba2DOTYR0KdLDUSQxPHX0bHsW2Fqo6TkmJZO0YGB1INOaQlwtqQki6eC405o1qmZl7z1Jo/UQZ9n24Ed2q8IpCuB9hE7rqvAhDmPy2ZWyZF+6minPwMD8ViqkPiYcOr1nCgUztMCxwzjOKpCqr0cmFR3W2dps47Drm6E0ieJM0UBOh9OyN2SukUN2OTq4iFStxa7Xc5a5QrYM55Y5C55u1a7hUvq7ufmePp/9z1abcbgB6CTJvWQk9UglplISH3Fz6HjsaQkMABXQpfNVr0SzTLX6uM6aT5OtBGdYrXf5h13MG6yLqu25Z+kbw7t4umwUP3etXeLUAGW2/F+StKEvgsdSXU1sNKZDprYbRTSZn1jtDe5XlN4rsuWqkDnCCzfocNNx0dGYGQkVtnsabdTJC70VuDpQ0Hsx4DdtvIixF1idhtn1vJVpLNBLXqIG1tMb0wugE4IZvLYQlwWbcWk+zq7q59DQvBRb/K53Esc81ge8FYeQNkK2o3szcCukys0uYd0+BaI9w89pDXe/NJBh9wrJN2ZtpNNZjJlCNqRQsqDg7rdnunRWic5ClMiR9UJKFQi4qnKMbtFh8riRmANDqaChoDU0+EBETpXoY/NSvRxmiBtlXdr4wG0Q7oD9PlpnODuKWST41JWbwW4AkqrYfm1Wkq0aODEv7vjLvC0oCD2pwhrKdnrZXTtWl00CaOdX00coZEpn4tpLfdVWLoN7fQzoQiHzIXWjHSXItAhgCkJxr3iIEzo7ueQpALeZy+TKJWh5MKf5iznyjQWRyKYCsg/1hR0pm1ynbpmSGvQhcnqZdh+HvRF6YYgJRvqWNe9VKfQWnvXWmhshLoTe93KcXcgUP74AqV02OCGDVpqaBP32FDDAI6W3G2xMhsiaL+7E1LW8dTgIGW0zRCXmriadMhoj5lGcp7ZBiulg2YZqxOuNzkDflE2gNco+Jo5XpPAFuKKmRb2ackf96g9zk930IbQZinU+YHTTvQATjbsFmGnCKMijN54Y5zkAfqiHL3xRkad0LAqhq8HBmDduuy5Yy2kQNWpQ+fD9OuJ5ReLRjUhdRepmi0+ebvx6DjzjkRaGolWayvbfcUhMpYiIpjuSv53p6eW9RG4YlySPxLpl1oB+6P0WH3pyIZMGlJ3rfXxDhg7l7izU7zpCrAlORQH9gA7dPx/SFYdXbuWu0XYbV7TRi5o2jZ77353HFkSx6oHrO+Q1BLHsLdIQrPftedTHvZ643KjdZYrReeHdBVPG3yUIvUbtPEQ9581mD4LHr/Qa1/oho5CQu6r0c7oPrQ2dj2a1HHeN8HejclPaX/Gen9/htRBH8/QGTSXICLjIvJjERkTkUx4n2j8tYjsEJH/EJEXevMXi8huEbn5eI3xmC12EXkWcBv6VHgS3SOwjcCikw++XhnMkBsYgA0bmBoaoom+mLqBpvnca5e7FW0duY+vJuNoehGptnPR4cQKrS+GHZUkCgS8krl5pBiy4PdHid4dw7fKzfu0U1XF/T9UZL1kN+KsqzwAFbP3QVmGRLLZvyUZpH9PMJ/XkFWmbBLT5OlJjZlmZCQEJwqmAXQa3f2CVXDP4GBmE+5mbdAQtVrQmr703e/WN3R/xvBwpo57I5SFTH4NFcgmO42aPIhQdM5eEbpsOYUV+tV5BSy3EtQVetqhF+iPnQ+QaCIVaDzH+GwWJ92y/LaIbverjseAW+GejVpyvORs0X1/JoCbdERSvO9mM81Ag3EXy09EmeTTPN/MseNXlVKTOfNeDZxvXr8MfNy8W/wZcNesjsbDbEgxTeDdSql/F5EzgFER+aZS6rh0BjlesHq5RQ9kHGbx4ySAqbo4RZLUET8WA9Wt+oKL9U2D6bOShhdxa7kmsDipgW5JPSZ0l7jd8EUXLrHb0MMMqVtE6ffpTmgeZchA6QAZPV21AAAgAElEQVRpM8+KUA3zNNBCoy/3QmN7Mu7A8PL69MT14xfrImHlBnAWdCxLk/X0nqTUwEWrnBk7ksJhFnX0hZCXSRpX3/SyjKOREe4bGUknKbVB6naaS+5xp6wNG7Qm7TndU7gd7Si2WIE+z7rTzcQ7Hst+tdzQx7DWic74jfRne6NMdb/aBx2T2iJvGOfnzqEhei7Lrteemq6efsILf504/FfgNqWThL4nIlUReaZS6mERuRTdrGMz+lc7LjhmYldKPQw8bP7fLyL3Aks5Ti2fjheW2p6P5rG3UavREElZF+78iimQFIfFmUgBNgN7zOPwCuKqixZTZycNLCC5oNymFuCRuv+cC2kt2iX0cm/iFAXPMdoGobci+HI9WaZZ0Z9LB5J1WZT2JoMLhUo2x2GJMxyXZc2d0m/WYeHeDO3NFKCjXzfwXnozidPa+Dfu2eOE6a2Crj4TCmk2XW+RIj968cWxhV8fHk7i4Net08MeG0t9L5RZ6evqbty8bRQSuxecJ8GuVcBWcy6aBiwRpB2gJiPokCm1YHvUVpqJaGblFvukWNkHLwFuPRMoJTX57fG152DvYb3urrUwYcZVRauKZdKyu0/qof2eR1DAN0REAX+nlLrFm78UeND5vAtYKiK/AG4E1qLb4h03zKrzVER6gRcA3w/MexvwNoDzzjtvNjc7a4h7Pg4M6GJcIyPcZ8LHQB8smy3Kh0ield0Suquhcwuxw8k+/jWqUHuGjmqpAdeYTkXdT+iL0TpIwauBbhGbRAaWbRpo7RqyhA7Geq46XyIhYUviLpmHiP3wQljo9dBrViCayC4LWqeP5RrnvdwLjS1wRiMs1dtFp3TJAauvu8RjkSIfNLGVXw8lawNVgGth5x5o3HEHjI2xfWiI3q3QuUq/pkw4H+vWxYSU0c1N/Z4yZJuaONa5S2J55J6CIXXrtbAGwlI73zxt2NwuWzsofiox0ooldWswxLPLMHFeumFKXGDtMU32V0ZwQ4e+gdpzb9JEXY3gyDNr4KLN+klnAmO9j43RHBmJs0zzIl9OFnJXp6XzQlqg29PNbwkQ98uUUg+JyNnAN0Vku1LKbWYdSPFFAX8MfE0p9aAEnPWziVkjdhGpAF8G1iul9vnzzcG5BXRJgdna7qzDRj+YuF17/YDm8Qh9gnfdjNYa0fJK/PjbTdz5x06fOlvr5jeXHMIuwbVlqCxMl9dNkbprqbv/u5bu/kgTuh/xAg6pt0nohxdm1xFb0wvJuCIbPdpqt5a8fQdD7nvTK3K1d0+ij2FuWP2ke6v66D8APY9qooojOiD1hNRxBZnmxHaTOw2pl0Fr5T4pjY1p4l63jvrAAJGxrOvr1wcJLLdRh1sYzkF5ZIQm+vxKPfXZ8Ngt6OShVbB0h5GW3GV6tPO9GaWLv/momaI7NXRnLHs8G1XtjO438/qdJ8VPmHLKzUgbJqUKcDN0bQc+ALtMQ/BlpFs2po7HSUDmTxGTM5UUUEo9ZN4fEZGvAC8GXGLfBTzL+bwMeAhYBfwXEflj9MW0UETqSqlrZ3MHYJaIXUQWoEn9M0qpv5+NdT7tOFvYvUcf7YaxzGya9JSxrECf+M0NGxgfGuLS3wa+lG5u0XwmYELJbG30HZU0adv13GBaxeXKLr5zFO+zdY66zskYnovQJfUQsbuk7jthW1k6rizjw8o0fjikf4Ny12/MVzcaJnaalvT/fXWoPmI05DiujvjJyTbMjmvHGJTR5LhzK0xZB9/goC7Y5Vnjrn9lFG2lloeG9HK2WJuBG9/uyiyRIUCf+JqYSJd162iOjTFpkqECvKzRB6Wbgd5ETqk75ZlvN92x1gS+6p5vby/B9Uug17RJ9JfrR4eWWiu/sVD7MeJxrYCuz0PXTbpuERDncbiYw6Q+I0SkEzjNyM6dwCuBD3uL/SNwjYhsQjtNHzeS9e8761kHrDwepA6zExUjwCeBe5VSf3XsQzoxuG9PoFiXseTKmNK5O6C8R6dPRxDr6m74XcOxwO3jrU/qkLbMa5CE+vmkHiJVS/7B+HSXLW0IY2eY0C2Z+w5M/0YSwyH/hYc9C98hd0v0052mcLzVjLxtWM+zO2yz7VrAgdp7CPoeSUdwtMR2U3/MKe1w3x6ou7LBHXdQN9UJl9tmHd/NOs1Hh4djK3WC5FxxL6C4wNjYGGXHGPBhpZXm2BiYqoglm9Bm5b0+5321DlG0yWpuSQlrFDClu171O8fNXcbeSK8D6Ag/CfnT6mVTadM9B3vh0Kd1Z6vpN+r9L6OvidmqPX88oKRtKWYmnAN8xUgpZeCzSqnNIvJ2AKXUJ4CvAa9Be0OeAN4yK1s+CsyGxf4ytDPgxyJiszj+VCn1tVlY99OGup/yPDREtVajWqvR9SFi//VF2+DAB0Z0nReT7BFNaXK3kS62wQXk9Ay1aDhWuusI9bNKISu/lHsdaz0nhBHySd0bR0uE5jccUg9dMLEVHyD1I860CbLkHiXJSv3o/qtXOiUFQFuu04ug1O2s18gw5YbpvrRDr7JsLOd6oIzs6J13wrp1TA0McM/QEL17sgGeAFVL/I8olmNI3FjctkyzW5LZyiz+0a7i1Aj6+oi+claTxI5b+aOPOKLqwHPSpG7hkjoQn08195C775GZF5nzMfC72Rr5vQTkGBcVfTOqGst9+SnSM1Yp9TPgksD0Tzj/K+C/z7CeYWB4locXYzaiYr5N2FkwZ3C3SKb4fxNTqOuraL38MqNxfgk6AxGojYX6wtviRbUE5RQ3LtK++xEiDZIol+Y4qSqKQVKPaEnqedEu7k3EH8MR876AMLkfrQXklx6w65wgnQdlicoQUS1KR/cBlO5FM2cfWK/i9CJD/OPEBa4idITJFCT9Y4eG4szHyuAgdVNGwMItXgU6FBG0I3a5nWijY8bGMpZ7HnrQkTushOlnZS3IZpRILNbRCck0FxlSt3Ajp0LGgS9/uTdTx7expsOQu5VjrNxV1yGQbEsii+YCjsJ5Oi9wypcU2GtjlE0SS5m4aKDuq1mHva+F8f5+ncBiEo8ev1B/v7EwraFbnXI1QMlrS+f+707zMykzF+RqTYi28FbGUerJLpAm9DynaIOkKQckRO7Dne5eHO44F5Jz8/B2xo79yHgyzdZ4dxNdzbastHB7B1TOheWP6VVOv9EkIn0S3VlpkanPbgtZmZBAa1VWQecdGKf4XpGk4YSHOsYit2TvlAfI1IAxstwyYFetRt2k9lcgrrVvj0LJxtPXtcPXNu22cOsC1cvQY342NzwW4HUQrhXU6t2FPccCNoF90tyEJveehpG+bk+vYu9GfZhtbfbRtWvzHcgFnnaccsSeCjsbHIw51UYsViEpnlQBNhmrZN06GBtjeuMIvCcJXdxEIrP0R9qB5YYupuATejvwHaOh6BcAIjh0TlpyObww66h0x2EJe3eL7S91/rfL+5ZPA4i8qJmQMzU0fmvBB6oX9Hvb2dQJ1zwHzvzXhNc6t6CdiwdJOihtBmwkiQltBC2xVLcagr1C15YJRXSUIS4R0Q4q6G0s25o8fCxdBdwApQngJiPv9GW/6/tkrH7eX8qeQ5MlQ+rj/kpIn1Oh38eeA77wb5Z19fXrzXv3tM6KjqOOdiTvcamDnJDPAicWpxSxZ2KJ77iDqeFhGBnhUqWYFkknG1WANbo06X1Gpy19Hn5xnteKzqA2Ba8zVmYq/R/SFpQrvbgXoVsuQK3QYYy+tQ6O1W6vVIfU8yQX34JzST1QxwbQZqhP+ku9z66F7cKOJc5QjbJSTFzm19kV83JJ3XU6b1kCg6vQfo8dwDW6QmHnw+jf6/0Kbhb2Gku5y/nubkO8la3QtzX//rp8FeweGdE1yAPdlu5xYpCNLM6BrVo7X25LCF8Djz8Hys+Bzh4SecQ9RJG2xsc70mWK4+QsT09nCk3qIWe3hbXG8yKbpkgyi5xIWLvNEeOgtuMDvR9nXk0SffQaBaE47LGxuOzCyUbwT0raLTTfMe+JPc4WDXTOAWDduvgRsqQUfNg7YbuB6+CCTehH/gv1xeheiCk0oDZB2AEK6QvOlV0sobshjMHYdDvdXJXTXTqeHLLhi6EwSdCkbgl7l/cOSWTGLtIFS+wyS9G6e2ZsTlSMfS/XTWSMOwA0qVsJxvf9RmkL0j/OO86GvrckYY3RlHaadpxvdPFHNPE0ICb/6Q/oyEg2bKA+PMx4rUYPaYu9aUbYaQqIVbfW2D44mClmdZFS8AXRN5U9yWHu3Ky3xRX6ZtNYCPXTobIYqs80Dl1ziKYXaenFhirWvPMl3mf3KS9PfnHh+0lCN11L7vaGatZ/faRJvfNhvdihs5JQXis9nrlEH+PlNpkPUqGd9ZxaOQWeXsx7Ym8HNiU8tthvQGeNmjTt8lnQcY2+GMdNy7Y1mNAxSHuQpghaZzHypJFWpO42wPBJPS+EMbRdmJnU3c/LSMjdX9Yn9wawMPC04FvtkJC6K73Y4zKh3zc51qQl+evRiUndTySr73wAncjTA9MNoXRx4Oa8TMsv0VZdFqBcq6WUn+bgII3+fppo8o/2AMbiXwGMDw7SRJfKXa6ckrU3Q2cdOm/VY7AVFS0ZVvZpx+P4YpMFei4se0h/dft5Oo9hE2Rr//jk3SDsJJ3JoW0PuXuMQ98375si6H8GLDOf3TDeyoRxWN8ucJvSAQfW9zA8TMQJKu7VJp6UHHl0nuIU2tXWuEeEHnQdEUvqtlgXAOfqN3tyxBalISIga1WFdO0GmsBBp9ZHzrRcp6iBm0kainjZN8OzZh6ph6x1F3lSjY99C/Vw/fIDMRqJrHQGadKBhLysUN1LLBtcj45jrzSTaoQd9wPbAFMfprQeeCwnQKtPa+C7TGp+l5HcLrgtIaP7jLzgcmjVDKPTxri/SZKKnXUzvhv0soecdoaRs5IeEzE1EcHEc51mKSGfi3/+5Pll/M9uz9nQgq78Zyd70VibqrCmEypOyejKPmPBjxM3Cd+1UfSTjhNJdDKT+qmIeU/scY0Pp5a1nxJeNsWVupagnabmcdmGMLqFuSyhx5YWhB1X9gL1NU/XKj+4BfaPm0H0eglHPutB7PnyY9NBk2qeaOw6SiHfUvfhSzE+jhCO4llIUmqg4xek4hdznb9kHag9yawa0O2creUG+nfarOWQaI/pGXtNsqoGaH/JFj3NOjg7lyTL8SmBt2pSWkZYvbI3Aa50xrSDuDTz9IuyES629G19caKjX4fxyewnbQy4Vrl7V3HPnf1RIl25cAu/hfwXdhvuOl3Yc7RX3zy7p5PInErTcYdv19mmE6TzAYBM4bQCxwYRuUYpdfOxrGPeE3sFUmVVd4vQMCnhVl+tYqIYHO4vPQaRRyxx7DCeJhrSNatk66ZDmthsWd14mutR9byrbiVGl9CfeEb+Y3ke0UPWSs8jeVeWse+uA3U/aUnGDjsm9QCr2LK97i46jlOiJD3eDb+jBD1lE1dtomCmtyaGZ3WPSb2/QseLd+7Q7QftA0DXKpMc1EfS2QrgWmHaVOTsNHq8q6ZN25uGxQRM/4FeZukOKH0Jms7TXfzIvzhxjG6CcLkIf5pvGOx3iLuxhRTKvVBZ7RzTLckTUSu452sV6DEyl7dYvWxuVpPAZkcxyvNXncSYY1LMW9HPoU8Zc2dXnyLqmJhlQ+496ESj6T2JcdS1Cm1OdZPUHjEn/6QT0hhnj7pk3YpUcyUWA7fPKJAhc2idOeo7SEPjcC31mbR1O62Vle7Dknpg6ElJAYspTTqyPduJrw1c/biuEdOoJtmQbvaj/V1LNjKlT0+L9phN9aHvGBPAzYnjc9Ksf9ke7QAtrYIe54ZRsjIMWqIroRuSVzeiwybP0jq6WyExZQj4N/iG93Kf7nCWUSs0cdsnOEvekCSpuSztNxe3sDdOH4bU1wRmAfTugzPvRctdfWQKqhU4eTH/iX3DBupDQ0yK0I0+xztJE0Ic3ghxS6RD5mK9Dq8cgG9h5WmXbsii7wgNxgi2QeouDudIL6EkI5fQ27HUXXL3Sd53noZ2BZIxx47ThialxhZdi90SjU9yPsEZ/MMBOOfrwO3QYRPEeqByIXTeC9wE7EisaJahTdAV0HU7cc0VImBCk/ounGbLwPahISpbkxZzVZyKiquBXiO59AA3QvQ+2Hmujq+3T3IWcVRLKJLF+mXsy9XH95OfWVxeAZUVzkEK6GCuo918LdXR0C5uSD1UM6bS1E7ejvtJHl1uUyy9DSacMMeTLaSxFabnlsX+SyKSqZCLzvBXSqnFM61g7uzqUSKOWff0dJcz4kxAH5PAWfpEiLsXQVgPjXUAs8wE2kqvrAvUcbEjCGk3/ujQ5GhDGd2EnxCpuzq6H3vuE3o7OFrLPQ+x1Z6jC+VY7f2kC6fFdcV3APXE/1E5HXqq0FkB1jkrqJA8fdniWpGx9Nscuk1kss3Gbfjf5Om6A9F1S7yKnKH98WUXGzVltXEwVbTGiXvTZornkPPZ3YiDPMvdXU2ULo285oA+xtFh7fjteNiMcxc6N8Agr8tUgVnFj5VSLziWFcxLYg81GU7hCtJZgCtIvEQ1YAt0rIGec9HmfYhEXallPzDV0BfMwRWJ7pmBf8FaiysKdzNq9Djhi044Y8jD50a8QJjIH5yGJ++Cyj3QcRE8dhlIC5rz9fWjgVunPZ7Wm+jrIYeew1t+NczvnakTkzp7dcIMJO3byg1dm7xjGJZuQjs5twOb0xp854egdAXQr/+/YDOwAw6YNnRx5It9LUNrd906Lj3jCPX32X/ywPt/gvR+h8Jag+dIKzge1vKKZLKbJ+CfLw6xrwFe8jh0/1zPihtbTxA7nQvMPcw7Yp+J1EurSCpKWd7pMdrpYyR1RuragrGWY4zQxVvu1f1RjvRmpZcYTphiCjl10v149FCJ3RChQ5jU1TQcXgdPjsGjB4FFcNoALBwOk7ufmDQTuedZrJneqKSPXUT6IcbAJ80a0Hs29JiH0Mq+dEghoDWUNejf1ThH3U11bgZW6CSz0hXE1RNTv4hb250kl8E2s0iVkQjtk2vBh+QWi0zJZd9hHjpXyLYitOdMuW4amzSS9U5vT48vInO8ew/pY1ly+6PW0VEwG/WudN2WHcZcQ5NwI5KTFF8EEJFnKKUedWeIyHOUUg/MtIJ5Rew+qVsNcCemUuNaTDGX9PfimPXnQuUSqLxFZzde5TUjCOrBsRMrBO9KatU02i3YtfBwut65T+b2Pc9KD+noB+7SpI7N7nkC1BicdRdMvSK9vI9l3v9+WQGf1FOWujlQze1J1IatFmkPje9IbZjysgZWLrCoPmKkAnebDTQh2U0H6rLYG3bpoOmy5HS+soimPJIjsdRTZSRCN3hL6G7YYqixeLnXWNctiNy90bvH0+96ZZ+M4qYmZn3lXmiuIG4a7h5r53jXy6Yx+FnOfk/oyKDtZnfcsgwFNESkhHYr71ZKXenNezu6bO80+ox8m1LqJyLya+ish4XAYeA9Sql/8detlPoL8+8/icirbUc6EXk+8AXyG4vFmNPE3so6dx07y5XS8cpGro4LZRvNtd6jHaW3225GZwaOnHvBulZYruwCyQXb1bqOC7Suwhh69zFTFumhe4CD6XnqIBz+CSx7BRnMRPDWcWphicMmJ8VkY4l9PByt0UJpsL/BGvNeaWrrsuMxEqeerd7mkjpoea1PR8SAc9h2EZP+rnOz26wshmqU3DiaUdK2cBMmzDUUtuha6/v9u8XqbCkId8dDdX7ccyGU0WsRk76nt1iDoz4OkZme83RUL0OlCp0H0YlI1+p7YGMelQc4SE4JkKeOdwH3AiFH5mdtfXYR+Q3gr9AC8CTw66Zfaj/wdbJmkou/QJP7a4ELgNtwujC1wpwl9jzrPBdvNfU97NPp6qTN2PhiTRpXYnpENswFbK9Be8GGNFL3PYZjHrkVF6F1hIu1YnEs9ZkIfSnZMEYXlozrF4EsAuXk48siWPj85HMrycVOd6NhfCw87FmQgUH7pQTwFjPTLJlba70GsWbSA3S633HlE7e7xXrtKK1uNas18ovNEO15NLHG41WVgbN156Byw8gwkckUdc+BPItdrUhfVamwRE+Os2RucxF82PPBr17lHvuysdqDBzMA73ebLEGfyealAdwOu/ek12KvtbkUBXM8ISLLgNcCfw78iT/f6/nciRZqUUr9yJl+DxCJSIdS6lBoO0qpr5q2o99A52q/Til1fztjnJPE3i6pj954o55vW3a9QcGbhOmNUPok8CItuVSauv5IZSFcUzbNHWzqoiX1cbKhjRm4sksLK92XWHypJSS9QH69dItWiUYXXAZTA9AY05Y6iyAagM7L9HzX4ndrxECY1PMCNmKnqUcwfn2YgOEa73s1sdb70dmQ/TY2vBOuLGuju/NhEgkGEuu94nxeD53WGbpSO1khKXRVbgCm2bjFRASVsj4vdlR0r9D49w9p6TbCJVS4za/C6ZeC8COc8v7Hm+4eu7gej3cwW0TF2ONrS/N23A9crytgxpsNNO0+BdAtItucz7copW7xlrkJ+J9osg1CRP47mvQXAoFHYn4L+FGI1EXkbzA3A4PFwM+Ad4gISql3zrQTc5LYU53f160L1oEedYoUjYpwqVLsFWES7Ui56Foo3Z7uTVovw80lJ4TNXjz2Qt4fQaORDh/OdDGqZp1eR5tU5KIVmbeqo55ZtgRLh7XWvucnsOT5mtSt4zSPzNuJiLESTNBan0q09SWkrXVXW/cIze15ahN+rNOy1gHXGkdq9RGIXSEOr00v0lZ5R0TsV5k+K6kIecj53y8bQQn6OzTp3VwyT292jL5z9EhEuBpn4JzwHeRubR//fGj1u/s31cyN1DupIvRjjnnZLFPbGPzMf4Xdbwxs5yTvY3o0OAopZlIptTJvpohcCTyilBoVkdV5yymlPgZ8TER+D/3A92ZnHRcBH0E3wg5hm/d5NLhUC8xJYr9048a4JnbD1MwevfNOME2JY2PGdMsBXeCpgen4MjDAfUNDXDAOXKov6i0mySQuuevXLbGO0iPj+sORSEfA5AmXIbik3orQQ9mi7SIvoQg0uS97BVRCBoT5bojQrQoYKtXrkno0kdbVaRAnJdmCZ761npOoVIuS2jypyKSG+Z0i3eHnygh6q9rp6Ta5blSNU7CqHYPlhnEOTuhKkIeeqZfZcXa6HrrdTj/6KWGTXad9gmspx0XpV6us4X0BuS30u/vHPnSaxY5TZ2DWWl/QSJVruB7d4avS1K9zvqNJfcK2hly/nkqtRhm4pOiIFMLLgN8Qkdegf43FIjKilLoqZ/lNwMftByPjfAV4k1JqZ+gLSqlPH+sg5ySxg66JPSoSdO7kuZrKQH1khObYGN2gy/bhWIPuBevKLinrHE3uKe3UIfZQiBqki3SFSN230EK6eUhq8Uvqhuqn02KZEHxS92vBuMZoxlI35mxze1gKcInM3jDd4+wRVz+aiPptm8Eome5i6mwtKVgrfPJ0TVyNhUBVT+88SHxyNKMAqTvWuO21mhqzG+0SW+u9JElFDqE3erIO8TzJrVXIqv297O8QOEbZ0g1GhnGzfL3v2OxSNnmn4003UV+/XhfOm0eYLeepUuq9wHsBjMU+5JO6iJzvaOGvBe4306vAV4H3KqW+MwvDycWcJXYg04Aa4BKTGTcd6vBiUDqjxqGfwneem0zrJ+ngnnKYQvLo3dHrWWjexWzhZozmSS/tWOV5ES55ZQDy5kG6eUaI3PP0dMhGv9h3NwImruI4lZCKC/cYTJAcY9dB7TCMTZ6xpXp7gGsiHZ3iygjVR3SUTGURTJ4HLEyaQbtoRjrJqHxWQuo2CiozxilvmhvG6OrpHb3EpZR9n4qVWkIKSah8civ4ZRwsguWRG05p5EbGl7EJXR2h+wmjq28O2BmmyFfhNG0fIvJhYJtS6h+Ba0RkEP1LP0Yiw1yD9va8T0TeZ6a9Uin1yGyPZ04T+6WXX55ypIaa6abKBlwB06+Hx3t0z1IXlkg2+bq6q6V2kP/IbWEttX0tYtBnakvXqoZLu8izzNvRzPMCsFxr3SV1P8N0puqCkD0udt09sCZKF6bK7SGbA5/UrUTTjJIoGNu5CDzr3ydyW2Y5Wg0d44mgH8ehm8gn90burgeyv3er3zHv93HJvaXq5xC7hTm2/ZGuubN8FLgVDmzUYY3NgIFUYGYopbZg8nOVUu93pr8rZ/nrSVrK5kJE3gts9qJojgqzQuwi8il0tOAjSqkZg+dnEyFrYq8IESZF/GozsQL0JHWyexpJ0okbebEJ0hc35v9UnXQHrULXWmnp3w9Mm8nqdhHSwp9KTRj3u36cut+82rXSwSN1Z0fdIlRuNIyFmyTj6e5rSCz1ellLJZbUt3iZgztMBExf3cRimzh3SH5n8IjdrBP0725lmNQuRKR/t7i6oiVzsxNWQ3dLP6SMAbNoXhG2VgXYQgQfaknoJ4O5pO4+AUVwx+Nwzpdhryk57LaMnO9WeYNZj2M/nngAeJeIXALcDfwz8A2l1GOtv5Zgtiz2YXT94JMi+dg+4cfZh+ZWc8jLNHRTjOMiU5Ats7o/MolIvqlkPrsXt0vq+9EXtXtBHqF9Ug/p5X4WaB6OJqrFRShO3R5Q10q3724qu5VhQqTuk7hbcRDiEDx7k4XkN7FRMaun9Q3Ztd7rZU3wPY1EY3eRcqiawmHuBR47ZkNVGK3D3H9KO3SO2b9AyOIU4WJsebLaQf/OX4f7TfZcq986IhCF5BD7gkYSCdOrzcTqI8Ctzq4a/9R8J/W5BqXUJoyNKSIvQGdg/L3Jdr0Dbc3/oNU6ZoXYlVLfEpHe2VjXbKAHLcFMbwW26pj16RclpD7xDB0P3XtI+09v6ADso/n9BJofrA5Y64GohzxL3b/A25VZ/GkHG7ArSqQa18p2if9oi3flJR+575ZILEKk7uvq7nddUu8hW0Z2KnFY9peypH7lIVjxc62nHzpLO0utlUZkjCUAACAASURBVA6G7E/XuQjxZ9I3bBuTXvN/J9856kswtvOznyEaIvSQdu7/5jGZexKWO22XOfju72slMtdp7ZO6jUAyhG5JvddGS1+ta80DTJniZ5wKFRsVrSPRTlIYOeZHwP8SkcXAr6F1iONP7O1ARN4GvA3gvPPOO67bKq0irs0NsHQLlPp0HW1bOqCf5KLvN9X62E62IULpCtINDRxdPRSfnodWMolrvS0KCKiuVWfJ3f1uniPU/9/dljsvpKn7kSqWSCB5z7PUIXHehSx1V6u32zDvm9CygdsGb/W0bvpgS8l2GJ8gTm9OSHISxjvSYYvd0/oG8TpIEo3cbbvEbgm9ujoZ9HRXNsrFvXnbJ7PQU9cuzO/XgsRD8H/nl3jzU4lgDR2FVB9OSN3cPO1TUI8da4+umdS7UR+KAnMHJqP1y+bVEk8bsZvsrVsAVq5ceXxNhO/qTvVLbzafjScumgIWJ/HK1ziRAqlUccy7WtE6PhkSC879HmR7jLrE7lvZPnG3C3sjCK23FaG7y+aFM7rEu/CwjlHv+IWZMIMTIUTqjo6eIla7TWc1m9BWu0WlqW/IvRdCZOLPQUsvfra9q5/HHYyM1c8ETn83bxfs7x13LLKZoobU83wnVloL3bhTlvlMxO4F6S6K0nJaq6zf5hgc+kQ6AcyBzS4tN8xmblOUNgvVPa1tkQJzF3MyKsZNRjJqJABL3UfK9ytYaUIebZG8hrHSjUX3duM4rbnSC97/PrGHSH3fwuzFDvlRELn6qosWF739nm/dhySYvKibkCPOJWGrp7uNqTNE7um6oNvenUG2FZurqVvr2CUhh9wtIa9BSwjW0T2+GFic6OkQJncLS+5+jHrmd44bQq9OR7u4Upt7A/B/3++R81vWySf10LJ+TQSyjtSMDDOhLXU/q9e81mCO15Sptf5yc408ophqERI87/Akp9RdbE4SOwDr1tEcGGBqaCi+FHabE3Wp0jVhuBIOmT4ktqCT1V7to3omIsL98f0U8RCp+xJMO6FtQY21EpiGN61Flb92nKpuyd12SN1KLtGEI7tYeLqu6yjtRUsBDZKYdfuVkATjj4EkFHGNscBtuYdNwPUduuOPCyvDuNmj8XqmvFeG0HtJFepyY9KfeIZTE4ZshMsu4OAkWdTj9zP8qpoO9rMod14MexOuAovdp6epcGip2Y01kfZNdD9hKmJeCXv36P6/o6cSqc9RiMiXgU8B/6yUevJovjtb4Y6fw+Q9iMgu4ANKqU/OxrpDsF78UYA77ogz5ezO7BRhuQl1nDI6bHRYZyOO+0kpIdiLvp3C/O4NoRWpZ6z0uvNeIU3oeUSeY8W3qu/iWvGtKjO2jHpxywQ4yCP1XhJij8hWRHQtdbvKSGvrrpW9ybyv6XBklQZcFwGd6Qga0OGQ1/nDDJE6pBON4sxRJyb98MIsobsaeko7z78ZW1LvpZmaOz7Tpec6ys3xyZC63RnblcqRvPqrcM10Yq0zDnv3wHh/P+N33hlM7pvXmJsW+8eBtwB/LSJfBIaVUtvb+eJsRcX87mys52gRJygF0qDv2wPLLoNzvgo7L9OP8taiS8F/RLevUKU+v1JjqH56K1KPUQ+8uxZ7CxnGR6vIF5/U85DnHM0tv9tICN0l9Zy6L5mYdbx350bj/z4uwbuwFn1/ydHTG2Q9gi6p27FlpJdqWku354Ndl/t7xjfoSfJI3bXQfUJ3p49T5gwOela7+b2txr6M9M04VfDLPDHJ9iTSaAC+PQ0VEwVT2QcdPzL1YDZsOFUrNs5JKKXuAO4QkTOB3wW+KSIPAv8bGFFK5ZaKm7tSjIFL7i7q69czXqvR81pY/nn4xavg1jMDKeTjZCNh1ArdjBrImLeuFOPrrqHSAK3ilYOf86ZDqh6te+HnbdOHlV98qzkim0Ga6cpjkOqE1EhIx9XQJ5zt2PUHiN1GbFhscv6301PSSiM7P/7OFOkCXb7k4xJ6ZXU4c9T1lYyTJfSMMzQ5XmdwEMWTTPMjFvEzFtHLGVxCL9kn6Hxr3SP1l5C+IYeenhpb4FlAL6zphZse1/HqzUi/Ord6Rb4KzCmISBdwFbAWHfb4GeBX0GUKVud9b84TO6DL+OZ0enHpNEjqlgys8y7Fvw4L+THrIWsdWpCBRbvOtBBcy74F3KiYkATjv/s1RzKWunMHC3VBCoUwus7RUCSHgVu98XqyGnkMx/pf432vZkm91ffiiBfPOepGvLix6CkN3UouFmlCB1A8ieJ6GtzPQY5wGgs4iz56eQ/PaYvcAw7TkLbu/ib1YVjSiOUv9yZZthLMNdkgJDgFk5LmYBy7iPw9uhfYRnTnJdsQ8vNezfgM5jyx3+3UXfcRAUuXAL3pbMVUGVgLn9Sb4w4JGFYq19MO06M+UWYi9ZkIP0DqodBGi2Xoiz4ke4RKA8Tx6b6mbv63CUhuyQBH180kHbmk7x4rx1q3OvlqM2s1Xm9Rd30Oqa+22rr1g0Rko15AE/7+iHSTcXMQXOnFjUXP1dDzdfOIH/BD7sX2R3iSIzzKdhT/Dswkf3iRMAC/TBJdZMNN4yeoKU3qi7bH/ozrzXGZiPRXyg0ofUXr6r4YdMqR+tzFzaGeqACtasbDHCf23SKUyZ64AOVaTbc4vUFnKs5YRMolBKsdx/proBSvTzozIhT6VvHmh6b78OblFfrydfV2SN2vpQ5kIl8s3DrfvjUestB9cvf3yvkR4xh2b/l+793+7za+Tg3fVmOMM0gdHDonTepuLHpM6PoxwNfMLaG7xP51vk+66Q0oFLv5Ab9CPw+ZS8231DP6+qJISzD2CSuVQ9DQY5r6oLbUnRuqX6Cp9EN0hS8HpzShz03naVVEftOb9jjw45kqQs5pYgdjJI6MUDdNNcoYS30tuixZJSn+BDkFn9zPkO84DYU3QjZqAu//lE7gRsCELHTXevPJvg0ZBvLj1CHtJIUcUncySf0SsC5cmcUSqS+9+Ba0kb3ceuf9Jf2yES4pZ6mz/tqUbllYI3Gabmok68xsM5Vw5DxaHDoHHjlPL9OC1H3L3CV1gHOd/xeRn3Pnk3qa3J3YdUvqS0nkFzeHoPk9nYj0LGdfo+Spp6eh+7h22Ergtym6boOu3JEVOMnxB8Aq4F/N59XorInniciHlVK5nVDmNLG7QYHLgM4PoRWpnvRytkZMrN+6ZONbkkci76g0iLsoh/qX+nA19jhywo7Sjjoh93RURJvE3Qqtko9CpB7D09Hdl1vIy1+nby37ZA/ZG6Gz3CYn7b3fDy+NdAik/b9mDNZa1bPS/TGBll+q69IDmu6Cep92kob0dI/UQ4TuknnT+f9iXsh/8EOUQ/CC0M8LgaylPk7Z/O4eqS9DSzDWoR3nEExo+eVZpH0aRoaJa9NbFfatp0D9lxMAEYmAb2GKeANfUkp9ILDcG4APoh/j7lZK/Z6Z/mZMBRPg+hm6JT0JXKiU+oX57jnoEMhfNmOYp8Terx9A68BErUb1A7B8LVqI7c4u35a2Dvklest1OOwUcnejYXILPR0NQjHtofmAmob6XTB+D1QugoOXwbMMM/oNMtxH+pZo03EQIvN2lo+8acbKjxucGLi/Uc13bFsd3R+DG9a4JzJRTU5ZAFuB85GFye9mLfX7XS29HpO6fZ3riX3NgPh3Ps+nlz4e5AGaNClT5jyey4X8Uup5bZyyQ+o9ZEi9l/RN2HalssW93MqYPfAPwEseh+6fQ+lB9MPhLuA1mSGe2pg9KeYQ8AqlVF1EFgDfFpF/Vkp9zy4gIuejuyy9TCn1mIicbaY/A/gAsBJN+KMi8o8tyvH2WlI3eAR4nlLqURFp2dp+ThO7H+I4tX49ezfW6OpDH7pu3dS4sVBruG7Mc4zQj53qjkQ67G9hRa8w7yQJOt0q3ntC3NnMw7zIF2d6tAAOr4OJMVAHQRZBNADLhmFpKWxZu85R/8nDyjBuOzu/UYP/vy97hCSYFt+NtXFDxjXncKdkGJfQQ85Rl+yncOSXAZhemp9B+n1MKYAkHj1kpZ9LM0jk7rQyZU7jNNbxDu7nHh5mN89kKRdyCadxWrxcmtQdS936RHrd42h19Qn9u0xv1tZ6D3EC0hqg/4C21G0/V27VVU1LcduHArMJpZQiuYAXmJf/ePTfgI9Zwnb08FcB31RKPQogIt9El+T9XM7m/k1Ebge+aD7/FvAtEekkHOwUY04T+6WXX55Kja4AXUvQckwl6VQPXrq57x90EWvrDlu5oY6uxp4Xu56Ba4nb1zTw78DPgOcCLyQJ8fBvCM7/iyKo/oshddMqSD0BR8ag+y5Y4DSqDjlLU/taT0oFNMf0o76LBY1EFnBhnXaQH6vuf8eztD9hKy5WSZ2iqRtvg/BvZf93yd4WbCv3Jk5vt5uVXUeG1MOyi7XSm957CO685VzIBVxMmXJKW7dSTJDUra5uj9vpjxpSH0+iX6zEaEgdc5zioIA6sAV2b9WH5KLc0RaYAd1eKOEtpoBhDFMXfRTd8eFjSim/w8LzzHLfQV/UH1RKbUb/yg86y+2idergfwd+Ex23Luh+F182N5dfbbUTc5rYAS5VKu5vWnJkGEvq1nFqQ+iC6eUWblRMuZcUE/lVHI8Kfgz0NFp++yn6ya4DfS58kGwdA4/clwHNe7Sl7qJ5ECZ/ApEhdp9YQ9q6T+pnBHbOd44ackk5UhvOZ+sQdWb3o7V0S+DXYxx9wPWdpkRA5CUsuZZ4CPY33B95pQHMXcfejA8HSD2u7zKRkV2AjJXuk7pvrbeCq60n/xtid6OX7A3UWutu9EsP4ePuog57N8JEf39uTscpjfbj2CdnCiVUSk0DA6Y59VdEpF8p5QoBZeB8tLNzGdry7keTc2hkGZibx9eVUoO0UabXx5wndoCSUvCF9DGzpN6MvK45eZKBS+r2vWzij2dKSvKxKCJb98kl939Hk7q78Z+a6S/KWWklqeZYuQhOWwRPOs09FyyCc56f/Zorw6QySVs9tpCOT7eSSWjRqmNl56yqhiHsajoBKS7xYL5Xc38TS9wh/d69McelAXqTAdvyALbeix+jbjT1kJbuCy8zWexWT3dhrfVxz2LP6OrLyMar26eo5nb9xJSDNehCaL0/06GN039gXDwDAzA2lvu9ArMHpdSUiGxByykuse8CvmdS/h8QkfvQRL+LdLboMkzP1MC6p0XkCRE5Uyn1+NGObV4QOwBvMDe+rwk0dIJGMyJb0tUlqwny46vjOPZGkpjkoqXrghxyB03wP0Nb6i4OsZCf0mFs1mB8s7Xw1GVw9gBMjmlLfcEiWDoAfZcl+2jffSdcRhw3Ox8oJJVpjDHlfBUtCXxi2tRL79CNwDeZRWIy99CPqbfeaeq7OLp5HAI5Qfg3CVnxsaO7J+0ktfJLkNQnOIM9QVIHWlrrIbjk7pN8QupLyI2AqaLll1Ro43iykilSlvr1wOo9cOYDwGaY/oBOoi4DlwQauhdg1pynIrIEOGJIfREwCHzEW+wf0LVdhkWkG/04/jNgJ/AXInKWWe6VaCdrHhrAj40WH1tlSql3zjTOOUvsd5vkJMs9nbYWu1UZSCisD7imop2n1/krcnFkhvCOhU7xb7fZc1553iAqaE29g/SZ1sFhnsdhV4MNYRewrASXD8NDd0HzJ9pSv+AyOL2UkHCoqBeQObtTFRoJ6+MWefp5DqwJ4xN8XArAH1LIIYs3LUTs9MChZ9OyKuNRknqepd5KeimbP2utB0n9/CjR1ONyAT93EpCmoLEZFmzWNdbtCd4L10faSl/2kBPWaFAFuk6FFncnHs8EPm2kktOALyilbheRDwPblFL/CHwdeKWI/AStu75HKbUXQET+DPihWdeHrSM1B181r6PGnCT20RtvhP5+yk41R84W/UDUh3Y0ocm986C23vuASgS1jqT8aybaArJV/3wWcz+2stqDNbotXo7+ve4lrbG/kLYTkR4qwbJXwFKjqZ9GltRbEbobBWNj1X391lro9qnGI/baFPxKVd8wM1KMq7s76Me4QVznqVl3f9XR443MnNH4Y708Mr/TAKkfxdXTU40w0pr6ahop6cUlcp/MH+JBzuVZ2Z0xKDt/WVJ3btTnRzquwS2cFjtKJ5ICa4u2J5p6b/Jk1L8HzryXTOBUaS103VaQ+tMBpdR/AC8ITH+/878C/sS8/OU+ha6x3s62Pm2eCs5TSt13NOOck8TOwADcdBON9eup12pEmBQiN4XakLu95qPDxA66YFy0hRtV0Q5ClnpLUgd9y/kIuh/tTmA58GIzPYfU/frcLvwa6ylCdxEILQmV3YWsVZ5npU85seYObORGyFrvnnb6zdqbQuQsa9flZrLa38o6S92CXq6j1CYeZWLUNamvpuHp6VlCfygVuJCGr6kfFam/hCypd+7UO2ZJXbanrPT+apKAdOYD6HDGjfD/t/f+UXJd1Z3vZ6tK6pK7LJWtll22flBGMhJxa6w8aQCFAHojZaIVnGTyhsRingQK9niYxCR6i87CiT0IiLOGPOR5ZsZ5cTwGFCQGkZcEwhBGxPKMCAwy2A3tcRlLuG03tiQ3+mE3dgmV2tU67499zr3nnrpVXbJaUnfrfr3KVXXr3FunSl3fu+8+3/3dufloIHMr2i0sQ3tMTROwXwW2o33MrhORlWiU/2vj7TsliT30Yc9Xq5SOacNqBu1tOVBUn5h6KW5iHTViaKu6SMsFkFw8HfePxJcsphFtDq0WXjPegRSzW5B6GhrFJLk3mXo5pc+QtQsgjhB9+GTfaq3VkZAnZ9xYgjtOJ/159uVib/XduZjQqy6fPmwjdZdjx87Jve5uiSuqctKhcZgUUh8Gaqzg1ZRF0iSxtyN0h/FIHcKq0nL8b+ckjY7UI7fGESJfdfd9lmOteuW0+qqzBw7u1L+mwjEo7YQFWaQ+nfExNOLbB2CMGRCR6zrZccb4QyYnIkOjLVvStR1efVD1SvhkF9zlX8q7HVrm1Vswmb/Jb5OWgOf/kYo0QzC3X8px/IXTdg2Ow48SEbz94JFNgI0O06oZC+Pc3HdAvK3Xe+wcGMv1+FZsxG6MLre+ewSqQyS/3hF0FdA/vn9i8bseOc+X2lJ4qTWpX86xtqTu/jsbUs8H//mkvo9CM6m/jdbqF2cVMHOPnsgqeusta/pl40movALdL+oiaa23F/bupd7b275CJUMSbvF0vNvkQiNFEdPRmXxKRuwOq9at45n16zUVsxlNpLt0jA2Uu34AS68C5trtfrReQHXCqb+QNiuE7iW/tMAn93HljmkEnvI4jdDT5uKTrp9fj1qo+ZMmvuS/nNbEXfceh9vKybFVP7Vln1fs6y5q35cLisMcabvptaorcGsgCUMvb7H0lVbpl+Qi6S9Qp0GD+lmmXoBEfj0k99dF6oVh6PoxEalfcSChVXdXPEuP6vpQvg4cn4y8k+E8oyoi/wrIWZuC3wO+08mOU4bY+x9+OPHcRexLjFGJIyh3OgI8AOzRhz1vho0rVI4XkQUouTiScuQQveBQ11RGvgZ4PjHuWK0wu+D5xaS5Nfrbg+dOrx6S+gLiqwTf6MstmEJM6lELNW+S+UpcyVjw9vU/bhiVuxOgT+o2Ko+ULcG5YzfoerCFS8GkZYOAZJQeYgTIbYDChniyvgKmTnKh1Eu/QFxJGv7nEJL6Ia+33kIqEan7Ebp7np5TL5OQNHZC6hXvey1or9KlR6H7OW9iw/Yrsq3tVj3xRIsvLMM0woeAO1GFxRdRtc0fd7LjlCH2EP0PP8yqdes47FkKlIHcl+yT5SghFqFW1k29hTaNrF3+Nqo4hQTL5WuaF+XKJEGFUftC77Ej6FN1ztq5MS1SD0ndRb1p0kY/r+400S5SD8k8LTp39wXv3tsvIVcskSDmiMi954m+o/54/woKkieROh6p2zxR2MbOt931SN0vPrrKRuutSN0nc4eFVKLHE0LqoU7dqV/C9BdJb/oIHzBwS1rhYoaOMAX92I0xP0OJ/c6z3XfKEjvYKN7r5Ti8dSuVm6vMc/a9S+Gnb9YcO2iUeZe7/A8JBDyfGLwX7K/NRe3OBMxhJkq0YZ59obctit7TIvd2bo4B3EnEz6uH0sZo0dQj9YSvOklSdwhrl6A5qg/TL3i9SYMTRaKK1E+p+CeKtB/aCPG6R9Qgw0XqV8Mri5PNpiM5YzOpr40kjXRM6BCT+rUsShA6tCP1oPgojdSdVr0xEOfUvbUNJwUt19G2dk7W+E5Nq84zhqHNm1PnnGH6QUTeBPSh13TRH6Ex5p+12sdhQohdRDYAn0alHg8aYz45Ecc9WxSqVV1rWw5U4KfXqRoGbBMHV7Xk53Yd8Zjl7Q/ud1EaLw3TCk0E7xDk1aG5C5Lvr+4TbZq0MdGvNAXe+Sq6dwTj+BFvexjde+Tc9goobWEqXM/wX3NSRtfGLnLZLMWR+gjJatJxST1d+eJIPUy9tELritIOIvXIU30EavfHOXWrfgFUTQRseFrXhbgXTuxvLjpalVWWvj5MQbkj6up4P/AgWujUMc6Z2G0F1p8Bv4T+1B61HsM/PNdjd4yBAUp9fbEXu5M52sC62NBFPGcyFfl/Oz7IbWjxTQTs59sK+BFpJw6PjrBnF9qPg+aF0rSK0DCv7u6dsZfPmP6VSP1Acv4QLdz1ljyiHva2k+KK6Uf14dzcWJ/E/ZOhr04aJl4cLVSSEbpv5uXSL6mkHtvurqWeWChNqyANSd0ndD9aDxFG6y091ZtI/XBstLao3hSpg2rV+56Gri/CiW368Rrbt2fdjy5tNIwxf/56dpyIiP0twKAx5lkAEdkN/DowscTuGxvZBSS3vdzXx4KPAxvU1TFnq02Lr1gpdEg6YdTuctAtUzF1oi5KaVhgby794lIwfr7dR9r2hSn3oZwxRMt+pf6YStLY7PLgo3lIRN8eETcZqPlv4x3nbrvpLn9MmN5xGCGpeAmrfU9fHTtqQqx+6YDU69Q5wxkOUuUIL3A11/JGljHM4QSphxF6SOppssZ9FFoXHzmXxhKeTYD9kPV9SafGctyAumcMVr4AXd+Cw9vsBZOXYsxwyeK/isjvAF/GM5cax4YAmBhiT/MYfms4SERuA24DWLx48Vm/yaoPf5jHRTTu2rQpYU1aAFhuSf0FoAa5cmwnMGzf7i4CknKXwiMHNA3gEFr2Rttr6VG7Q6uF1BALU8a1kzSmvZ8frbcidcfCzg8mzZYXktF1SNpe7t21qau67454+0bgJluUtDFni43c/uHU/O0RqVfiN0rzv3ca90fwmpnE6Ze11CP1yxnOsJP/l8M8z2uMMpNZzGM+y7kRQcYl9TS9+ncsoT/B5STqFBypV4gj8DnPW+WL/ZBu0bpMopq09yT0/Ezz6V0/gJO3WNtd20Dmkm4+PdEwAq91jT9ucuVr3m/v/8DbZlCzqbaYCGLvyGPYmtU/ALB69erXVS53o8039q9YEW9cuZKhTZso3bwrVsSALc/Th8WGVp02+ZkkCKoOrw7FEW5T9G6Rr9lv7cpkOqJOvJDqonfojOBDUk9bJA1TL229YCxLuwXTdqTu7x6mTtyhbNog0QTDG+f7qLtOVYnjDsdjo2PW0ZNpoeK9aG+O1P0GGUN4pH4cF6mHxUd1GhykGpE6wGuMcpRh5nMNPVx91qTup1/akrprPu0i9fqeuKK0QtQhqRetJo1I/Wngdv2ImZd6BgdjTEdVpmmYCGI/BImE5ELgyAQctzWClnhs2cLAli2U1q9nyRpgq24emw31UpyKaUVMcYqh3oLUgxx7vgaXQaRr9zM2bpGzQrpJWFpaJiR1X5/ubgm5HBqlj6Wlh+xlv2/udXkwJEyPpAX6eK/ZW7XQXGR0N1pV+sEc0OVJG8M8fFrUnlC9lNrb7iZIvblBhpK6vsEwhyNSdzjDGDVeoYero21pqRd3n0bqkUvj7J6481GFJKm7Sl9X3eunX+y/pTsJ5uvWpfFueOYYTQbBTtKbYSLQ1TpYS+DA+Z5IxxCRy1AjscXGmNtskdIyY8zXxtt3IiwFHgWuF5HrRGQWelX+1Qk4bir677mn5Wsjvb0c3k/i3ybqd4pX+t4OjhAjpJC6uw9bzhWCm4vgIRnJhzf3eoXxST13Mm6YET12rDmSlDY6cy9/fiXv3t2gmXi99Uv3etV7GxeFrx2Lv9/qCDBAvCCalh1y7+VSMC4f5ppjvA5Sv8rm1UEXSedTZiZJI35BKDInet4Jqbv0SyKnnkbqV9pF0qhP6UhyodQj9Y32O6u8At0Pwcl3wcH9ev3R2LQpuX6U4VLH54BR4Bfs80PES1ltcc4RuzGmISK3o1VROeCzxpgnz/W4aegXAdf6q4MfQO5lKBWgUALmwEabYtudFpH6qO/T+8La+EWXX3f+K47oWy0O1okrIv0ovEJzr1R/kdQhTL/4uf1EJySPQRPe6vV4Hj5Jhycf957Dwdi0E6Af4dvHx23apakIKUy7+PufWg7FLURnFmfk5adfWpD6Cl4FSJB6KGl8I8uYx3yOc5QGrzGDHHMpMY+rgNapF3d/1CP0OKde7kDOaL+Y2o44UnffQ0HXIm4fg5VHoPtOeHKn/Upa/E1n0frkg4gsQnuPltGypweMMZ8OxswFdgGLUY7dboz5nH3t/wbejQbVDwG/b21+07DEGHOziLwXwBhzSkTSUt9NmBAduzHm68DXJ+JYrfCMiKq9q1XqfX00tm9vT+62uKPrZX1anmUlj7lA8ugQErOLevPL9XnupP6Aw05KafCPFZK2n4uH1qoXR+qhTt29v9NFu4P6VxppLdV8wk6pdmyZYw8/U5BeqRaCnLqDH7G7fZyHenFtPAl/obSDSB1oyqlHX419PIMZrOPXOMLzDPE0ReYwj6sQpG0+HYhIfR+B78v1heZWdr6c0X3Q0E/dft8bC7H/S0Tqbf6GM1KfYMisiUrFNIAPG2O+LyKXA/0i8lAg7/5d4IfGmF+1HZcOisgXgNXA9ZLQQQAAIABJREFU24F/Ysd9G3gXLdrjAaPWj90AiMgSmluvpWLKVJ4uMYaT9mQ1DNT6+vRn3Bs7fheqVRbMJ+7oXoibWjszqkQpvCM5P0eeCst6OZpz26NXNg8vBPdpCFMu4b5hmscRfCSh8xgzkTpqg/BzuueQ3jUvbX8H+73dhfWNcdF6GNXXiQm9q+IVHRXU86VRVMtdd2JpYeYFMaE7nbqj9DRDrxk2y1jh+mjKaRr1tHx6TOplWsoZE17qA8lF6gqR/0svsOs0LH1BHRo5rk2nGwB9fazKuh5NKRhjXgRetI9fFZGn0L8Mn9gNcLmNrovAS+g/ucGGbKjoZCbwkzZvtw11vFpkTwxvB7Z0Ms8pQ+xg29+9T6jY4ruc/VE4v5gF84FPoinbK5TQAWpzdAE14TDo4JM7aLTrStqjhdSA/fwIdVYRXD7XT2X47xO+n0+oIbH7hO7n8xN9S0fak/lrheaoPZxLGnH78x4OtqWNQ8fttvdtST1qYGK/cGcP4Ma1qCZ1Do1Ay2rS0HY3tArwVTCtNOrN+fRycz7d/a2EpF7boYTuWtnZvPpXgLf9FHqetzLcx4D7dB43ZoR+YSGdLp7SIyKPec8fsIq+5kOKVNBuSt8NXroPXWc8goYCNxtjzgD7ReR/oCcGAe4zxjzVaiLGmIdE5PvoX6CgaZvxuvgAU4zYnxGhgDX78iwzyig3jB2D3IPArZDrVUWMI/elNjVTzXk2sw4+uYdFS00qmYDp8vbAs4CCJfhZo6q/bpWn9tMgYfVo07HpjNT9bfkKvDbUvHjq3j8N4ckozUogjPTr49xakfrJJRqlh/YAKR2PWilf2tkEQLo9QCdyxgSpu3x61JvUuTPaqya/mtQtRJdjnXrlNNq5y2EQDurFR1ZROnlx3BizerxBIlIE/gbYaox5JXj5l1EZwT9D26M9JCLfAq4C3kwsmXhIRN5pjPnHVu9je6Wedd/TKUPs/Q8/DJs2kd+1C9DuMYWdGqm7DEABKO23AfWtUJgNxxfHP66lNW1qHZE7JHPGTgEyEkTtjtzTCpd8Qk58mzYX7hO8v+uctL6kKUiQepBLD+GfjNIQnsz8e39+nm7doWlx1H/s9vUzRH5TjLCi1JF6nRaVpOmkHtoDjDLKsxxkkKc4w5koj54G39ArdGncR6E1qTvTtASpD6uc0ZF6hSin7uSfxYaaeeXrGlzkhuHwTqht357+b5NhykBEZqKk/gVjzN+mDPlt4JN2UXRQRJ5DE8TvAh4xxtTscf4b+pfWkthfL6YEsYde7DXgODr5UtoOAA9C7hBcvUHNwHzcjSpjdruoOVURUteoF4KI3eYpcvadwzZ0oSzSN/hyCrzxyLzrx/ZBmrQkQBqZu0g9bUHW3x6+7oqRLKp45B5G7GEePkqpFAIjL0vqTs74M8/2+DXaGnmF+XQ/BTPKKF/iQQ7xY84whjCDuZRYydsS5B46NIaLpKkadZ/US3g+6kHTaV+jXtK1BlBSX3oUuvcT1zPtsx6eAwOZkdfFQMeLp+McRvPmnwGeMsb8hxbDngfWAd8SkauBZcCzwHXAvxaRf4+mVt4F3NviGOeESU/sIalH1aciNLZvZ3hggOKuXRTR32DObyG6BxiEubcrubtCpcppuN0qZHYXoOrSMH70HvJolGuvxxF8jpjYnQTRl0W654njeFG6P869luh65CYTzoNmvfrMYHgaebcjdQc/Ii94C81hbt5PX6UttrYjdd8eIEq9xJ4v7dwZ3eNneIrneQ5X5Gw4wwgvcZyfMN9aJrYj9WSUHihffpHmoqMwUg9aCm4krpVYehS6Pwcnt+nXlNsMDOq7LMtIfarj7cBm4AkRGbDb/giVNmKMuR9thrFDRJ5ACfwjxpjjIvLXaHrmCfQPd48x5r+2ezMR+UXgemPM56zCpmiMea7dPjBJiT0k8wheufUqY7RYaeVK2LWrmdQdBoEhmFuH8o3NhmBN0sfE4qaN2n0jLfBMwwq2U5B3UF9jnvCVGQ7053idmdxz2wMzgXr8nuEcfL36eEROsL0dRlK2tTpWylQ1r76ShN3u6Kykj3qC1JNSxrXUbdFRsjGGn0//AY+Q1v7xGC/y87ztLEm9GLchfCtx6qU46BWAjcSk7iL1ivYm3UicfnGFRye2eevPO5X/F2QLphcRMxj/D3h8GGO+TbqNij/mCPDPU7aPAf+m0/cSkW2oRHIZWqw0E9XHv328fSclsa9at64luftl1qs+/GHd6O59fCL47gfh6hrU3qXSx/vStNcO/r//a3hEXgny1y78DFYaQwKP7HRT3iCsHIV0Anck7nLX4XxDQk/bHi5+kvJaK7RKvbhtI3iLpXah1GnUW5J6c+rF16enkXq8QJr+27qMYmKRFOKF0nEtd1365apRmPs4iRRYfU+cfqkAFe1NuhG9AizXoXTUWgTcYbN7mzZFX9XIwEC2YJrhbPEbqOrm+6AnDKudHxeTktihc3JviY+qNJK1wNfQyP1WKL8E93rONlEOOVTFRKirOZiPUIGS6PlmkavHmvdEMVEYWnurjmmE7h5f7rPyAX3o/onbRet+uiSN5P3nfgqqlaInvDmk+ak7jbpTv4zTwi5cJIV4oTTsTdrdZIADIPw8b2urfGnSqKeRenGQaOHF6dNnH1A5o43U7y6pm2W5rvbQ3S+ieas74MljtC0+ypChQ4waY4yIuAKlNt7hSUxaYoc25D4wQL83Jg0nReheA9yEkvtSoKwKhV60sCZCnTjoHqY5GxKSexOxO4YbpollcyELtrgcdFYAfg7dv/dJ1J0/0uSIhZSbP8bfzz9XuHyx49w04nfnIKeTf807eETq9kAuBVMnKWn0curOmRFIkHrYnzRNztjD1ZSYxyuMcIYxZpBjIRWW0ZuqfIlJ3S6SOlJ/D0kjr+KgbYxhF0jlAFxBouPRRkvqS49aQh8CvqaFR4eAhvVSz6pHJxOEiUjFXGD8lYj8BVASkX8NfAD4z53sOKmJHZrJPb9+PQWgZj3Z+++5J07JACdEIp6KUETFRj36tGdMF06dhK/XXlJXC4GPOOjfwghQr2tKN9S3p0kLo1RJsOIYEXfJ214nsndNy+M7ODJNI2wHv3o0baE0JGt/PcG3GfAuIhLHfrXFDyNyaSzB2AISvi8jJO0BUuSMjtDrLQgdmouOBGElb0UQXuI4S1nOm7iBWcxKTb00KV/8wqMyzUVHvoe6Kzoq6d/KHadh+fPWbvcxYNAjdStnzEg9w7nCGLNdRH4JeAXNs3/UGPNQJ/tOemKH+EfS//DDNPbupb5+PYVdu5R3BgZ4sq+PG4xhzFagFoDu+cAGYrWhrSUrXAHFOdoM4i43OEQHATbQeSm/I/qI8EeSufOQ0EOVS51kdJyGkJxDpOXQ/Si/Xf7c3999BreoHPYmDUl9iCYfdb8xxtlE6SEWoTrWt/COpkrSdFIvJ1vYVXRT3BjDntEcqXtFR47U7x/TSL3rW2gnyg3AUv0K8wB9fVlV6aTExCyeXgT8CDDGmL0icpmIXG6MeXW8naYEsTu46L2xdy+NrVsp7NpFHv0QJ21VahSpb0CjdKek26ObuypQvkozNHSp3HEjyWIcIEmQJeIIvpUWNm2B1d/WSqboHCP89/Xv3eN6PZlqSUPaa2kplfCxi9BbRfZpJ4V8hWTno1Jz8dEQTeqXFbwaSBnjPLpDJ6TuEBYetW5hV04aec0kIHU7aZcSc/YATs5Y8pwZ98PYLfrxSvth3nwNJJYfi20uIPNTz3BusOmX29DGD0tQT5r7UY18W0wpYgdL7ps3Q7VKATUHczjoWuehlanlnVZDvBbNsS9UUzBQadrGBtyUhyG/Y5a/4BimPYbtoqUjOrN8/KjdkXhhbfw8DWnRdihn9zM4/knHn1M4Lg3tovdWj/0rhqZIvRDn1M+K1JOLpM3Kl9Z4C+/QqbSJ1FM9X1L7krqU2D4iI6+A1O84baWMLwL36seLCm2PxRm7JXZ+LSW7GTJ0jt9Fe0p/F8AY87SIXNXJjlOO2AEKu3Y1kTrAMmOU9FeuZHjHDkaqVUo7YYFVxGD9Y0DVMXlLXuWrYNBaDegbkCTIENFC4oF0+SEQtURz0Xh9Xzq5u/fxid3PfQ97Y9x2R7qByrIpHTPi3VqRfNpCakjsCXuAihel2wO4AiRnu+sqSiMMRwul45l4+c2mfTjvl4VUUjXq43q+uNRLWHjkN5t2V1HuOyyrTv0Of6G0Cif223Pu3r3Ud+zguP17rPX2xrUVGSYZpmQq5rQxZtRZsItInrTijRRMSWK/oU0Oc9XOnVGXpQZqPbBgECjD6St0TPEVVcfk08guLZ0RSgf9bfU6cMAjxKF4Pz/FMrMOY3tgdvAePhk7+PNyi3edwM+Zp20rBNta7e8/byL15d7ELanXliZJ3Tf1ChZKW5F6K0KHdFJ3cI0xgBQf9WLcpapC/D3na3ahdLhZ/eJOssvh7gJsPKlBQPdzOrypteGWLTRWrqS2Y0fWhDrDROObIvJHwGy7iPo7QNtKVYcpSezjIW+92p1C4URfX1QcUhiB3Ckl+XoJ6rPgkblBL9Q0NUkYTfvwo1yXC3ek6FeGlrx7/7389xvxjlWimdTrwRgfoQbfjUtT0oTv5W9vlVOPrkwCUh9t1yCjFqhfFI7Uv8e3gNb59JDUE1Min+L5kuKjXiH+7iP1y1DL1AvL4dtj0HsMisPWcredvc/KlRmpT3aYGS16BE9q3AHcgloQ/Bu0mdGDnew4LYm9QPrvMF/XdnnUIG9z7YPFWP4IyYbXzvtjt1s89Qk4jQxHvMd1SGXI8KQREq5/VeCTej0YE5K4H6m3Imo/xePuQ/IPx6aioE0yfKuAJlLX3LqL1kP1y/M81zZKh3Tr3TSbgHFJvUyKkddQXPhVJkq9OB/13pdg7nNEGnUG9XbimG7CatV9ZKSeYaIgIjngL40xm+hQu+5jWhK7UzjW+voAVS3wKwaekEglk3sZ8iVdRAWtIKzkodplXQ0L6gLZMwbkPCdIaM51+0TsEyfB661SPK22h2hF7i66D4k9TMkMB887WTBt8qIvp/u/QJx+iWSNx1hL3XNoVD/18VIvRIcbampl5xO7y6kDyfTL9QX4l7ToSzqi6ZdoUZuErHFjAXpPWlKvAg/Cyf3eUoWtnwiRkXqGiYQxZkxE5ovILGPM6Ph7JDFliN1XGXTyIyrifbijmpPPrTDwKVHp432QvyYmdne/scv28URJvdiwfVIL1gUyRBjthvfQnAJptzgbjnXvkTbGRdytThDheoBjp7Rjum2JvLrdHi32ltJJ3eXVH6GJ1P28+tmQuo+0yD1MwSQkjW8jbo4RpV7sJZUjdd/3xUbqG1FJ48IjQBVO3qLnwgb2CjDFJiAj9KmCGXTUr3hyYQj4nyLyVSAyoGpjFxxhyhB7O++YEAuM4bDVtYO1F7AFTIPoT3zhzTDv76GnBMcvU2Mw1xf1brTd276ctfbFpmUKcZqm6ogSmgnc3YdyxVbpFx9+1NxqsTNcCPXHh1cRPkJljT//Vnl1dx91QAr2aUHqLq9+lY3W25F6q85H/uNQ1ugI/QkuJyJ159C4wE41agbukbpbKF1E7NBYijseVV6Brh/AYUvqUU/dLIeeARCRRcDn0b+eM2jrvE8HY5ajboz/G3CnMWa791oJzZP3ogqXDxhj9rd4uyP2NgNSzZFaYsoQO7w+co9wlTCEStK4916G1q+n9G6Y+10o2IXU4qx4+MYuJXRH6hvtdtcMuxeN4Ks+QTqEKZu0KLwd0o7p9g317n5OPkyvpF0V+KQeppRaITL2shj1FDDfJZXU11JPdWkMEUbtLv3SSt7o9ydtIvW3kWyQURi2xUfDSTMvG6E735e7T0LPz3RhvetFOHkzDNu/Ex8ZqWdAL+A+bIz5vnVa7BeRh4wxfjPrl4DfA/5Fyv6fRn3Y3yMis4DLwgEistMYsxkYCU8aneKciF1EfhP4GNrH7y3GmMfa73HuOJsfVxO5++jtZbBaZdlW6NoAXcthbhlK11iSvww2dttGzXiRuneIu4GemfDBMlRd/jok01YRdDuCd2Xsfv7aLzoKc+3+9pD03ethfj3tcRp8Y6+xbhLGXs6tMWppFypgmvuSgpJ3u/6kPrGHmvUjrSL1t5FskOHMvJxGvb5P29g5Ui+pnPHWn8LVjxNXKH/NflXValwM55AR+9SFyU1IKsYY8yLajBpjzKsi8hR6jfhDb8xR4KiIvNvfV0TmAO8Etthxo0Ba/nyViLwB+ICIfJ7Ao9oY89J48zzXiL0K/B/AX5zjcc4rCqB+Hse8jStX0qhWObxfy8K75+uYrq1qO1C4QtM0vZfB7u5mpYzz4C6Mwq458Mky7K5rusa3J3Dpm2ohkFSmwB3bXR3sLsUnluiKoQBVlysveemhMLXiqibdscrEro1pkX3LVIw9UzhS9429Iq26wterQ7LjEcSE/npJ3U/BJBwafVK/0jPzYkQbY7jUS4VEb9KtL9j2dX5gPqiHccVvq3bu5EkXHGTdjzJ4EJEK6pf+3Q53eSNwDPiciNwI9AO/b4wJGjhwP7oS+EY7xid2Y7e3xTkRuzHmKQBpFRVPJmyAgp/JGhgA4gD35DEl/tyt8ZC5T8HcYdi6BgauheM5q5IhmZN3947UHZlWTuviqxt7PGju4dI6Do7Y0xClf9wGS8ZNTbkdbKolet3XrIdj/edhQVJaaztXWTpOtB4WIIWk3g6tPGAS/i9+f1JH6oVhIodGv4WdJXSs2unO7wFbVfESGcYBDOpXdVyEZZbc2xXEZZgiMB0vnvaIiJ95eMAY80A4SESKaEPrrcaYVzqcRR7Nu3/IGPNdEfk0qlX/d4mpGvMfgf8oIn9ujPm3HR676Y0uCETkNtTQhsWLF1+otwUgNx9dVAtQwEohgRHr9zFvHxq5f1F7Vo4ACzbD22+HZ26I1TO1fNxmz3Vj2kisplla0wpX16W+UYC6l8Ov5VVe2ev71LTCiI3U/QVb7zWgWa5oiTzq59pJMZJP6oW1xALv4D1ctB6ReroPTJiCSYMvaWwHZ+wVp2CKsaHX5Wj6pTBsvdQH4hZ2fuFRQdMvdzwBY2/VC5gSRKopPiHqKXQstgdYldadK8N0xnFjzOp2A0RkJkrqXzDG/O1ZHPsQcMgY4yL8v0aJPRWvl9ShA2IXkb1guwMncacx5u86fSN71nsAYPXq1RcsBDoOcAwWHGqR1l4KDHqp6aUkBkYPj6sMruE2zFGSH+pSQneRPMRRehgfFEaT5O6jZ8x242noCQGg51pV5dzlNOutIm1/+0jwOI34076IJlK3TOja27lox4/W25C6U8FAegrGwV8kdfdhtO5LGxOyRl+rHvWUrZOoJnXKl4L+O936U/037JqvV2jd873JfNTAR2MjrwwZQoimJz4DPNWJ7NCHMWZYRF4QkWXGmIOoS+MPx9vv9WBcYjfGrD8fb3yh0Ni+nXpfH3zU0P1RYMUKfWHlSuo2x75gPsxbg5L6anWAbLwXupdD9wHdRk0bK3QVYcx6zpRRAi/XVVUBStxRimYOFCyJFizB1q0323DByilJRvjdL2LPRrC6DuVrgW6P3EOkLYT6uvq03Hl4QijYJ68NeYPsxB2ph7n1oLIUSE2/tMN4pF6gkNL9qJzUqhdQrboz9GociL3ULanfXdBm025NpF6C/MMw71FgH/BZgQ9k6ZZpDSOq5jp3vB3YDDwhIgN22x8BiwGMMfeLSBntADEHOCMiW4GfsymbDwFfsIqYZ4HfnohJhZhScsfXhZUrky6QnoTNmYSVjkH3BtSkvabeIIl0+G60pHyp3nLLobuiLy1FyWLuc/q8cAUULSemmYyVjlrVTQNuClSAeUfEtfh5saGkdHcuhdxDnbsvrQwll+3y6q3govW03LqXVwciHxhH7OPl09NIHZJ6daA1qS8nvaq0vk/NvCpEpL7Ryhnz9fjfpF6C/DugaylK7p8QjdgzZGgDY8y3adVJPR4zTGryF4wxA2ioeF5xrnLH3wD+E9qW4O9FZMAY88sTMrOJgtcfNcLWreSrVRq2i/zwrl0s2YMS+wEiX5Cx/UEGY7/yZG4NcCs4S6Hu51ByQGWTXS417UjaRvkNL3pf6hFrbY5G8LVrVYlTtFcExxfHefxe1MPkLqedb5dWgeaK2HDhNLx3qZjiFqI0jIvWW+TWL+cUb2AUww94iGe5mmupcH2LCcVoReppenUldbdYWkiaeiVIfVjz6rMPNJF6+aX4iil3Su+7QM/qu+GkFbt03ye6iPr5jOCnHc7QWTAzTXCuqpgvA1+eoLmcF5T6+hhxxSYDA9DXF68Z2vLw+q5dPHMMKjfH+6UJSJyKsHuQKHLtXg7cqx7dAPM+jqaBnS7aImdvp6/R5340X7Lp4JGrYGiOaughTun46IVY7eKT+0jwPCyKSlPP1IkbaPh9S/3c+uisltH6GxhlgE9xlGcZ5TSzmMW1LOY3+UDzxGldVQrNpJ5UwJRjBcwCmkm98UhT8VGvTb8UG/F3nTsVfP5hYI93EXQMuvekTj1DhimFaZ+KSat5LKAfvNbXB729+hiVEhZp7nkBRJYEgOZfAPaocuYQMYeXtkHuS2jU2ENMJPZk0LUaxhZ5i7AWXS/C1U9DzyKovkm3Fb3JH88l5ZEbiaWMTZbD7orB162Hi6/R+9e97khOs+7JG91+KdH6ZXyfH1lSBxhllMM8z2P8TzpBq0jdb22XqoApoDn13EmiSH1+PVmZGyAidVeoVUTz6ncEV9VHs2g9w9THJUnsEJuE1atVtRmwjn0j1hHSKeRy84lNxIwni9ujUboW0qMGUcBwXx8L7gXuVc/3yCp4H0ruCyFXhIYXuedeIMrj57bC8gIcutbO036Anrx2eNrY4vNEskZbtBRtGyEmeIcwD089LkYKNet1UtrcaW5deC4idYfXGGWIp9umZFoVIbnFUkfqqYulLlrPqyKHkY8lZY1losIspzCK1i6G0VTbctTtE/Tf1p6wuzO9+vRFloqZXvCLS/q3bk18YBeRj1SrNGxzDj9Sr4Nenr9PNO/6PVEmt8siJWCetS1wJ4QFa+zB74Oum4ijZ7swy9I4HQMaued6vNcrmpLx0zDu8e1jqqSB5sKmqovUR2wO3n2AEkpkQ8TRapiOMcutyVepOa+ekoJxSpgrWMAP6EqQe56ZFJlDKzhSd00z0jxgUkm9Yj9fIgVzIPZT94j9K6j17sIj9sR5yn72r8HYzmTDaVDriQwZphOmPbH7WPXEExz0qmQX2IrDBd5i2QnrCumyFwCVnZBbK8qgNjPArZD7ju63wBgW2CjeeXcX9sO8PcDtaDPtStxI2+V8GwUlca6C4jW6PST14YJGnpHGvaDaed/iIKpKjc5URGmJjejj3RW7w5C9d/n1wloobCDyWXd5dTemKQUTOzeW6WUBiznEj3mNUUvqlyMIx/kJBoN4AgKf1MN+pU6rnvCACUk9aj7t2QUsIiL03hLcb6WjPc8HOfUDSurDaJo+wyUGQxaxT2csM4Z+Z+l7DEo7YWSnRFGbC2obxGmcIWDJHcDtMLZNx5SBXK/AW2wkP6ipmaa/nYXxQXIuDWBPDrkaXP0DoAI/WQbMabYqOJ6L0zGFUSjmNS3j7IQT5mQpksaE5YCvhIlIfS2RwRfEhUgtonXAsw04wyb+LQepcoTnOUiVlznO/+JRZpBjLiVu5K0RubdqRH3Uql9SfdUrxOkX59ToGk/Pr8dXRJbUVx6xtQBDxAsf5fgjZchwKeCSI3Yg0V1+BKj39jLs7IBTWp6NrF/P4WPAttjOdXj9elbdBHxN4D44sTOFODYQhdQnb9E4swzkNgODcXRfBq7+e/jJ25O+Mo60e2wRFNiovRG38qtCelUquq3qtg8TV6JGi6XQtNKYr8ErV7awDiDSrTvMYAbX8SbGGOMRvhmZfp1hjJ8ywgmOspK3Au39XxIeMCGpzxmF7idV/VLboX1kXeNpS+rfHoPVz2oRGcPEnSE3oIukvwJsmwKeRhkyTAAuTWKHqLt8Y2AgtdVZAnv3Mrx1q8ojg7FjsyG3tHmXEqh65j500TTEUsBKJHO2rL34SuzdULZeMmmo5ZsNxIBkgZJ77qti/OKmppZ3FqnWAYrLURF4mnvjTzhCg9cS0znDGDNsqVdnpF6G63tiUi+gpD73caJIfWawUFqBrxRg5QuW1B9DffEAvpPMnbt/hiwVcwnCQPDnOa1xSRJ71LBj5cqmVmfhOLBt+YKmCwAcNUpbL4taEvgEvkHTNi7d3W29SXJr0Jz7IeheYx0nN+guhZE4/14iaR42bPPntbwuoLq8uhuc2vAD0kndIbLlpbV1AOAL8n1Sjw5Dnqu5lpnM4jXPXjrPTJayPJXUv2PfN/aAsVp1vwPSVdZX3ZG6swoISH39C7ZIrA2p999zT6RcypBhuuOSJHbooBvT1q3wxBMA5NevpxGmaHp748fvNJzcbyVzn7HbHgyqVh2p29ciewLs42FbxNQTL7KCJfqS5taHC5qiCa19qxDr130Cb0fqEEft7awDTrlPUEs7AnnyNGjwJm7gWhZxhBd4jVFmMotrWcQbWdYyUgdiBYxvweuMvYqDtqq0rjl1Idaql7Wy1DdfixCSujuJZ7h0MYoGK5cILklijwjdVqICrLKLpz7ZPyOS2m+aHTvIV6s8I8KSNXB4v829A5Vbqsz7ErBV/WfCSkZnU8B+jdhZikaph1Ct+1qb2nGBdMF6xsxKNt6uAHTFKZmNKOHvLtlGHP5ZJQ1RKsazDugAQ+Sp0IgWP13kPoMZ/Ba38CwHOcoRruJa3sQNzEIVNm3TL2mkPud5W4A04jXTJiL2XvsPczxnzdaugNwGUv1eEldeGTJcArgkiR2AgQHyfX0UIWES1kkkX6hqdrsBHNwPtU2botz78Pr1lG6G3Mdpzq0PJmWU9f3xOCd9AAAZkElEQVQwbzD5OktRYUgl9jUBdYks2n+taAG1K5mScY06IruBVsTuk7ojdtqMD6L1IfJcm1L6NYtZLGcFy1mRUL24+2anxmJrUu86qL7qLgVzBfFVSSXpfV86aovA3tlej342PXMzZJjKuCSJfdW6ddq/slUThfXrIw7pNoZDTvs+MADVKjcYw+Mi1DdtasrTu0LN0jatZaqjPL1kPk2oA4ePwYI9aJ7dRe/lOB3jyL0wog4FzjAM4rZ5oAS/222L8j8kyTqSOVZINKjO15IRe4oaxkeojMkHf0YhqR/x0i6dkfrjMHKHFh9dQ0L94pN65TT0/ghyj8LYLTBI3PWoFbKG1JcoXHrxEsGMiz2ByYaDIpHdgCsxv9GP6MNtYe520yZqwFBvL7W9e2ns3cvI9u08fgxOHEsPik8eIxmto4TuSN3dOxtfiLs2+UiQuoNveANK6GEvUx9u/+BHEPYzdQgLjcJcuh+lJ+13e+A9aK9SR+qXvaSRen1PrFF3twoRqW9ESX358z6p65VT/8MPJ24ZMkw0ROSzInJURNq2MRaRfyoiYyLynmD7HBE5LCL3na85XpIReyv0r1iRWBTtX7GCVXYBNYryPazaubOZPKyM0hH+qnXr6Acae/cytH59pOLzg+mId10+3hUyVfSpi94bBa+/KnHEHhUnhUVITv5Y8ra9us9a8wYSx9DFMUJgU+nBLZy2i9DdfcuKUje/y7wm1H5OHRKNuR2pl+vQZT14Ihn/wIB+35s3634DA9ECeIZLHKNMZMS+AxUyf77VABHJAX8KfCPl5T8Gvjlhs0lBRuwWrSSN/Q8/3NHluz+mP9gW5Xb37mV4/XpN8ayx9r8Qu0UOEpP7HmArUEnvmery6bshvVG1I3X/LFIA6lZh0sobpgPlgCPsaz1Sh5jQ3ZimRdKw+MiRemQVMByTun9SsqR++1hsr1A6ClTVJsBJSpe5k/DOnTw5FRqsZ5iSMMb8o4hUxhn2IbQv6j/1N4rIKuBq9Bd+3hpuZKkYi1bkPR6pr1q3LknqLoLfsYMTllwOi8COHdGYEsCtqI+MJfXD++HJY96Bb0XtW66IlTE9z0PllaSdbwSPBCPUaSZ811DDh4vWfUT9X+Lce6w518dHgtuQl3ZpSeq+Tj0idWsVEMkaPb16RU297rBRetl+D11fhBM2BTMCTbn1G4yhDvRnBJ/h7NAjIo95t9vO9gAisgD4DeD+YPsM4B7gDyZmqq2RReweQtXE2S60hWmZOqjFb4DIXfDr1mNmp5VLbtnCM319LFkDp9+R7LiUewGowdxTULkOhvzFWBuVO0MwoDmKd3r2RI7dYtYoFFLIfXYBTmka5lVmczmnImIPMRRE65E9gLs5Um9aKP0xUa/SxhCM7Unk1e9GnRqLDduC8CngQThsI/W6Z7ncBFt70OlVV4ZpjM5TMceNMecaSd8LfMQYMybJwOJ3gK8bY16Q8xxwZMQe4JwIYGAgfrxyJcd37aK0Lfla3St0OvFuvS+AtuobGNCA+5Pqx97zM0vqLxOnues2eh8jaMzaAj7Bm+VetD6e0N2HawmVhB/BO8SEbvdzyhcXqScWSn+SnEd9nzZZtAVIGwva2q7yrF1AHgJuVyVRhgyTGKuB3Za8e4BfEZEGsAZ4h4j8DvojmSUiNWPMHRM9gYzYJxB5W+zkFDOP9/UxjBJ3sVqltmMHbNmikf3AQFTinu/r04zMrl10fxNOXqekXhwOSN2DbwSWQLh46t/yleZoHZILpz5OObEmiUk4MlcSBz9dEz/38ukL0CYZ/kJpYZjEiaW+L/ZWt7c7TltSfxR4MDZN85GvVrmxhWw1061niGCMV0V9vt/KXOcei8gO4GvGmK+gWUW3fQuw+nyQOmTEPqG4Mcjz3mibcIBSXWPXLuoDA7pI68kkXYOPnLUj6H5u/PcqNmBtHu7Owe5Ci0Fn+3c8M3g+u+D9GMq8Gilk0qpU/Sjdpl38RdICcKXffNrmh2o74IoD8GaivHpvGXY5jfrL9rCDmoX3z3F5mr/zEFkKJsNEQ0S+iDo+9YjIIWAb9tdjjLm/za4XDBmxn2e4tcxuY1iA6uRr69cD+uUX0PXT3BpiDxns/XJiu8cATva4dizwZnck7+vXHTc3hjxFjB00npXA7AJEFbDpKZnoNUfo7lYhtt2NfF9sbqgxBKfvh58n7oBU0Jx6lH55OT76yWM2WveMvBp+6itDhgsEY8x7z2Lslhbbd6CyyfOCcyJ2EfkU8Kvo0sQzwG8bY1pZTl2S6F5DwpRqmTE8Y6N4J3sEUu0HWGjvA+51hUqO3HssuTdVS/gyxwLw6lBS6hgWJ0G6p+0hP3K3k5ltzxqn6s2E7o7jbHejKH04PrmM7dET13LNpbvqgZtOa0u73KPe2x2ypO4qfTNkOGuM0ToomX4414j9IeAPjTENEflT4A+Bj5z7tKYRfKfB9wljO5sy3Np0w1alzoPI+RGIV/IXEgXM+brm4MOmHE1NN3xFTFpaxjf/8ie1wN6cORkouTss9F8rJLeF+fTCsKd8GUoqXyrq0HiTbZlay6ucs+sHWk0K9kpmUH+ShV27uGHnzngeWZolQ4ZUnBOxG2P+wXv6CFoknqEFTtqem06GXgdKQTu9k8eg2xH7PmKrAYjJHSX3CnD8MlTq6BZSQ5mju70Wnk48n5jRK1tPeqF3f4hkVO5H5xATurtFckY7ifo+mH1A+5TaBdK1Y3F3qEjOeLt+7DxQ3q9e9ktpbkKdIUPnOEMWsb8+fAD40gQeb9rhgLUrKFarNIBaby8j1Wqi6h9QVvMJ3ZmDeT4yecvL5Tr0FKB8GawtwAcrUB0mLsdsibqSOsCsIuA1sJ5JUiETErmDH537uf3L/EVSq1Gv7YBF9bhHaVmrSct1JfTCCHS9CNwBzxzT7wb73RSOZaSeIcPZYFxiF5G9pC/h3WmM+Ts75k40EP1Cm+PcBtwGsHjx4tc12SmPe+9NVqmuWKH9VwcGGLYWwhGR+XYDHqkDcBy6QK1s0ei9XIfiHLi/CB8se71OfXKP5I4Adet3jiV4S+4lYk1hux5yM4ldF2eNxsfJ16w+3VZE1ffpbVE9XkwtqJ/67WOerPMp4F54Zr+Nq2zR0UhfHyUgZTUgQ4YMLTAusRtj1rd7XUTeD9wErDOmdVhljHkAeABg9erVl2z45VdBOoOxfoC9e6lt3UreRqnd0Ezq7v4xYBDyfxD7yNTmePn2sPaogEbgoZWAj5Dc044R3s8ajcncRf9+lO681F3qxfm+FKxG/RWvAOtrSup19I/SKV4KaComQ4ZzQ7Z42jFEZAO6WPouY8zPJmZK0x9hiXtUSHPvvYxYL/goQrWGiBEG4cQ2JcAFS6HxXhi+Uj3a9+XgrjqR7UrSl91KDBPkXtd3ynt/8Hk0NRNaDDhCTyXzkyRWauv7dIF0NnF/Unu7G10sXXpUrzS6XiayVRhx9gA7dlDYtQtINkHJkCFDZzjXHPt9aFbgIVs++4gx5oPnPKtpinbFMqvWraP/nnvIo0mM2n5YNogahblg2vJolF1ZGDe/vi8Hu30FTJhfLwH1A6pJL2AJvqSk7JQxPsFT1LS72+aPKQwHZO7JcBype7n0KFpHZY3luo3UT9lde2GeMaoIgtYNUDJkeJ2YgaGbU+OOe/UCzOVC4FxVMUvHH5WhU+T7+mgADVuE82RfH5VtntYdYAN0f9xa/lrxd/kl2NjtWfg6Unf58kS+/YDq2YtbvEIlAlInjshdHt41vM7XkoVGEFvthiZeNlrv9U4yPWO6WBqR+jjt7DJkyHD2yCpPJxEKdJgFXI4WNAfoLUDVb6qR9gagYYkj4zx2JdZD7kTzvjm32OoqRw/Y99mnx5pp3zAgdVd8tHYM3t4P7IGxbcBm4PMZqWfIcD6QEfskgqtKHenrI08sOomwgaiz0tjseOF02MrQe4FqmVjVEqZjXKqmKc/urMoK3iAHf7XUpl2cxa4z7vK1685H3ZL6RtR2d8k98OQ2tdktUWVJRuoZLiC6ME1tHdMwXfptZcQ+ybDEsxwoATmnCLGNN+hRUq+XYiXMsEfgvdio3aVkfL8YV6jk/6tHC6ouMV+PK0TBk0h62+v7NEK37UsT72NJ/Ssooff8DOY+ahd8t2+HlStZklWMZshwXpER+ySCc4JsqitailZ9OrnjIuheZKj9VCK/GIj9VsCSu+8TE6EOrw0pQUdNN7zXwn6jjszdbWY99lQPFkZBo/S7T6rfS75upYy7m+12M2S4kJiVRewZzhfaeoNv3dpcgergipX2AN8xUX+Nq+caSk8Lj72xeZdeoOoOGNoLjNh8eH1fsmgpX0kuhLr7mXWVLoYpF0vsrv8qaNFR+SVbRWpdfk/spIOfVIYMGSYKGbFfZPga9vqOHRzftYuetIGDMLYfBkVooD09T4gw75tACrFXIY7YyzSbgdXd/w5AfSh9ci46L6TcvCbTjtTdfW1OrMM/cbNm8PMAfX3Q2wtPTJe4KMOlCFu/82m0h9mDxphPBq//P8D/bp9eBlxljCmJyErgz4E5aMXUnxhjzosNS0bsFxBOx+4i96gC1ZH7li00BgYYqVZTI/c6GvnWt2+n/557qHiv+da9VWheD00QOsH6aD3dVyYkc39bSVU4G4HK6eYG26evga63awapsX27Fh1Vq+SrTebCGTKcd3S6eDoeRCQH/BnwS+if96Mi8lVjzA/dGGPM/+WN/xDadQDgZ8D7jDFPi8i1QL+IfON8WJ1nxH4RkFaotGrdOvo3b267X/cadTs8ZFvwzVujbfQcqZbrcHsBqjmbY/fJeNg+DqtJ/TGl4PUSzTl6O24jmnZZeSQoNnJNlu7V3qQFoDYwANUqeVT5kyHDFMZbgEFjzLMAIrIb+HXghy3GvxftsIQx5kduozHmiIgcRSUIGbFPZ6yyXuNP2kXUbp8Ef0G3zVtDVKF5+i817eEcEkFJ/u65cJcl415sFB+SdlpqxVfQeNtciqVKnHpxtgDdT6DeNXtI+Mq7t3ItAQtk9gAZpgR6ROQx7/kD1ufKYQHwgvf8EPDWtAOJyBuA64D/nvLaW9Da7mfOecYpyIh9EsKpFY+LRD09D+63nuSfQR8UVHVSOqrSR4daPk7LJOCvzIZpmVL6670F7T3qrgh8BU65rr1Zx94NQ3abu9ANm+2VgAUZqWe4iJiJ4drOUjHHjTGr27wuKdta/XFvBP7aGDOWOIDINcBO4P3GmDOdTOpskRH7JEajt5f+e+6BgQGKQG4zsUd7XVducqfUSOv0FUmCd5F6Ao7AQ2K3UXtvIe6b2gt8xfUePaXaeYiLogojwHtgABJ9SAFq1oK4SEboGaYdDqF+pQ4LgSMtxm4EftffICJzgL8H7jLGPHJeZkhG7JMSC4xhePPmyJOcHTu0x8VaGHtXcmxuM2oU5qFch1sbMFhUx0fwUjJNuvZYrujr4KPeoy8ANchVYnJ3+vSxY6jKJexDai2Il2XqlwyTBAZDY2JEt48C14vIdcBhlLz/VThIRJahHRP2e9tmAV8GPm+M+f8mYjKtkBH7ZIUl9VXr1iXkgcO3CMO2E1OpGpfmj9hipWIDrj4IHIfS9VBcTOQF40fxVft8I5q6gTjlUmx4WnTr/c5NSu4RatqPtLi/Sm29WvavstG5k29myDDdYPs73w58A71o/qwx5kkR+QTwmDHmq3boe4HdQY+K3wLeCcwTkS122xZjzMBEzzMj9kmMNPVMFM2vXMkSz962+Aowx5K6Ze+uIlRKsPZKbcAByai8cjq58AoajTsrYGrEbfpqwHHIFVGFzdd0ew/6R+QvjEbyzQwZpiGMMV8Hvh5s+2jw/GMp++0Cdp3XyVlkxD5J0da73apnIvyC0H0rdPeiK5n77PayknSxEd98lF/SBVAACppqyZ3SfH3X9YaxupC7CSV11+xjELhD1S+hmKbT+WfIcKFhYKJSMVMCGbFPcvQ//DDs2EFx1y5VxQSLkQdF6MFKIG/1XhgEhqF7GG7sjd0gIV78zJ0i1rcXIVcDDkDXR/U9cisMPCq6XHRIjzm2U3epQeRAOS9rXZchw6RCRuxTCD6p969YAdUqRaxKcSkaqd9kb9bY/fAter9gDeRuBYqaoqEIY1dArkyM3errcnybsMwYDkqs7OohrnOqoYodbBXpvKOZ8iXD5MYELp5OCWTEPkXgV2z233OPLk6uX08Pau075rIzOyH3cWAPHN4Pw1aKONLXR2m/ZlRya4ANkFsNY4usA+Og7jMC1Hp79Uph71495sBAJGFsYEnd9iZdkilfMmSYdBBzEXTGq1evNo899tj4AzM0oV9EdeNOYmgbP7ty/YPWJMyZLzoVTA9x3XIZmOdbD9XQVMt98MwxbSq9KiPsDJMEItI/TtHQuLhGFplb2DruuD+h75zfazJgxsWeQIbOEaVGfN34li3UN22KIvplxlDftInhvXtxbVEdqdf37qXe28s8Y+C37M2RurUGLkBG6hmmITQVM95/0wVZKmaKoF9Ei4GcPnzHDr1fuRK2bKH/4YcTLpFRtepnYME+qLlUjacvPywSKVtc7hz0BJKZdWXIMHWRpWImOdL04E5KeNJG8Ad8wodIRbPMaxh9UoQDkEzjbN1KwS6A1l3e3L02MMAqTyefIcPFwkSkYubJVebd/Na443byZ9MiFXNOEbuI/DFqWXkGOIpWUbXyTcgwAfD14d3GwPuE0s4qI1u3Nld7eg2ju41hFfC4CI1Nm5TEbXMPILYvcAhtAjJkyDBlcK6pmE8ZY/4dgIj8HvBR4IPnPKsMEcYr9Dm803p5BXn3WgtivtFX17i0TYfvlSFDhqmBcyJ2Y8wr3tNuWttXZjhPGAHqLgL30UHE7Zf+Z6SeYTrjNUY5FBlMT3+c8+KpiPwJ8D7gp8R9/tLG3QbcBrB48eJzfdsMFj6pvx6PlozQM2SYfhh38VRE9hI7hfi40xjzd964PwQKxpht471ptng6cUiLuLMoPMN0wkQsnhZktlnMknHHPc2Tl8biqTFmfYfH+i+ogfy4xJ5h4tCqf2qGDBkuXZxTgZKIXO89/TVQRV2GDBkyTFeIyAYROSgigyJyR8rrXSLyJfv6d0Wk4r32h3b7QRH55fM1x3PNsX/Sdgo5A/yYTBGTIUOGSYjRCVo8FZEc8GfAL6E124+KyFeNMT/0ht0CvGyMWSoiG4E/BW4WkZ9De9vcAFwL7BWRN4U9UScC56qK+ZcTNZEMGTJkmAJ4CzBojHkWQER2o7U8PrH/OvAx+/ivgftEROz23caY08BzIjJoj7efCcZFsRTo7+8/LiI/9jb1AMcvxlzOAlNhjjA15jkV5ghTY55TYY5wbvN8w7m+ueHMN05xsqeDoQUR8ZUdDxhjHvCeLwBe8J4fAt4aHCMaY1vp/RRtmbAAeCTYd0GHH+GscFGI3RiTaM0gIo9N9pXoqTBHmBrznApzhKkxz6kwR7j48zTGbJigQ0nKtlBa2GpMJ/tOCDJ3xwwZMmToHIeARd7zhUBooxKNEZE8MBd4qcN9JwQZsWfIkCFD53gUuF5ErhORWehi6FeDMV8F3m8fvwf470YLhr4KbLSqmeuA64HvnY9JThbb3gfGH3LRMRXmCFNjnlNhjjA15jkV5ghTZ55tYXPmtwPfAHLAZ40xT4rIJ4DHjDFfBT4D7LSLoy+h5I8d91foQmsD+N3zoYiBi2TbmyFDhgwZzh+yVEyGDBkyTDNkxJ4hQ4YM0wyTgthF5I9F5H+JyICI/IOIXHux55QGEfmUiBywc/2yiJQu9pxCiMhvisiTInJGRCadDG68cuzJABH5rIgcFZHq+KMvDkRkkYj8DxF5yv57//7FnlMIESmIyPdE5HE7x49f7DldKpgUOXYRmeO83W3Djp8zxkw6ewIR+efoCndDRP4UwBjzkYs8rQRE5M2oxcNfAH3GmEljo2nLsX+EV44NvDcox77oEJF3oi1gP2+M6b3Y80mDiFwDXGOM+b6IXA70A/9iMn2Xttqy2xhTE5GZwLeB3zfGPDLOrhnOEZMiYp8qDTuMMf9gjHGtzB9BdaiTCsaYp4wxBy/2PFogKsc2xowCrhx7UsEY84+ommHSwhjzojHm+/bxq8BTnKcqxtcLo3A90mfa26T8bU83TApiB23YISIvAP8n2mJvsuMDwH+72JOYYkgrx55UZDQVYd0Dfx747sWdSTNEJCciA2hP5IeMMZNujtMRF4zYRWSviFRTbr8OYIy50xizCPgCcPuFmtfZztOOuRPVoX5hss5xkuKClVRfKhCRIvA3wNbgyndSwBgzZoxZiV7dvkVEJmVqa7rhghUoTZWGHePNU0TeD9wErDMXaYHiLL7LyYYLVlJ9KcDmrf8G+IIx5m8v9nzawRgzIiL7gA3ApF2Uni6YFKmYqdKwQ0Q2AB8Bfs0Y87OLPZ8piE7KsTN0ALsw+RngKWPMf7jY80mDiMx3yjERmQ2sZ5L+tqcbJosq5m+ARMMOY8zhizurZtgS4S7ghN30yGRT74jIbwD/CZgPjAADxpjz1qnlbCEivwLcS1yO/ScXeUpNEJEvAmtRq9mfANuMMZ+5qJMKICK/CHwLeAL93QD8kTHm6xdvVkmIyD8B/hL9t54B/JUx5hMXd1aXBiYFsWfIkCFDhonDpEjFZMiQIUOGiUNG7BkyZMgwzZARe4YMGTJMM2TEniFDhgzTDBmxZ8iQIcM0Q0bsGTJkyDDNkBF7hgwZMkwz/P+g1LObS+U4swAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plots free energy diagram\n",
    "fig, axes= mplt.plot_free_energy(np.vstack(Y)[:,0], np.vstack(Y)[:,1])\n",
    "cc_x = cl.clustercenters[:,0]\n",
    "cc_y = cl.clustercenters[:,1]\n",
    "plot(cc_x,cc_y, linewidth=0, marker='o', markersize=5, color='black')\n",
    "title(\"CatS TICA Cluster Centroids\" )\n",
    "\n",
    "fig.savefig(\"tica_FE.png\", dpi = 100)\n",
    "\n",
    "print(cc_x)\n",
    "print(cc_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-7.87387408e-02, -3.68212266e-02, -2.35290772e-03],\n",
       "       [-3.76031836e-02,  5.99394478e-04,  3.38061107e-02],\n",
       "       [-6.46345438e-03, -1.96376026e-03, -2.17023188e-03],\n",
       "       [-8.00493419e-01,  5.78973509e-01,  6.27073917e-02],\n",
       "       [ 8.10931749e-02,  6.79184138e-02,  2.74853845e-02],\n",
       "       [-1.11191129e-02,  1.14119456e-02,  3.37291436e-02],\n",
       "       [ 6.22111594e-01,  7.36973140e-01,  4.36155321e-02],\n",
       "       [ 1.47632245e-01, -2.81385849e-01, -1.03837994e-01],\n",
       "       [-9.18043855e-02, -9.94803241e-02,  3.66813226e-02],\n",
       "       [-1.95983013e-01, -1.36016416e-01,  1.30610116e-02],\n",
       "       [ 7.38214829e-02, -9.03526272e-03,  9.86100963e-01]])"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#This tells us how much each feature contributes to each TIC\n",
    "tica_obj.feature_TIC_correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.62313897, 0.86850045, 0.97029445, 0.98760017, 0.99341148,\n",
       "       0.99699391, 0.99849689, 0.99953202, 0.99997108, 0.99999433,\n",
       "       1.        ])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#This tells us the number of TICs TICA produces to retain the variance we want, which is .95\n",
    "tica_obj.cumvar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([0, 0, 0, ..., 8, 8, 5], dtype=int32), array([0, 0, 0, ..., 7, 7, 1], dtype=int32), array([0, 0, 0, ..., 6, 6, 6], dtype=int32), array([0, 0, 0, ..., 8, 2, 2], dtype=int32), array([0, 0, 0, ..., 6, 5, 9], dtype=int32)]\n"
     ]
    }
   ],
   "source": [
    "#This prints the discrete trajectories we saved out before, then prints the maximum of them...\n",
    "print(dtrajs)\n",
    "#dtraj=max(dtrajs)\n",
    "#print(dtraj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From here on it's clustering and extracting centroids, as well as saving out the centroid coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#makes dictionary for each cluster and frames in each cluster\n",
    "Clusters={}\n",
    "\n",
    "for count in range(0,10):\n",
    "    frames=[]\n",
    "    indices = [d for d, x in enumerate(dtraj) if x == count]\n",
    "    frames.append(indices)\n",
    "    Clusters[count]=frames\n",
    "#print(Clusters)\n",
    "\n",
    "traj=md.load(\"/net/jam-amaro-shared/bccgc4/Strided_Traj/joined_traj_all.xtc\", top=\"/net/jam-amaro-shared/bccgc4/Strided_Traj/system.prmtop\")\n",
    "\n",
    "#makes .xtc's out of each cluster\n",
    "for key in Clusters:\n",
    "    print(key)\n",
    "    a=traj[Clusters[key]]\n",
    "    print(a)\n",
    "    a.save_xtc('cluster_trajs/Trajcluster_all'+str(key)+\".xtc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#saving cluster centroids as pdbs for 10 clusters\n",
    "avg = [100, 100, 100, 100, 100, 100, 100, 100, 100, 100]\n",
    "indices = np.zeros(10)\n",
    "for i in range(len(Y)):\n",
    "    c = cl.clustercenters[cl.dtrajs[0][i]]\n",
    "    \n",
    "    v = Y[i]\n",
    "    \n",
    "    newavg = np.linalg.norm(c-v)\n",
    "   \n",
    "    if avg[cl.dtrajs[0][i]] > newavg:\n",
    "        avg[cl.dtrajs[0][i]] = newavg\n",
    "        indices[cl.dtrajs[0][i]] = i\n",
    "        \n",
    "print(avg)\n",
    "print(indices)\n",
    "\n",
    "traj=md.load(\"/net/jam-amaro-shared/bccgc4/Strided_Traj/joined_traj_all.xtc\", top=\"/net/jam-amaro-shared/bccgc4/Strided_Traj/system.prmtop\")\n",
    "\n",
    "for i in range(len(avg)):\n",
    "    traj[i].save_pdb('rep_struct/testing_pdb_res%s.pdb' % i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "\n",
    "path=\"/home/jegan/Documents/TICA/*.xtc\"\n",
    "\n",
    "files=glob.glob(path)\n",
    "\n",
    "print(files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "#writing a string to select the features in vmd\n",
    "poc = [[146,187],[116,136],[142,143],[71,212],[27,29],[117,187],[64,160],[119,20],[216,191],[24,146],[114,119]]\n",
    "\n",
    "pock = []\n",
    "for i in poc:\n",
    "    for k in i:\n",
    "        phrase = \"resid \"+str(k)+\" or \"\n",
    "        pock.append(phrase)\n",
    "\n",
    "f= open(\"vmd_poc_string.txt\", \"w+\")\n",
    "\n",
    "for y in pock:\n",
    "    f.write(y)\n",
    "    \n",
    "f.close()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}