StudyC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "from sklearn import linear_model\n",
    "from sklearn.feature_selection import chi2\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import PyGnuplot as gp\n",
    "import glob\n",
    "import os\n",
    "import random\n",
    "import time \n",
    "sigma=3.504              #Parameters for LJ of Argon\n",
    "kbT=0.010142612711"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Descriptor functions ###\n",
    "\n",
    "# Lennard-Jones function\n",
    "def LJ(r,sigma):\n",
    "    return 4.0*((r/sigma)**(-12.0) - (r/sigma)**(-6.0))\n",
    "\n",
    "# Gaussian Function\n",
    "rcut=9\n",
    "def gaussian(r,p,q):\n",
    "    return np.exp(-p*(r-q)*(r-q))*0.5*(1+np.cos(np.pi*(r/rcut)))\n",
    "\n",
    "list_Gauss=[]\n",
    "for a in np.arange(0.5,2.5+0.001,0.25):\n",
    "    for b in np.arange(2,10+0.001,0.5):\n",
    "        list_Gauss=np.concatenate((list_Gauss,[a,b]))\n",
    "list_Gauss=np.reshape(list_Gauss,(int(np.size(list_Gauss)/2),2))\n",
    "\n",
    "list_LJ=[]\n",
    "for a in np.arange(0.7,2+0.00001,0.02):\n",
    "        list_LJ=np.concatenate((list_LJ,[a,1]))\n",
    "list_LJ=np.reshape(list_LJ,(int(np.size(list_LJ)/2),2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Read functions  #####\n",
    "def readFile(input_file):\n",
    "    Y=np.genfromtxt(input_file,invalid_raise=False, usecols=[5,6,7],dtype=float,skip_header=9)\n",
    "    Y=Y.reshape(1, 3*np.shape(Y)[0])\n",
    "    Y=Y[0]\n",
    "    \n",
    "    list_bin=[]\n",
    "    for i_type in list_itype: \n",
    "        search=input_file[:-10]+\"_f\"+str(i_type)+\"_*bin\"\n",
    "        search=search.replace(\"Input\",\"Input_\"+str(i_type))\n",
    "        list_bin=np.append(list_bin,sorted(glob.glob(search)))\n",
    "    \n",
    "    XMAT=np.genfromtxt(list_bin[0])\n",
    "    XMAT=np.transpose(XMAT)\n",
    "    for input_bin in list_bin[1:]: \n",
    "        XMAT_tmp=np.genfromtxt(input_bin)\n",
    "        XMAT_tmp=np.transpose(XMAT_tmp)\n",
    "        XMAT=np.concatenate((XMAT,XMAT_tmp),axis=1)\n",
    "        \n",
    "    return(XMAT,Y) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Select the database\n",
    "\n",
    "list_itype=['1','2']  # 1=Gaussian 2=Lennard-Jones\n",
    "str_type=''\n",
    "for i in range(np.size(list_itype)):\n",
    "    str_type=str_type+list_itype[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training is done with\n",
      "Input/1_Liq.18.lammpstrj\n",
      "Input/1_Liq.22.lammpstrj\n",
      "Input/1_Liq.27.lammpstrj\n",
      "Input/1_Liq.31.lammpstrj\n",
      "Input/1_Liq.36.lammpstrj\n",
      "Input/1_Liq.40.lammpstrj\n",
      "Input/1_Liq.45.lammpstrj\n",
      "Input/2_Crys.0.lammpstrj\n",
      "Input/2_Crys.18.lammpstrj\n",
      "Input/2_Crys.22.lammpstrj\n",
      "Input/2_Crys.27.lammpstrj\n",
      "Input/2_Crys.31.lammpstrj\n",
      "Input/2_Crys.36.lammpstrj\n",
      "Input/2_Crys.40.lammpstrj\n",
      "Input/2_Crys.45.lammpstrj\n",
      "Input/2_Crys.5.lammpstrj\n",
      "Input/3_Clus.13.lammpstrj\n",
      "Input/3_Clus.18.lammpstrj\n",
      "Input/3_Clus.22.lammpstrj\n",
      "Input/3_Clus.27.lammpstrj\n",
      "Input/3_Clus.31.lammpstrj\n",
      "Input/3_Clus.36.lammpstrj\n",
      "Input/3_Clus.40.lammpstrj\n",
      "Input/3_Clus.45.lammpstrj\n",
      "Input/3_Clus.5.lammpstrj\n",
      "Input/4_Def.0.lammpstrj\n",
      "Input/4_Def.13.lammpstrj\n",
      "Input/4_Def.18.lammpstrj\n",
      "Input/4_Def.22.lammpstrj\n",
      "Input/4_Def.27.lammpstrj\n",
      "Input/4_Def.31.lammpstrj\n",
      "Input/4_Def.36.lammpstrj\n",
      "Input/4_Def.40.lammpstrj\n",
      "Input/4_Def.45.lammpstrj\n",
      "Input/4_Def.5.lammpstrj\n",
      "\n",
      "Testing  is done with\n",
      "Input/1_Liq.5.lammpstrj\n",
      "Input/2_Crys.13.lammpstrj\n",
      "Input/3_Clus.0.lammpstrj\n"
     ]
    }
   ],
   "source": [
    "#### Data reading #### \n",
    "list_files=sorted(glob.glob(\"Input/*.lammpstrj\"))[::5]\n",
    "N_files=np.size(list_files)\n",
    "list_Train=sorted(random.sample(range(N_files),int(0.9*N_files)))   #Spliting between Train and Test\n",
    "list_Test=sorted(np.delete(range(N_files),list_Train))\n",
    "list_Train=[list_files[i] for i in list_Train]\n",
    "list_Test=[list_files[i] for i in list_Test]\n",
    "\n",
    "#Reading Train samples\n",
    "print (\"Training is done with\")\n",
    "input_file=list_Train[0]\n",
    "XMAT_Train,Y_Train=readFile(input_file)\n",
    "for input_file in list_Train[1:]: \n",
    "    print (input_file)\n",
    "    XMAT_tmp,Y_tmp=readFile(input_file)\n",
    "    XMAT_Train=np.concatenate((XMAT_Train,XMAT_tmp),axis=0)\n",
    "    Y_Train=np.append(Y_Train,Y_tmp)\n",
    "\n",
    "#Reading Test samples\n",
    "print (\"\")\n",
    "print (\"Testing  is done with\")\n",
    "input_file=list_Test[0]\n",
    "XMAT_Test,Y_Test=readFile(input_file)\n",
    "for input_file in list_Test[1:]: \n",
    "    print (input_file)\n",
    "    XMAT_tmp,Y_tmp=readFile(input_file)\n",
    "    XMAT_Test=np.concatenate((XMAT_Test,XMAT_tmp),axis=0)\n",
    "    Y_Test=np.append(Y_Test,Y_tmp)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-07 0.9999941549973493 0.9999828553561407 0.00013492933758579433 0.00019542019270199357 6124 6165\n",
      "1e-07 0.9999927370268697 0.9999856082261952 0.00015040810887267134 0.00017904506927344282 1585 6165\n",
      "1e-07 0.9999909565410328 0.9999838013490427 0.00016783453058889745 0.00018995233080151353 1196 6165\n"
     ]
    }
   ],
   "source": [
    "#### Fitting from Lasso to Ridge #### \n",
    "alpha=1e-7\n",
    "for ratio in (0.00,0.50,1.00):\n",
    "    clf =linear_model.ElasticNet(alpha=alpha,l1_ratio=ratio)\n",
    "    clf.fit(XMAT_Train, Y_Train)\n",
    "    coeff=np.copy(clf.coef_)\n",
    "    \n",
    "    # Score/RMSE/MAS\n",
    "    Y_OUT_Train=clf.predict(XMAT_Train)\n",
    "    Y_OUT_Test=clf.predict(XMAT_Test)\n",
    "    scoreTrain=clf.score(XMAT_Train, Y_Train)\n",
    "    scoreTest=clf.score(XMAT_Test, Y_Test)\n",
    "    RMSE_Train=np.sqrt(mean_squared_error(Y_OUT_Train,Y_Train))\n",
    "    RMSE_Test=np.sqrt(mean_squared_error(Y_OUT_Test,Y_Test))\n",
    "    MAS_Train=mean_absolute_error(Y_OUT_Train,Y_Train)\n",
    "    MAS_Test=mean_absolute_error(Y_OUT_Test,Y_Test)\n",
    "    \n",
    "    print (alpha, scoreTrain,scoreTest,RMSE_Train,RMSE_Test, sum(abs(coeff)>1E-10), np.size(coeff)) \n",
    "\n",
    "#    output_file=\"Coeff\"+str(ratio)+\"_\"+str_type+\".txt\"\n",
    "#    with open(output_file,'w') as fw:    \n",
    "#        for i in np.arange(np.size(coeff)):\n",
    "#            fw.write(\"%g %g\\n\"%(i,coeff[i]))\n",
    "#\n",
    "#    output_file=\"ResultsTest\"+str(ratio)+\"_\"+str_type+\".txt\"\n",
    "#    with open(output_file,'w') as fw:    \n",
    "#        for i in np.arange(np.size(Y_Test)):\n",
    "#            fw.write(\"%g %g %g \\n\"%(i,Y_Test[i],Y_OUT_Test[i]))\n",
    "#    \n",
    "#    output_file=\"ResultsTrain\"+str(ratio)+\"_\"+str_type+\".txt\"\n",
    "#    with open(output_file,'w') as fw:    \n",
    "#        for i in np.arange(np.size(Y_Train)):\n",
    "#            fw.write(\"%g %g %g \\n\"%(i,Y_Train[i],Y_OUT_Train[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-07 0.9999868801090451 0.9999945880034787 0.0001915434751010273 0.00016878860343828248 63 6165\n"
     ]
    }
   ],
   "source": [
    "alpha=1e-7\n",
    "#### Fitting with LassoLars #### \n",
    "clf =linear_model.LassoLars(alpha)\n",
    "clf.fit(XMAT_Train, Y_Train)\n",
    "coeff=np.copy(clf.coef_)\n",
    "\n",
    "# Score/RMSE/MAS\n",
    "Y_OUT_Train=clf.predict(XMAT_Train)\n",
    "Y_OUT_Test=clf.predict(XMAT_Test)\n",
    "scoreTrain=clf.score(XMAT_Train, Y_Train)\n",
    "scoreTest=clf.score(XMAT_Test, Y_Test)\n",
    "RMSE_Train=np.sqrt(mean_squared_error(Y_OUT_Train,Y_Train))\n",
    "RMSE_Test=np.sqrt(mean_squared_error(Y_OUT_Test,Y_Test))\n",
    "MAS_Train=mean_absolute_error(Y_OUT_Train,Y_Train)\n",
    "MAS_Test=mean_absolute_error(Y_OUT_Test,Y_Test)\n",
    "\n",
    "print (alpha, scoreTrain,scoreTest,RMSE_Train,RMSE_Test, sum(abs(coeff)>1E-10), np.size(coeff)) \n",
    "\n",
    "\n",
    "output_file=\"CoeffLasso\"+\"_\"+str_type+\".txt\"\n",
    "with open(output_file,'w') as fw:    \n",
    "    for i in np.arange(np.size(coeff)):\n",
    "        fw.write(\"%g %g\\n\"%(i,coeff[i]/kbT))\n",
    "\n",
    "output_file=\"ResultsLasso_Test\"+\"_\"+str_type+\".txt\"\n",
    "with open(output_file,'w') as fw:    \n",
    "    for i in np.arange(np.size(Y_Test)):\n",
    "        fw.write(\"%g %g %g \\n\"%(i,Y_Test[i],Y_OUT_Test[i]))\n",
    "\n",
    "output_file=\"ResultsLasso_Train\"+\"_\"+str_type+\".txt\"\n",
    "with open(output_file,'w') as fw:    \n",
    "    for i in np.arange(np.size(Y_Train)):\n",
    "        fw.write(\"%g %g %g \\n\"%(i,Y_Train[i],Y_OUT_Train[i]))\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "459\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEQCAYAAAC5oaP8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5hdZbn38e89LdNSJ70TSEJJhMCEEpQSWihBLPEYUcRCVFRUzlEReEFEjmIXBBTFgyggRVBEiBAVIkICKZBCCqSRPqmTTC/7fv94ZpKQTKbt2bN2sn6f65prZfZe+3nuJDPrXuup5u6IiEj8ZEQdgIiIREMJQEQkppQARERiSglARCSmlABERGJKCUBEJKYiTQBm9lszKzGzRVHGISISR1E/AdwPTIo4BhGRWIo0Abj7TGB7lDGIiMRV1E8AIiISkayoA2iJmU0DpgEUFBScdPTRR0cc0aFt+3ZYtQr694dBg979XsITzN84n8KcQkb3Ht3+SnbuhBUroF8/GDw4uYBFJGlz587d6u599n897ROAu98L3AtQXFzsc+bMiTiiQ9srr8CECXDttfD1rx/4fs/bezJh8AT+dvnf2l/JsmVw9NFw1VVw663tL0dEOoSZrWnqdTUBxcywYeHYtWvT7x/X5zgq6iqSq2TkSMjJgZqa5MoRkZSKehjow8ArwGgzW2dmn4kynjjo3j0cN2xo+v1+Bf1YunVpcpVkZMCAAfDmm8mVIyIpFfUooKnuPsDds919sLvfF2U8cZCfH47Tpzf9fkl5CZvKNlFRm+RTwM6dMHNmcmWISEqpCShmzMLX7t1Nv39kryMBmLdxXnIV9ewJFUkmERFJKSWAGMrKgvLypt87ts+xAMzZkGRn+4ABUFcHtbXJlSMiKaMEEEPZ2VBZ2fR7Jw04CYDFJYuTq6Sxt3nhwuTKEZGUUQKIoS5doLq66ffGDxoPwNs73k6ukpEjw3H+/OTKEZGUUQKIoaOOgqKipt/r1qUbPXJ70Cf/gDkjbXP22eGYSCRXjoikjBJADI0aFTqCD2Z00Wh2VO1IrpITTwzH0tLkyhGRlFECiKHs7LAkxMH0yuvFopIkV+ju3j2MOZ2X5GgiEUkZJYAYWrEi3JgfbKLujsodbCrbRE19kjN56+rgxReTK0NEUkYJIIYaZwNv2tT0+yN6jgDgjU1vJFdR165hQpiIpCUlgBjq1SscN25s+v3GuQCvbXgtuYqKiqCqKrkyRCRllABiqDEBHGw9oHEDxgEk3w8waFAYBaSOYJG0pAQQQ717h+PBmoBOGXQKAG9tfyu5io44Ihznzk2uHBFJCSWAGBo7NhwbnwT2V5RfRPcu3RnUdVDTJ7TWBReE48FmnYlIpJQAYmj48HDMamY7oKN6HcXm8s3JVTQuNCWxZUty5YhISigBxFDjktBvN7PaQ2FOIfM3JrmMQ+Oeky+/nFw5IpISSgAxVFgYjs0t119ZW8nm8s2U1xxk2dDWyM8PU461L4BIWlICiKE+Dcv8NDdE/7QhpwHw5NInk6ssNxe2bk2uDBFJCSWAGMrMbH5TGICPHPcRAJ5a9lRylXXvDrt2JVeGiKSEEkBMZWQcfFMYgNMGn0aGZSQ/GaxvX40CEklTSgAxlZ3d/I6NZkb/wv5s3H2Q6cKtNXhwOK5dm1w5ItLhlABiasQIGDq0+XM+O+6zVNdXs72ymaVDWzJ5cjg2t/yoiERCCSCmhgxp+Zyzhp8FwOx1s9tf0XveE44HW3dCRCKjBBBT7rB+ffPnNG4PecfsO9pfUeOQoxdeaH8ZIpISSgAxtW5dyzflhTmF5GflM3t9Ek8AjY8a//lP+8sQkZRQAoipbt3CU4B78+eN6DWCHVU72r85TG5uGHd6sLWnRSQySgAx1a1bOO5oYevf9w55LwB/XfbX9leWn69OYJE0pAQQUz16hGNL/QAfPvbDAPx56Z/bX1nPnlBW1v7Pi0hKRJoAzGySmS0zs7fN7LooY4mbxqWgW0oAZw0/i0zLZNXOVe2vrF+/sD9wfX37yxCRDhdZAjCzTOAu4ELgWGCqmR0bVTxxc0rY82XP/sAHk5mRyXlHnkdpdRK7en04PEWwbVv7yxCRDhflE8DJwNvuvtLda4A/Au+PMJ5YaRyc05pVGk4bfBqLShaxs6qdG7yPGhWOa9a07/MikhLNbAmScoOAfdcHWAecElEssdPYBPTKK3DWWc2f27+wPwB3vXoXN5xxQ9sra9yB5qGHYPz4tn9eIpHwBDX1NdTU19A1pysAJeUlbKnYQmVtJZV1lVTVVVFTX8MxvY+hNlHLki1LWL9rPdX11VTXV1NbX4vjjB84nrpEHfM3zWf97vXU1tdRn6ijutrx6q6M6zeemhpn/vpFbC0rJafXRsiqoqq0G75jOMcUjaW2FpZsXs6uyiryjpiPZVdRsWkIbDyBIbnHUV+Xwapdy6iqrSb/+L9j2ZVUvjOGjHUT6GOjSdRnsL52MbWJWvJPfQDLqqHqrdPJXHMm3etH4AmjxJaQoJ68iT8AnOpFF5O5ZiKFNYMhAdtyl0FmNbnn3YIDtXOn0mX1WXSt7AMJp6Tbaiy3lMzzb8aBxCtfIG/VBLpV9CBBgs09N5JRsAkuuAkAf+EbFK45mcLyQuozEpT02kZGjxUkzr8JMJj+v3RbN5aheWu54edd+dCUszv0/zjKBGBNvHbAoEQzmwZMAxja0toF0mrDhoXjjBnwrW81f+7FIy8O566c0b4E0LgH5fPPt/2z8i6VtZWsLV3LprJNlFSUsL1yO9srt3NkzyNxnCVbljB/03zKa8qpqK2gqr6KqroqTux/ImbGW9vf4u3tb1OXqKMuUUd9nVNf2Y2heceRqM5jy45KdpcloGg5dF8H5b1h8RSsfCBenb83kPc8CAPnwebj4F+3wvYKqM+B7OFgQ+CCa2H4v2HFufDE76GqJ3gG2CSoz4ZPToQjXoQFH4Mn7uHx/f+i004K5b/2efjbzzlgJspVxdBzNfbqJPyFm9h/b6OMHn8lUbCLzLnjqH/1C+wdhHwWAFXDfk9ddiY5C4upWfBRSkmQQYK6jOOxjHrqj/sFteRSsGIMdcsmUFuXwHDI609m9m4qxw6klhy6rTuKzLVHkVMd3s8qO46c3M2UbhyNeyY9Nw8kf1tPulc5bsYu709eTRWb150IQNGOnvTcDUVVu6nLcCq75FCQkcnatcWA0bc8k77VuxheWsW6hdUwJekfoXcxb2kgeIqY2WnAt939gobvvwXg7t872GeKi4t9zpw5nRTh4S8jI7TOLF3a8rl5t+VRkF3A1m+0c23/Hj1CR3DMRwMlPMGmsk28vf1tVu1Yxbpd69hYtpGh3YeSaZks2bqEmWtmUlZTRmVtZbiLTtRSlFdEVV0Vu6p34fveJzmwaxBUFkFFEVT2hKpe0GcxDJkFFb3gmbvI3HE0VtWTRAISiQwyxv2OrOOewrcdRe0jDx4Qpx39JN5/AewcDK9/Zt+/AWTUwaDXoOvGUOeGk6CmYVxxVjkZWZVkFb4DuTXk1xhdd+fSvdrIp4LsrJ10sXJyc7ZSYkX0ri1jZM1GetVXkkcVXaycPCop9W6stiMZ4Bs4iTfIppZsasmknizqWJV3LJsLjqRP7XrGlb5AJvV7vrKoY9kxl7Gr95H03raMo1Y8CxmZZGSCZWZiWRmsOeMKavoMpue6hQxY+i/IysKzssI+qdnZlFzwCRLde1G4ZhGFb82HrGw8KwvLyYKsbMrPvQwKCuiyZhk561aRkZMFWZlYTjaWnUX9GROxnGyyNqwha1sJ1iUby8oko0s2GTnZ2Pji8Au4aVP4nWisOzMzHIcNC2u279wZfm/y86F37+b3cW2Gmc119+IDXo8wAWQBy4FzgPXAa8DH3H3xwT6jBNCx8vOhoKB1W/YW31vM3I1zWfD5BYztN7btlZ10EsybF5Ygzctr++fTXG19LRt2b2D5tuVMXzGd1TtWs2H3BrZUbKG0upSivCJq6mvYVLaJyrrKZsuyRDbZ9d3pUlBNblYudbOn0b3sZNhxBOU78ynfnUW3oWvoe+yb1Jb2ZskTH26iFKfph+xwjSksDNeTvLywLHh+/t6fh64FCfp3r6Bf7i66VO4gf+ta+pStpG/5KnqVraWwfBPVdVlUVRuFZRsYwUqyqSODejLwvbVmZ4cJJ/t3/pvBOefA6adDTQ388Y8hoIIC6No1/HnKlDBSYdeu8Jjao0d4PS8vnHfCCeEvUF4OJSXh9dzc8JWTEy6uskfaJQAAM7sI+BmQCfzW3W9r7nwlgI41YECYCFZV1fK5z771LBc9dBHnjTiP5z7xXNsr++IX4e674fHH4UMfavvnI+bubNi9gUUli3jszcdYvm05a3etZVvFNipqK0h44t135g0yLZMuWV0YVTSKMX3HkJ+Vz9rSDfRIDGfE4O4MKhzMP387kQ1LB7JpXR5bt2Swe7fRr18YPbthw8ETdO/eYaml+vow1aJPn/CZAQPCKtz9+oXrZo8eYbRX9+7h+ppdUQorV8Lq1aFj/t//hlWrYPPmcMfZ3DrhjUaNggsuCIU+8ggUFYWvPn2gf//w3plnQm0tvP763gC6dQsXa2s6OUlqHCwBRNkHgLs/AzwTZQxx1qdPeAJ1b/n38cKRFzK8x3CWbF1CwhNkWBvvsC6+OCSAV15J6wSQ8ASz183muRXPMXv9bJZuXUpJeQk19TXUJmrfdW6mZdI1pyvDewyneGAx5xxxDgMKB1BZV8mYvmMYVTSKzIxMnnoK/v53mD8fVqwIO2R27Rou0KtXh5vgfTXeyA4ZAqeeGi7uI0bAoEHh4t6/f/i/a7Y1oLQUli2D5cvh/hmwaFFYAGr79nBR3pfZ3jVB8vJCJaecAp/+dKhs5sywdvigQWGDn759w1+g8Yfm1lsPHkdODpx8cqv+7aXzRZoAJFoXXAALF4ahoLm5LZ9/28TbuPyJy5mxcgbnH3l+2yo7//xwcdn/4hOhRCLB8yuf56llT7Fy50oqaiuYv3E+u2v27pWZYRn0yO3B+IHjmXLcFEYVjSIrI4vigcUU5hTuOW/7dnjtNfjPX+Cll+Cdd8LoqkWLYO7c0IzbKDs73AyPHQuXXRaae0eMCIOlhgwJLRytVlkJb7wBzz0Hs2aFi/7GjeH1/ZmFivv1g+Ji+NrXQuXr14dhYQMGNJ1VNHLrsKUEEGPHHBOOGzaEC1BLPnTMh/hS7pe49u/XsujqRW2rLCsLxo0L/QARqaqrYta6Wfxs1s94ee3LbK3YuqfZxjBOGXwKVxx/BQlPMKbvGC486kKO6HnEAeUkErBkCYwcCQsWwP/7fzB9+oH17dwZtkO48spwHDsWRo8ON9XtagEpLw8X+b/+FWbP3nuXn0i8+7zCQjj+eLjlllBhdXW4wA8a1HTbeO/e7QhGDgdKADHWOAv4hRdalwC6ZHVhbN+xzHxnJi+ufpEzh5/ZtgoLC0OHXnU1dOnS5njbqrymnPvm38djbz7Gm1vepKy6jJpEaG/pktmFkUUjOXXQqUwePZmLR15MXvbBO6fffhuefTZ0Ybz2WrjBzsvbe6NdWBgu8OecA6edFvooBwxIsqn7nXfCMtqPPw4vvnhgZ+oFF4RZ1u7hP3PixJDVD8NOdkkNJYAY69cvHF98MTT3tsaPzv8RJ//mZL4545vM+uystlU4bFi4W50+Hd6fmknfG3dv5Mcv/5g/LPwDm8s373m9ILuAz5z4GS4aeREThkygV16vZsupqQl5qqwMfvQj+MlP3v3+wIHh+nv++eGCP3RoB/Rrrl0LDzwQ7vAXLtzbGZuTE/7dRowIzTGTJoWKBw5MskKJOyWAGDv++HBcvbr1nxk/aDwDCwfy6vpXKa0qpXtuC4sJ7evii+HXv4ZnnumwBJDwBI8seoRfzf0V2yq3sagkNE1lZWRxbO9juXT0pVw9/mqGdG95D8zq6nCXf//9odO2T5+9e9l37w7nnhtGJ06cuHejs6RUV4c7/IcfhsceC006jfLz4ctfDu1HRx8dOmk0tFE6mBJAjHXtGq4pbd2r5ZpTruG6f1zHdTOu455L7mn9BydNCsdXX21bhftJeIKHFz7Mna/eydyNc6lLhB7W4T2Gc9vE25g8ajJj+o7BWnlL7g5XXRVWqti37zQ3F7773bCv/dixHTRysbQUfvWrcKe/YkUYg5uVFf4jTjwxJMZPfAKOOLDvQaSjRToPoK00D6DjFRaG5vi2LNRZV19HwfcK6JHbgzevfpOi/KLWf7gx65S2bXVRd2dRySKmvz2dBxY8sOdOvyiviEtGXsI33/tNjulzTKvL27kztESZhTv+xk7ckSNDMpgyZe8SRkkrK4M77oD77gvj7xuNGQO33RYeKQoKNDZeUiYt5wFI9Lp3D+PS2yIrM4vpl09n0oOT+NCjH+K5TzxHTmZO6z48dizMmRPatFvRpLFw80JufuFm/r7i71TUhjbxUwadwjWnXMMXx3+RUUWj2hT74sXhrv7xx/cOzRw6NKyHdMUVobWlQ9TXh2adp56Cv/xlb3v+8OEhu3zlK2FUjkiElABibtIk+HM7Nvs6+4iz+e2lv+XjT36co39xNMu/tJyszFb8OH3mM2Ey2IoV4Xa7CfWJer7/0vf5+eyfs6UiTIPNzcpl8qjJ/Oj8H7X5og+hn+PKK0OHd6P3vQ+uvz70p3ZY8/rSpSGbPPtsaOPv0SNkllNPhYsu6qDOA5GOoQQQcyNGhElMFRWh37EtLn/P5Ty08CGeefsZzvrdWbz06Zda/tBJJ4XjK6+8KwFU1lZy56t3snL7Sv761l/ZsHsDGZbByQNP5rr3XsdlR1/W6jb9fZWXh2Gbt9wSLv45OfC5z8G113ZgE4976Mi98cawpAKEoZhTp8LPfhZmzoqkISWAmMtpaLmZOXNvH21bPP2xpzn6rqP5z9r/UHxvMb+e/GvGDRh38A80zj674w52TJnMPa/dw8OLH+bNLW/uWWJi8qjJTB0zlUtHX9rs2PzmvPVWuNDPmQO7d4dr8He/C1/6Usu7oLVaeTn88pehbX/JkvDacceFRPCRj2jUjqQ9JYCYa9wcftas9iUAM2PetHmMvWcsczfO5cR7T+Sz4z7LjWfcSL/CfuRk5pBhGbg7y7ctZ8HmBZzZqwvrShZw0g/2jsXvX9Cfj475KNe/73r6FLS/mWT79tC8/uCD4cY8Pz+M4f/85ztwftT27XD11aEjob4+PNU88ACcd16Y5ityiFACiLkTTgjHt95qfxkFOQWs/MpKnl/xPN9+4dvc/8b93Df/PjIsg3qvxxoWCG5cduGOo+DK12FcvxO4eNQlfPnkL9O3MPlmktmzw4CaioowsvKaa+Db3w4DjzrEpk2hD2P69NCJXVAQss13v6sRPHJIUgKIuTFjwrEjtus978jzOO/I81hbupYH3niAf63+F9srt1NZV0ldoo6RvUZy5rAz+fDaeXR99VHmZX8JJn6m5YJbUFcX1kL7whfCxf/97w8tMx12M757d2jLv+220LHbs2foPb72WjXzyCFN8wBkzwZEK1Z0UoVr1oQe2Pe8J6xk2U7u8OMfw3e+E67RxxwTLvxnnNFBcVZXh7ajJ58M8xYuvBAuuSRkGt3xyyFE8wDkoPLzO3mnxmHDwsqgq1eH2/d2bHO3dWu403/55XAtvu66MNInp5XTEZrlDj/8Idx8c5ip269fWBvilFM6oHCR9KHnV+HssyMYqXjjjWE7shkz2vzR6dPDA8TLL4dr89y58L3vddDFf/r0sIznN78ZktNXvxrWy9fFXw5DSgDCkCF7Fz3rNBdfHHpnv/3tNn1szZowkba8PKyEvHJleJhIWmlp6NC96KKwNeKkSWGv2Z/+NGyiK3IYUgIQ6uvD9S+ZkUBt1qVLaAqaPTtsVdiCysowubZxd8GHHgorLbR18toB3EP70fDhcOedYfLAggWhsp49kyxcJL0pAcieYZKvvdbJFX/lK+F4003NnrZuHRx7bHho6NYtxDl1agfU//LLYajQ7beH5p7Zs+Gee8J6RSIxoAQgHHdcOC5e3MkVf/rTYYPcJ5886CkLFoQF2lavDvsXvPpqByzYVlMDH/wgnH56aOaZPDk0+2jvW4kZJQDZ04beqU1AEMbQn356WJu5iT0CZs4M1+Tycrj88nDnn3SrzLx5IeM9+WTYC3fmzLBiZ9JtSSKHHiUA2XNH/c47EVR+443h+MtfvuvlsrLQH1tTA1//Ovz+9+0aLbrXzp0hi5x8cij8ppvCXf/73pdEoSKHNs0DEHJywte+u2F1mnPOgeLi0P7eMCegtBQuvTTc+d9+O3zjG0nWcffd8LWvhWwydSrcdZc6eEXQE4A0OOGEvZvEd7qvfx3efBOuvJL774fRo/dulZvUxX/VqtDc88UvhqFOt90Whg/p4i8CKAFIg0jmAjT6yEfgtNP484NlfPpTCTZvDhf/j360neW5w29+A0cdFRLLuHFhKNH113do2CKHOiUAAUKz+PLlYfmbKDzxub/zER7lJOby6C2LmTKlnQU9/3xYlvmqq2Dw4LAu9Lx5WqZZpAlKAAKEfQESCXj99c6v+5ln4MOf6koGCZ7hQqb8eEJYg6ctZs2CUaPC/o6zZoV2/pUr4WMfS03QIoeBSBKAmU0xs8VmljCzA1aok8539tnh+OijnVtvfT3cemtotfnBz3Pp8/kpsGsXnHkmbNvWcgH/+EcYK3raaWEc60knhSx29dVawkGkBVE9ASwCPgjMjKh+2c9//Vc4tmNttnarqQlrrc2aFZZ1vuYawkzcE08M8wKGDQtNObNn7/1QfX3orLj33jCH4Nxzw76PRx4Zeo7nzAlt/yLSokiGgbr7EqBdm3xLavToEb6WLQt346n+r9myJWxGU1IC//3fYW+VPebOhYULw5LM990XOnSzssLFf9/9K445JuzQ9bnPaRavSDuk/TwAM5sGTAMYOnRoxNEc3iZMCMve79rVgRunN6GiIrTYlJSEG/gf/KCJk8aODfvsfv7zYc2gbdugsDAsBtSzZxg6+r73aWMWkSSkLAGY2QygqaEXN7j7X1pbjrvfC9wLYUewDgpPmnDVVaFDduFCeO97U1NHIhFWWl6xIgzRf/rpFnZVnDAhglXqROIhZQnA3c9NVdmSGo2rIjz9dOoSwNVXw7//HSadvfRSWBVaRKKhYaCyR1FR6Ae4447UlL9tW1jDPzc3rMTco0dq6hGR1olqGOgHzGwdcBrwNzP7exRxyIHGjAlrAq1c2bHlVlWFCb9lZfDcczBiRMeWLyJtF0kCcPcn3X2wu3dx937ufkEUcciBJk8Ox9//vuPKXLEibLP7z3/Cr3+tBThF0oWagORdPv7xcHzmmY4pr7R075L/n/scXHFFx5QrIslTApB3GTgQCgo6Znewurqw2vPmzaFT+e67ky9TRDqOEoAc4Oyzw3DNurrkyvnUp8KcriFDYPr0FoZ7ikin06+kHGDq1NARvGBB+8vYuBGeeCLstPjyy+GpQkTSixKAHOCMM8LxV79q3+crK8OIn/r6sLbQ4MEdF5uIdBwlADnA4MGQlxc2ZfE2zr1+4YXQj/DSS/C734UlH0QkPSkBSJOOPx52727bU8C8eWE5/p07w+6LjSuMikh6UgKQJjXOBr7mGli9uuXzly8P4/tra8MSz9p9UST9KQFIk8aPD8s019aGxdvq6w9+7qJF4fyKCrj8cvjJTzovThFpv2YXgzOzXq0oI+HuOzsoHkkjt98OTz0VngDWrGl6+YZdu+Dmm8PxwgvDCs5aoVnk0NDSaqAbGr6a+5XOBLRQ/2EoMzPsD3D88WFM/2OPhU7hRCKs6XPLLWE3r1WrQpv/dddprL/IoaSlBLDE3cc1d4KZze/AeCTNHHEE3HknXHklDBp04OSw3r3DyB+t7yNy6Gnpfu3UVpShgX6HuSuugClTwsXfLDwZZGfDxImwZIku/iKHqpaeAP5kZg8Df3b38qZOcPeqjg9L0okZPPpo1FGISEdr6Qng18AlwCoze8TMLjOznE6IS0REUqzZBODuf3H3qcAw4Angk8A7ZvZbMzuvMwIUEZHUaNWYDXevdPdH3P0DwPnAOGB6SiMTEZGUalUCMLN+ZvZlM/sP8GfgOeCklEYmIiIp1dJEsKuAqcBoQhPQN9z9P50RmIiIpFZLo4AmAN8HZrh7ohPiERGRTtJsAnD3TwFY8HFghLt/x8yGAv3d/dXOCFJERDpeayfu302Y8DW14fvdwF0piUhERDpFS01AjU5x9xMbl31w9x2aDyAicmhr7RNArZllAg5gZn0A9QmIiBzCWpsA7gCeBPqa2W3AS8D/piwqERFJuVY1Abn7g2Y2FziHsDT0Ze6+JKWRiYhISrW2DwB3XwosTWEsIiLSiZptAjKzeS0V0JpzmvjMD81sqZktMLMnzaxHW8sQEZHktPQEcIyZLWjmfQO6t6Pe54FvuXudmd0OfAv4ZjvKERGRdmopARzdijKa2S68ae7+3D7fzgI+3NYyREQkOS3NBF7TCTF8GnjkYG+a2TRgGsDQodp6WESko7S6E7itzGwG0L+Jt25w9780nHMDUAc8eLBy3P1e4F6A4uJiT0GoIiKxlLIE4O7nNve+mX2SsNvYOe6uC7uISCdrUwIwswKgyt3b3O6/XzmTCJ2+Z7p7RTJliYhI+7Q0DDTDzD5mZn8zsxLCPICNZra4YSjnyHbW+wugK/C8mb1uZr9sZzkiItJOLT0B/AuYQRimuahxTwAz6wWcDXzfzJ509z+0pVJ3P6o9wYqISMdpKQGc6+61+7/o7tuBPwF/MrPslEQmIiIp1WwTkLvXNjQDLWrunI4PS0REUq3F1UAbmn3eaNgFTEREDhOtHQU0AFhsZq8C5Y0vuvulKYlKRERSrrUJ4JaURiEiIp2u2QRgZubBiy2d0/GhiYhIKrXUB/AvM/vy/u3/ZpZjZhPN7HfAJ1MXnoiIpEpLTUCTCIu1PWxmRwA7gVwgE3gO+Km7v57aEEVEJBVaWg20CrgbuLthvH9voNLdd3ZGcCIikjpt2RKyFtiYwlhERKQTtTgPQEREDk9KACIiMdXSaqC/MLMJnRWMiIh0npaeAN4CfuV3dj4AAAmgSURBVGxmq83sdjM7oTOCEhGR1GtpMbifu/tpwJnAduD/zGyJmd1kZqM6JUIREUmJVvUBuPsad7/d3ccBHwM+ACxJaWQiIpJSrUoAZpZtZpPN7EHgWWA58KGURiYiIinV0lpA5wFTgYuBV4E/AtPcvby5z4mISPpraSLY9cBDwP807AImIiKHiZaWgji7swIREZHOpYlgIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU0oAIiIxFUkCMLNbzWyBmb1uZs+Z2cAo4hARibOongB+6O7vcfcTgKeBmyKKQ0QktiJJAO6+a59vCwCPIg4RkThr9abwHc3MbgOuAEoBLTkhItLJUvYEYGYzzGxRE1/vB3D3G9x9CPAg8KVmyplmZnPMbM6WLVtSFa6ISOyYe7StL2Y2DPibu49p6dzi4mKfM2dOJ0QlInL4MLO57l68/+tRjQIauc+3lwJLo4hDRCTOouoD+L6ZjQYSwBrg8xHFISISW5EkAHfXdpIiIhHTTGARkZhSAhARiSklABGRmFICEBGJKSUAEZGYUgIQEYkpJQARkZhSAhARiSklABGRmFICEBGJKSUAEZGYUgIQEYkpJQARkZhSAhARiSklABGRmFICEBGJKSUAEZGYUgIQEYkpJQARkZhSAhARiSklABGRmFICEBGJKSUAEZGYUgIQEYkpJQARkZhSAhARiSklABGRmIo0AZjZ/5iZm1nvKOMQEYmjyBKAmQ0BzgPeiSoGEZE4i/IJ4KfANwCPMAYRkdiKJAGY2aXAend/I4r6RUQEslJVsJnNAPo38dYNwPXA+a0sZxowDWDo0KEdFp+ISNyZe+e2wJjZWOAfQEXDS4OBDcDJ7r6puc8WFxf7nDlzUhyhiMjhxczmunvx/q+n7AngYNx9IdC38XszWw0Uu/vWzo5FRCTONA9ARCSmOv0JYH/uPjzqGERE4khPACIiMaUEICISU0oAIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU0oAIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU0oAIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU0oAIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU0oAIiIxpQQgIhJTkSQAM/u2ma03s9cbvi6KIg4RkTjLirDun7r7jyKsX0Qk1tQEJCISU1E+AXzJzK4A5gD/7e47mjrJzKYB0xq+LTOzZe2srzewtZ2fTbV0jS1d44L0jS1d44L0jS1d44L0ja2tcQ1r6kVz944JZ/+CzWYA/Zt46wZgFiF4B24FBrj7p1MSyN545rh7cSrraK90jS1d44L0jS1d44L0jS1d44L0ja2j4krZE4C7n9ua88zs18DTqYpDRESaFtUooAH7fPsBYFEUcYiIxFlUfQA/MLMTCE1Aq4HPdUKd93ZCHe2VrrGla1yQvrGla1yQvrGla1yQvrF1SFwp6wMQEZH0pmGgIiIxpQQgIhJTh30CMLPfmlmJmaVVR7OZDTGzf5nZEjNbbGZfiTqmRmaWa2avmtkbDbHdEnVM+zKzTDObb2ZpNXrMzFab2cKG5U3mRB1PIzPrYWaPm9nShp+306KOCcDMRu+zHMzrZrbLzL4adVwAZva1hp/9RWb2sJnlRh1TIzP7SkNci5P99zrs+wDM7AygDHjA3cdEHU+jhpFQA9x9npl1BeYCl7n7mxGHhpkZUODuZWaWDbwEfMXdZ0UcGgBmdi1QDHRz90uijqeRma0Git09rSYOmdnvgH+7+2/MLAfId/edUce1LzPLBNYDp7j7mohjGUT4mT/W3SvN7FHgGXe/P8q4AMxsDPBH4GSgBpgOfMHd32pPeYf9E4C7zwS2Rx3H/tx9o7vPa/jzbmAJMCjaqAIPyhq+zW74Sos7BTMbDFwM/CbqWA4FZtYNOAO4D8Dda9Lt4t/gHGBF1Bf/fWQBeWaWBeQDGyKOp9ExwCx3r3D3OuBFwlD6djnsE8ChwMyGA+OA2dFGsldDM8vrQAnwvLunS2w/A74BJKIOpAkOPGdmcxuWMEkHI4AtwP81NJv9xswKog6qCR8FHo46CAB3Xw/8CHgH2AiUuvtz0Ua1xyLgDDMrMrN84CJgSHsLUwKImJkVAn8Cvuruu6KOp5G717v7CcBg4OSGR89ImdklQIm7z406loM43d1PBC4EvtjQ/Bi1LOBE4B53HweUA9dFG9K7NTRLXQo8FnUsAGbWE3g/cAQwECgws49HG1Xg7kuA24HnCc0/bwB17S1PCSBCDe3rfwIedPcnoo6nKQ3NBS8AkyIOBeB04NKGtvY/AhPN7A/RhrSXu29oOJYATxLaaaO2Dli3zxPc44SEkE4uBOa5++aoA2lwLrDK3be4ey3wBDAh4pj2cPf73P1Edz+D0LzdrvZ/UAKITENH633AEnf/SdTx7MvM+phZj4Y/5xF+IZZGGxW4+7fcfbC7Dyc0GfzT3dPizszMCho682loYjmfNFjixN03AWvNbHTDS+cAkQ802M9U0qT5p8E7wKlmlt/we3oOoY8uLZhZ34bjUOCDJPFvF+Vy0J3CzB4GzgJ6m9k64GZ3vy/aqIBwN/sJYGFDWzvA9e7+TIQxNRoA/K5hZEYG8Ki7p9WQyzTUD3gyXC/IAh5y9+nRhrTHl4EHG5paVgKfijiePRrasc+jc5aDaRV3n21mjwPzCM0r80mvJSH+ZGZFQC3wxYMtpd8ah/0wUBERaZqagEREYkoJQEQkppQARERiSglARCSmlABERGJKCUBEJKaUAERSwMzuNLN5ZjY+6lhEDkYJQKSDNcwE7kuY3JQ2y1WL7E8JQKSdzGy4mVXuM5MbAHcvJ8ymfgG4o+HcvIZNT2rMrHfnRytyICUAkVayYP/fmRUNq6bue14RYQ353UA9gLtXNpyXLuvKiygBiDSn4S5/iZndTVgbpjVrr99IWE9+MXBsKuMTSYYSgEjLRhO2FB3X0o5VDZv7TAAeIawgeVzKoxNpJyUAkZatacN+yN8FvuNhlUUlAElrh/1y0CIdoLw1J5nZCYT12d9rZncBucDCVAYmkgwlAJGOczsw2d3/AWBm/QhryYukJTUBiXQAM5sIFDRe/AEatjgsMLNe0UUmcnB6AhBphruvBsa04rx/Av9s4vXuKQhLpEPoCUCk/eqB7vtPBGtK40QwIBtIpDwykVbQlpAiIjGlJwARkZhSAhARiSklABGRmFICEBGJKSUAEZGYUgIQEYkpJQARkZhSAhARian/D0knC561V3mLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "459\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEQCAYAAAC5oaP8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hU1dbA4d+e9IQUSEgIEAIRpDcJVbGACijY9aooKCp2Rb1XL+pnR+x6FUS9iqh4RREQRUFEBEVUegmCSA+hEwLpyczs74+dkAghM5nMZEbOep9nnkkypywgnHXOLmsrrTVCCCGsx+bvAIQQQviHJAAhhLAoSQBCCGFRkgCEEMKiJAEIIYRFSQIQQgiL8msCUEpNVErtU0pl+DMOIYSwIn8/AUwCBvo5BiGEsCS/JgCt9Y9Atj9jEEIIq/L3E4AQQgg/CfZ3AK4opUYCIwGioqK6tWnTxs8R/b0VFsLvv5uvTzkF4uIqPtNoVuxeQaN6jThUeIjIkEjS6qd5fpImTSArC1q2hNhY7/wBhBA1tnz58gNa64bH/lz5uxaQUqo5MEtr3cHVtunp6XrZsmU+j+lklp0NTzwBb7wBn34KV13118/fXvY2XRp14e7Zd5MQmcA3Q7/x7CRvvgndusEFF8DkyTB0qFfiF0LUnFJqudY6/difSxOQxTRoAPfea74uLj7+81vTb6Vn057EhceRU5Tj+UkefRS6djXfHzni2XGEED7l72GgnwC/AK2VUjuVUjf5Mx6ryC7rdi8qOv6zPXl7WLRjUe0SAMCePXDwoPn68GHPjyOE8Bm/9gFora/x5/mtqn9/817VE8CEpRN4+senuaHzDbVLAFdeCcHB5iUJQIiAJE1AFhQTY96rSgDN45qj0QTZgjhUdAiP+4hSUmDHDtP5K01AQgQkSQAWVD4gp6omoNS4VACCbcEU2Ys4XOzh3XuzZpCZabKNPAEIEZAkAVhQ+dDPEz0BAKDM284jOz07SUoKlJZCZKQkACEClCQAC4qJAaWqTgBNY5qiUBTbzYe1SgAAoaHSBCREgJIEYEF33gkREVU3AYUGhTLjHzO4ocsNQC0SQPfu8MEHkJgoTwBCBChJABY0eDDUq1f1EwDAxW0upnfT3tiUjczDmZ6dJDkZhg2Dhg0lAQgRoAK+FITwvn37TtwEBLB+/3oy9mXQqF4jz58AAJYvN/0A0gQkRECSJwALeucd2LsXCgqq/nzymslcM+0amkQ3YWduLRLAddfB6tXmCcDPJUeEEMeTBGBB5fMA8vKq/rx5XHMc2kF8ZHztngBSUiA/HxyOE2cbIYTfSAKwoPJ5ACe6JpcPBY0Kjqp9Aihv/5dmICECjiQACyp/AsjPr/rzo5PBgoI5UnyEI8UeXrybNau48EtHsBABRxKABZU/ARQWVv15s9hmADi1E/DCXACQBCBEAJIEYEEdOkDnzmA7wb9+eHA4v938G8M7DwdqkQAGDoTXXjNfSxOQEAFHEoAFJSZCWho4nSfepkeTHrRPbA/UIgE0bgz9+pmv5QlAiIAj8wAsyG431+MTjQICWLRjEct2mdXXPE4AWsPPP5uvJQEIEXAkAViQwwHz51e/TO/XG7/m5V9eJjEy0fMEoJRZGQykCUiIACRNQBYUFmba/0tLT7xNeuN0Sp2lNIhoQOYRD8tBAKSaEUXyBCBE4JEEYFFhYaYp6ER6Ne0FQEhQSO3nAthskgCECECSACzKVQJoEtOElJgUiuxFtUsAzZqZvgBpAhIi4EgCsKiICDMKqLoSPb2a9iKvJI+cohzySqrpMa5OSoo5yYEDnu0vhPAZSQAWdeGF5r26foD/Dvkvz5/7PABZR7I8O9GwYWbSwYmmHQsh/EYSgEW1aWPeq1oUplxseOzRWcEedwQnJUGjRpCb69n+QgifkWGgFnXwoHk/0ZoA5d5d8S5Qi7kARUWwZw8cOuTZ/kIIn5EEYFG//WbeXSWADQc3ALVIAMHBsGaNWRxeCBFQpAnIourVM+8nKghXrk/TPgDsOLzDsxMFB5uTVdfWJITwC0kAFlVeEtpVy0zvlN6AWSbSY/Xrm+nH1fU4CyHqnCQAi4qONu+uRmeWTwjbmrPV85MlJJh3mQsgREDxawJQSg1USv2hlNqklPq3P2Oxmrg4856dXf12qbGpJEUlcaioFp24SUnmPSfH82MIIbzObwlAKRUEjAcGAe2Aa5RS7fwVj9X06GHeGzSofjulFPf0vIeC0gIKSj1c1/faa827DAUVIqD48wmgB7BJa71Fa10CTAEu9mM8ltK0qXl3OFxvmxJjVvbKPOzhXIAmTcy7q8cNIUSd8ucw0CZA5SvKTqCnn2KxnIgI875qFQwZUv22pU7TeTt9/XRG9x1d85OV152eM6digRgR8JzaSYmjhBJHCdGh0SilyC7M5lDhIUqdpZQ6SrE77Ti0g/TG6QCs27eOnUd2Ynfaj76CbEFc0uYSAOZvnc/m7M3YnQ6c2kFxicZWHMfV7a+jtBSmZXzJ1gM7iU7KxhZSQm52BOS04Io2V2O3w9S1M8g6vJdG7TYRFFZE9s4E1N4uDG55CQ4HTF07nQP5OTTv+wu20CL2/dGC0N196df8PBwOmL5uBkeK82l54ZfooGL2rOpCxO7+dE/uhcOhmfXHNxSXlpB29Vs4tZO9v5xD1J4BtEtoh9Pu5PutC9FBRaRc/SJONPt/uIyYXf1oEZmCw+7gl+y1BIUfJmnoGFMB5ZubqL/7DBoHNaDU6WB14Q6CY/YQd/UTaA05X/yTBrt7EO+MpsTp4E+1n+D4LURd8TgayPtsDA32duKcbrt4+M3OtCm/mfISfyYAVcXPjqtMo5QaCYwEaNasma9jsoz69c37okWut+2W3A2AX3f+6tnJWrRw/2SiWkX2Ig4UHCC3OJe8kjzyS/PJL8mnb2pfYsJiWLVnFXM3zz3aZFdYWkihvZCx/cfSMKohn2Z8yjsr3qHIXmReJSUU5IYx49J5hDji+M+P7/HB0s+x11+PM2Y75CfAuit59qz/4CgJYWbGUpZtXwedPobGK2Bve9SCZ+if0onCIicb9gRxMDceBtwPzX+CzeeiZnxEVFAxdruiuPh0tP1sGN4PWiyENdfC9I+57+if8CLzNrKbOf6yW2HWaF47+vmlAATd1Q4SNqF/HoXzu0v439HPLwNgie1ZHNH7CVr4EPafz2PWMfv/Ef8o9tASQn84n5KVZ7AIJzYc2G2DUDYHu9NGU6pCiVqcjHPDqeywl6DQHI7oSVBILmtaOSjVocSsSaZgSzJFxebz0ujWqPC9bFiShNbB1N+USEFmFNlFRWilCYpNJDyvgO1LOgGK+F2xBGc5CS7KAZsmsn4oUdpG5oqeoBWJh4KJyz5M7pdFfN5jBY8+4t0EoHR11cB8SCnVG3hCaz2g7PvRAFrrsSfaJz09XS9btqyOIjy5FRaauVldusDKldVvW1BaQNSzUaTVT2PzPZtrfjKtzXyAxETYvduzgE8STu0kuzCbAwUH2J+/n4OFBzlYcJD+af1pHtecFbtX8MLPL5BTlMOhokMcKT7CkeIjTL1yKn1S+jB5zWSun3F9xQE1cKQJ/xs0l6Sgdny27DveXvQZNFxHeNpKwkqTKPnyVXrEn09JQTjb9x7iQE4Jjc6bTKN+0yje14zVj049Ls7GV4+hQd8p5GW2ZNtzMyo+sJVCSBFhl9yFs/0nlO7sAF9MguAiCC6GoGJsqpSgHq9RmrSBsMwu1F96LdElEOFwEKxKCNOlHG66jt9j40jZH8X520poUFpClLOYcAqo5yxkbnQ6s0MH0TN/LaMLxxNCKSGUEoSDYOw8yWN8yyDOZw6v8ABBOI6+grEzlMks4kyGMJMxPEopwTgIKjtKCLfZxrMlOI2Beja3OCdSqoKxE4RDBWEniBcj7uKQLY4+9l84s3QxGoUDGw5t/spnhJxPiS2EDo4NtHJuR1Nx96rQLA1uj1KQ6swiwZmNQmNDo8pemUGNCMFOA2cO9XQBymZD2RS2IMBmIz8khmAbRFJMiM2BIyiE8165lhaXn+nR751SarnWOv24n/sxAQQDG4H+QBawFLhWa73uRPtIAvAerU2Z/mbNYPt219tHjImg2F7Mjvt20DSmac1PGBdnOoGLiiAkpOb7/00cKT7C3M1zyTqSxa7cXezO282evD2M6jWKC1pdwM87fuaM9884br/PrviMK9tfyaIdixgxYyQxtkY0iAsmLjyOrPkX0SFiAI4jDdmalc+WzHy6nrmbi0aupKC4lDt73XLc8VoMmkHyFS9xIKeIzWOn4gg9BGGHICwXwo5A+8/g1G+guB6sugFC8iE0D1tYERGRTuo1PEhkeCgheU1J2J1EYq6iYZ6dhnl24nJLWR/WkUWRZxORY+edTbcRaz9CnCOH+jqbMEp4iOd4gYdozla2kvaX2OwEcSfj+ShiJB3CN/FW/jAKQ6IpCalHSWg9SsOi+LH5MHYldKJR4VZ67poONtNYoXCitJPcsIbgdBJZsJ+4vCyCivMJLcoluOAwYaW5hFFCGMVHX6GUHH0F4ahofggLg6goczcUGWnaRsvfK7/Cw49/hYW5/woNNa/KX4eGmv8LQUFe/i08XsAlAACl1AXAa0AQMFFrPaa67SUBeFdYmPk9dmetlrbj27LhwAYeO/MxnjznyZqfrHVr2LjRdDp07lzz/QNEkb2Iz3//nC2Hthx97Ti8g3t73st9ve9j66GtpL1uLnihQaEk10umUb1G/PuMf3NJm0vYn7+fj1dPIcKRTFrjWBIiE5jyVnO2bYwhc3sQO3bA7t2a/oPyeej138g8ksld/a6g4HAkoTE52KL34YjYTckp06HnOBPUihshPAcisgmpl0d8vCIxIYiE2CjiI+JpENGA+uH1iQ2rjy2/EfbDjSjKTqDwUCy5B+rRbO0vRO7eRdShPcQcyaJ+4S5+cJzJc4wmCDslhGI7pnX23ah7ea35ayTFFvHqH4MoikqgtF597LEN0HH1OdTpbIo69yQmvITkvasIT4olPD6KyKBiogoPEJWThW3vbvNEuGcP7NsHe/ea9/37q68eGxFhbigqv2JiTF9TbKz5Ojr6r6969SpeUVEVrzq4+AaCEyUAv9YC0lp/A3zjzxisLCLC3JRrbZbvrc4p9U/hcNFhhrR20WN8Ii1awObN5j98gCeArCNZZOzLIGNfBuv2r2PDgQ2cm3YuT53zFDZlY9iMYWg0TWOa0iKuBX1T+9KyQUsAUmJTWH3bahpHNyY+Ih6lFF99BT9PgvfXwx9/NGTLlrvp0MnOuBlLWLtvLZ9Mi+fQgRKC6mdS0ngzzlP+5Lt6y/juo5kmoFseIL5BMClxjWkS3YTG0Y1JrpdMcvQEkqKSaDSiEUn1kkiMSqQkL4pNmxQbN8KWLbBtG4SuWoRtewYhOZtJdS6hBVtZR3uG8yE2G+zkTpKdu3CiOBKeRG58E0K7QNoVkJAQzJ+/vUdEchyRLZKIPiWRsJREbo6O5mYFEA78YOIsLTUn3LwZtq6A1TPM99u3Q2am+bd3Ov/6l22zmabBpCTzOvVU833DhmYCYUKCGascH2/e69c3dy3CK6QYnIU1bQrJye5t2zGxI3M3z6VDYgfPTtaqFSxeDAMHera/D2it2ZS9iWW7luHUToZ2GgpAt3e6sTd/LwBJUUm0bdiWxKhEwNzVb7x7I01jmhIeXHEhys6Gb7+F5cuDWbWqE3v3woIFmt25u3lhvI1fvk8guvEebAkbCTp9DSsbLOH0iWVdl1comsWlkFY/jRZxLUiNTSU17hJSY+8lJTbluHOBuVleswZ+/xnmZtjJW7GRsD+/oUnuetrxOxDC40wmORm+y32A9nlLsAeFkpvUgpImLUjq045dD5vrbPDvs6F+fWzJycQFBxMHpABnl5/s0hv/+he3f7/pONqwAdavN+9//mku9JXHFYeEmDWhU1Ph/PPNL1zTpmZYcOPG5pevYUPL3IUHIkkAFhYfDyUlru/+oWKR+Onrp7M7dzcP9HmgZidr0sQ8buTnm0dvP/rv8v/yxR9f8OvOX8kuNHMTOiZ2PJoA3hr8FvXD69M+sT0JkQnH7Z8W15L166FtW3MD+9hj8PTTFZ/HJWcT2ngDic9dyYGSXdC5PvTIJTSmAW0T2tI6vjWtE7rROv5aWsW3IjU2lbDgsCpj1dpcV5cuheXLYf3yAhzLV5F0aD0TuQmAL0Ku4eLSz832SlGQ2AJ7tx4UTiu7Wf79fYiOJrhJE+rbqpj606lT1X9RDoe5uC9fbpru1q41r717K7aJjDSLS3Tvbib8nXKKeaWlmQt8VecTAUMSgIVFRMDq1eblqlWme5PuAEzJmMJXG79i8KmDaZ3Q2v2TlT9qdO9ubl2Dff+rV2wvZnHmYuZtmceqvauYdc0slFIsyVrC9pztXNrmUno26Un3Jt1p37D90f3Kx6xXtmkTzJ0L330HP/5o7vi/WLCVrMjZLA46RNzgYHIafAuNl5MbkU/7xPaclnw+XRt1pWNiR9ontj/6FFGd/HxTqnvxYvNauhRaH1jEcD7gOn6lLb8ThGlGGT7jUlr3bkDiqltg30XQsSOqdWuiyid5lGvn5gT7zEz45Rfz+u0384tRUDb7Ozwc2reHCy6Ajh3N123bmsQuF/m/LUkAFhYVZfrffvrJdQJIiUkhMSqRiOAIQmwhvLXsLV4d+Kr7JytPAOvXm1fHjp4H7sKCbQt4cfGLLNi2gILSAoJUEOmN0zlYeJCEyATeHvI2NlX9RaukxKyVEB1tLvwDBpifJzTOo0GXFTgbfsEl33wAkdkk10vm3GGn06fpEHqnPEvnpM5EhERUe/xyubnm7/+HH2DhD05KV2bQ17mAs1nAz61eZMiQUxhWsoEzvpqGrVdPbL0ug27dID2dMxuX1fEYcH7N/5K0hj/+gIULzeunn2Bn2ZoP4eHmHLfcYt5PO83c5UtTzUlHEoCFlddoy3SjwoNSivTG6azbv47L213OpNWTGNN/DJEhbi70UrmzYcUKryUAp3by287f+PKPL7m247V0TOpIXkkefx78kxu73MiAUwZwdvOziQ6LPrrPiS7+xcUwezZ8/jnMmgUPPAD/uGMjS9Q3pA11sjX+bQ7U34iOjOf8tP6c2+I5+rXoR1r9NJQ77WiYPtCVK815Zs82N9opjq2MtT3CI7a5xDnNUm2O5mlc9m4WnHkKlAyD4BG1v9Pes8c8wsydC/Pmme/BtMefeSacfjr06mXuBk7iobqigiQACyuv0OBOAgDo3rg7czbN4eUBLzMlYwqfrP2Em067yb2dGzc276Ghpk15+PCaB1zGqZ0szlzM1HVTmbZ+Glm5WQTbgjk1/lQ6JnXkwlYXMvjUwW4fT2sYNQo++sisjxBX30Gr03/no8P/4bFx7wHQpW8X/u/UqxnSeginJZ/m8gmispIS+P57k1i+nqVpuC+DS/iCM9LactaDVzCwWz363vkD6vwLoH9/OPtsglJTKw4QGur2uf7C6YRly+Crr8xr9Wrz84YN4dxzTVmOs882bfZuJjBxcpEEYGHli8K4MxEMTAJwaifhQeFHy0O4rUEDcyFr2NAkgBrSWpNdmE18ZDyFpYUMmDwAh9PBoFaDeL7t81x46oXEhZsa1+7cjefkmJvgK64w1749e+206Z1JXpsJrI16hWVBDno06cGrHV7lsraX0Sy2ZmVI7HZzoz1lCnz5JbQ4vJIbQ//HsyHTSWILWinUkHvg2SuAhnDZLu9chO1206QzdSrMnGnu8m02OOMMGDvWtGV17izt9gKQBGBp5U8AWVnubV9e8GvZrmX8POLnE45cqZJS0KiRGdfdpYvbu209tJXJaybz8dqPiQiJYOWtK4kKjeLb676lc1LnvzTtuGPdOnjtNfj4Y1MO47vlf/Ll3jf4tstHHC7OoWWDljzd+Qmu6XANpzQ4pUbHBsjIgEmTzPFL9xxAN0jgssvg5cV3ELd5GarvuXDpQ6iLLjJ/H+Vqc/F3Ok2dpU8+gWnTzDDNqCgYNAguvti8x8d7fnxx0pIEYGHlTwCzZ7u3fVK9JFJiUli2e9nRi/83f37DR2s+YvKlkwmyuegkTE42varjx7s81+w/ZzN20Vh+2vETAGelnsXQjkNxaic2ZeOMZseXU6jOtm2mmWfmTIiI0Jw1JIvDHZ/jvK/GExoUyuVtL+eW027hrOZn1ah5B8z8p2nT4I03YMXiQq6yTWNug/doF/orjt93EZpUHzL+a/783rwQb9kCH35oXlu3miGZQ4bAlVeai36km/0zwrIkAVhY+RNAdbPuj9W9SXeWZi09+v2OwzuYkjGFpKgkXhv4WjV7YvoBNm40je5FRRU1qTElFmZtnEXfZn1JqpfEvvx97C/Yz7P9nmVop6E1boIpV1BgroP16sHy5ZpLbl3D+rQ7mFO4mEYRjXim+zPcmn5rleP9XTl4ECZMMC92ZfFM7EtcE/4+4UWHISYNRj1KUFjZnX0HDyfQHau01GSxCRNg/nzz5NCvHzz5JFx2md/nWIi/Ga313+bVrVs3LbxnwQKtQevBg7XesMG9fcb+NFbzBPpgwcGjPxs1e5TmCfQtX96iN2dvPvHOd9yhdVyc1k2baj1qlM4pzNFT1k7RQ6cN1TFjYzRPoMf9Nk5rrbXdYddOp9PjP9vGjVpffrnWPXpoXVxaqj9Y9YE+5dXWmifQHd/sqCetnKSLSos8Ova+fVr/+99a14ty6lgO6QEDtF7w9gbtDA3V+uqrtZ4/X2uHw+PYq7R7t9aPPaZ1crL5R2vWTOtnntF6xw7vnkeclIBluoprqjwBWFj5E8CsWXDDDaZemyuV+wHOP8WMP3/p/JcodZbyzvJ3mLhyIjvu20Hj6MYU2YsIDQrFpmxordneKJz8kBzaN2nNkVW/kfRSEsWOYuIj4rm87eVc2/Fazml+DoDr5qQTyM6Gxx+Ht96CsDDNoOEZtHvjH2w+sp7OSZ2ZMXAGF7W+qMbNPOXHHjsW3hrvYHDhVNbGPUf981oSO+dzoDVcuadioQVv2bABXn7ZNPOUlppSGu+8Y5p4ZFy+qCVJABZW3gcAFXOAXKkqAQTZghh3wTge7vswczfPpXG0GfJ5xsQzWL57ORHBESilKHAWcP4A+NbWmpjp03npuRfpmnwavZr28viCX9maNXDeeaZpZsg1u9ne5SY+z5tNp4hOzBjk+YW/pMR0W4x9soRLjnzIn/WepxGboFEbuL7SKqbevPgvX27qS8ycaSZm3XQT3H8/tGzpvXMIy5MEYGHlTwAhIe7PBYgLj6NVg1Ys3bX0uM8aRzfmhi43HP3+9vTbyTySSX5JPnannTZ77HR/exzcmgZ5edxlPw2anV7rP4fdbipLnHoq9OqbR36PJ/ii8GUaq8a8f/H7XN/peo8SjNbm+vuvf5lSEO+3fJEbDj8Kp3aDh6fBJZd4fzjlqlXwxBPmxPXrm0JDd91lhs8K4WWSACwsumwEZUyM+wkATEfwwm0LXW533CSxVatg9zgz8SgiAiZPNrNPPaQ1vP02/Oc/8POvJby95mW+6/o0NruNp85+igf6POD+TOVjZGXBnXdogr6cTq/URN6Y3ZeB3W+D37qa5hdvT5zatAlGjzazxWJj4amn4N57//qYJoSXSQKwsNBQ07oQHW3KILgrPTmd/639H7tzd5Mc7WY9aagoB5GTA6+/bgqKeejAAbj5ZnOj3P3MbHpNGMKfxYu5rO1lvDbgNVJiUzw6rtNpmtg/e+A3xhTeT28W4+x9DbaBfYF4UwzNm3Jy4JlnzN9HaCj83/+Zpp64OO+eR4gqyHRAi4uNNe3mX3zh/j7llUGX7arh6mzltd937zZX7969a7Z/me+/NxWMZ8/W9L75M5aenYA9YjdfX/s1066a5vHFPzMTru6zg+jbhzK/oBfdE7bAe+9hm/yRR8erltNpeqpbtYJXXoHrrjM19Z96Si7+os5IArC4mBg4cqRm+3Rt1BWbslXZD1Atm81UoCtfGH7lSnjvvRodQmvTNxoSmU/83YP5tenVPHD6/WTckcEFrTy/O58xw1RISF35Bf8ImY5+5FGCN2+EESO8P9omI8OUZrj9dvMUtGwZTJzo/uo8QniJJACLi401s2QvvLDiuuxKVGgUnZM6M2fTnJqfsHHjihNNmgR33GGaQVwoLDSF2grtBaSNfJAd/0gkuvkmFo1YxEvnv+RxW39hIbwx+Fs+uOwL0tLg1tV3ELx5I+qZpys6SbylqAgefRS6djV3+x9+aOpAn3aad88jhJskAVhc+RPAN9+YagLuGtF1BEt3LWVJ1pKanTA5GXbtMl9ff70ZYzl1arW77NxpqhUPuvgIXd86jUl/vsR9Z97KqltX0SelT83OX0nW0l0sanwld389kFeavMzinzUt2wRDimdNSNVavtxc+MeMgaFDzZoI118vVTiFX0kCsLjYWDO/CNyfCwAwrPMw6oXWY/xS13V9/iI5ueIJoFs3M/ts8uQTbp6RAT17atb+XsyylOHkleYyb9g8XhnwituLrhzH6WTLg28R3bMtfXO+4o/rnyFt8zxCw3xwMXY4zOyxXr3M6i9z55onn4Sal54QwtskAVhcTIxpmYCaDQWNCYtheOfhTMmYwr78fe7vmJxsqlWWlpq73+uuM2ssVlGTevFi6NvXycHCgxQP78ZFFynW3LaGfi36uX++Ksz7v4WkvXg7GWHp7PxmLa0/fATCalDZ1F3bt5t6+w8/bOr0lM9UEyJASAKwuPImoHr1apYAAO7qcRcljhLeXfGu+zuVLwxTvrD40KHmMWTt2r9s5nDA1dfnkRu0DUacwTs33cu0q6YRH+lhNc3iYvSPP/Hkk3Des+dwf6d5tM6cR8tBrTw7nitz5pi2/TVrzEozU6aYNRGECCCSACwuNta0TPTsWfPqwW0S2nBu2rlMWDYBu9Pu3k7lI13K+wFatIB9+2BwxQpepY5SHp7/EJkXnEarf41gxUPTuaXbLW4vu3icX35Bn3Ya9n7n8dYTuxk+HJ5b2p/4BB80+TidZpjSBRdA06am7f+666StXwQkSQAWFxNjhlbOmAHPPlvz/e/qfhc7j+xk5oaZ7u1QnhoMRMIAABwESURBVAAqDzkKDTUXzvnzeX3iHlIGTOOFn1/gtgH9WXH/bNo1bFfzwMCMLrrzTvTpp3Noey4XO6Zz1T3JTJzo+SqLLs938cWmfMPQofDLL1K7RwQ0SQAWV14PqKZzAcoNPnUwqbGpjFs6zr0dqkoAgJ4wgS/6v86omxM4uDmVyRdNY8LgCZ539ObnQ4cO6Lfe4ttWd5Oav46uD1/Aa6/5aDXErVvNxLY5c2DcODPEUxZkEQFOEoDFlZea+fhjs1JjSUnN9g+yBXFH9ztYsG0Ba/eudb1DUpJpDqmUAA4WHOTsrTlcyWekq2Ws+yyEoaddVrNAyu3YYd6jonA++BCPnr+EQRv/wyNjoxkzxkctMUuWmFE+e/bAd9/BnXdKk4/4W5AEYHHlTwC5ubB6tVllsKZu6noTkSGRDPtiGHvy9lS/cXAwJCbCrl1orZm0ahJpo0by42v/pEnaXr6NuIJT77/PlPisicxMUzK5RQv4+We0hlGb7ubZOd149ln4979r/udyy4wZZqRPVJQZtnT22T46kRDe55cEoJS6Uim1TinlVEql+yMGYZQ/AaSmmvdFi2p+jPjIeKZdNY2NBzdy+sTT2Zy9ufodkpNZd2gjZ006ixtn3khyaCs6drazYmkK9d95wQTx0kvunXzzZlOv+dRTzXyCUaOgTRueecas0Xv//T68+L/1Flx+uSlM9Ouv0Latj04khI9UtUyYr19AW6A1sABId3c/WRLS+9auNSsMfvqp1klJWl93nefH+jXzVx3/fLxOfDFRL9+1/LjP7Q67nvPnHP2Pu5N18GPo+s8m6neXv6sdTsdfV1C8/36zrKLWWu/Zo3VmZsVnTqfWhYVlB7Rr3bCh1kppPWyY1tu2aa21fvNN82caNsz7KzMe9eKL5iQXXqh1fr6PTiKEdxBIS0JqrdcDng/rE15T/gSQm2vKLSxcaEYFefJP07NpTxaNWMSAyQM4a9JZnNP8HBIiE0iITMDhdPDZ75+x88hOGsSFMWJ5Kgu3bCK8WTC204DK53v55Yqvx441y3E1a2ZG2Rw+bEo1bN5sirR9/LG5827aFIAvvzRN8EOGwLvv+qDDV2uzAPuTT8KVV5qnDp8MKRLC9wJ+PQCl1EhgJECzZs38HM3JpzwBHD5sJqvGxpq1AcLDPTtem4Q2LB6xmHvm3MPm7M2s2L2CAwUHKHWWMrDlQF4d8Cr9P1zDoK8vYHt4EC1auDjgffeZC/2ePaZMcvmroMDMXqs0s3bDBjPkvls3+PRTs9KZV2kNDz5omqduuMFkGFmXV/ydVfVY4I0XMA/IqOJ1caVtFiBNQH5lt5uWjMcf9905nE6nLnWUaq1Nk8xlXTZrG3Y94519XjtHTo7WrVubFqEdO7x22ApOp9b/+pf5y7rzTh+2LQnhfdR1E5DW+lxfHVt4T1CQuZE+fNh873SaKg3eLE2vlCJYmV+1Rx+F6avSeI17uaTJAKD2K2w5naaw5ubNMG+eb4p58tRT8OKLpnz1G2/IME9xUpBhoILY2IqJYMOH12qZ3mppbWqu3TaimHt43dTJ8YKnn4avvjILa511llcO+VcvvWQWar/hBrn4i5OKv4aBXqqU2gn0Br5WSn3rjziEERNT8QTQrZuZ1FrTwnCuOBzmuvn44/Dmu2Go1FQz8aCW5s0z1+bhw+Guu2of53EmTDDDTK+6yke9ykL4j19+m7XWM7TWTbXWYVrrJK31AH/EIYzKTwDld9A//ui942/ZAh07mqHyUHYD3alTrZ8ADh0yN+Vt2sCbb/rgxnzq1IohRZMnS4evOOnI7Yz4y7rAnTqZhLBwoXeOffiwuX7u2QPxlSs5d+oEf/xRsRiBB+64w/RXfPSRD8ruLFpkOhb69PHRkCIh/E8SgCA2tqIJKCjIrFfujScAux2uvho2boRp06BV5dL7nTubdqH16z069iefmBL7jz8O6d6eS75hA1x0ETRvDjNnQoSHBemECHABPw9A+F7lJwCAf/4T8vI8nxBW+Thz5sA778A55xzzYadO5n31arNWbg1kZpq7/969fVDmYc8eGDjQ3PHPnn3MY4sQJxdJAOIvncDgnXpmdru5UI8aBbfcUsUGLVua2WY17AfQGm680awo+eGHprac1xQUmIVp9u83bWAuZ6kJ8fcmCUAQG2vK5zscFf2cS5aYTtYBHnTPa20uzFOnmq+rFBQEHTrUOAF89BF8/70ZnOPVtVa0hhEjYMUK0+zj9XYlIQKP9AGIo+UgKjcD/d//mUqaJ7yAn8BPP0GPHmbFR5vNxcCZzp1NE5CbJ8nONs1KvXrByJE1i8ulsWNNZ+/YsabXWggLkAQgqlwV7OKL4fffTa01d61caVpQcnPdHDTTqRMcOFCxQLwLjzwCBw+au3+vDsefOdMc/NprTa0fISxCEoCo8gng1lvNjOA774Tt210f488/Td9pbCzMnQsNG7px4sodwS4sWQJvvw333GNWLvOajAxTQa57dzPRS2b5CguRBCCOPgFU7ggOCjLt7VrDsGGmf+BE1q2Dfv1MTZ7vvjOVm93SsaN5d9EP4HDA7beb+kRPPunmsd2RkwOXXALR0WZlLxnuKSym2k5gpVQDN47h1FrneCke4QdVPQGAGQTz5puunwCaNDFj/F95BVq3rsGJ4+PNzi4SwIQJpm/2008rYq218gpyO3aYET9NmnjpwEL8fbgaBbSr7FXdc3EQIIX6/8YqrwlwrOuuq/j64EEzvNPpNPMEJkwwfaZxcTB/vocn79y52gSQnW06pM8916y/4jVjxsCsWTBunJlQIIQFuUoA67XW1c7SUUqt9GI8wg+q6gQ+1ty5ZvnbvLyKn4WFmUVkzjijFifv1Mm0G5WUVLmy1tixJjG98ooXm+dnzzZTiK+/3swoE8KiXCWAXm4cQ26f/uaqewIod8op8NxzZvSNzWYuxmeeaQqx1UqnTmZW14YNFZ3CZXbsMNWXhw2r6C6ota1bYehQc8C33pJOX2FprhLANKXUJ8AXWuv8qjbQWntezUsEhKgoc1Gv7gnglFPMiCCvK7/or1lzXAJ47DHz/tRTXjpXUZF5jHE6Yfp0H1SQE+LvxdUooP8Cg4GtSqlPlVKXKKVkBeyTjFLHl4OoM61bm6afY/oB1q41pR7uvrsGo4pcuftuM1nho49MRhPC4qpNAFrrmVrra4BUYDowHNihlJqolDqvun3F30vlNQHqVHAwtG9/3FyA0aNNUho92kvnef99M85/9GiZ6StEGbfmAWitC7XWn2qtLwXOB7oCc3wamahTx1YErVOdOpk7c6cTMKMyv/7aXKsbuDMQ2ZVVq0xnb79+XmxPEuLvz60EoJRKUkrdrZT6GfgCmAt082lkok5VXhOgzg0caCpwLliA1vDww2ZY/j33eOHYOTmm3b9BA7OIgFfLhwrx9+ZqItgtwDVAa0wT0INa65/rIjBRt2Ji3C7J430XX2wy0KRJLFD9WLwYxo/3wsRch8OM+Cmf7JWY6JVwhThZuLod6gM8B8zTWjvrIB7hJzExZuUuv4iIgH/8AyZP5pkdE2nUKJgRI7xw3Mcfh2++MdmkTx8vHFCIk4urTuAbtdZzAa2Uuk4p9RiAUqqZUqpHnUQo6oTfOoHLDR/O4oLOzF8YzL/+ZdaKqZVp08xs35tuMoWEhBDHcbcY3JuYCV/XlH2fC4z3SUTCL/w2DLRc7948EzmWhJAcbr21lsfKyIDhw83CAePHy2QvIU7A3QTQU2t9J1AEoLU+BMh8gJNIbCwUF5uXPyxfoZhdcBb3lz5P1L6tnh/o4EHTpxAdbZ4CwsK8F6QQJxl3E0CpUioI0ABKqYaA9AmcRE5UEbSujBkDcTFO7uRNMwPME3l5cOGFkJVlZvo2buzdIIU4ybibAF4HZgCJSqkxwCLgWZ9FJepcfLx537+/7s+dkWHK8d8zykZM/+7wwQdH5wS4rbgYLr0Uli6FKVOkwqcQbnB3ItjHwIPAWGA3cInWeqovAxN1KzXVvLuz+pe3Pf881KsH994L3HCDKdi2aJH7B3A4TN3qefPgvffMIi9CCJfcnhWjtd4AbPBhLMKP/JUAsrLMDftdd5XN+r30UpMNJkyAvn1dd+A6HHDbbfD556Zm9A031EXYQpwUqn0CUEqtcHUAd7apYp8XlVIblFJrlFIzlFJxNT2G8K7kZDNJtq4TwLhxprXn6KzfqChTdnTKFLjvvuqbgvbsgQEDTI2fRx4x2wsh3ObqCaCtUqq69foUEOvBeb8DRmut7Uqp54HRwEMeHEd4SVCQqbq5bVvdnTMvz5Tkv/RSs/zkUc8+a0o3/+c/Zmzqf/97fAmHuXPNgi65uSYBeGXmmBDW4ioBuLPcRzXLhVetbHJZuV+BK2p6DOF9qal1+wTwwQemVM/99x/zgc0Gr75qeqYfe8xs9PzzsHOnyVBLlsDbb5sqovPnm3chRI1VmwC01nVxORgBfHqiD5VSI4GRAM28VhheVCU11dxY1wWHA157DXr2PMGAHaXMYsBxcaZ96IsvKj6z2WDkSJMkZFEXITzms9KISql5QKMqPnpEaz2zbJtHADvw8YmOo7V+B3gHID09XfsgVFGmeXPYtcuMqPT1/KlZs2DTJjP+v9p+3rvvNuWit2wxAaamQtOmVa4fLISoGZ8lAK31udV9rpQajlltrL/WWi7sAaB8JFBmJrRs6dtzvfKK6XO47DI3Nj7rLPMSQniVuxPBAFBKRZXNCK4VpdRATKfvRVrrgtoeT3hHXQ0FXb4cfvzRjPuX8vxC+I+rYaA2pdS1SqmvlVL7MPMAdiul1pUN5Wzl4XnHAdHAd0qpVUqptzw8jvCi5s3Nu69HAr3+uhnqf9NNvj2PEKJ6ru6/fgDmYYZpZpSvCaCUagCcAzynlJqhtZ5ck5NqrX3cwCA80bSp6V/15RPA/v1miP/NN5sCdEII/3GVAM7VWpce+0OtdTYwDZimlArxSWSizoWEmPppvkwA770HJSVmrpcQwr9cLQhTWtYMlFHdNt4PS/hL8+a+awJyOEyFh3POgXbtfHMOIYT7XHYClzX7rFZKySB8C/DlZLBZs8zyvHfd5ZvjCyFqxt0xGMnAOqXUEiC//Ida64t8EpXwm9RU00Zvt3t/hM64caaf4SL5rREiILj7X/xJn0YhAkbz5qapJiurYlioN2zYYKo1jxkjQz+FCBTV/ldUSiltLHS1jfdDE/5QeS6ANxPAm2+aybs33+y9YwohasdVH8APSqm7j23/V0qFKqX6KaU+AIb7LjxR13wxGSw3FyZNgquugsRE7x1XCFE7rh7GB2KKtX2ilGoB5ADhQBAwF3hVa73KtyGKulReb8+bI4EmTzZJQIZ+ChFYXFUDLQLeBN4sG++fABRqrXPqIjhR9yIiICnJe08AWpvO327dTOVPIUTgqMmSkKWY9YDFSc6bQ0EXLoTff4eJE12v7iiEqFs1KgYnrMGbk8HGjzdr/V59tXeOJ4TwHkkA4jipqWbCVnXL8bojKwtmzDCrNUZEeCc2IYT3uKoGOk4p1aeughGBITXV1OvZu7d2x3nnHZNEbr/dO3EJIbzL1RPAn8DLSqltSqnnlVJd6iIo4V/eKAtdUmKW7R00CNLSvBGVEMLbXBWD+4/WujdwFpANvK+UWq+UekwpdWqdRCjqnDfmAkyfbp4gZOinEIHLrT4ArfV2rfXzWuuuwLXApcB6n0Ym/MYbCWD8eHPnP3Cgd2ISQnifWwlAKRWilBqilPoYmA1sBC73aWTCb6KjzcgdT5uAVqyARYtM279NhhkIEbBc1QI6D7gGuBBYAkwBRmqt86vbT/z91WYuwAsvmCRyyy3ejUkI4V2uJoI9DPwP+GfZKmDCIlJTYePGmu+3ZQtMnQoPPCBLPgoR6FyVgjinrgIRgaVFC/j2WyguhrAw9/d7+WUICoJRo3wXmxDCO6SFVlTpvPOgsNAkAXft329KPlx/vVlbWAgR2CQBiCqdey4kJMD//uf+PuPGQVER/OtfvotLCOE9kgBElUJC4Mor4csvIS/P9fb5+SYBXHwxtGnj+/iEELUnCUCc0LXXmmagmTNdb/vee5CdDQ895Pu4hBDeIQlAnFCfPmaBGFfNQPn58NJLcMYZ0Lt33cQmhKg9SQDihGw2U8Z57lw4cODE240eDZmZZsF3IcTfhyQAUa1rrwW7HT7/vOrPFy6EN96Ae+6BM8+s29iEELXjlwSglHpaKbVGKbVKKTVXKSWDBgNUp07Qrl3VzUB5eXDjjdCyJTz7bN3HJoSoHX89Abyote6kte4CzAIe81McwgWl4Jpr4KefzCIxlf3736Ze0PvvQ1SUX8ITQtSCXxKA1vpIpW+jAO2POIR7rrnGvE+caJLA5s2mSWj8eDPj94wz/BufEMIzSmv/XHuVUmOAYcBh4Byt9X5X+6Snp+tly5b5PDZxvF694Lff/vqzVq1g1SqIjPRPTEII9yillmut04/7ua8SgFJqHtCoio8e0VrPrLTdaCBca/34CY4zEhgJ0KxZs27ba1OkXnjsjz/ghx/MBLHQUPPerx8kJvo7MiGEK3WeANyllEoFvtZad3C1rTwBCCFEzZ0oAfhrFFCrSt9eBGzwRxxCCGFlrtYD8JXnlFKtASewHbjNT3EIIYRl+SUBaK1lOUkhhPAzmQkshBAWJQlACCEsShKAEEJYlCQAIYSwKEkAQghhUZIAhBDCoiQBCCGERUkCEEIIi5IEIIQQFiUJQAghLEoSgBBCWJQkACGEsChJAEIIYVGSAIQQwqIkAQghhEVJAhBCCIuSBCCEEBYlCUAIISxKEoAQQliUJAAhhLAoSQBCCGFRkgCEEMKiJAEIIYRFSQIQQgiLkgQghBAWJQlACCEsShKAEEJYlF8TgFLqn0oprZRK8GccQghhRX5LAEqpFOA8YIe/YhBCCCvz5xPAq8CDgPZjDEIIYVl+SQBKqYuALK31an+cXwghBAT76sBKqXlAoyo+egR4GDjfzeOMBEYCNGvWzGvxCSGE1Smt67YFRinVEfgeKCj7UVNgF9BDa72nun3T09P1smXLfByhEEKcXJRSy7XW6cf+3GdPACeitV4LJJZ/r5TaBqRrrQ/UdSxCCGFlMg9ACCEsqs6fAI6ltW7u7xiEEMKK5AlACCEsShKAEEJYlCQAIYSwKEkAQghhUZIAhBDCoiQBCCGERUkCEEIIi5IEIIQQFiUJQAghLEoSgBBCWJQkACGEsChJAEIIYVGSAIQQwqIkAQghhEVJAhBCCIuSBCCEEBYlCUAIISxKEoAQQliUJAAhhLAoSQBCCGFRkgCEEMKiJAEIIYRFSQIQQgiLkgQghBAWJQlACCEsShKAEEJYlCQAIYSwKEkAQghhUX5JAEqpJ5RSWUqpVWWvC/wRhxBCWFmwH8/9qtb6JT+eXwghLE2agIQQwqL8+QRwl1JqGLAMeEBrfaiqjZRSI4GRZd/mKaX+8PB8CcABD/f1tUCNLVDjgsCNLVDjgsCNLVDjgsCNraZxpVb1Q6W19k44xx5YqXlAoyo+egT4FRO8Bp4GkrXWI3wSSEU8y7TW6b48h6cCNbZAjQsCN7ZAjQsCN7ZAjQsCNzZvxeWzJwCt9bnubKeU+i8wy1dxCCGEqJq/RgElV/r2UiDDH3EIIYSV+asP4AWlVBdME9A24NY6OOc7dXAOTwVqbIEaFwRubIEaFwRubIEaFwRubF6Jy2d9AEIIIQKbDAMVQgiLkgQghBAWddInAKXURKXUPqVUQHU0K6VSlFI/KKXWK6XWKaXu9XdM5ZRS4UqpJUqp1WWxPenvmCpTSgUppVYqpQJq9JhSaptSam1ZeZNl/o6nnFIqTin1uVJqQ9nvW29/xwSglGpdqRzMKqXUEaXUKH/HBaCUuq/sdz9DKfWJUirc3zGVU0rdWxbXutr+fZ30fQBKqTOBPOBDrXUHf8dTrmwkVLLWeoVSKhpYDlyitf7dz6GhlFJAlNY6TykVAiwC7tVa/+rn0ABQSt0PpAMxWuvB/o6nnFJqG5CutQ6oiUNKqQ+An7TW7yqlQoFIrXWOv+OqTCkVBGQBPbXW2/0cSxPM73w7rXWhUuoz4But9SR/xgWglOoATAF6ACXAHOB2rfWfnhzvpH8C0Fr/CGT7O45jaa13a61XlH2dC6wHmvg3KkMbeWXfhpS9AuJOQSnVFLgQeNffsfwdKKVigDOB9wC01iWBdvEv0x/Y7O+LfyXBQIRSKhiIBHb5OZ5ybYFftdYFWms7sBAzlN4jJ30C+DtQSjUHugK/+TeSCmXNLKuAfcB3WutAie014EHA6e9AqqCBuUqp5WUlTAJBGrAfeL+s2exdpVSUv4OqwtXAJ/4OAkBrnQW8BOwAdgOHtdZz/RvVURnAmUqpeKVUJHABkOLpwSQB+JlSqh4wDRiltT7i73jKaa0dWusuQFOgR9mjp18ppQYD+7TWy/0dywmcrrU+DRgE3FnW/OhvwcBpwAStdVcgH/i3f0P6q7JmqYuAqf6OBUApVR+4GGgBNAailFLX+TcqQ2u9Hnge+A7T/LMasHt6PEkAflTWvj4N+FhrPd3f8VSlrLlgATDQz6EAnA5cVNbWPgXop5Sa7N+QKmitd5W97wNmYNpp/W0nsLPSE9znmIQQSAYBK7TWe/0dSJlzga1a6/1a61JgOtDHzzEdpbV+T2t9mtb6TEzztkft/yAJwG/KOlrfA9ZrrV/xdzyVKaUaKqXiyr6OwPyH2ODfqEBrPVpr3VRr3RzTZDBfax0Qd2ZKqaiyznzKmljOJwBKnGit9wCZSqnWZT/qD/h9oMExriFAmn/K7AB6KaUiy/6f9sf00QUEpVRi2Xsz4DJq8Xfnz3LQdUIp9QlwNpCglNoJPK61fs+/UQHmbvZ6YG1ZWzvAw1rrb/wYU7lk4IOykRk24DOtdUANuQxAScAMc70gGPif1nqOf0M66m7g47Kmli3AjX6O56iyduzzqJtyMG7RWv+mlPocWIFpXllJYJWEmKaUigdKgTtPVErfHSf9MFAhhBBVkyYgIYSwKEkAQghhUZIAhBDCoiQBCCGERUkCEEIIi5IEIIQQFiUJQAgfUEq9oZRaoZTq7u9YhDgRSQBCeFnZTOBEzOSmgClXLcSxJAEI4SGlVHOlVGGlmdwAaK3zMbOpFwCvl20bUbboSYlSKqHuoxXieJIAhHCTMo79P7O5rGpq5e3iMTXkcwEHgNa6sGy7QKkrL4QkACGqU3aXv14p9SamNow7tdcfxdSTXwe082V8QtSGJAAhXGuNWVK0q6sVq8oW9+kDfIqpINne59EJ4SFJAEK4tr0G6yE/AzylTZVFSQAioJ305aCF8IJ8dzZSSnXB1Gc/Qyk1HggH1voyMCFqQxKAEN7zPDBEa/09gFIqCVNLXoiAJE1AQniBUqofEFV+8QcoW+IwSinVwH+RCXFi8gQgRDW01tuADm5sNx+YX8XPY30QlhBeIU8AQnjOAcQeOxGsKuUTwYAQwOnzyIRwgywJKYQQFiVPAEIIYVGSAIQQwqIkAQghhEVJAhBCCIuSBCCEEBYlCUAIISxKEoAQQliUJAAhhLCo/wcXYNyhC6dGMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "459\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEQCAYAAAC5oaP8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3iUZfbw8e+dRgqkEEIvoQgCoiBFQVQURVRAFFHsHf2tffW1rK51dV0b6rqKqIgKiiIqiAiICgiiUqR3lIQeSICE9Jm53z9OQgIkmclkJjMw53Ndc02SmXmekyjPee52bmOtRSmlVOgJC3QASimlAkMTgFJKhShNAEopFaI0ASilVIjSBKCUUiFKE4BSSoWogCYAY8xYY0yGMWZVIONQSqlQFOgWwDhgYIBjUEqpkBTQBGCtnQdkBTIGpZQKVYFuASillAqQiEAH4I4xZiQwEiAuLq77iSeeGOCIjm35+bBmDUREQHw8tG5d9prFsnTnUhrXbUyzes1YlbGK2MhY2iS18e4kbdrAwYOwdy906+bbX0Qp5bElS5bstdamHPlzE+haQMaYVGCatfYkd+/t0aOHXbx4sd9jOp5lZcFbb8EHH8Cpp8KkSYe//s7id+jauCunNT+NPu/3ITYyltnXz/buJMOHw8cfwwsvQHExGOO7X0Qp5TFjzBJrbY8jfx70LQDlW/Xrw+OPy4W/uPjo12/vcfuhr5Njk9mWvc37kwDExoLTKSeLivIyaqWUPwR6GuinwEKggzFmmzHmlkDGEyq2lVzTK0oAuw7uYn76fACSY5LJzMv07iS7dsHmzZIAAPLyvDuOUspvAj0L6CprbRNrbaS1trm19v1AxhMqOneGjIyKE8Dbi97mrA/OoshZJAkg38sEMHw43HorxMTI95oAlAo6OgsoBMXHg7UVJ4DUxFQslm3Z20iOTSavOI8CR0H1T9KiBaSnl7UA8vNrFrRSyuc0AYSghARwuSpOAK0SWwGwZf8W6sfUByAr34ulGi1bwtatEB0t32sLQKmgowkgBFWVAFITUwFI259GckwygHfjAC1ayAmKiuR7TQBKBR1NACEoPh4cjooTQPP45hgMW/ZvITm2JAF4Mw7QooU8Z2fLsyYApYKOTgMNQXfeKVP1K7omR4VH8dWVX9G5YWdyi3IBL1sAPXvChx9Cw4byvY4BKBV0tAUQggYNglatpBVQkUtOvIR29dsdagF4NQbQpAlcfz00by7fawtAqaCjCSAEZWRIz0xFXUAAa/esZdLqSYcGgb2eCrpkCaSlydeaAJQKOpoAQtCYMTBzZtn47JHGrxjPVZOvIio8iuiIaO8Xg117LYweLV9rAlAq6GgCCEHx8fJcWQJITUzFaZ3syNlRs8VgLVrAzp3ytY4BKBV0NAGEoIQEea4qAQCHZgLVKAFs3y5fawtAqaCjCSAElbYAKhsDKF0MVroWwKtBYJDFYLt2QXi4JgClgpAmgBBU2gKoLAG0TGgJlGsBeDsGULoWICZGE4BSQUgTQAg66SQYPFiqNFckOiKa3279jb/1/Bv1o+t73wU0cCD88IPUA9IEoFTQ0QQQgho2hK5dZR1AZfsB9WrWi+TYZJJjpQvIq42DmjaFc8+FuDgdBFYqCOlK4BDkcJTtCeB0yvaQR5qfPp9lu5aRHJOMw+UguzCbhOiE6p3IWvjyS/laWwBKBR1tAYQgp1O2hITKVwN/u+Fb/j7z78TXkRFjr7qBjIE77pB9gTUBKBV0NAGEoDp1ZGIOVD4Q3KNpD4pdxWQXSjG3Gs0Ecjo1ASgVhDQBhKjSjboqSwCnNz8dgPQD6YCXBeFAZgIVFekYgFJBSBNAiHKXAJrFN6NFfAs27dsE1KAeUMuWUFioLQClgpAmgBBVulNjZQkApBWQtl+KudWoBVBcLOMASqmgogkgRF19tTxXlQDeHfwuS29fCtSgBXD99XDlldoFpFQQ0mmgIapzZ3muKgGUTvtMjE70fhC4USNo3BgKvNhYXinlV9oCCFGlRTqrSgAAd02/i3AT7n0LoKAAVq3SLiClgpAmgBD100/y7C4BLNqxiCJnkfdjABERcjKn0/3JlFK1ShNAiKpbV54rKwld6vRmp5NbnMvevL3enSgioqz6nI4DKBVUNAGEqHr15Dk7u+r39W7RG5d1sTNnp/cnqy9bS+pUUKWCiyaAEFWaAPbvr/p9pQvCsgq8HAQGSEmRZ20BKBVUApoAjDEDjTHrjTGbjDGPBDKWUFPaK+MuAbRKaEXz+OYUOApwuCopHOROw4byrAPBSgWVgCUAY0w48D/gQqATcJUxplOg4gk1ffrIc2Ji1e8zxvDwGQ8DNagHdM018qxTQZUKKoFsAfQCNllr/7TWFgETgUsCGE9IadZMnl0u9+9NjkkGYE/uHu9OVtoCyMnx7vNKKb8I5EKwZsDWct9vA04LUCwhp3QPgLVr3b83p0gu3N9u/JbODTtX/2TR0fL8ww+yQYw6JrisiyJnEUXOIupF1cMYQ1Z+Fvvy91HsKqbYWYzD5cBpnfRo2gOA1Rmr2Za9DYfLcegRHhbO0BOHAvDjXz+yOWszDpcTl3VSWGQJK0xkROdrKS6Gyaum8tfebdRrlEVYZBE5WTGwvzWXnzgChwMmrfyK7Qd207jTJsLrFJC1rQFmd1cGtRuK0wmTVn7J3tz9pJ65kLCoAjLWtyZq55mcm3o+Tid8uforsgtzaXfxVAgvYveyrsTsOo/ezfricsHXa6ZS4Cii/VXvYq1lxy9nE58xkG6NTsXldPHNhu+wpoC2172BxbJj9mCSdg2gY3xbHA4X3+/8hfDobFJveFVen3otybv70aZOE4qcDhYcWE1E/G6a3/ASADs+/xsNM/rQLCyZfGcRSwvTiEr+i0bXvwDAzo8ep9GebtxxiYO/vdSBOhF1fPrfOJAJwFTws6O2nTLGjARGArRs2dLfMYWM0mmgCxe6f2+nBtIzt3jHYu9O1qSJPC9Y4N3n1SEFjgL25u0lpzCHg0UHyS3OJbcolzNbnUl8nXiW7VrGrM2zyCvOI684j/zifPId+fy7/79JiUvhs1WfMWbpGAocBfIoKiIvpw5fXTqbSGcir897nw8XfYEjaS2u+DTIbQCrh/P82a/jLIpkyqpFLE5bDSdPgKZLYXdnzJx/MTC1B4WFsGpHFBkHkuGCv0Pqz7D5PMzXH5MUKctAcgv64iruBzecC63nwoqr4csJ3H/oNxwiTyO7Q9OlmMV3YKc9ymuHXr8UgPi/9yCy0Rby5/4fedOH8smh1y8D4K/ENwhPziT7+5PJ/vF8ph3x+eyOT0PdQrLnnc+BhX2Yh4NwnDjCB4BxUtznLhyRERQuHsqfy1LZ6NiLwXIgpjvhkTms6fsnxSac8FXx7F+fwL7CDAyW/HptiYrezcoNe3ESRt0/Y8n802ALdmKNxSYkE5Gdw4qVLqw1JG+PpDA9lwMFB3GEWeokRRDncrFuaX2wYTTMdBG1aw/zRm8i7Iw87r38dJ/+/xTIBLANaFHu++bAjiPfZK0dA4wB6NGjhxf7EqqKlM4C2uvB9P7mCc0BWJWxyruTJSXJ859/evf544jLusjKz2Jv3l725O4hMz+TzLxM+rfpT2piKkt3LuXFBS+yv2A/+wr2kV2YTXZhNpOGT6JPiz58seYLrvvqurIDWiC7GZ9cOItG4Z34fPEe3pm/EVJWE93mD+oUN6Jo6ig2/TcWRz6kZwwgI+ssWl00iRMGzqJ4Tyqz/vEm3f5ResBbgFu44O6p9Om3jD2bm/HmS7fwj+nyanjE+UTXOYfbLjuNsy7eSdraZN78oTt798o+E20bNqVNoySuu+A/nNDzIBuWwLdpW0myaUQ4cnE4sjGOg3TpcDZFffqTE+2AxU+RUJhPVHEBkc58YooK2WoS+K3NebTYsYnLuYhIiomkmHCcRODghbUOpmefxLn7FvAqnQnHeegRgYNrvr6Y+XGdGLx7E8+Ft6XYROAkvOQoEdw26hU2mo5ckLeQO8IvpshGU2zrUOyMpogoHn98BXtoyFnMZQCv4iBCPpsfiSM/gtHPriePOLqzmC7MwYEc35ETgTMnnCkvr8dJBCmsoznbJDIbTsJ+gysnhayxvxEe5iTJpJMYsQ9nApgwS1trsTkn4PxuHOHGQaLdR3RiAX+Z9tyUU9E9c80Yr/Z69cWJjYkANgD9ge3AIuBqa+3qyj7To0cPu3ixl3eh6jDFxRAVBa1bu78uZxdmk/BCAmEmjLT70mge37x6J7NWdqAJC5OpoJGR3gce5LILs5m1eRbbs7ezI2cHOw/uZNfBXdx3+n1cdMJFLEhfQN8P+h71uc8v/5zhnYczP30+N381kviwxtRPjCAxOpHtPw7hpJgLcGan8Nf2XP7cmkv3s3Zzw71biDb1GNCh31HHu/9+y6uvGg4ehG7dZNZXYqIk/vh4uOIKuPhiGZYZN062ba5bF+pEO3BGZtOgRSaRCXvZm72fvWt3YvalE7FvG1GZuwnfm8mGFrEs7BCLKyuTZ19bTt2DxcTnOkjMc1HHCQ+fBy/2hdR98Nfrh8fmMHBn95MZ07wrbfcaPl7xAwfD6pJHXXKJJ9fG83Hc5Sy2fWiUd5AhebPILk4knzjyiSGXOJbRlUwaEEsuDckgnxgKiKYAuYBHRFnqRDupE+MkOsZFdLQlOsYSHQ0xsRATDTExhuhoQ2x0GDHlH2EO6tgCohz51HHlE+XII7JYHlHFuUQW5RJZeJDIwoNEFR0ksiBHvi/IIaogm8i8A0TlHyAyP5vI/GyiinJK0o6DcDwYdCsVGSlZNTpaHh99BOecU73/IUsYY5ZYa3sc+fOAtQCstQ5jzF3ATCAcGFvVxV/5VmSk7Ni4b5/799aLqke4Ccdpnby75F2ePufp6p3MGLnCHDwIa9bAKad4F3QQKHAU8MWaL/hz35+HHukH0rn3tHu5v/f9ZOZlMnzScACiwqNoUrcJjes2psgpS67bJ7dn1PlvEONsQpumCTSIbcDE0al88a94RqVBenpfdu5cwyWXlG2n3Ogu+HWvjKU3bhxH+2Zx9OvUkMEdugDw/vtyca9fv+yRnCx3i3Fxlt9X7Gd7znZ2Hdx16DHn4G4+/zqDPbl7aLN5E/E7s4jfm0ODfUU0y4HvWsELZ0K4E4qePXq2yNi+9dkU0Zu6OZ1xFuwkLaY+2XH1OWBS2OdqyOZtZ9NsQl/y9odzRuQ69hbX5wAJZBNPvo2BxQYWw2agD5KYEhIkOZU+Lo4vaanW60KzevJ13bry/Le68nVcXBx167Y+lMBiY+UREW5hfzbs2SPN3MzMsuesrLLHrn0yF3pfyfP+/VK2xJ3w8LJgJBB5NKoLcfUhrkXZz+LiygIr/3VMzOFfx8TIhb70uXTbPj8KaDVQa+10YHogYwhlkZFyB2itXKMrY4yhQWwDUhNTGdxhsHcnq19fmh07dwZ9AtievZ1VGatYlbGK1XtWs27vOs5rcx7PnPMMYSaM67+6HouleXxzWie25sxWZ9KufjsAWiS0YPkdy2larynJMckYY/jmG1gwDj5YC+vXp/Dnn3fTtSv8/ruc767ZsGsXpKbC+edD8+bQvXtZPGvWyAW+ouuB0+XkvGHb+XPfn/y17y/mZKWR9lca6QfS2XpgK9tzttNtUx4nZUDbfdB6H5yzH9Y1DOPRG5qREpfCRx+l02BfIS4DuUl1yU1OJjqlH1GN7iZ3b30+6j+NjMIU0vOasXF/Y1ZlNGLH/HowX/6n6VkunqQkWffXoD50bwDJyZCc3Jvk5LLklJRU9khMlIu9R9c6a+HAAdixQ/4/2rlT/nC7dsHu3WWPjAy52Fe24XVExOHBNGwI7dtLMImJZZkoIUEepc2m0ue6deXOvKp/NMcILQcdwuLj4YILPHtvcmwyzeKbHZrtUW0NG0LHjjBwoHef9wNrLZuyNrF4x2Jc1sU1J8t6he5jurM7dzcAjeIa0TGlIw3jZCprVHgUG+7eQPP45kRHRB86VlYWzJwJS5ZEsGzZyezeDXPnymsffghTp8o1pnNnGDoUunQpi2P+/KqvJfXrW3bn7mbtnrVsyNzA+sz1bMjcwMasjWzZvwVnURHtM+GkDOi0Fy7dF0OdmDg+eKA/g9sP5sEJX9BkTTquqEgcLZtjTmpLu+7n0rTro2zcCB8NWcGanUn8sasJm7ZEkL0R2IhMzAYiI++mWTOZOty0LQxvKhW+mzSRat+ljwYNpFvRawUFsGULpKVBenrZY+tW2L4dtm2ruJxITExZEKmp0KuXZKGGDSWolBTJRA0ayEW/Xr3j4uLtC5oAQljpvwNP/i0kxySTlZ/Fyt0rmbV5Fg/0eaB6J0tKcr/suJa8u+Rdvl7/Nb9u+/XQ4rYuDbscSgCjB40mKTqJzg070yC2wVGfb5PYjrVrJZ+FhcETT8Czz5a93rq1NHKKiuSCOGYMfPpp5UMf5f/+OYU5rNi9guW7l7Ni9wpW71nNmj1rDsUZUwSn7Y3i7LyGxAw5nUtPvJSRL3xPmx9k4x5rDKZ1E+jYiwsu/5TsbNjw8C1M/7Mev21rxqo1Yaz5HQ58D8hMQ+rUOZk2baBNG+hzlsSfmiqbubVoIdfRMF+tGMrNhfXrYcMG2LQJNm6U57/+KqtRXiosTLJO8+byB7344pIs1FSyT+lDL+he0wQQwqKi4I8/YPly970yybHJbMjcwDcbvuGxHx9jUPtBdGjQwfOTJSXJvgDdusGiRWULEfyo0FHIL1t/Yfafs1m2exnTrpqGMYbft/9O2v40Lj3xUk5rdho9m/Wkc0rZ+obSOevlbdoEs2bB99/DvHlyx798OZx8sozLxcbKjWf37mVlNkqV1sI7Uk5hDkt3LmXRjkUs3rGYxTsWs3nf5kOvJ0UncVLDk3jE9uXCpXtIXbeTuE3pGFcRsA1ef0cO7joXbthNUYcuLMvvwMJlMfz6Kyw+QeKWhfZyA9y5syzM7tQJOnSAE06Q66vPu5sLC6XvasUKeaxZI4/09MPf16wZtGsnTdHWrcuyT6tWcqGvhf9PQpn+dUNYdLRcxH7+2X0CSE1IZdbmWdzc9WaemvMUoxePZtTAUZ6fLClJ7v6WLZPVZ+X7QHxszpY5vPTLS8zZMoe84jzCTTg9mvYgMz+TBrENeGfwO4SZqm9pi4rkGlavnlz4S7vKUlPhkkvgzDPLVlOfc45nkzO2Z2/n5/SfWZC+gPlb57Ni9wpcVmaFtExoSc/G3Xm43oX03lxIm+VbiRn1X0y7dvDee/DzQ3DaaTDieskyPXqQW6c+v3wPc+cNYM4cyaul5b1btICePeHGGyXndutWthzD54qKYOVKCWDRIli8WC72pX3w0dHSXOrbV547dpT+sDZtZFBUBYwmgBAWHy8t561b3b+3bf225BXnYbEM6zSMccvH8Vz/54iNjPXsZElJZf23S5f6LAG4rIvftv3G1PVTubrL1XRp1IWDRQfZmLmRm7rexAVtL6Bfaj/q1al36DOVXfwLC+G77+CLL2DaNHjgAfjnP6Vu0ptvwoABcrPqaW/DgYID/LTlJ2b/OZsf/vqBdXvXARAXGcdpzU/jsTMf4/Tmp9OrKIUGz42STJNZsvFOmzYy2NmuneyrfPPNuAhj6VIZa5j5siziczjk7r17d7jnHujdW/JEaXLyi5wc+OUXGbz4+Wf47beyOk/JyZJ5Bg2S5tEpp0gzoxZmtKjq0wQQwmJjpRvIkwRQOstl877N/K3H35i4aiKfrvyUW069xbOTJSXJ1SouDpYsgRtu8Dpul3Xxy9ZfmLR6EpPXTmZ7znYiwiJon9yeLo26cPEJFzOo/SCPj2ct3HcffPyxzAasXx8uvRTOPlter1sX7rzTs2NtztrMNxu+YdqGacxNm4vD5SAuMo6zWp3Frd1upV+rszllbzgRU6dBbEc45yKZqvjTT3DRRdC/P/TrJ10gSKNp1qwopkyB6dPlrSB39A88IC2PPn3KFvb5hcMhU5Zmz5Y+sF9/lZ+FhUkgd9whmadnT2kiaX/8MUMTQAiLjZUbM49aAEltAbnAXX/K9XRv0t3NJ45QWna0UydJANVkrSUrP4vk2GTyi/O5YPwFOF1OLjzhQv7T8T9c3P5iEqPlHMaDC9D+/XI9u/xyuV7t2QMXXgjXXQfnnVe9rue0/WlMXDWRiasnsmzXMvk1UzrxQO8HuOiEizi9+elErVgNn3wCX14pK++MkVv2yy+XWSo7dhy6cGZnw5SPYdIkud4WFMif76KLJMbzz5cJL3514ADMmAHffCPNoqwsia97d3jwQanpdPrpfs48yt80AYSw2Fi5idtxVAGOo7VKbEWYCWPzvs0YY1hw84LqFaYqLQdx5pnVKgv9176/GL9iPBNWTiAmMoY/bv+DuKg4Zl47k1ManXJY144nVq+G116DCRNkUfKWLXKzPWFC9W5cDxYdZOKqiYz9YywLt0lBpdOancarA15lSIchtK3fVuaiNyiZRfS3v0nf+HnnwcMPw5AhMpeyRGGRYepUmDgRvv1WuqNatICRI8vGHPy+gHr/fpgyBT77TDKPwyFdOoMGyQyc/v3le3Xc0AQQwmJipEdmwwb3740Kj6JlQks2ZW0COHTxn75xOh+v+Jjxl44nPKyKft7SBDB0qFzN3Phu43f8e/6/+Tn9ZwDObnU213S5Bpd1EWbC6Nvy6HIKVdmyRbp5pkyR3/uaa+D//u9QT4vHF/+lO5cyZskYJqycwMGig3RK6cTz5z7PiJNG0DqptWSVyZPh/Vulq2THDvnd331XRmGPuID+8QeMHSuNg6wsyQkjR8KIEXKD7bPpl5UpLJS7/I8+ksGFoiL5o9x3n/y3Ov107b8/jmkCCGGxsXIz7um/77ZJbQ+bpgiQfiCdiasm0iiuEa8NfK2ST1KWAPbtk073ggK5EpcocBQwbcM0zmx5Jo3qNiIjN4M9eXt4/tznuebka2iZ4F0l2Lw8+T3r1pWepyefhLvvrt6NrMu6+HbDt7y88GXmpc0jJiKGK0+6kpGnjuT05qdLl9P27fDM/fDBB9J90qYNPP54WWY56aSy37VA7vT/+18ZD69TR8YcbrpJbrJr5Xq7dKnEWpp5mjaVgY4rr5T5rNqPHxqstcfMo3v37lb5ziOPWBsZae3tt1u7bp3794+cOtIm/yf5qJ/f9919lqewt029zW7O2lzxhzdtshas/eADa5s3t/a+++z+/P124sqJ9prJ19j4f8dbnsK++dub1lprHU6HdblcXv9uGzZYO2yYtb16Wet0ys+Ki6t3jCJHkR27dKzt+GZHy1PYFq+2sK/+8qrdl79P3uByWbuv5Ot166yNirJ2xAhrf/yx7KTlbN9u7WOPWduggfwpOne29r//tTYz0+tfs3ry860dN07+KGBtnTrWXnmltTNmWOtw1FIQKhCAxbaCa6q2AEJYTIyU53nnHRlY7OBmXVe7+u3IzM/kQMEBEqLLVju9POBlil3FjFkyhrF/jCX9/nSa1mtKgaOAqPAowkwYNjGRtETIzVxL5xYtyF72G41ebkShs5DkmGSGdRzG1V2u5pxUmVBfZXdSFbKy5C5/9Gi5s37oIenKjoryfGDXZV1MXDWRJ+c8yaasTXRt3JUJl01geKfhRIZHSrGwiRPhhRdkmuYXX8gfb9euspZOOWlp8taxY+XvPWSIjP+ec04t3Wjv2AFvvCHrCTIz4cQT5ftrr60wXhVCKsoKwfrQFoBvvfSS3AiCta+95v79k9dMtjyFXbJjSYWvb8/ebj/444ND33d/p7vlKWzMv2Js7HOxlqewA55sY+1dd1lbt679769v2Plp863D6Zu7z+XLrW3Y0NrwcGvvuMPanTur93mXy2W/Wf+N7fJWF8tT2JPfPtlOXTe1rCVSWGjtu+9a266d/NFOPNHajz6q9Hh//mntzTdbGxFR1tLaXEkDyS9Wr7b2ppvk5GFh1l52mbU//CAtFxVS0BaAOlJpF3xMTPWmgm7K2sSpTU496vWm9ZpyY9cbD33/fz3+j63ZW8ktysXhcnDic6Pp2bOXzBl/803ucpwKLc+o8e/hcMjdffv20of+0EPQtWv1jrF+73runXEvMzfP5IT6J/DpsE+5ovMVhy8ae+kl6dfv3l0GeocOrXCUdt8+eO456eM3RqbJP/SQzOqpFcuWwVNPlY1433473H+/jEsoVV5FWSFYH9oC8K2xY+VGNjXV2iuucP/+7IJsy1PY5+c9790JU1Otve46a7OzrY2Jkdv0GnC5rH37bbkRz8nx7hg5hTn24e8ftpHPRNr4f8fbUQtH2SJHUdkJvvjC2nnz5Pu9e6399ttK76CLiqx9/XVr69e31hi5+d62zbu4vLJihdzlg7WJidY++aS1e/bUYgAqWKEtAHWk2JIqDs2ayWxAd+rVqUfDuIZHzQTyWGKi3B7Xqyd90J292GC+xN69cOutcpN7wQUy+7J0n2NPfb/5e2795lbSD6RzY9cbeaH/CzSqW7LC6rff4O9/l5IHV10lU1eTk2U1VgV+/13iWblSWiEvv1z9VojX0tLgH/+QGT3x8TIIct99ZYvvlKpMRVkhWB/aAvCtKVPkZnHRIs8/0/u93rbfuH7enfCcc6w94wzvPlvO7NnWNmkik25Gjapwwk2VDhQcsCOnjrQ8he3w3w52QfqCshfT0qy9+mr5wzRubO3771c5QyYnx9r77pMu9mbNrP3qq1rsYt+/39qHH5bZPDEx1v7jH9ZmZdXSydWxBG0BqCOVtgDy8z3/TLv67ZizZY53J0xKgnXryr7/4w+Zj36Lh/WEkCHrZ5+VksvffiulaKrjp79+4sYpN7L1wFYe7P0gz5zzDDGRZesR+Ppr2Yvx8cel476KUgfz5kmdtrQ0Wej773/LDbjfuVwyh//RR6WGxfXXy6BD82ru1axCnr/XGaogVpoAFi2Slf5H7sdRkbZJbdmWvY0Ch+flHA45clOYcePkyunBRjH5+dJ7ZIzMuly8uHoXf6fLydNznqb/R/2Jjohmwc0LeGnAS3LxnzlTLvwg8WzYIFmmkou/0wlPPy3TOCMjpSDm//5XSxf/1aulSt2tt8rU00WLZMsxvfgrL2gCCGGls83xddAAACAASURBVID27ZNKk3/95f4zbeu3xWL5a58Hbz5SUtLhu9Bfd52UHpg0qcqPbdsGZ50FV1whLYAGDapXRn7XwV0MGD+Ap+Y+xTUnX8OSkUvo3aK3zI8fPly2qXzlFTl4RESV03W2bZM+/qeeknISS5dKmXu/y8+XO/6uXWU/hbFjpQnSw8stOpVCE0BIK20BlN7obtvm/jOHqoJ6MxCcmCgXstIR5+7d5S52/PhKP7JqldS3X79eFk9Vd+HUvLR5dB3dlYVbF/L+kPf5aOhH1I2IlZViHTtKHZx//UtKg7o5+Jw5cv1dvFhuuj/6qJaKYS5aBKeeKqvJrrtOutFuuknLNaga0wQQwsqvAwDPN4YBKQtdbeXrAYFcwK69Vu5k09KOevsvv8idf+nXgwdX73TvLX2P/h/1JzE6kd9v+52bu90sdXvmzpVKcD16yLSdxx6TZcNVHes9WS3dsKHUFLr++urF4pXiYmlq9O4NBw9Khc6xY8sqjCpVQ5oAQlhpC8BamULpSQJIiU2hblRd71oARyYAkH6UhAS5EJfjdMrYcHIyLFhwWC01t5wuJ/fPuJ/bvrmNc1ufy6+3/spJCSdIZz1I5/3s2fI44YSqj+WU8ve33SYl8BcudF8ywyc2b5adXp5+WqahrlwppaSV8iGdBRTCShNAQYF0s8R6sLujMYZ29dvVLAGUH/Rt3RoyMqRYTznh4TLHPyGhepufZBdmM+KLEXy36Tvu6XUPr1zwChG/LZJB082bZaCjSRPpyHcjP1/KMk+dCnfdBaNG1dIe5V99JZv5hofL+Mjll9fCSVUo0gQQwqKj5TkvT26GPdU2qS0rM1a6f+ORKmoBgFz8XS6YM4fP957LL7/IxbZ9++odPiM3gwsnXMjyXcsZffFobm93Jdx9L7z9tsyS+fJLj3dGP3hQirbNmSMlHe66q3qxeKW4WBZ0vfyybK/4+eeyxaJSfqIJIISFhUkSqM46AJAEMHX9VJwuZ/WqdpauTD0yAQCMHs3Xd87i6rB+9O4TduR2AW6l7U9jwPgBbD2wlalXTeWipmdLX83OnbIBwL/+5fGI7YEDsvXi77/LPsHXXON5HF7bsweGDZNuqjvvlFlJbsYllKopHQMIcbGx0gL45BOZ4VJU5P4zbeu3pdhVzLZsD6YNlVdZCwD4tsmtXGEm0dMsYfq726t18V+dsZozxp5BRm4Gc8/9mItOuEjmiT78sFzFX3/d44t/Zqb0Di1eLDsj1srFf+1a6YNbtEj2pnzzTb34q1qhCSDExcRIAigqguXLZb9yd7yeClpJApg9Gy4bEcUpnZ18F30p9W4bISU+PbB051LOGncWjfYVsXl5P3r2vUJGjUHu/Lt7vnn9/v0y0LtqlawLGzbM449674cfZJZPXp7MTrr66lo4qVIiIAnAGDPcGLPaGOMyxuhKlgCKjZUuoNNPl+/nz3f/Ga+ngkZGyp35EQng4EFZ1TtzbjSJY16UIF5+2e3hlu5cyq2jzuH5GcUsejmb+pOnSxG0E0+sXlzI32DwYLkZnzq10ppvvvX++7IIrUULKT7Xq1ctnFSpMoEaA1gFXAa8E6DzqxKlLYAOHWS2zdy5MmGmKi3iWxAZFnlog/hqKVcOorhYcsLQoTLgGhaG3AEvWSJdIgC7d8sbS0sdWAuFhSzdt4YB4/qz7u2DJOdazHXXwTPPlO3yXg0Oh2yFu2CBbPQ1YED1f61qe+klqTV0wQUy2FsrdSSUOlxAWgDW2rXW2vWBOLc6XOkYgDGy6GruXLnGViU8LJw2SW1Yu3dt9U9YUhJ6zx4Zc5gwQX582L4qr7wic/VBKqy1bg1t28qigMhICtu35fxx/YmLiaf4ww8w6emyNNeLi7+1Msf/m2+kns8VV1T/V6r2CZ94Qi7+V14pJ9aLvwqQoJ8FZIwZCYwEaNmyZYCjOf6UdgEBXHaZzLsvLCybIlqZs1qdxWerP6PYWSz75HoqKYm8vXkMHizjDa1bu3n//ffLfPhduyAxkYzIYv63cTwNw5KZfsMcmiS5O0DVHntMatI99ZQsDvYra+GBB2SO6y23yGbM4d7tfayUT1RUI9oXD2A20tVz5OOScu+ZA/Tw9Ji6H4DvXXyxtaeeWv3Ple4PPHfL3Gp9zjn4EntZ/Pc2LExq51fH5qzNtsnLTWzTV5razVk131z300+l7P/IkbVQw9/lkk2Bwdp7763+JgZK1QC1vR+AtVbXrR8DyrcAQNZj7d7tfr1U/9b9iQiLYMamGZzV6iyPz/f4llv4Mvs8XntN+v49tSNnB+d9dB6FzkLm3TiPNkk12992+XK4+Wap5Fm6d6/flN75v/MOPPIIPP+8FnJTQUGngYa40jGAUjfcAGd4sE97QnQCfVr0YcamGR6fy1qoExfJHZHvcc89nseYmZfJ+R+fz568Pcy4ZgadG3q/lSTIXP+hQ6F+fam0cEQVCt97+mnp9rnnHr34q6ASqGmglxpjtgG9gW+NMTMDEYcqmwVUqnt3KZfjSWG4gW0H8seuP9h1cJfb9zqdct17cuBvvFV8G8bp2Tz/g0UHuXDChWzO2szUEVPp2aynR5+rjMMh9X127IDJk6Fx4xodzr1XXpEEcNNNkgT04q+CSKBmAX1lrW1ura1jrW1krb0gEHGoo7uAzj5bnufNc//Zge0GAjBr86wq3/fnn9ClC/z6K5CUhAGPdgErchYx7PNhLNm5hM8u/4xzWp/jPig3nnhCFp69/XbZTFO/ee89KSU6fDi8++4RU52UCjz9PzLElbYASqd+nnyyzASaO9f9Z09pfAqN4hrx3abvKn3PgQOywGrXLpnFWVU5iPJc1sXNU25m1uZZjBk0hktOvMTD36hy8+bJniq33CL9/341axbccYfM8x8/Xmf7qKCkCSDExcbKwG9pDaDwcBkY9aQFEGbCGNhuILM2z8Lpch71eml3y4YN0t1ywgl4lACstTw460EmrJzA8+c+zy2ner5pfGUOHJBNXNq0gddeq/HhqrZ6tdz1d+oki7z8PsiglHc0AYS40j0AyncDPfigVGJwtyAMpBsoKz+LxTsWH/Xagw/CjBnw1ltl67o8SQAvLniRUb+O4p5e9/BI30c8/E2qds89suXl+PGy+Y3f7N4NF18sf9hp03SRlwpqQb8QTPlXadXNvLyyas39+nn++fPbnI/BMGPTDE5rXtap7nDIQPJ998lK20NKT1LJGMCYJWN45IdHGHHSCEYNHCVbONbQpEmyf++TT5bVPPKL/Hy45BLZ4GbePNCFiyrIaQsgxFXUAgCpojzTg7lZybHJ9GrWixmby6aDWis7Z02aVEFNtypaAJ+t+ow7pt3Bhe0u5MOhHxJmav6/5/btcPvtUmftscdqfLjKWSsn+v13qW/RQ2scquCnCSDElSaA8lNBAf75T/j73z3vBvp9++9k5mXy889ysd2xQya9HDX2WUkCmL5xOtd+dS19W/bliyu+ICrcN/3md98tW16OHy+F5/zmnXdk95gnn4RLL/XjiZTyHU0AIa58F1B5l1wCa9aUFWurysB2A3FZF2O+XcSgQZCTU8XFNiZGNjsplwDmbpnLsM+HcXKjk/nmqm+IjfRgc2IPfPutbK/7xBNu936vmUWL4N57ZRuxf/7TjydSyrc0AYS4yrqAbr9dVgTfeSekpVV9jJ5Ne9KosC//vKUH8QkuZs2ClJQqPpCUdCgBTNswjYETBpKamMqMa2aQEJ3g/S9TTn6+3P137CgtGb/JzJRN25s0kRaAzvVXxxD9vzXEVdYCCA+X65m1Mn3SefQsz0PWrQ3HNe57nE7LaY887n7sMzER9u9n/IrxDJ04lM4pnZl34zxS4qrKGtXz/POyovmtt/w4C9PplD0jd+2CL74oWeig1LGjyllAxpj6HhzDZa11v6xTBaXKxgBASjW/9Zb7FkCzZnDSidG0HfE+7+38N5+sPImru1SxtWFSEm/ErODer77gnNRzmDJiCvXqeLZnryfWr4cXX4TrrqvejKZqe+klGSkfPVoHfdUxyd000B0lj6rm4oUDOt/tGFVZF1Cpa68t+zozU6Z3ulyyjePbb8t+LYmJ8OOP4HDdztpxn3LHtDvo1awX7eq3O+p4Ww9s5ZHuaXzSYAdDTxzKp8M+JTrCzeYD1WCtdFvFxMj12W+WLpXBheHDYeRIP55IKf9x1wW01lrbxlrburIHkFkbgSr/qKwL6EizZkFqqhRPa9oU2reX1sGiRWXviQiL4JNhnxAeFs6IL0awMXMj+cWSWfKK83h6ztN0eLMDk+vv4p/LE5k0fJJPL/4AX34p+6w//7xscekX+fmSGVNS5O5fC7ypY5S7FoAny2Z6+yIQFRjuWgCl2raVOjphYfIo3ULyyP3XWya0ZOyQsVz2+WW0f7M9AA1iG2CtJTM/k+GdhvPi/BhSf5oKYb5dh+hwyFz/Tp1kENtvHn5Ydo+fNUtqSit1jHL3L3CyMeZT4GtrbW5Fb7DWFvg+LFVbqhoDKK9tW+la8cSlHS9l6cilrNi9gq3ZW9l6YCvZRdn8X4//k81jVvxTivO4XD6dNTNunPT/f/21H2uvzZwpO8jcey+cf76fTqJU7XCXAN4FRgCjjDE/AZ8C0621RX6PTNWKyEi5BrtLANXVrUk3ujXpVvGLSUnSWX/gQNnCsBrKz5d9fXv3hiFDfHLIo2VlSV3/zp1l8EOpY1yVt1/W2inW2quAVsCXwA1AujFmrDFGb3+OA8YcvSeA35XuwrJzp88O+eabUvbhhRf82CX/wAOwZ48sKy4dPFHqGOZR+9tam2+t/cxaeykwAOgGeL4XoApqR24L6XetWsmzu/mlHtq/X27IL7xQxiX8YvZs6WP6f/8Punb100mUql0eJQBjTCNjzN3GmAXA18AsoLtfI1O15shtIf3OxwngxRdlYfHzz/vkcEfLy5OpniecoKUe1HHF3UKw24CrgA5IF9BD1toFtRGYqj213gXUpImUC/VBAti9WzZ4ufpqP96YP/mkLCueM0e7ftRxxd0gcB/gBWC2tdZVC/GoAKj1FkB4ODRvDunpNT7U669Ltc8nn/RBXBVZuhRefVU2NSjdMFmp40SVCcBaexOAEdcCbay1zxhjWgKNrbW/10aQyr9qvQUA0g1UwxZAdrYsRrvsMlmY5nMOB9x6KzRsKP1MSh1nPJ2E/Ray4Ouqku9zgP/5JSJV62p9EBh8kgDGjJGZpA8/7KOYjvTOO/DHH/DGG2U7mSl1HPE0AZxmrb0TKACw1u4DdKfr40StdwGBJIAdO6C42KuPFxbCqFGy13DPnj6ODWDvXhnwPfdcKfes1HHI0wRQbIwJByyAMSYF0DGB40RAuoBatpSVwNu3e/XxCRMkf/jt7v+xx6SP6Y03tNaPOm55mgDeAL4CGhpjngPmA/6adKdqWcC6gMCrbiCXS7rku3aFAQN8HBfIwO+778Jdd8mqX6WOUx5V47LWTjDGLAH6I6Whh1pr1/o1MlVrAtYFBF4lgKlTpebPJ5/44ebcWtlKrEEDqS2h1HHM43KM1tp1wDo/xqICJCBdQC1ayLMXCeDFF2WzmuHDfRwTSN/SL7/A++/rwK867lXZBWSMWeruAJ68p4LPvGSMWWeMWWGM+coYo//SAigmRubSu2pzVCcmRqZXVnMtwOLFsHAh3H+/rCXzqdxceOghGVW+8UYfH1yp4OPun1BHY8yKKl43gDe7eH8PPGqtdRhj/gM8CvhrOE+5UVoSuqCg7Ota4cVU0NGjJcbrr/dDPK++KgXqJk3Szd1VSHCXAE508zpAFduFV8xaO6vct78COs8ugMrvCVDrCWBFVfcXh9u/X/r9r70WEry57ajK7t3St3TppXDGGT4+uFLByd1KYN9U66razcBnlb1ojBkJjARo2VK3HvYHT7eF9LlWrWDaNBl49WA09+OPZazCL7t9Pf20HPyFF/xwcKWCk9/aucaY2caYVRU8Lin3nscABzChsuNYa8dYa3tYa3ukpKT4K9yQ5um2kD7XqpX0O+3Z4/at1som9D17Qndf16Fdv16WFd9+u59qSigVnHw9jHaItfa8ql43xtwADAL6W2utv+JQ7nm6LaTPlbbo0tJkQLgKP/8s2/COHeuHOB55RJpBfqsop1RwqlYLwBgTV7IiuEaMMQORQd8h1travuyoIwS0Cwg8Ggh++22ZlXnllT6OYf582UT44YfdJiGljjfupoGGGWOuNsZ8a4zJQNYB7DTGrC6ZynmCl+d9E6gHfG+MWWaMGe3lcZQPBLQLCNwmgIwMmDwZbrjBx4PU1sq0zyZNZF6pUiHGXRfQT8BsZJrmqtI9AYwx9YFzgBeMMV9Za8dX56TW2nbeBKv8I2AtgMREqFfPbQL44AOpGefzwd9vv5VFBe+8A3FxPj64UsHPXQI4z1p7VLlGa20WMBmYbIyJ9EtkqtYErAVgjIwDVLEYzFopy3P22dCxow/P7XLB449D27Zw000+PLBSx44qu4CstcUl3UCrqnqP78NStSlgg8DgdjHYwoWwebMfrtGTJsHy5fDMMxCp9zAqNLkdBC7p9llesguYOg4FrAsI3CaAjz+W+C67zIfndDik1n+XLjBihA8PrNSxxdNpoE2A1caY34Hc0h9aa4f4JSpVqwLWBQSSAPbtg5wcGQ8op6gIPvsMhg496qWa+fBD2LgRpkzRkg8qpHmaAJ72axQqoALeAgAZBzii9v706ZIbrrvOh+crLJRVv6edBoMH+/DASh17qkwAxhhjxVx37/F9aKq2hIdDVFSAWgDlF4MdkQDGj5ep+eef78PzvfMObN0K48bpTl8q5Llr//5kjLn7yP5/Y0yUMeZcY8yHwA3+C0/VloDsCgaVrgXYvx+++QauusqHZZ9zc+H552Uj4XPP9dFBlTp2ufunNRAp1vapMaY1sB+IBsKBWcAoa+0y/4aoakNAdgUDWYQVGXlUApg0ScYArr3Wh+f63/+k6ueXX/rwoEodu9xVAy0A3gLeKpnv3wDIt9bur43gVO0JyK5gIIOwLVoctRZg/Hg48UQfFn7Lzob//AcuvBD69PHRQZU6tnk8BcJaW2yt3akX/+NTwLqAQMYByrUAtmyBefPk7t9n3fSvvQZZWfDssz46oFLHPp0Dp4AAdgGBjANs2XLo208+kedrrvHR8bOy4JVXZLMXn9eSVurYpQlAAQHsAgJZkLVjh2zHCEycKJtypab66PgvvyzrDJ7W2cxKleeuGuibxhjtMA0BAe0C6tdPnufOZeNGWLkShg/30bEzMuD116WOdJcuPjqoUscHdy2AjcArxpgtxpj/GGO61kZQqvbFxASwBdC1K8THw5w5TJ4sP/JZ6Yfnn5ddx556ykcHVOr44a4Y3OvW2t7A2UAW8IExZq0x5gljjO6ddxwJaAsgPBzOPBPmzmXyZOjVSyYG1diWLbKTzE03QYcOPjigUscXj8YArLVp1tr/WGu7AVcDlwJr/RqZqlUBHQQGOPts0tblsXgxDBvmo2M+8YRMM9W7f6Uq5FECMMZEGmMGG2MmAN8BGwBf/TNVQSCgg8AA/frxJdLv45MEsGKFLCa45x5o3twHB1Tq+OOuFtD5wFXAxcDvwERgpLU2t6rPqWNPQLuAALp1Y3K4g1MSt9K2rQ/6fx59FBISZMN3pVSF3LUA/gEsBDpaawdbayfoxf/4FBMjZfKLA7S9z849EfziPI1h5quaH2zuXCkl+uijkJRU8+MpdZxyVwrinNoKRAVW+T0BArFB1ldfgSWMYXtHQ8YIKQPqDWvh4YehWTO4+27fBqnUcUYXgikgwNtCApMnw4mt8unEWrmD99akSfDbbzLwW7rRgVKqQpoAFFC241Z2du2fe+9eueYPuyoK6taFOXO8O1B2Ntx3H3TrBjfe6MsQlTouaQJQADRoIM979tT+uadMAacThl0RLjUgvG0BPP447Nolm774bBMBpY5fmgAUACkp8hyoBJCaKguC6dcPVq+ufiCLF0u9/7/9DXr29EOUSh1/NAEoIHAJID8fZs+W7XmNAc4+W16YN8/zgzidcMcdMnD83HN+iVOp45EmAAUELgH89JMkgUGDSn7Qo4eMSFdnHODtt2HJEqn5n5DgjzCVOi5pAlAAREfL+GttJ4Bp0yAuruzGn8hI+ebzzz0LZuNGeOwxGDAArrjCr7EqdbwJSAIwxjxrjFlhjFlmjJlljGkaiDjU4VJSZEZObbFWEsCAAVCnTrkXXnhBdoW/+WZ5U2V27YILLoCoKBg92ofbhykVGgLVAnjJWnuytbYrMA14IkBxqHJSUmq3BbBiBWzdWq77p9TJJ8v+vdOmSfdORbKzZX/fjAxZ9du6td/jVep4E5AEYK0tP9s8DqjiNk/VltpOANOmyfNFF1Xw4j33yN39Aw/IrKDyCgtle8dVq2QFmc76UcorARsDMMY8Z4zZClyDtgCCQoMGtZ8AevaExo0reDEsDMaNkxVqV18NBw7A0qXysyFD4Mcf4YMPJEkopbzit9UyxpjZQEX/tB+z1k6x1j4GPGaMeRS4C3iykuOMBEYCtGzZ0l/hKspaANb6vzs9I6OsYkOlGjeGsWNljmhiYtnPY2PhjTfg2mv9G6RSxzm/JQBr7XkevvUT4FsqSQDW2jHAGIAePXpoV5EfpaTI7om5uTIjyJ+++04SzVH9/0caNEhW9qanwymnyKNtW9lFTClVIwFZL2+MOcFau7Hk2yHAukDEoQ5Xfi2AvxPAtGnQtKmU7XFr5Ej/BqNUiApUwZQXjDEdABeQBtwRoDhUOeUTgD8n1RQVwcyZMGKEztxUKpACkgCstbqdZBCqrdXAP/8MOTnSta+UChxdCawOqa0EMH26rN0691z/nkcpVTVNAOqQ0gTg79XAM2ZItYe4OP+eRylVNU0A6pC6daUkgz9bAFu3wpo1On1fqWCgCUAdYoz/VwPPnCnPAwf67xxKKc9oAlCH8fdq4BkzoHlz6NTJf+dQSnlGE4A6jD9bAMXF8P33cvev0z+VCjxNAOow/kwAv/0mRTy1+0ep4KAJQB3Gnwlgxgyp4NC/v3+Or5SqHk0A6jApKbJIq7DQ98eeMQN69z68rptSKnA0AajD+GsxWEaGbNur0z+VCh6aANRh/LUY7Pvv5Vn7/5UKHpoA1GH81QKYMUOmmJ56qm+Pq5TyniYAdRh/JACXSxaAXXCBbPSllAoO+s9RHaZBA3n2ZQL44w85nnb/KBVcNAGowyQlyVRNXyaAGTPk+fzzfXdMpVTNaQJQhwkLg+Rk3yaAmTOl779RI98dUylVc5oA1FF8uRjswAH45Red/qlUMNIEoI7iywTw44/gdGr/v1LBSBOAOoovE8DMmVCvnqwAVkoFF00A6igpKb5ZCGatDAD37w+RkTU/nlLKtzQBqKOkpEBWlnTd1MSGDZCWpv3/SgUrTQDqKCkpcveemVmz45RO/9QEoFRw0gSgjuKrxWAzZ0L79tC6dc1jUkr5niYAdRRflIMoKIA5c3T2j1LBTBOAOoovEsDPP0N+vnb/KBXMNAGoo/giAcyYAXXqwNln+yYmpZTvaQJQR0lOlueaJICZM+HMMyEuzjcxKaV8L6AJwBjzoDHGGmMaBDIOdbjISCkK520C2LoVVq/W7h+lgl3AEoAxpgVwPpAeqBhU5WqyGnjqVHkePNh38SilfC+QLYBRwEOADWAMqhI1WQ08ZQp06CAPpVTwCkgCMMYMAbZba5cH4vzKPW9bAAcOyPTPSy7xeUhKKR+L8NeBjTGzgcYVvPQY8A9ggIfHGQmMBGjZsqXP4lNVS0mBX3+t/ue++w6Ki2HIEN/HpJTyLb8lAGvteRX93BjTBWgNLDfGADQHlhpjellrd1VwnDHAGIAePXpod1EtadoUMjIgN7d6M3mmTJHkcfrp/otNKeUbtd4FZK1daa1taK1NtdamAtuAUyu6+KvAOe002cz99989/0xREUyfLoO/4eH+i00p5Ru6DkBVqHdvMAbmz/f8M3PnQna29v8rdazwWxeQp0paASrIJCbCSSdVLwFMmQIxMXBehZ1/Sqlgoy0AVam+fWHhQs/2BbBW5v8PGACxsf6PTSlVc5oAVKX69oWcHFi50v17//hDVgBr949Sxw5NAKpSZ5whz550A02ZAmFhMGiQf2NSSvmOJgBVqZYtoXlzzxNAnz5llUSVUsFPE4CqlDHSDTR/vvTxV2bRIli+HIYPr73YlFI1pwlAValvX9i+HdKrKNn36qsQHw833lhrYSmlfEATgKqSu3GA9HSYNAlGjpQkoJQ6dmgCUFXq0gXq1YMFCyp+/Y035Pnuu2svJqWUb2gCUFUKD5fB3YpaANnZ8O67cMUVMmCslDq2aAJQbvXtC6tWwb59h//8/fclCfz974GJSylVM5oAlFtnnCGzgBYuLPuZwwGvvw5nnQU9egQuNqWU9zQBKLd69YKIiMPHAb78EtLS9O5fqWOZJgDlVlwcnHoqvPIKNGgAdevCiBHQrp2u/FXqWBbwaqDq2PDsszLds04diI6W50su0br/Sh3LNAEojwwYIA+l1PFDu4CUUipEaQJQSqkQpQlAKaVClCYApZQKUZoAlFIqRGkCUEqpEKUJQCmlQpQmAKWUClGaAJRSKkRpAlBKqRClCUAppUKUJgCllApRmgCUUipEBSQBGGOeMsZsN8YsK3lcFIg4lFIqlAWyHPQoa+3LATy/UkqFNO0CUkqpEBXIFsBdxpjrgcXAA9bafRW9yRgzEhhZ8u1BY8x6L8/XANjr5Wf9LVhjC9a4IHhjC9a4IHhjC9a4IHhjq25crSr6obHW+iacIw9szGygcQUvPQb8igRvgWeBJtbam/0SSFk8i621Pfx5Dm8Fa2zBGhcEb2zBGhcEb2zBGhcEb2y+istvLQBr7XmevM8Y8y4wzV9xKKWUqligZgE1KfftpcCqQMShlFKhj/ZdSQAABahJREFULFBjAC8aY7oiXUBbgNtr4ZxjauEc3grW2II1Lgje2II1Lgje2II1Lgje2HwSl9/GAJRSSgU3nQaqlFIhShOAUkqFqOM+ARhjxhpjMowxQTXQbIxpYYz5yRiz1hiz2hhzb6BjKmWMiTbG/G6MWV4S29OBjqk8Y0y4MeYPY0xQzR4zxmwxxqwsKW+yONDxlDLGJBpjvjDGrCv5/613oGMCMMZ0KFcOZpkxJtsYc1+g4wIwxtxf8v/+KmPMp8aY6EDHVMoYc29JXKtr+vc67scAjDFnAQeBj6y1JwU6nlIlM6GaWGuXGmPqAUuAodbaNQEODWOMAeKstQeNMZHAfOBea+2vAQ4NAGPM34EeQLy1dlCg4ylljNkC9LDWBtXCIWPMh8DP1tr3jDFRQKy1dn+g4yrPGBMObAdOs9amBTiWZsj/852stfnGmM+B6dbacYGMC8AYcxIwEegFFAEzgP+z1m705njHfQvAWjsPyAp0HEey1u601i4t+ToHWAs0C2xUwoqDJd9GljyC4k7BGNMcuBh4L9CxHAuMMfHAWcD7ANbaomC7+JfoD2wO9MW/nAggxhgTAcQCOwIcT6mOwK/W2jxrrQOYi0yl98pxnwCOBcaYVKAb8FtgIylT0s2yDMgAvrfWBktsrwEPAa5AB1IBC8wyxiwpKWESDNoAe4APSrrN3jPGxAU6qAqMAD4NdBAA1trtwMtAOrATOGCtnRXYqA5ZBZxljEk2xsQCFwEtvD2YJoAAM8bUBSYD91lrswMdTylrrdNa2xVoDvQqaXoGlDFmEJBhrV0S6FgqcYa19lTgQuDOku7HQIsATgXettZ2A3KBRwIb0uFKuqWGAJMCHQuAMSYJuARoDTQF4owx1wY2KmGtXQv8B/ge6f5ZDji8PZ4mgAAq6V+fDEyw1n4Z6HgqUtJdMAcYGOBQAM4AhpT0tU8EzjXGjA9sSGWstTtKnjOAr5B+2kDbBmwr14L7AkkIweRCYKm1dnegAylxHvCXtXaPtbYY+BLoE+CYDrHWvm+tPdVaexbSve1V/z9oAgiYkoHW94G11tpXAx1PecaYFGNMYsnXMcg/iHWBjQqstY9aa5tba1ORLoMfrbVBcWdmjIkrGcynpItlAEFQ4sRauwvYaozpUPKj/kDAJxoc4SqCpPunRDpwujEmtuTfaX9kjC4oGGMaljy3BC6jBn+7QJaDrhXGmE+BfkADY8w24Elr7fuBjQqQu9nrgJUlfe0A/7DWTg9gTKWaAB+WzMwIAz631gbVlMsg1Aj4Sq4XRACfWGtnBDakQ+4GJpR0tfwJ3BTgeA4p6cc+n9opB+MRa+1vxpgvgKVI98ofBFdJiMnGmGSgGLizslL6njjup4EqpZSqmHYBKaVUiNIEoJRSIUoTgFJKhShNAEopFaI0ASilVIjSBKCUUiFKE4BSfmCM+a8xZqkxpmegY1GqMpoAlPKxkpXADZHFTUFTrlqpI2kCUMpLxphUY0x+uZXcAFhrc5HV1HOAN0reG1Oy6UmRMaZB7Uer1NE0ASjlISOO/DezuaRqavn3JSM15HMAJ4C1Nr/kfcFSV14pTQBKVaXkLn+tMeYtpDaMJ7XXH0fqya8GOvkzPqVqQhOAUu51QLYU7eZux6qSzX36AJ8hFSQ7+z06pbykCUAp99KqsR/yv4BnrFRZ1ASggtpxXw5aKR/I9eRNxpiuSH32vsaY/wHRwEp/BqZUTWgCUMp3/gMMttb+AGCMaYTUklcqKGkXkFI+YIw5F4grvfgDlGxxGGeMqR+4yJSqnLYAlKqCtXYLcJIH7/sR+LGCnyf4ISylfEJbAEp5zwkkHLkQrCKlC8GASMDl98iU8oBuCamUUiFKWwBKKRWiNAEopVSI0gSglFIhShOAUkqFKE0ASikVojQBKKVUiNIEoJRSIUoTgFJKhaj/D0mWqHUDJ1KdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "r=np.arange(1,9,0.1)\n",
    "V_AA=r*0\n",
    "V_BB=r*0\n",
    "V_AB=r*0\n",
    "\n",
    "\n",
    "N_coeff=np.shape(list_LJ)[0]\n",
    "for ratio in (0.00,0.50,1.00):\n",
    "    list_coeff=np.genfromtxt(\"Coeff\"+str(ratio)+\"_\"+str_type+\".txt\",usecols=[1])/kbT\n",
    "    \n",
    "    \n",
    "    N_coeff=np.shape(list_coeff)[0]\n",
    "    N_coeffLJ=np.shape(list_LJ)[0]\n",
    "    N_coeffGauss=np.shape(list_Gauss)[0]\n",
    "    \n",
    "    print (N_coeffGauss*3)\n",
    "\n",
    "    for i_coeff in np.arange(N_coeffLJ):\n",
    "        a=list_LJ[i_coeff,0]*sigma\n",
    "        b=list_LJ[i_coeff,1]\n",
    "        V_AA=V_AA+LJ(r,a)*list_coeff[i_coeff+(N_coeff-3*N_coeffLJ)]\n",
    "\n",
    "    for i_coeff in np.arange(N_coeffGauss):\n",
    "        a=list_Gauss[i_coeff,0]\n",
    "        b=list_Gauss[i_coeff,1]\n",
    "        V_AA=V_AA+gaussian(r,a,b)*list_coeff[i_coeff]\n",
    "\n",
    "    plt.plot(r,V_AA,\"red\")\n",
    "\n",
    "\n",
    "\n",
    "    for i_coeff in np.arange(N_coeffLJ):\n",
    "        a=list_LJ[i_coeff,0]*sigma\n",
    "        b=list_LJ[i_coeff,1]\n",
    "        V_BB=V_BB+LJ(r,a)*list_coeff[i_coeff+N_coeffLJ+(N_coeff-3*N_coeffLJ)]\n",
    "\n",
    "    for i_coeff in np.arange(N_coeffGauss):\n",
    "        a=list_Gauss[i_coeff,0]\n",
    "        b=list_Gauss[i_coeff,1]\n",
    "        V_BB=V_BB+gaussian(r,a,b)*list_coeff[i_coeff+N_coeffGauss]\n",
    "    \n",
    "    plt.plot(r,V_BB,\"green\")\n",
    "\n",
    "\n",
    "    for i_coeff in np.arange(N_coeffLJ):\n",
    "        a=list_LJ[i_coeff,0]*sigma\n",
    "        b=list_LJ[i_coeff,1]\n",
    "        V_AB=V_AB+LJ(r,a)*list_coeff[i_coeff+N_coeffLJ*2+(N_coeff-3*N_coeffLJ)]\n",
    "\n",
    "    for i_coeff in np.arange(N_coeffGauss):\n",
    "        a=list_Gauss[i_coeff,0]\n",
    "        b=list_Gauss[i_coeff,1]\n",
    "        V_AB=V_AB+gaussian(r,a,b)*list_coeff[i_coeff+N_coeffGauss*2]\n",
    "    \n",
    "    plt.plot(r,V_AB,\"blue\")\n",
    "\n",
    "\n",
    "\n",
    "    plt.ylim([-5,1 ])\n",
    "    plt.xlabel('r [$\\AA$]')\n",
    "    plt.ylabel('V (r) [eV]')\n",
    "    \n",
    "    plt.plot(r,LJ(r,sigma),'r--')\n",
    "    plt.plot(r,LJ(r,sigma*0.88)*0.5,'g--')\n",
    "    plt.plot(r,LJ(r,sigma*0.80)*1.5,'b--')\n",
    "    plt.show()\n",
    "   \n",
    "    gp.s([r,V_AA,V_BB,V_AB])\n",
    "    os.system ('mv tmp.dat OutLasso.txt')\n",
    "\n",
    "\n",
    "    gp.s([r,LJ(r,sigma)*kbT,LJ(r,sigma*0.88)*kbT*0.5,LJ(r,sigma*0.80)*kbT*1.5])\n",
    "    os.system(\"mv tmp.dat Out\"+str(ratio)+\"_\"+str_type+\".txt\")\n",
    "    \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AA\n",
      "LJ 1.0000000000000002 1.0 0.992174\n",
      "Gauss 2.25 2.0 1.20451\n",
      "Gauss 2.5 5.5 -0.00189808\n",
      "Gauss 2.5 10.0 0.998957\n",
      "BB\n",
      "LJ 0.8600000000000002 1.0 0.00862156\n",
      "LJ 0.8800000000000001 1.0 0.493029\n",
      "AB\n",
      "LJ 0.8000000000000002 1.0 1.49876\n",
      "Gauss 1.75 9.0 -0.0268035\n",
      "Gauss 2.5 8.5 -0.00263855\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEQCAYAAAC5oaP8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxV1bn/8c+TERIChBmRMCgCCioQZ+usVevc2mona22xt62/69XWau11qFO9bfW2td7r2NpbqrYoaqtVwTorIqAyKzLJPE8hCcnJeX5/rBNByJyc7AP7+3698trk7L3XeghhP3uvvQZzd0REJH6yog5ARESioQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU5EmADN72MzWmNmsKOMQEYmjqJ8A/gicHnEMIiKxFGkCcPfXgA1RxiAiEldRPwGIiEhEcqIOoDFmNhYYC1BYWDhm2LBhEUe0Z9u0CRYsgN69Yd99P7sv6UneW/kenfI6MbTH0JZXsnEjLFxYdyUi0u6mTZu2zt177vp5xicAd78fuB+gtLTUp06dGnFEe7YZM+CQQ+D734cbbth9f5dfdOG4kuP4+1f/3vJK5syBgw6Cyy+Hm29ueTki0ibMbEldn6sJKGYGDAjbTp3q3n9gzwMpT5S3rpJhwyA3F6qqWleOiKRV1N1AHwXeBoaa2TIzuyzKeOKgqChsV66se3/vwt7MWzevdZVkZUGfPjB7duvKEZG0irQJyN0vjrL+OMpKpfxJk+rev7Z8LSu2rqAyUUmHnA4tr2jzZnjttZafLyJppyagGMrOhm3b6t63X/F+AExfMb11lXTrVn8lIpIRlABiKCen/mvz8B7DAZi6opUv2/v2hUQCqqtbV46IpI0SQAzl5UFlZd37xvQdA8Dsta1svy8pCdtZmuVDJFMpAcRQfj5s3173vsP3PRyA+Rvmt66SAw4I2+mtbEoSkbRRAoih4cOhS5e693Xt0JXOeZ3pWbDbmJHmOfHEsE0mW1eOiKSNEkAMDRvW8HV5aI+hbNq+qXWVjB4dtps3t64cEUkbJYAYyskJszXUp7hjMTNXz2xdJV26QGGhmoBEMpgSQAx99FF4B1Bez4DfzZWbWVm2kuqaVvbgqa6GV19tXRkikjZKADFU2/6/bFnd+wcVDwLgg9UftK6ioqIw+5yIZCQlgBjq1i1sV6yoe/+BPQ4EYOryVo4F6NGj/v6mIhI5JYAY6t49bOubD2hU31EAzFzTyvcA/fqFt816ChDJSEoAMdQz1cNz9eq69x/R7wigDcYCDApNSUyb1rpyRCQtlABi6OCDw7a4uO79PQt70jm/M/t2buViLqenlntWM5BIRlICiKHaG/OsBv719yvej1Vlq1pXUe1YgDVrWleOiKSFEkAMFRSE7cKF9R/TKa8T7616r3UV7bNP2L79duvKEZG0UAKIodrVwN58s/5jKhOVrCpbRXlVK1YH69ABzOD111tehoikjRJADHXuHLYNzdJwZL8jAXj6w6dbV1nHjrBuXevKEJG0UAKIKTPYurX+/RcedCEAT89rZQLo0qXhikQkMkoAMZWTU/9UEADHlhxLlmUxZcWU1lXUu3eYd8K9deWISJtTAoip3NyGE4CZ0buwNyu21jNcuKn69w/bJUtaV46ItDklgJgaOnRHJ536XDbqMrbXbGdTZStG8p57btiuX9/yMkQkLZQAYmrAgMbXajlh4AkAvLPsnZZXVDvqrL6Jh0QkMkoAMZVM1j8XUK3D+4XlIX/zzm9aXlHtvBOvvNLyMkQkLZQAYmrx4sZ7ZxblF9ExpyOTl01ueUUaDCaSsZQAYqp2LEBjzUD7Fe/HxsqNVNVUtayivLzQ5aixxw0RaXdKADFVuyjMqkam+zmm5BgA/v7h31teWUEBbNjQ8vNFJC2UAGKqa9ewXb684eMuPLANBoR169Zwn1MRiUSkCcDMTjezD83sYzO7NspY4qZ2UZjGOuecMPAEsi2bRZsWtbyyvn0hkQhrBItIxogsAZhZNvB74AzgQOBiMzswqnji5qijwrb2XUB9srOyOXW/U9m8vYGJgxrzla+EbX0r0IhIJKJ8Ajgc+NjdF7p7FfAYcG6E8cRK7ZoAZWWNH3tkvyOZuWYmmytbmASGDg3bRa14ihCRNpcTYd39gKU7fb8MOCKiWGKnd++wfestOPvsho/t06kPAPdMuYfrj7u++ZUNGRK248bB5z7X/PMlEolkgm1V2yirKiM/O5/qZDVry9eyrnwdlYlKKqorqKqpoqqmiiHdh1BdU83HGz5mVdkqqmuqqUpWUV0Tmv1G9R1FIplg1upZrCpbRVWihupEDVXbjUR5F4YXHcH2imzmrJ/NxvLN5HRZhWdXUbmlM75+EPva4dQkslhSOYNtia3k958FOeVsX7cvrDqEbpVjSCazWGPvU0UF+cMnQk4l1SuGw9Jj6FQ+Ek9msanjNGqoIm/MOMiuJrHoSGzJiXQoG4a7sa1oOkmrIfdzdwNOYt7p2OJTyNs6GEtCRbfpkL2d7BPvwHH8gy+Tu/AkOm7pCzhbes3D8jbjJ92KG9g736Hjos9RsLU7bjVs6LWErMJVJE6+BYDsN66iYNERFJZ1piarhvU9V5PVdSHVJ/887J/4c4qWHkJP38RVv6ni8ktPadN/4ygTgNXx2W4zhpnZWGAsQElJSbpjio3a7vmvvtr4sWcdcBY8C5MWTmpZAhg0KEw/+vLLzT9XPlWTrGF9xXqWbFrCmm1rWFexjo0VG9lUsYmBxQNxdz5a/xEzVs9gW/U2KqorqExUUllTycheI3GchRsXsmTTEhLJBIlkgpoE1JR3pbi8lOqtxWxLbCGRtRW6fwRdl8K2HjD7Qvj481BdCEXLoWA9HDwO9pkOqw+Cf90KS/eHZA506gm55fD5q2Dg67DgFHjy/2BbL8Agqzocd8lJMOhVmPFVePJ/eX7Xv+zYMaH8qd+DZ+9g7qc7Uo0E3y2F4sXYnM/hr/yYHZ3ZTgCgfJ8nSBZuJPujL1Ez5dts+XT/GACSQ/9AIjdB3oJDqJpxPlUkySJJImswllVDzaF3UW35FC4dQvLj0XRMJDGcqk3Hkp27lZothSTIpfOqErJX96fr9rC/smYYeR1WU7a5B5BN0YaeFGwupFtlBW5O5cZedKwuZ/2GsNxql82FdC1P0L1qPYksp7o8l8JsWL1uAGB0rYBe1Rs5MDGfglVd2/YXCjCPaJZGMzsKuMndP5/6/joAd7+jvnNKS0t96tSp7RTh3i8nB0pKGl4ZrFaHWzvQKa8T665p4dz+PXqEnkAx7w2USCZYuGEhCzYuYOHGhSzbsoxVZavoW9QXgI/Wf8S7K96lvLqciuoKqpPVJJIJ8rPzqUhU7F6gA1v6wYrRsGEIuEHHzdBzNvR7N1x4//kbspYfDds7k8ypAKuGEX+FwRNh4yD45z27l3vAM9D3fdjcD96/bKcdSchKkDv4bTr0WgoVPSj/6Ehykp3IzgbL20Z2djklQz6kqPs22JDD1o+70a3K6OjlZNtm8mwbg/t8QrKokJz1FWQtL6dbdSX5yQpyqaCDV9KH1Wz3jhhJ8qkil2pyqSabGnJIUEAZNeSSBDqynWxqPv3KIUFOrmEdO4Inydm2mSyD7KwkWVlGTraTPWx/sou7YJs2kL1iWfjPkJ0dtjk5cPTRoa/0ypXwySdh9sSdv447LnRvXr48HJObG87Lywt/Puqo8OcVK2Djxh37a88/5JDw/Zo1UFGx4/PaMkpKQjxlZaEDRadO0KdPOKYFzGyau5fu+nmUTwDvAkPMbBCwHLgI+GqE8cROQUH43WyKEb1GMG3lNGaunsnI3iObX9mQITB5cliFpnYQwl6kJlnDqrJVzFs3j39+/E8WbVzE8q3LWVu+ls2Vm+mc35ny6nLWbltLkrpH3+Vk5dAptxNlldvJrepBbl5H8hPFVL1zCR1Xn0jHjf2prHSqEgmyOq8i2XkJia1dYckJdZTm5OYaXbqEH3luLuQXQIcOnSgogOE9rmTkgCvJGwKzCkJngC5dwlfXrjB06DkMGHAOuYkKqmdOp3jNRxSvm0+HVYvDBa+gIHytXAk2GRLboTIBZakbyrWN/MCqi6BXr9RAwWVQkBdWkOvQISwiNGpU+J1JJODDD6GoKFwEO3UKQY4YEXqXJZOwbVsIujb4goJw8ZRGRZYA3D1hZj8EXgCygYfdfXZU8cRRcXHT52i75cRbOPMvZ3L1i1fz4jdebH5lRx8dEsDTT8M3v9n88zPAhooNzFw9k8dmPcaH6z/kk82fsLZ8LduqtuE4Sd/9wp5t2XTI6cDg4sEc2udQCnMLWbR+GXmbh9LJ+mKbBzH9H6NYt7QHmzbksKXcSCYhAex8v18GZGVBfj50KYDBvfty7LGj6NkT5s0L73T22SfMvl1SAoMGGd27h5a3OrmHuUCWLgWeCYV8uDT01Nq4MVRWVQVbttR9fk4O9OsXxnjk5oZfpqKicAEuLg4X8DFjoLAwjDbs2zcE2K9fyDb1BibtKconANz9OeC5KGOIs969w9NtMhn+vzfkjCFnMLDrQOasnUPSk2RZMzuQnXsu3HVXSAIZnADcnfdWvscLC15g8rLJzF03l1Vlq0gkE7s1wWRZFp3yOtG/S39K+5Zy8uCT2afTPmyv2c7IXiPZp8MQPpyXzbhxMPkZeGZxGBC963AIs3DDWlQEgweHJ/2RI+Gss8I1Mz8f9tsv3Cw3S1VV6Hn11FMwdSosWBDu2DdtCkHU1Ox+TnZ2uAsfNAhOOikEM39+CGTwYNh//9Crq1cvXcT3ApEmAInWuefCu+/C2rU7egU15PaTbuerT36VSQsncdp+pzWvsmOPDe8BKitbFmwaJJNJXl7yMk/Pe5pFGxexvWY701dOZ33FjrULsiyLorwiRvUZxQXDL+CA7geQm53LEf2OoLhj8afHLVgAjz0GT7wSbqbXrw9/1Z1fseXnh8lRS0pCC8bXvx4u7H36hBvqFtu6Ff75zzDj6owZ4aK/ZUtoW971Ip+fH+7ShwyBCy8MwWzfHi7qBx3U4jZm2TMpAcTYgalhd8uWNS0BXDD8Arp17MZVL1zFrO/Pal5lWVlQWgrTpjU/0DZSmajknWXv8N+T/5u3lr3F2m1r8Z06no3uO5rzhp2HYRzU6yBO3+90hvUctls5iQT84x/h5/bWW/D88599l5KTEy703/lOaAnp1y/c0Xfs2AZ/iWXLYPx4eP31UNGMGSHj7Cw3N2SVK68MF/qcnHBxHzFCbePyGUoAMVY7Cvi110JzbWPyc/IZ2Wskry55lZcXvcyJg05sXoX5+eGCtXVraO9Is63bt/LQew8xfs545q6bS9n2MqqSYVbTvOw89u+2P4f3O5wvDPkC5ww9h8K8wjrLcYeXXoL77w/X3dWrd9zZ9+sXrquFhXDqqeGmunYVzFarroYPPoA//xn+/vfQXr9z+1FJCYweHf7xunULTTYnn9wuP1vZOygBxNjOg8H+4z+ads6vTv0Vhz14GNdOupZ3vtvMlcIGDgzbZ56Br32teec20eJNi7n77bt5bPZjrNm25tPPC3MLuWz0ZZyx/xkc1f8oehT0aLCcsjJ48UWYOBEmTNgxi0VOTmgGP+44+P73wx1+mzWFV1aGdqTx42HKlJAoa5vMzEK7+4gRofIvfjHc1Yu0ghJAjA1LtW588knTzyntV8o+Rfvw7op32VS5ia4dmjE45eyz4Te/gRdeaLMEUJOs4Y/v/5GH33uYteVrmb9hPgC5WbmM6DmCc4aew7+V/hv7dtm30bK2bIHbboPHH9+xhn1hIRx/fOh9OHZsuMluswu+O3z0EfzlL/DAA59dMyEvL7wF/vKXQ9PZoEGNv6kXaSYlgBirHfPS3DnarjziSq6ZdA0/mfgT7jv7vqafePzxYdvK9wA1yRoenP4g9027jxmrZ1Dj4UXn0O5D+fVpv+asIWdxQI8DmlSWe2jW+cEPYPbs8L1ZaDr/9rdDM3qHDq0K97Oqq0OTzgMPhDfHa1JPKdnZ4UXsaaeFig89tA0rFalbZCOBW0Ijgdte167hore5GfO81SRrKLy9kKK8Iub8YA49C3s2/eTi4nARbMosdDtxd9745A2enf8s42aOY9mWZQD0LuzNuUPP5SfH/ITB3QY3ubzFi+Gmm0Lz1/z54frbv394cXv11W180U8kwgX/3nthzpwdy7ANHAg//jGcfvqO6TJE0iATRwJLBiguDh1LmiM7K5uJ35jIaX8+jfMfP5+XvvkS+Tn5TTt51KjQXbGysklX2Tc/eZObXr2J15a8RlVNFdmWzWn7ncYlh1zC5WMup3+X5r1xnTABrr02tLxAeH/6yCNw/vlt/O7UPXQP+sc/4G9/C31tIQyIOvvscOHff/82rFCkBdx9j/kaM2aMS9u6+mr37Gz3RKL55z4+63HnJnzA3QO8OlHdtJPGj3cH9ylT6j1k6/atftULV3nXO7o6N+HchHe+vbNfPP5iX7V1VfMDdfcnn3Tv3TtUDe7Dhrk/+miLimrYokXuF13kXlAQKurQwf3LX3a/+273ZcvSUKFI44CpXsc1VU8AMTd4cBgrtHr1jhlCm+rLB32ZcTPG8cxHz/C5P3yOt7/zduMn1fY3ffttOOywTz/eun0rt79+O/PWz+OFj1+gIlFBblYuJw86mRuPv5HPDWjZNNKTJsHtt4eJSLOy4IQT4MEHwwCsNuMOTz4ZHi0+/jh8lp8P550Hv/1tG/YLFWlbSgAxVzsu6I03QoeT5nrqoqc46N6DmLx8MqP+dxT3n30/h/U7rP4TSkpCW/cDD7Dy0gu5a/JdTJg7gYUbF+I4hbmFfOvQb3HRQRdxbMmxZLWw58vzz4d3qStXhjFRd90Fl166Yy3kNlFREbLJQw+F/voQXuRedx184xvqtSMZTwkg5moviO++27IEYGZM/e5UDrnvEN5f/T6HP3g4XxvxNf7z+P9kQNcB5GXnkWVZYY6dVe8xbcU0zuyRz7LNczjyrvDIYRgDuw7k26O+zdVHXU3H3JYPmZ0/Hy64AGalBiqffHLoVt+mF/4NG0K3ofHjwwvegw8OL3nPOCOMDBPZQygBxNzo0WFb23LREgV5Bcy/Yj6vLXmNG16+gfFzxzNu1rh6j//tILj0PTiu3zGcf9CXuLz08lZd9CG0wtx2G9xwQ/jziBHw6KNh22aWLw/dhF58MfTkKSyEH/0otDGpB4/sgZQAYq62LXzp0oaPa4rjBhzHK996hTXb1vCXmX9h4oKJrK9YT0V1BTVew7Aewzhl8Cl8bescOk35Ha+u/QIcdWWr6507F374Q/jXv6B7d/if/wlTMrSZTZvg17+G//qvMMNmcTH89Kdw1VVq5pE9msYBCHl5YVqItkgCTbJxY5i7ZtCgpi1HVo9kEi6+OPSyLCqCO+8Mo3Xb7JpcVhYKfOaZsOjIueeGPvuXX647ftmjaByA1KuoKFzf2k1xMZxyShgPsHFj+L6Z5s6FE08MvZd69Ai9fNqsuSeZhBtvDBmlujq067/5ZljGT2QvoudX4eyzwyp67eoXvwgvUP/2txadOmJEuPhfcEFYcKrNLv6PPRYyyq23hu+vvTZMDKSLv+yFlACE/v1Dd8lEoh0rHT06NAH98pdNPsU9XJevuy7MYfT44/DEE200xf2yZaEb1MUXhzb/884Lo3fvuENz6MteSwlA2L49tHq89147VmoWOuh//PGOPvQNWLcuTNfwn/8Z7vprr9etVlUVBgzst1+Yc/+GG8Ii5BMm7JWL14vsTAlA6NQpbKdPb+eKr7subG+8scHDJk8O48eefjo0y48fH1bcarUnnwzNPX/4Q/ghzJ4NN98cpgIViQElAGHkyLCdO7edKz7rrLBO4osvfnbx3J2MHw/HHBMG3d54I1xzTRt0wFm3LhT6xS+GRVe+8Q1YsSLMiyESI0oA8ulgsAUL2rliszD/fUVF6Gq5i3vv3dGf/+mnw/TNreIO48aFWTjfeis8VnzwAfzpT2HuHpGYUQIQSkrCtrnTQreJm28O2wcf/MzHEyaEwV25uWHBlnPOaWU9H34IRx4JX/96WArt/vvDogAHH9zKgkX2XEoAglmY1aC8PILKDzkkNMW8+WZohgH+7//Cnf+IEfD++3D00a0oP5GA//f/YPjwsM7uLbeEur77XQ3mkthTAhAAjj0WOneOqPJbbgkLxJx5Jlf8MMk3vxnWPX/zTTjwwFaUO2lSWIDld78LTTx//CP87Gfq1imSogQgAOy7bztOBbGr4cPhS1/iqQ8G8vvfh96XtdM7tMjGjXDZZXDqqeGF7xlnhFFjl1zSpmGL7OmUAAQI18zVq5u9VG+buWGfB7mQv1HKVD68cwLdu7egkOrqMCXo0KHhbv+EE+DVV+G55yJ8vBHJXEoAAuyYjufVV9u/7ttvh1vuzKOoY4IXOJ3eV3zl0/cBTZJMwu9/H6YC/dnPQrPPtGlhgqDjjktf4CJ7uEgSgJldaGazzSxpZrvNUCft74tfDNvHH2/feqdMCR2BCgpg2uyOFN9xTbiTHz0annoqXNzrs2FDeH/Qp0/oMlRWFoYLv/QSHHpo+/0lRPZQUc0GOgu4ALgvovplF5//fJhG+fXX26/OOXPgC18I7x/eeCPcuHPttWHo75Qp4WK+//7wve+FQWObN4e2quXLQ7POM8+EZAFhYNef/qTBXCLNEEkCcPe5EJYTlMyQlRVmPV66NNx0p3udk9mzYdSo0DnnrbdSF/9aTz0Vum+OHw9XXBFW3frRjz5bQI8e8P3vwwEHhOlMtfC6SLNl/HoAZjYWGAtQUjtiSdLilFPCtDjz5rWy+2UjVq+Gww8PN+/33lvP1Ds5OXDRRWHw1o9/HN4JFBeHhWR69gxt/S1YR0BEdkhbAjCzSUCfOnZd7+5PN7Ucd78fuB/CimBtFJ7U4eqrQwKYPDl9CaC6OjTPl5eHa/h3vtPICQMHtmjNABFpXNoSgLufkq6yJT0OPDB0pHn++TBDcjoce2xYwOVLXwrvb0UkOuoGKp8yC71xJkxIT/n33Rfe7Y4cCX/9a3rqEJGmi6ob6Plmtgw4CnjWzF6IIg7ZXWlpeP86eXLbljtxYuipefrpoYu+3v+LRC+SBODuE9x9X3fPd/fe7v75KOKQ3dWusvXII21X5ksvhQt///5hyd3c3LYrW0RaTk1A8hkXXBC2L7/cNuUtXgxnnhmm4n/wQa2yKJJJMr4bqLSvvDzo1QsWLmx9WRUVMGZMWHb3t7+Fk05qfZki0nb0BCC7ueCC0F1z+fKWl+Ee+vpv2ACXXhrGc4lIZlECkN1cdlnYtmZaiLvvhlmzwkvlhx9um7hEpG0pAchuDj00dAf9859bdv4TT4TBu+eeG6Z5EJHMpAQgu8nJgY4dQ9fNhibjrMvdd4flHEeNCuuvq8ePSOZSApA6nXZaeHk7dmzTz/nTn+Cqq8KL5EcfDesMi0jmUgKQOj3ySOiy+dBDTVsk5rnn4FvfCsvtvv56PRO8iUhGUQKQOuXmwrPPhj+ffTZs317/sQ8/HI6BMI/QYYelPz4Rab0GE4CZdWvCV9f2Clba1zHHwA9+AFu3wh137L6/pgZuvDH0GurQIUzff4qmABTZYzQ2EGxF6quhmVuyAU3Uv5e6556wENett4bZQgcMCBf+tWtDcli+HC65JBzXqVPU0YpIczSWAOa6+6iGDjCz99owHslA99wDL7wAX/nK7vt++lO47bb2j0lEWq+xBHBkE8o4qi0CkczVpQs8+WSYzjkrK7zozcoKC3apyUdkz9VYAnjCzB4FnnL3bXUd4O6VbR+WZJpjjw1fIrL3aKwX0APAWcAiM3vczM4zs7x2iEtERNKswQTg7k+7+8XAAOBJ4BLgEzN72MxObY8ARUQkPZo0DsDdK9z9cXc/HzgNGAU8n9bIREQkrZqUAMyst5ldYWZvAk8BLwJj0hqZiIikVYMvgc3su8DFwFBCE9A17v5mewQmIiLp1VgvoKOBXwCT3L2Z80KKiEgmazABuPulABZ8HRjs7j83sxKgj7tPaY8gRUSk7TV1Mrh7CQO+Lk59vxX4fVoiEhGRdtHUReGPcPfRtdM+uPtGjQcQEdmzNfUJoNrMsgEHMLOegN4JiIjswZqaAH4LTAB6mdltwBvA7WmLSkRE0q5JTUDuPs7MpgEnE6aGPs/d56Y1MhERSaumvgPA3ecB89IYi4iItKPGVgSb3lgBTTmmjnN+aWbzzGyGmU3QqmIiIu2vsSeA4WY2o4H9BnRpQb0TgevcPWFmdwLXAT9pQTkiItJCjSWAYU0oo6a5lbr7izt9Oxn4UnPLEBGR1mlsJPCSdojh28Dj9e00s7HAWICSEi09LCLSVpr8Eri5zGwS0KeOXde7+9OpY64HEsC4+spx9/uB+wFKS0s9DaGKiMRS2hKAuze4WqyZXUJYbexkd9eFXUSknTUrAZhZIVDp7s1u99+lnNMJL32Pd/fy1pQlIiIt01g30Cwz+6qZPWtmawjjAFaa2exUV84hLaz3HqAImGhm75vZ/7awHBERaaHGngBeBiYRumnOql0TwMy6AScCvzCzCe7+5+ZU6u77tyRYERFpO40lgFPcvXrXD919A/AE8ISZ5aYlMhERSasGm4DcvTrVDDSroWPaPiwREUm3RmcDTTX7fJBaBUxERPYSTe0F1BeYbWZTgG21H7r7OWmJSkRE0q6pCeDmtEYhIiLtrsEEYGbmwauNHdP2oYmISDo19g7gZTO7Ytf2fzPLM7OTzOwR4JL0hSciIunSWBPQ6YTJ2h41s0HAJqADkA28CNzt7u+nN0QREUmHxmYDrQTuBe5N9ffvAVS4+6b2CE5ERNKnOUtCVgMr0xiLiIi0o0bHAYiIyN5JCUBEJKYamw30HjM7ur2CERGR9tPYE8B84NdmttjM7jSzQ9sjKBERSb/GJoP7jbsfBRwPbAD+YGZzzewGMzugXSIUEZG0aNI7AHdf4u53uvso4KvA+cDctEYmIiJp1aQEYGa5Zna2mY0D/gl8BHwxrZGJiEhaNTYX0KnAxcAXgCnAY8BYd9/W0HkiIpL5GhsI9lPgL8CPUquAiYjIXqKxqSBObJTd72wAAAkSSURBVK9ARESkfWkgmIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITEWSAMzsFjObYWbvm9mLZrZPFHGIiMRZVE8Av3T3g939UOAfwA0RxSEiEluRJAB337LTt4WARxGHiEicNXlR+LZmZrcB3wQ2A5pyQkSknaXtCcDMJpnZrDq+zgVw9+vdvT8wDvhhA+WMNbOpZjZ17dq16QpXRCR2zD3a1hczGwA86+4jGju2tLTUp06d2g5RiYjsPcxsmruX7vp5VL2Ahuz07TnAvCjiEBGJs6jeAfzCzIYCSWAJ8L2I4hARia1IEoC7azlJEZGIaSSwiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISExFmgDM7Edm5mbWI8o4RETiKLIEYGb9gVOBT6KKQUQkzqJ8ArgbuAbwCGMQEYmtSBKAmZ0DLHf3D6KoX0REICddBZvZJKBPHbuuB34KnNbEcsYCYwFKSkraLD4Rkbgz9/ZtgTGzkcBLQHnqo32BFcDh7r6qoXNLS0t96tSpaY5QRGTvYmbT3L1018/T9gRQH3efCfSq/d7MFgOl7r6uvWMREYkzjQMQEYmpdn8C2JW7D4w6BhGRONITgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMRUJAnAzG4ys+Vm9n7q68wo4hARibOcCOu+291/FWH9IiKxpiYgEZGYivIJ4Idm9k1gKnC1u2+s6yAzGwuMTX1bZmYftrC+HsC6Fp6bbpkaW6bGBZkbW6bGBZkbW6bGBZkbW3PjGlDXh+bubRPOrgWbTQL61LHremAyIXgHbgH6uvu30xLIjnimuntpOutoqUyNLVPjgsyNLVPjgsyNLVPjgsyNra3iStsTgLuf0pTjzOwB4B/pikNEROoWVS+gvjt9ez4wK4o4RETiLKp3AP9lZocSmoAWA5e3Q533t0MdLZWpsWVqXJC5sWVqXJC5sWVqXJC5sbVJXGl7ByAiIplN3UBFRGJKCUBEJKb2+gRgZg+b2Rozy6gXzWbW38xeNrO5ZjbbzP496phqmVkHM5tiZh+kYrs56ph2ZmbZZvaemWVU7zEzW2xmM1PTm0yNOp5aZtbVzMab2bzU79tRUccEYGZDd5oO5n0z22JmV0YdF4CZ/Ufqd3+WmT1qZh2ijqmWmf17Kq7Zrf157fXvAMzsOKAM+JO7j4g6nlqpnlB93X26mRUB04Dz3H1OxKFhZgYUunuZmeUCbwD/7u6TIw4NADO7CigFOrv7WVHHU8vMFgOl7p5RA4fM7BHgdXd/0MzygAJ33xR1XDszs2xgOXCEuy+JOJZ+hN/5A929wsz+Cjzn7n+MMi4AMxsBPAYcDlQBzwP/5u7zW1LeXv8E4O6vARuijmNX7r7S3aen/rwVmAv0izaqwIOy1Le5qa+MuFMws32BLwAPRh3LnsDMOgPHAQ8BuHtVpl38U04GFkR98d9JDtDRzHKAAmBFxPHUGg5Mdvdyd08ArxK60rfIXp8A9gRmNhAYBbwTbSQ7pJpZ3gfWABPdPVNi+2/gGiAZdSB1cOBFM5uWmsIkEwwG1gJ/SDWbPWhmhVEHVYeLgEejDgLA3ZcDvwI+AVYCm939xWij+tQs4Dgz625mBcCZQP+WFqYEEDEz6wQ8AVzp7luijqeWu9e4+6HAvsDhqUfPSJnZWcAad58WdSz1OMbdRwNnAD9INT9GLQcYDfyPu48CtgHXRhvSZ6Wapc4B/hZ1LABmVgycCwwC9gEKzezr0UYVuPtc4E5gIqH55wMg0dLylAAilGpffwIY5+5PRh1PXVLNBa8Ap0ccCsAxwDmptvbHgJPM7M/RhrSDu69IbdcAEwjttFFbBizb6QluPCEhZJIzgOnuvjrqQFJOARa5+1p3rwaeBI6OOKZPuftD7j7a3Y8jNG+3qP0flAAik3rR+hAw193vijqenZlZTzPrmvpzR8J/iHnRRgXufp277+vuAwlNBv9y94y4MzOzwtTLfFJNLKeRAVOcuPsqYKmZDU19dDIQeUeDXVxMhjT/pHwCHGlmBan/pycT3tFlBDPrldqWABfQip9dlNNBtwszexQ4AehhZsuAG939oWijAsLd7DeAmam2doCfuvtzEcZUqy/wSKpnRhbwV3fPqC6XGag3MCFcL8gB/uLuz0cb0qeuAMalmloWApdGHM+nUu3Yp9I+08E0ibu/Y2bjgemE5pX3yKwpIZ4ws+5ANfCD+qbSb4q9vhuoiIjUTU1AIiIxpQQgIhJTSgAiIjGlBCAiElNKACIiMaUEICISU0oAImlgZr8zs+lmdljUsYjURwlApI2lRgL3IgxuypjpqkV2pQQg0kJmNtDMKnYayQ2Au28jjKZ+Bfht6tiOqUVPqsysR/tHK7I7JQCRJrJg1/8zC1Kzpu58XHfCHPJbgRoAd69IHZcp88qLKAGINCR1lz/XzO4lzA3TlLnXf0aYT342cGA64xNpDSUAkcYNJSwpOqqxFatSi/scDTxOmEHyoLRHJ9JCSgAijVvSjPWQbwV+7mGWRSUAyWh7/XTQIm1gW1MOMrNDCfOzH2tmvwc6ADPTGZhIaygBiLSdO4Gz3f0lADPrTZhLXiQjqQlIpA2Y2UlAYe3FHyC1xGGhmXWLLjKR+ukJQKQB7r4YGNGE4/4F/KuOz7ukISyRNqEnAJGWqwG67DoQrC61A8GAXCCZ9shEmkBLQoqIxJSeAEREYkoJQEQkppQARERiSglARCSmlABERGJKCUBEJKaUAEREYkoJQEQkpv4/5hYLOqmXrrcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "r=np.arange(1,9,0.1)\n",
    "V_AA=r*0\n",
    "V_BB=r*0\n",
    "V_AB=r*0\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "list_coeff=np.genfromtxt(\"CoeffLasso\"+\"_\"+str_type+\".txt\",usecols=[1])\n",
    "N_coeff=np.shape(list_coeff)[0]\n",
    "N_coeffLJ=np.shape(list_LJ)[0]\n",
    "N_coeffGauss=np.shape(list_Gauss)[0]\n",
    "\n",
    "print ('AA')\n",
    "for i_coeff in np.arange(N_coeffLJ):\n",
    "    a=list_LJ[i_coeff,0]*sigma\n",
    "    b=list_LJ[i_coeff,1]\n",
    "    V_AA=V_AA+LJ(r,a)*list_coeff[i_coeff+(N_coeff-3*N_coeffLJ)]\n",
    "    if (abs(list_coeff[i_coeff+(N_coeff-3*N_coeffLJ)])>1e-10):\n",
    "        print (\"LJ\",a/sigma,b,list_coeff[i_coeff+(N_coeff-3*N_coeffLJ)])\n",
    "\n",
    "for i_coeff in np.arange(N_coeffGauss):\n",
    "    a=list_Gauss[i_coeff,0]\n",
    "    b=list_Gauss[i_coeff,1]\n",
    "    V_AA=V_AA+gaussian(r,a,b)*list_coeff[i_coeff]\n",
    "    if (abs(list_coeff[i_coeff])>1e-10):\n",
    "        print (\"Gauss\",a,b,list_coeff[i_coeff])\n",
    "plt.plot(r,V_AA,\"red\")\n",
    "\n",
    "\n",
    "print ('BB')\n",
    "for i_coeff in np.arange(N_coeffLJ):\n",
    "    a=list_LJ[i_coeff,0]*sigma\n",
    "    b=list_LJ[i_coeff,1]\n",
    "    V_BB=V_BB+LJ(r,a)*list_coeff[i_coeff+N_coeffLJ+(N_coeff-3*N_coeffLJ)]\n",
    "    if (abs(list_coeff[i_coeff+N_coeffLJ+(N_coeff-3*N_coeffLJ)])>1e-10):\n",
    "        print (\"LJ\",a/sigma,b,list_coeff[i_coeff+N_coeffLJ+(N_coeff-3*N_coeffLJ)])    \n",
    "\n",
    "for i_coeff in np.arange(N_coeffGauss):\n",
    "    a=list_Gauss[i_coeff,0]\n",
    "    b=list_Gauss[i_coeff,1]\n",
    "    V_BB=V_BB+gaussian(r,a,b)*list_coeff[i_coeff+N_coeffGauss]\n",
    "    if (abs(list_coeff[i_coeff+N_coeffGauss])>1E-10):\n",
    "        print (\"Gauss\",a,b,list_coeff[i_coeff+N_coeffGauss])    \n",
    "plt.plot(r,V_BB,\"green\")\n",
    "\n",
    "\n",
    "print ('AB')\n",
    "for i_coeff in np.arange(N_coeffLJ):\n",
    "    a=list_LJ[i_coeff,0]*sigma\n",
    "    b=list_LJ[i_coeff,1]\n",
    "    V_AB=V_AB+LJ(r,a)*list_coeff[i_coeff+N_coeffLJ*2+(N_coeff-3*N_coeffLJ)]\n",
    "    if (abs(list_coeff[i_coeff+N_coeffLJ*2+(N_coeff-3*N_coeffLJ)])>1e-10):\n",
    "        print (\"LJ\",a/sigma,b,list_coeff[i_coeff+N_coeffLJ*2+(N_coeff-3*N_coeffLJ)])    \n",
    "\n",
    "for i_coeff in np.arange(N_coeffGauss):\n",
    "    a=list_Gauss[i_coeff,0]\n",
    "    b=list_Gauss[i_coeff,1]\n",
    "    V_AB=V_AB+gaussian(r,a,b)*list_coeff[i_coeff+N_coeffGauss*2]\n",
    "    if (abs(list_coeff[i_coeff+N_coeffGauss*2])>1e-10):\n",
    "        print (\"Gauss\",a,b,list_coeff[i_coeff+N_coeffGauss*2])    \n",
    "plt.plot(r,V_AB,\"blue\")\n",
    "\n",
    "\n",
    "\n",
    "plt.ylim([-5,1 ])\n",
    "plt.xlabel('r [$\\AA$]')\n",
    "plt.ylabel('V (r) [eV]')\n",
    "\n",
    "plt.plot(r,LJ(r,sigma),'r--')\n",
    "plt.plot(r,LJ(r,sigma*0.88)*0.5,'g--')\n",
    "plt.plot(r,LJ(r,sigma*0.80)*1.5,'b--')\n",
    "plt.show()"
   ]
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}