{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning atomic charges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "import quippy\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load a database of atomic charges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "atAll = quippy.AtomsList(\"data_GDB9.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a SOAP descriptor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Z = 1\n",
    "desc = quippy.Descriptor(\n",
    "        \"soap atom_sigma=0.3 n_max=10 l_max=10 cutoff=2.5 cutoff_transition_width=0.5 Z={:d} n_species=5 species_Z='1 6 7 8 9'\".format(Z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the descriptor for all input configurations. This will take some time to run..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d=[]\n",
    "c=[]\n",
    "for at in atAll:\n",
    "    at.set_cutoff(desc.cutoff())\n",
    "    at.calc_connect()\n",
    "    \n",
    "    descAt = quippy.fzeros((desc.dimensions(),desc.descriptor_sizes(at)[0]))\n",
    "    desc.calc(at,descriptor_out=descAt)\n",
    "    \n",
    "    for dd in descAt:\n",
    "        d.append(dd)\n",
    "    for i in quippy.frange(at.n):\n",
    "        if at.Z[i] == Z:\n",
    "            c.append(at.charge[i])\n",
    "d = np.array(d)\n",
    "c = np.array(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train a GP on a subset of the charge data using the SOAP descriptor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nTeach = 1000\n",
    "epsTeach = 0.01\n",
    "zeta = 1\n",
    "xTeach = d[:nTeach,:]\n",
    "yTeach = c[:nTeach]\n",
    "\n",
    "cTeach = np.dot(xTeach,xTeach.T)**zeta + epsTeach**2 * np.eye(nTeach)\n",
    "alpha = np.linalg.solve(cTeach,yTeach)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xTest = d[nTeach:,:]\n",
    "yTest = c[nTeach:]\n",
    "cTest = np.dot(xTest,xTeach.T)**zeta\n",
    "yTestPredict = np.dot(cTest,alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0143970580438\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAF5CAYAAAC83HEwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt8XHWd//HXJ2laKEhAqC1SbKm03FtobJXFqCxr0VWx\nuvhjg1bBFRdFZesFXbta0KV4QfjJrigW14JKVt2L6HqBxR+sFcWWpgWBQru2Dba0UKSES0ozTT6/\nP75nOpPpTDIzmcs5M+/n45FHM2fOmfnOoQ3vfL43c3dERERE4qyl3g0QERERGY0Ci4iIiMSeAouI\niIjEngKLiIiIxJ4Ci4iIiMSeAouIiIjEngKLiIiIxJ4Ci4iIiMSeAouIiIjEngKLiIiIxF4iA4uZ\nXWJmm81st5ndY2bzRjn/nWa2zsyeN7PHzOxbZvbiWrVXRERExiZxgcXMzgO+AiwFTgPuA24zsyMK\nnH8GcBOwHDgROBeYD3yzJg0WERGRMbOkbX5oZvcAv3P3S6PHBvwRuM7dv5Tn/I8BF7v7zKxjHwIu\nc/eX1ajZIiIiMgaJqrCYWRvQAfwyfcxD4roDOL3AZb8FjjazN0avMRl4B/DT6rZWREREKiVRgQU4\nAmgFHs85/jgwJd8F7v4b4F3A981sANgO7AI+VMV2ioiISAWNq3cDqs3MTgS+ClwO3A4cCVwN3AC8\nr8A1hwNnA1uAF2rRThERkQZxADAduM3d/1SpF01aYHkSGAQm5xyfDOwocM2ngLvd/Zro8QNm9kFg\npZktcffcag2EsPK9SjRYRESkSb0TuKVSL5aowOLuKTNbA5wF/Bj2Dbo9C7iuwGUTgYGcY0OAA1bg\nmi0A3/3udznhhBPG2OrkW7x4Mddee229m1F3ug8ZuheB7kOg+5ChewHr16/nXe96F0T/L62URAWW\nyDXAiii4rAIWE0LJCgAzuwp4qbu/Jzr/J8A3zexi4DbgpcC1hJlGhaoyLwCccMIJzJ07t1qfIzHa\n29t1H9B9yKZ7Eeg+BLoPGboXw1R0SEXiAou7/yBac+VzhK6gdcDZ7r4zOmUKcHTW+TeZ2cHAJYSx\nK08TZhl9qqYNFxERkbIlLrAAuPv1wPUFnrswz7GvAV+rdrtERESkOpI2rVlERESakAKLjKqrq6ve\nTYgF3YcM3YtA9yHQfcjQvaiexC3NXwtmNhdYs2bNGg2eEhERKUFPTw8dHR0AHe7eU6nXVYVFRERE\nYk+BRURERGJPgUVERERiT4FFREREYk+BRURERGJPgUVERERiT4FFREREYk+BRURERGJPgUVERERi\nT4FFREREYk+BRURERGJPgUVERERiT4FFREREYk+BRURERGJPgUVERERiT4FFREREYk+BRURERGJP\ngUVERERiT4FFREREYk+BRURERGJPgUVERERiT4FFREREYk+BRURERGJPgUVERERiT4FFREREYk+B\nRURERGJPgUVERERiT4FFREREYk+BRURERGJPgUVERERiL5GBxcwuMbPNZrbbzO4xs3mjnD/ezK40\nsy1m9oKZbTKzC2rUXBERERmjcfVuQKnM7DzgK8D7gVXAYuA2M5vl7k8WuOyHwCTgQuAPwJEkNKyJ\niIg0o8QFFkJAucHdbwYws4uBNwHvBb6Ue7KZvQHoBGa4+9PR4Udr1FYRERGpgERVGcysDegAfpk+\n5u4O3AGcXuCytwD3Ap80s61m9oiZfdnMDqh6g0VERKQiklZhOQJoBR7POf44cFyBa2YQKiwvAAuj\n1/g68GLgb6rTTBERaTYDgwO0WAvjWpL2v9ZkSFSFpUwtwBBwvrvf6+6/AD4KvMfMJtS3aSIi0gjW\n7VjH/OXzufo3V9e7KQ0raTHwSWAQmJxzfDKwo8A124Ft7v5c1rH1gAFTCYNw81q8eDHt7e3DjnV1\nddHV1VVis0VEpBENDA5w5a+uZNmvl3HSpJN4w7FvqHeTaqq7u5vu7u5hx/r6+qryXhaGgCSHmd0D\n/M7dL40eG2EQ7XXu/uU8518EXAu8xN37o2NvBf4NONjd9+S5Zi6wZs2aNcydO7d6H0ZERBJr3Y51\nXPCjC3hw54Ms6VzCpzs/zfjW8fVuVt319PTQ0dEB0OHuPZV63SR2CV0DXGRm7zaz44FvABOBFQBm\ndpWZ3ZR1/i3An4Bvm9kJZvYawmyib+ULKyIiIiMZGBxg6Z1Lmbc8LAG2+qLVXP66yxVWqixpXUK4\n+w/M7Ajgc4SuoHXA2e6+MzplCnB01vnPm9nrgX8CVhPCy/eBz9S04SIikniqqtRP4gILgLtfD1xf\n4LkL8xzbAJxd7XaJiEhjyh2rsvqi1Zw65dR6N6upJDKwiIiI1IqqKvGgwCIiIpKHqirxosAiIiKS\nQ1WV+FFgERERiaiqEl8KLCIiIqiqEncKLCIi0tRUVUkGBRYREWlaqqokhwKLiIg0HVVVkkeBRURE\nmoqqKsmkwCIiIk1BVZVkU2AREZGGp6pK8imwiIhIw1JVpXEosIiISENSVaWxKLCIiEhDUVWlMSmw\niIhIw1BVpXEpsIiISOKpqtL4FFhERCTRVFVpDgosIiKSSKqqNBcFFhERSRxVVZqPAouIiCSGqirN\nS4FFREQSQVWV5qbAIiIisaaqioACi4iIxJiqKpKmwCIiIrGjqorkUmAREZFYUVVF8lFgERGRWFBV\nRUaiwCIiInWnqoqMRoFFRETqRlUVKZYCi4iI1IWqKlIKBRYREakpVVWkHAosIiJSM6qqSLkUWERE\npOpUVZGxUmAREWkQCxYsore3r+Dz06a1c/vt36lhiwJVVaQSFFhERBpEb28fGzb8eIQzzqlZW0BV\nFamsRAYWM7sE+DgwBbgP+LC7ry7iujOAu4Dfu/vcqjZSRKQKcqsoW7duJ5UaYu/eZ3A/sI4tG05V\nFam0xAUWMzsP+ArwfmAVsBi4zcxmufuTI1zXDtwE3AFMrkVbRUQq7e6719LfPwNYD0zMeuZAoI1M\nFaUdGN79s3HjFo47LlNlqUYXkaoqUi2JCyyEgHKDu98MYGYXA28C3gt8aYTrvgF8DxgC3lrtRoqI\nVEMqNQH4MTAPGKmwvH/3j3tbTpdRZbuIVFWRamqpdwNKYWZtQAfwy/Qxd3dC1eT0Ea67EDgGuKLa\nbRQRqY2hejdgn4HBAZbeuZR5y+cBsPqi1Vz+ussVVqSiklZhOQJoBR7POf44cFy+C8xsJrAMeLW7\nD5lZdVsoIlIlCxYsIpXaS6iMDIxy9vYatEhVFamdpAWWkphZC6EbaKm7/yF9uNjrFy9eTHt7+7Bj\nXV1ddHV1Va6RIiJFCoNt74senTLK2Xuq2haNVRGA7u5uuru7hx3r6ys8tX4skhZYngQG2X/Q7GRg\nR57zXwS8AjjVzL4WHWsBzMwGgAXuflehN7v22muZO1eTiUQkjnJ/91oE5P6PoiPr+z7ggIq8s6oq\nkpbvl/ienh46OjoKXFG+RAUWd0+Z2RrgLMKoMyz08ZwFXJfnkmeAk3OOXQKcCfwVsKVqjRURqak+\noh+LBcwGpo/pHVRVkXpKVGCJXAOsiIJLelrzRGAFgJldBbzU3d8TDch9KPtiM3sCeMHd19e01SIi\ndeVjulpVFam3xAUWd/+BmR0BfI7QFbQOONvdd0anTAGOrlf7RERqo5/h05K3jXK+YbaeceNCqb6t\nbQ/Tpp026ruoqiJxkbjAAuDu1wPXF3juwlGuvQJNbxaRBNq6dROZkDKR4V1Ao6+p4n4CqVT4vq1t\n06iLxqmqInGSyMAiItKMMovGwfABtcWYTnbASaUKX6+qisSRAouISAIsWLCIvXtTWUeqs+6nqioS\nVwosIiIxt2DBIu644z7cs6cyH09pY1hGpqqKxJ0Ci4hIDKV3Zd66dRP9/a2ELp1thD2Espfl3wPM\niP6cTZgNNB5IMXwa8/BFMLOpqiJJoMAiIhJDvb190UaF5zB8cG3uAnHpysqM6M+HCFunbWPkdVlU\nVZFkUWAREYmhzIyg3pxnilkgbhuwl8zA3BbgSEKVJcwM8sn9zF8+X1UVSQwFFhGRGEl3BfX3txCC\nyejTlYebTuFAcw60DkDnleztfASYraqKJIYCi4hIjGS6gk4mhJXcwbTZa7Hks77w81PWMuG8Sexp\nf5YZW09m1dJVqqpIYiiwiIjESKYrqFCFJXstlnzmRH9mun/SVRU6f8LxL53DioUrVFWRxFFgERGJ\nkd27WwlhI3c35mKNY1jQmbIOFl4Akx6ElS9h1R2qqkgyVWflIRERKdmCBYsIe7b2MebfJ1uH4HVL\n4aJ54fHy1XDXZIUVSSxVWEREYmDBgkXcdddDZNZOSa9q206Y+dMGHEWY/TOKKetg4UqYdBusXAIr\nPw2DCiqSbGUFFjObSqg3voywQtE+7v7RCrRLRKSp9Pb2kUodlXVkiP3XXNnO8EXjcrQOQOd26JwH\nOyeGqsoOjVWRxlByYDGzswgdpJsIa0M/QPiVwICeSjZORKQZLFiwiI0bewk/kp8hhJKDgPsIP1qN\nTGXFCEEmZ6flfWNVdsDKz8LKe2EwN6x41T6DSLWVU2G5Crja3Zea2bPAXwFPAN8DflHJxomINIL0\n2ippW7duinZehr17n8H9QOD+6Nn0rKCRZgJlzRzaNwNoGew8CZZ3wo7LyT+1udyBvCL1V05gOQHo\nir7fCxzo7s+Z2WeBW4GvV6pxIiKNILO2Slr2cvvnENZOOZnMPIi2UV5xGzAHphwKC9fBpGdh5UxY\nORUGD8s65xxgC5lxMf1j+yAidVROYHmezLiV7cDLgQejx0dUolEiIo1rESFEnEro+hkCDmT4CrWj\nrG7beiR09kDn/bDzQFh+LOw4hPAjuSW6/kRCt9HsrNedXckPIlJT5QSWe4BXE34l+BnwFTM7BXh7\n9JyISFPL7LS8nVRqiFRqL5kQsh44ABgEfh8dO5kQNoqwbwbQs7DySFg5FwazV6go3JVkpi4hSa5y\nAstHgYOj75dG358HbIyeExFpWLnjUbJt3bod2E0qNYFUak10ND3TZxNhldoDCINfXwCOJVRXWoA9\nZEJN7nL85IxVmQjLj4MdExm+nFYvYcPDo6LvxxGK4rOB6Rx44GC5H1uk7koOLO6+Kev754GLK9oi\nEZEYu/vuh+nvXz3CGfMIOyNDCCQTCWNHJpKZ6dMSPU6vrZK1jD6wX5fQsNVql0QzgPJVZMYTgs/m\n6D3/jDDTCGbNgmnTTivuQ4rEUDnTmjcB89z9TznHDwV63H1GpRonIhI3qdQI66AAMECY8TObEEru\nJ3T5TB/hmrX5D+83Ayi9rso5hIXlstuyBdgNvIoQUk6K/uwHBnnkkZFmHYnEXzldQtOB1jzHJxB+\nVRARaRK5C7tB+LF6CKGKMkgILumNDAuZw/CqyqMw5QRY+FjODKDPRs+n9xqaQQgqc6LHJxAqNR3R\n+83m9a8/ndtvz1mzRSSBig4sZpb9r+lsM8v+V9oKnEX4lyMi0iT6yB9EZhN+LK6LHncU+Xrt0Pot\n6DwFOjfCzgOidVXaCeNa0r8T9pGZouyEkDInz+s5vb19LFiwSKFFEq+UCsuPoj8duCnnuRQhrHys\nAm0SEUmYfJWWPSVcH+2wPOW1sHA+TNoYrVZ7Kwy2R+ekpymnpUNQetCtE34Mp7uLTgZOY8OG77B1\n67wS2iIST0UHFndvATCzzYQxLE9WrVUiIomSW2mZw/DZO6NoHYLOpdC5EnbOzlqttofCXUkthKrL\nc1nvty7r+XNIB5zRx92IxF85s4SOqUZDRESSaRGZVWUBHiLMBjqkuMunrIOFj8CkZdFYlVUweG70\nZPbrbiGMUxlHqKCku4QOJHRB5Xb5pKc4QyqVYvz4Dl70okH+9Kd1iCRRubs1HwS8lvy7NV9XgXaJ\niMTOggWLSKWeIjNe5BnCWJV2QkAYJMwMcsKMnXTYyF44LtJ6MHTOjGYAtUUzgD5L+JGaDirZ3UC5\ny/nDyAN5pw17PpWCp56areAiiVXOtObTCCvcTiRsJ/oUYUn+fsImiAosItKQwoJxpwCPErpiXkSo\neuROkOwlVD4KBIop62Dha2FSf7Suyk/y7Ky8hbBOSweZikpYAC4EpPyL142sjVRqDbt25RugKxJv\nJXSy7nMt8BPgMDKT/qcBa4CPV65pIiLxElayhVBVOZGwxsp9hACRbVr+F2gdgNcthYvmASlY/mdw\nVw8Mnphz4lGEsLIm+gI4PTr2Y/bv/sknXaVJfy3a94x7EZeLxEw5XUKnAn/r7kNmNghMcPdNZnYZ\nYfbQf1S0hSIiMVF48Gq+6c3FrFb7XznXpENGehpz2nQy66sU66icNo2yoaJIzJUTWLKXV3yCMI5l\nPeFf7NEVapeISGMouFrt7OirLTqxheFjVmYTQsYWQjF7DmFsjEhzKiewrCVslrER+B/gc2Z2BKHe\n+EAF2yYiUlEjbVwIMG1ae2UXWNuvqvJpGEzPUxhi9B+ZWwgDew+P/pxAGEOT3jMou+KyNzqe3sco\nt5tKJNnKCSyfJow0A1gC3Ax8nRBg3luhdomIVFxvbx8bNow0s6ZC3SatA9D5CHTOy6mq5EpXTZww\neBdCKBkgBJoZZDZGTM8SSu/GnK8LSvsFSeMqZx2We7O+fwJ4Q0VbVAQzu4QwwHcKYcTbh9097/ap\nZvY24AOEsTcTgAeBy9399ho1V0Sayb6qSnq12uyqSrb0PrHbgRcIs4AOJFRIcndvroReMqFIJHnK\nWoelnszsPOArwPuBVcBi4DYzm1Vg9d3XALcDfw88TagC/cTM5rv7fTVqtog0gLa2PaRSm8g7wXK/\nsSpvhB09wLlkFn1Lc0IF5bTo8SZCWPl9Ea1oIQSa3GrQtjzn7tdIoIOWlueKOFckXooKLGa2liJH\ne7n73DG1aHSLgRvc/WYAM7sYeBMhiHwpT3sW5xxaYmZvBd5CqM6IiACwefM2jjsuEwRyx7RMnToj\n6lI6lrCi7ezwxJTd0R5A+caqQKYbJ9s2wlyFdkJYGXnvoba2bRxzzDls3bqbVOqhrOMtTJ16JJs3\nh8XhCl8/jmOOOSr6XLnTqEXir9gKy4+yvj8A+CDhX+tvo2OvAk4Crq9c0/ZnZulVlJalj7m7m9kd\nhEUKinkNI4zBeaoqjRSRxEqlYMOG7cBO4BA2bOjH7BTC72st0Z/HsG8YX+sQdO6EzieinZWnw45/\nA/496xqA5wkTLNNVlnFkAsp9mO3msMMO4ogjCo+hmTbtxBEHBIcBxeVfLxJ3RQUWd78i/b2Z3Qhc\n5+6fyT7HzK6g+tOajyDUNB/POf44cFyRr/EJwgq9P6hgu0SkIaQHs84hFGDzDWSdDdyfNVblKVi5\nNKqq/E10XRth1k56QO1B0Z/9hEXCU4QwswV4lpkzT+GRR8Y2YFZhRBpdOWNY3gG8Is/x7wL3EuOZ\nQmZ2PvAZ4JxidptevHgx7e3DpwZ2dXXR1dVVpRaKSOy1erSzcu66KhC6eO4f4eLZgNPSYhx77PR9\nR6dN0xRkSabu7m66u7uHHevrK2fbiNGVE1h2A2cQpjFnO4Mw1L2aniTsLjY55/hkYMdIF5rZXwPf\nBM519zuLebNrr72WuXOrPSRHRKppwYJF3H33WlKpCdHGhccSqhyFLCr81JR1sPB/o52V841VGY3h\nPlKgEUmWfL/E9/T00NFRyqrMxSknsPxf4OtmNpcwSwfglYTKyucr1bB83D1lZmuAs4jqtNGYlLMY\nYdNFM+sCbgTOc/dfVLONIhIvvb199PdD6II5hNxdjPeXZxzIsBlA40ZYV0VEqqWcdVi+YGabgEuB\nd0WH1wMXunstxoVcA6yIgkt6WvNEYAWAmV0FvNTd3xM9Pj967iPAajNLV2d2u/szNWiviNRAehXb\nrVu3D9vzJ5UaIISV3YSlmEazjTDeJKq07Lda7b/n2VlZRKqtrHVYomBSl0Gr7v6DaCuAzxG6gtYB\nZ7v7zuiUKQwf/HsRYaDu16KvtJuI8XgbEckv3/L6W7dup79/F2GKcAv7TyHOXjull/2rKNkLtUUD\nb1tPg84deVar/c8KfhoRKVbiFo4DcPfrKTCF2t0vzHl8Zk0aJSJlK2WPn8LL66dn9JzK/ouo7QYO\nJjNzJ/f535MJMdtgymth4cMwac8oq9WKSK0kMrCISGPJH0IWEWbdbGLDhlbGj+9g795ncHcy4WIT\noYCava5lK7Am57XSYSb9mrmiFddaD4bON0djVdqyVqs9g7AybYqwpkqh9U6KWW1WRMqhwCIiMdVH\nCBnnAO2kUvcR1q3MXuJ+AnBidO626Ny9RbxmHlPWwcLXwqT+aKzKrTD40zwnLiKsm5nudtpOJsz0\nE3ZTzl2G3wBn3LjdI7RNREaiwCIiMXYq4X/448i/vsk57B9CSpxOOWwG0MSssSr/QeGxLvkWlBt+\nbNasc8a8GJyIZCiwiEhdZM/q6e8fAE4mdOdA6HZJr5XSRggls0lXKsLX7qxzZkfnPUPYl6eDUPFo\ni57fG52TmT0E5JkB9AMYfBvw4ui9crt4egkh6mVj+uwiUrpiNz+8ptgXdPePlt8cEWkW+49bya5Q\n5FYw8lU0Tgam5xx7Jjq2icwU5qfYt/cPLcCcaA+gJ8MsoJ0H5FRVWth/DEy2yi+IJSKjK7bCclrO\n47nRtY9Ej2cRVqAd6V+5iMg+W7duH+MrTAAeZXg1Jj0DqDXrvBcxrDtpX1XlyWgPoFuz1lWZTtjf\nJ4ybyUx1FpF6K3bzw31Tg83so8CzwHvcfVd07DDg28DKajRSRBpP9uJu+ys2zLQCJ1C4MgP7KiLD\nxqpkr6vyk5zzp5MZ7Ju35YS1MudEj8cRqjLHF9lmESlHOWNYPgYsSIcVAHffZWb/ANwOfKVSjROR\nZjVSmMm1JfpzEaG6kj0WBsDzjFUZy7oqRghJaRpYK1IL5QSWQ4BJeY5PItNRLCJSAel1U9JTltM2\nE4JDP5nuoD5Cr3RWlaV1ADpPybNabbYUmXEpLYSwNMIGiDhhWjPA6WV8JhEpRzmB5T+Bb5vZxxi+\n+eGXCSPWRETKdCqhOpIOEXuB+0Y4fw7DF43Lsq+qsnGU1WrbyL/Q3EgsalsfZnOYOXPafmdMm9Y+\nymuISCnKCSwXA1cDtzB8zuC3gE9UqF0i0nQeBA5i+JTh0VaOnUaYapyldQg6l2aNVemEHZdXsqGE\nLqFtwI+ZOVPrrYjUQjm7NfcDHzSzTwAvjw7/wd2fr2jLRKShtbXtIZVKD1ydRggruYvDjVbp6GVY\nhWXKOli4EibdljVW5dwyW7iF/FOY9aNOpB7GsnDckdHXr9x9t5mZh00+RERGNXXqDDZs2EYo0P6Y\nsLBbqaKumNYB6LwHOufmrKsCYXpydvDpJfzo2xv9OVjgtY38KzXMBjZhtoeZM89R149IjZQcWMzs\ncOAHwJmEX21mElZp+paZ7XL3j1W2iSLSiKZNa2fDhnQ30GiVlBFMGYSFk2DSM9G6KvdmrasC+6+l\nkh6Um28KdNo2woDe/ZnBzJkzhu0gLSLVV06F5VrCiLiXERYjSPs+cA1h2rOIyD7pZfj/93+3MDSU\nXl4fQliBMOtmYv6LC2kdgM5HoHNTNFalPxqrMobws89RBZ+ZOXO6xqyI1EE5gWUBcLa7bzXL3pGU\njeyrz4qIZNx998P0968e4YwSu4PyzgB6RfRkugtoG/mDh7pwRJKonMByEPlrpS8G9oytOSLSiEZe\n1RbC2icteY5njz/phdajoXNj+Nr5Ilj+xqwZQOlJi+lumpG6fCAEmg4Kh5rRZiiJSC2VE1hWAu8G\nPhM9djNrAS4D7qxUw0SkWSwCXiB/l1DWGJEpM2HhVpi0aZR1VYqVIvzulb2rM4TBuADPEUJPLzAN\ns15mzpymQbYidVJOYLkM+KWZvQIYD3wJOIlQYTmjgm0TkZhKj0kpJHdAaiqVKnDmIsLCcIcxfMXZ\nLK1D0PkYdD4BO+cUWK22XAcRxtOMtG/rHGbNgmnTZmuQrUgdlbMOywNmNgv4EGETxIMJK9x+zd3H\nuv2qiMTcggWLuOuuh0ilCv9PfvPmDg466GRSqQmjvFofoboxBKzb/+l9Y1WegJWnwMqpMPjZ6Mle\nwu9M6e6m3FCUO515S9b3Q4SqTjuZikp+EyeO1yBbkRgoZ1rzy4A/uvuV+Z5z90cr0jIRiaXe3j5S\nqcKzaABSqaMIRZX0/+hLHFS7387Kx8OO3EXlTmV4WOmP3ic9GcAJOy+nB9/2Eyoq6eeNMG5m5B+D\nU6ceWVrbRaQqyukS2kxYMO6J7IPR+iybGb5NqohIafLurJxvk8HsisxsQrdSOsDsJVRuthHmAvQS\nisHpa04mTGyE0C0VKjFtbds45pjhYUxjVkTioZzAkr2IQraDCTVWEZEc/YRQsIlQ1ZjOfrNw9quq\nZI9VaSH/+ipbCD+OdhN+NKV/pA2S2WdoAiG0zCAEm+nAaVmvkRmXcswx2hdIJK6KDixmdk30rQOf\nN7Psqc2thB2b83RCi4gYIaCkpy7/mFDZeCg8nNIPC+fnVFWyZwAdSf4pyucA24HfF3jf3KnNJ5PZ\nLDF0U02cOGFft4+qKSLxVUqFJf0riQGnAANZzw0QhvpfXaF2iUhDOYQwjiR7s8LvQOuboXMNdO4Y\nZQZQL/tXWLYQKjf5uosKORDILGA3a5YqKiJJUXRgcfczAczs28Cl7v5M1VolIo1v387Kz8DKw3Jm\nAOXKXatyC/AMoSf6YTJhJnd129yKyQBtbWHqdFtbC9OmHT+WTyAiNVTOGJa/y3edmb0Y2KsgI9IM\ncqcMZ+sldLdkr9MySAgZu6N1VY6LVqudAMuPgh0HAH/MOt/JzOYZD7yK4ZsYnkwYVNtKGGib3vUZ\nRlrddtasaaqoiCRUOYHlX4FbgW/kHP8/hJ9gfznWRolIfIVxHoUXjXvySeOII/rYunUTqVSoZqRS\nDjwb7ay8AyalYOVLYOUkGNxDWOX2vqxXWZTzHg8xbIl+dhOmKD9HqLIMEaowLcAchs8LMFpanGOP\nna4xKiIJVk5geSWhypLrLmC/tVlEpLGUs9rrX5x9PncO3snQnz0JO8fB8lXRWJVTCTN4niezyu0e\n9t9XqJ+9C42hAAAaZ0lEQVRMBSW7iJteqfZZ4PDoujCAduLETTz//AMlt1VE4qmcwDKBUKPN1UYY\n0SYiss+6Het48m0PMbT9cfjVZ2HlChh8d/TsEJkfG+mVc3OrKxC6k4bIdBMdHP05h+FdRRlTpxbq\nshKRJMq3PepoVgHvz3P8YkbekENEmsjA4ABL71zKvOXzAHjV/W9k4qqfYkPthOCR7sIh+r4j+nqI\nMCYls9NHS8tu3H+P+/24b8b998yaNZ1CYUVEGk85FZZ/AO4wsznAL6NjZwHzgAWVapiIJNe6Heu4\n4EcX8ODOB1nSuYRPd36a8RcX3ll59M0UX1WNZopIgpSz+eHdZnY68AnCQNvdwP3A37j7xhEvFpGG\nNjA4wJW/upJlv17GSZNOYvVFqzl1yug7K2sXZBEZTTldQrj7Ond/p7uf5O6vcPf31jKsmNklZrbZ\nzHab2T1mNm+U819nZmvM7AUz22Bm76lVW0Waxbod65i/fD7Lfr2MJZ1LWHXRqqLCiohIMYqqsJjZ\nIen1VczskJHOrfY6LGZ2HvAVwjiaVcBi4DYzm+XuT+Y5fzrwX8D1wPnAXwA3mtlj7v7f1WyrSDMo\nt6oyVmGKcuGBtZrCLNJYiu0S2mVmR7r7E8DT5N/8ML0pYrV3a14M3ODuNwOY2cXAm4D3Al/Kc/4H\ngE3ufln0+BEze3X0OgosImOQd6xKa+GxKpWkbiSR5lJsYPlz4Kno+zOr1JZRmVkbYRrBsvQxd3cz\nu4PCG4q8Crgj59htwLVVaaRIE6hXVUVEmldRgcXd/yff93VwBKGC83jO8ceB4wpcM6XA+YeY2QR3\n31PZJoo0trXb13LBrRfw0M6Hal5VEZHmVewYltnFvqC7319+c0QkrlRVEZF6KrZLaB2Z3cjyjV/J\nVs0xLE8SdlGbnHN8MrCjwDU7Cpz/zGjVlcWLF9PePnzgXldXF11dXUU3WKQRqKoiIvl0d3fT3d09\n7FhfX+E1lcbC3EfLH2Bm07IengZcDXwZ+G107HTgY8Bl7v6jSjcypy33AL9z90ujxwY8Clzn7l/O\nc/4XgDe6+5ysY7cAh7p73o0azWwusGbNmjXMnTu3Gh9DJBFyqyorFq5QVUVERtTT00NHRwdAh7v3\nVOp1ix3D0pv+3sx+CHzE3X+Wdcr9ZvZH4PNAVQMLcA2wwszWkJnWPBFYEbXvKuCl7p5ea+UbwCVm\n9kXgXwir8p6LdpUWGZGqKiISJ+UszX8KsDnP8c3AiWNrzujc/QdmdgTwOULXzjrgbHffGZ0yBTg6\n6/wtZvYmwqygjwBbCavy5s4cEhE0VkVE4qmcwLIe+Hsze5+7DwCY2Xjg76Pnqs7drycsBJfvuQvz\nHPsVmb3rRaQAVVVEJK7KCSwXAz8BtppZekbQbMJg3LdUqmEiUjuqqohI3JWz+eEqM5sBvBM4Pjr8\nfeAWd3++ko0TkepTVUVEkqCcCgtRMPlmhdsiIjWkqoqIJElZgcXMFgF/C8wATnf3XjNbTNiz59ZK\nNlBEKk9VFRFJmpZSLzCzDxCmFv8cOIzMQnG7gL+rXNNEpNIGBgdYeudS5t84H8NYfdFqLn/d5Qor\nIhJ7JQcW4MPARe5+JbA36/i9hCnPIhJDa7evZd7yeSz79TKWdC5h1UWr1AUkIolRTpfQMcDaPMf3\nAAeNrTkiUmkaqyIijaCcwLIZOBXozTn+Bmq0DouIFEdjVUSkUZQTWK4BvmZmBxA2Q5xvZl2EhePe\nV8nGiUh5VFURkUZTzjosN5rZbuAfCXv43AI8Blzq7v9a4faJSIlUVRGRRlRSYIl2Rj4a+Hd3/56Z\nTQQOdvcnqtI6ESmaqioi0shKrbAY8L/AScBGd+8H+iveKhEpiaoqItLoSprW7O5DwEbg8Oo0R0RK\noXVVRKRZlDPo9lPAl83sA+7+QKUbJCLFUVVFRJpJOYHlZsJg2/vMbADYnf2ku7+4Eg0Tkfw0VkVE\nmlE5gWUx4JVuiIiMTlUVEWlW5UxrXlGFdojICFRVEZFmV3RgMbMW4OPAW4HxwC+BK9x994gXisiY\nqKoiIlLaLKElwDLgWWAbcCnwtWo0SkQ0A0hEJFspXULvBj7o7t8EMLO/AH5qZu+LpjuLSIWoqiIi\nMlwpFZaXAT9PP3D3OwiDb19a6UaJNCtVVURE8iulwjIOeCHnWApoq1xzRJqXqioiIoWVElgMWGFm\ne7KOHQB8w8yeTx9w97dXqnEizUAzgERERldKYLkpz7HvVqohIs1IVRURkeIUHVjc/cJqNkSkmaiq\nIiJSmnJWuhWRMVBVRUSkdAosIjWiqoqISPkUWERqQFUVEZGxUWARqSJVVUREKkOBRaRKVFUREakc\nBRaRClNVRUSk8hRYRCpIVRURkepQYBGpAFVVRESqq5TND+vOzA4zs++ZWZ+Z7TKzG83soBHOH2dm\nXzSz+83sOTPbZmY3mdmRtWy3NLa129cyb/k8lv16GUs6l7DqolUKKyIiFZaowALcApwAnAW8CXgN\ncMMI508ETgWuAE4D3gYcB9xa3WZKM9DOyiIitZOYLiEzOx44G+hw97XRsQ8DPzWzj7v7jtxr3P2Z\n6Jrs1/kQ8Dszm+ruW2vQdGlAGqsiIlJbSaqwnA7sSoeVyB2AA68s4XUOja55uoJtkyahqoqISH0k\npsICTAGeyD7g7oNm9lT03KjMbALwBeAWd3+u8k2URqaqiohI/dS9wmJmV5nZ0Ahfg2Y2qwLvMw74\nIaG68sExN1yahqoqIiL1F4cKy9XAt0c5ZxOwA3hJ9kEzawVeHD1XUFZYORr482KrK4sXL6a9vX3Y\nsa6uLrq6uoq5XBqAqioiIoV1d3fT3d097FhfX19V3svcvSovXGnRoNsHgVdkDbpdAPwMmJpv0G10\nTjqszADOdPeninivucCaNWvWMHfu3Ep9BEmQ3HVVVixcoanKIiJF6OnpoaOjA8IkmZ5KvW4cKixF\ncfeHzew2YLmZfQAYD/wT0J0dVszsYeCT7n5rFFb+nTC1+c1Am5lNjk59yt1Ttf0UkgSqqoiIxE9i\nAkvkfOCfCbODhoB/Ay7NOWcmkO7HOYoQVADWRX8aYRzLmcCvqtlYSRatVisiEl+JCizu/jTwrlHO\nac36vhdoHeF0EUBVFRGRuEtUYBGpNFVVRESSQYFFmpaqKiIiyaHAIk1HVRURkeRRYJGmoqqKiEgy\nKbBIU1BVRUQk2RRYpOGpqiIiknwKLNKwVFUREWkcCizSkFRVERFpLAos0lBUVRERaUwKLNIwVFUR\nEWlcCiySeKqqiIg0PgUWSTRVVUREmoMCiySSqioiIs1FgUUSR1UVEZHmo8AiiaGqiohI81JgkURQ\nVUVEpLkpsEisqaoiIiKgwCIxpqqKiIikKbBI7KiqIiIiuRRYJFZUVRERkXwUWCQWVFUREZGRKLBI\n3amqIiIio1FgkbpRVUVERIqlwCJ1oaqKiIiUQoFFakpVFRERKYcCi9SMqioiIlIuBRapOlVVRERk\nrBRYpKpUVRERkUpQYJGqUFVFREQqSYFFKk5VFRERqTQFFqkYVVVERKRaFFikIlRVERGRalJgkTFR\nVUVERGpBgUXKpqqKiIjUSku9G1AKMzvMzL5nZn1mtsvMbjSzg0q4/htmNmRmH6lmOxvdwOAAS+9c\nyvwb52MYqy9azeWvu1xhRUREqiZpFZZbgMnAWcB4YAVwA/Cu0S40s7cBrwS2VbF9DU9VFRERqYfE\nVFjM7HjgbOBv3P1ed/8N8GHgr81syijXHgV8FTgf2Fv1xjYgVVVERKSeklRhOR3Y5e5rs47dATih\ncnJrvovMzICbgS+5+/rwUEqhqoqIiNRbkgLLFOCJ7APuPmhmT0XPFfIpYMDd/7majWtEmgEkIiJx\nUffAYmZXAZ8c4RQHTijztTuAjwCnlXP94sWLaW9vH3asq6uLrq6ucl4uUVRVERGR0XR3d9Pd3T3s\nWF9fX1Xey9y9Ki9cdAPMDgcOH+W0TcAi4Gp333eumbUCLwDnuvt+XUJmdinwFULoSWsFhoBH3X1G\ngTbNBdasWbOGuXPnlvJxEi+3qrJi4QpVVUREpGg9PT10dHQAdLh7T6Vet+4VFnf/E/Cn0c4zs98C\nh5rZaVnjWM4CDPhdgctuBv4759jt0fFvl9fixqWqioiIxFXdA0ux3P1hM7sNWG5mHyBMa/4noNvd\nd6TPM7OHgU+6+63uvgvYlf06ZpYCdrj7xho2P9Y0VkVEROIuMYElcj7wz4TZQUPAvwGX5pwzE2in\nsPr2gcWMqioiIpIEiQos7v40oywS5+6tozyfd9xKs1FVRUREkiRRgUUqQ1UVERFJGgWWJqKqioiI\nJJUCS5NQVUVERJJMgaXBqaoiIiKNQIGlgamqIiIijUKBpQGpqiIiIo1GgaXBqKoiIiKNSIGlQaiq\nIiIijUyBpQGoqiIiIo1OgSXBVFUREZFmocCSUKqqiIhIM1FgSRhVVUREpBkpsCSIqioiItKsFFgS\nQFUVERFpdgosMaeqioiIiAJLbKmqIiIikqHAEkOqqoiIiAynwBIjqqqIiIjkp8ASE6qqiIiIFKbA\nUmeqqoiIiIxOgaWOVFUREREpjgJLHaiqIiIiUhoFlhpTVUVERKR0Ciw1oqqKiIhI+RRYakBVFRER\nkbFRYKkiVVVEREQqQ4GlSlRVERERqRwFlgpTVUVERKTyFFgqSFUVERGR6lBgqQBVVURERKpLgWWM\nVFURERGpPgWWMqmqIiIiUjsKLGVQVUVERKS2WurdgFKY2WFm9j0z6zOzXWZ2o5kdVMR1J5jZrWb2\ntJk9Z2a/M7Oppb7/wOAAS+9cyvwb52MYqy9azeWvu7zhw0p3d3e9mxALug8ZuheB7kOg+5Che1E9\niQoswC3ACcBZwJuA1wA3jHSBmb0cWAk8FJ1/CvB54IVS3njt9rXMWz6PZb9expLOJay6aFXTdAHp\nH2Cg+5ChexHoPgS6Dxm6F9WTmC4hMzseOBvocPe10bEPAz81s4+7+44Cl/4j8FN3//usY5uLfV+N\nVREREam/JFVYTgd2pcNK5A7AgVfmu8DMjFCJ2WhmvzCzx83sHjN7azFv+PDOh5u2qiIiIhInSQos\nU4Ansg+4+yDwVPRcPi8BDgY+CfwMeD3wn8B/mFnnaG/47h+9u6nGqoiIiMRV3buEzOwqQqAoxAnj\nVsqRDmQ/cvfrou/vN7M/Ay4mjG3J5wCAtx7+Vi7ruIyhx4boeaynzCYkX19fHz09zfv503QfMnQv\nAt2HQPchQ/cC1q9fn/72gEq+rrl7JV+v9AaYHQ4cPsppm4BFwNXuvu9cM2slDJ49191vzfPabcDz\nwOXuvizr+BeAM9w9b5XFzM4HvlfqZxEREZF93unut1TqxepeYXH3PwF/Gu08M/stcKiZnZY1juUs\nwIDfFXjtlJmtBo7LeWoW0DvC290GvBPYQomziURERJrcAcB0wv9LK6buFZZSmNnPCONSPgCMB/4F\nWOXui7LOeRj4ZLriYmYLgX8FPgTcCbwRuAZ4rbv/trafQERERMqRpEG3AOcDDxNmB/0X8Cvgb3PO\nmQm0px+4+48I41UuA+4H3gu8XWFFREQkORJVYREREZHmlLQKi4iIiDQhBRYRERGJPQWWSL03VoyL\ncu9D1vXfMLMhM/tINdtZC6XeCzMbZ2ZfNLP7o78L28zsJjM7spbtHiszu8TMNpvZ7mhl6HmjnP86\nM1tjZi+Y2QYze0+t2lptpdwLM3ubmd1uZk9Ef2d+Y2YLatneain170TWdWeYWcrMGmZhkjL+fYw3\nsyvNbEv0b2STmV1Qo+ZWTRn34Z1mts7Mnjezx8zsW2b24lLeU4Elo24bK8ZMyfchzczeRtgmYVvV\nWldbpd6LicCpwBXAacDbCFPq91sjKK7M7DzgK8BSwme4D7jNzI4ocP50wgD4XwJzgK8CN5rZ62vR\n3moq9V4Q/n7cTpiJOJcwK/EnZjanBs2tmjLuQ/q6duAmwiSJhlDmvfghcCZwIWFJjS7gkSo3tarK\n+DlxBuHvwnLgROBcYD7wzZLe2N2b/gs4HhgCTss6djawF5gywnXdwE31bn+970N03lHAo4T/wW8G\nPlLvz1Ove5HzOq8ABoGp9f5MRbb3HuCrWY8N2ApcVuD8LwL35xzrBn5W789S63tR4DUeAP6h3p+l\nHvch+ntwBeF/aj31/hz1uBfAGwjbxxxa77bX+T58DNiYc+xDwKOlvK8qLEHNN1aMqZLvA+y7FzcD\nX3L39YXOS5iy7kUeh0bXPF3BtlVFtDJ0B6FaAoCHnyx3EO5HPq9i/9+gbxvh/EQo817kvoYBLyL8\nDyuRyr0PZnYhcAwhsDSEMu/FW4B7gU+a2VYze8TMvmxmFV2yvpbKvA+/BY42szdGrzEZeAfw01Le\nW4ElqPnGijFVzn0A+BQw4O7/XMW21Vq592IfM5sAfAG4xd2fq3gLK+8IoBV4POf44xT+zFMKnH9I\n9PmTqpx7kesTwEHADyrYrlor+T6Y2UxgGWFZ9qHqNq+myvk7MQPoBE4CFgKXErpDvlalNtZCyffB\n3X8DvAv4vpkNANuBXYQqS9EaOrCY2VXRANBCX4NmNqvMlx+2saK73+/uXyT0519cmU9QGdW8D2bW\nAXyE0D8be1X+O5H9PuMIfdcOfHDMDZdEsbAf2WeAd7j7k/VuT62YWQthH7al7v6H9OE6NqneWghd\ny+e7+73u/gvgo8B7Eh7oS2JmJxLGt11OGN91NqECV9T4yLS67yVUZVcD3x7lnE3ADkLFZB8LGyu+\nOHounycJ4xlyu0DWA2eU3NLqquZ9eDUwCfhjqIADIX1fY2Z/5+4zym10lVTzXqTPS4eVo4E/T0h1\nBcLf6UFgcs7xyRT+zDsKnP+Mu++pbPNqqpx7AYCZ/TVhMOG57n5ndZpXM6XehxcRxm2dambpKkIL\noYdsAFjg7ndVqa3VVs7fie3AtpyfAesJIW4q8Ie8V8VbOffhU8Dd7n5N9PgBM/sgsNLMlrh7brUm\nr4YOLB7PjRVrrpr3gTB25b9zjt0eHR8tGNRcle9FdliZAZzp7rvG3uraiP5OryF8zh/DvnEYZwHX\nFbjst4RZMdkWRMcTq8x7gZl1ATcC50W/TSdaGffhGeDknGOXEGbJ/BVhQ9lEKvPvxN3AuWY20d37\no2PHEaouW6vc5Koo8z5MBAZyjg0RKtDFV+DqPdo4Ll+EcSj3AvMIFZJHgO/knPMw8NasxwsJU5jf\nB7yc0B83AJxe789Ty/uQ5zUSP0uonHtB+AXgVkJgPYXwG0f6q63en6fIz/x/gH7g3YSZUjcQAt6k\n6PmryJoZR9iR9VnCbKHjCN1fA8Bf1Puz1OFenB999otz/tsfUu/PUsv7kOf6RpolVOrfiYOinwff\nJ8ygfE30c+Qb9f4sNb4P7wH2RP82jol+nq4CflPS+9b7g8flizCb47tAH2Ew0HJgYs45g8C7c45d\nAGwAngd6gDfX+7PU4z7kPL+JxggsJd0LYFr0OPtrKPrzNfX+PCV87g8SfhPeTaiUvCLruW8D/y/n\n/NcAa6LzNwKL6v0Z6nEvCOuu5P73HwT+pd6fo9Z/J3KubZjAUs69IFTdbwOeI4SXLwET6v056nAf\nLgF+H92HrYR1WY4s5T21+aGIiIjEXkPPEhIREZHGoMAiIiIisafAIiIiIrGnwCIiIiKxp8AiIiIi\nsafAIiIiIrGnwCIiIiKxp8AiIiIisafAIiIiIrGnwCIiDcfMfmtmy+rdDhGpHAUWESmZmQ2Z2WD0\nZ+7XoJl9tgLv0W1mt1SivUW819lR28fX4v1EpHTj6t0AEUmkKVnf/zVwBWGTt/RW8c/VvEVjY5S6\n1b2I1JQqLCJSMnd/Iv1F2M3a3X1n1vF+ADObY2a3mdlzZvaYmX3LzA5Nv46ZdZnZA2a228x2mtkv\nzKzNzK4CzgPOy6razM/XFjM72Mxuid7jj2b2oTznXGhma8zs2agdN5nZ4dFzxwE/i07dHb3X9dFz\nbzazu83s6ah9PzKzaZW8lyJSHAUWEamKKBD8P+DXwKnAm4BjgO9Gz78MuBn4J0J15kzgJ9Hl/wjc\nGn1NBo4E1hR4q+uAecBfAm8E3gyclHNOK/BJ4GTg7cDxwA3RcxuB86PvXxa912XR4wOBL0Ttfz3Q\nBvyw2HsgIpWjLiERqZa/A37l7p9PHzCz9wMbzGwqcBShC+Y/o0rNH4EHolNTZvYCUeWm0BuY2WHA\nIuCt7v6r6Ni7gUezz3P3G7Me9prZx4A7zWycu+81s13Rc0+4+0DWdcPCiZn9LfComc1w903F3woR\nGStVWESkWuYAfxl1wzxrZs8CawljRV4OrAbuBh4xs381s/ea2SElvsdMws+xVekDUfgZFibM7JVm\n9l9m1mtmzwC/iK6bOtKLm9lxZvZ9M9sUXbc+av/LSmyniIyRAouIVMvBhO6T2YTwkv6aCfzO3fe6\n+2sJXTiPAIuBh83spZVshJm1Az8HdgBdQAdhoDDAaLOCfk7oFroQeAXwGkJVSLOJRGpMgUVEqqUH\nONndN7v7ppyvF9Inufvd7r4UOI0w1uSc6KmB6PFINgJDwCvTB8zsJcCMrHNOAtqBT7r7b9x9I8Nn\nOaXfi+z3i4LTdOAKd/8fd98AHE6osIhIjSmwiEi1fBWYambfNbMOM5thZm80s28BmNmrzewyM5sb\nDcB9B3Ao8FB0/RbgVDM71swON7P9wou77wK+A1xrZq8xs9nAt4E9WadtAfYCl5rZdDN7O2EALjnn\nALzFzI4ws4nATsIMqIujtr8e+OIY74mIlEmBRUSqwt3/CJwBHAT8N3A/8GXgyeiUp4GzCN0u64F/\nAD6YHjwLfJ0weHYt8AShKyefS4F7CVOTf04Yn/JgVjseA95HGJz7EPAR4GM5bd0MXEmYcbQDuNrd\nU4QupDMIg4Gvyr1ORGrH3FXdFBERkXhThUVERERiT4FFREREYk+BRURERGJPgUVERERiT4FFRERE\nYk+BRURERGJPgUVERERiT4FFREREYk+BRURERGJPgUVERERiT4FFREREYu//A+CpM9eiY1urAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdeafe1b110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(yTest,yTestPredict,\"s\")\n",
    "plt.plot([-0.6,0.6],[-0.6,0.6])\n",
    "plt.xlabel(\"Test data\")\n",
    "plt.ylabel(\"Predicted data\")\n",
    "print np.sqrt(np.mean((yTest-yTestPredict)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def likelihood(atIn,atom_sigma=0.4,n_max=10,l_max=10,cutoff=3.0,cutoff_transition_width=0.5,Z=1,zeta=1,eps=0.1):\n",
    "    desc = quippy.Descriptor(\"soap atom_sigma={:f} n_max={:d} l_max={:d} cutoff={:f} cutoff_transition_width={:f} Z={:d} n_species=5 species_Z='1 6 7 8 9'\".format(atom_sigma,n_max,l_max,cutoff,cutoff_transition_width,Z))\n",
    "\n",
    "    d=[]\n",
    "    c=[]\n",
    "    for at in atIn:\n",
    "        at.set_cutoff(desc.cutoff())\n",
    "        at.calc_connect()\n",
    "\n",
    "        descAt = quippy.fzeros((desc.dimensions(),desc.descriptor_sizes(at)[0]))\n",
    "        desc.calc(at,descriptor_out=descAt)\n",
    "\n",
    "        for dd in descAt:\n",
    "            d.append(dd)\n",
    "        for i in quippy.frange(at.n):\n",
    "            if at.Z[i] == Z:\n",
    "                c.append(at.charge[i])\n",
    "    d = np.array(d)\n",
    "    c = np.array(c)\n",
    "    \n",
    "    cTeach = np.dot(d,d.T)**zeta + eps**2 * np.eye(len(c))\n",
    "    alpha = np.linalg.solve(cTeach,c)\n",
    "    \n",
    "    l = -0.5 * np.dot(alpha,c) - np.sum(np.log(np.diag(np.linalg.cholesky(cTeach))))\n",
    "    return l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAFkCAYAAACNTikJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8VNX9//HXJws7BlxY3HEDqhVIWEVAZVcBEUQCFq1o\n3YqKWq1VxOVrLbZq/Vm3FlrRCQEBFUUBBaqg4JawuICKiloUBdmkihpyfn+cSQ0xQGYykzvL+/l4\n5KHMvXPnM2NM3pxz7ueYcw4RERGRWMsIugARERFJTQoZIiIiEhcKGSIiIhIXChkiIiISFwoZIiIi\nEhcKGSIiIhIXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiISFxGFDDO72MxWmNnW8NcSM+sXPpZl\nZhPMbKWZbTezdWY22cya7+F6c8ys1MwGVni8sZkVhF9js5lNNLP60b1FERERCUKkIxmfAdcBuUAe\nsBCYZWatgXpAW+AWoB0wGGgJzKrsQmY2FtgJVLZ5yhSgNdATOA3oDjwcYa0iIiISIKvuBmlm9jVw\njXPuX5Ucaw+8BhzmnPtPucfbAk8D7YH1wBnOuafDx1oB7wJ5zrll4cf6As8CBzvn1lerYBEREakR\nUa/JMLMMMxuOH8FYupvTGuFHKraUe15doAC41Dn3VSXP6QJsLgsYYfPD1+kUbb0iIiJSs7IifYKZ\nHYcPFXWAb4DBzrnVlZxXG/gTMMU5t73coXuAl51zs3fzEs2AXcKHc26nmW0KH9tdXfsBfYG1wI4q\nvyERERGpAxwOzHPOfR2ri0YcMoDVQBsgBxgKPGpm3csHDTPLAqbjRx8uLff4QOAU/NqNWOuLHyER\nERGR6IzEr4uMiYhDhnOuBPgo/MdlZtYRuAK4BHYJGIcAp1QYxTgZOALYamblL/uEmS1yzp2CX6PR\npPxBM8sE9g0f2521AKFQiNatW0f6tqSKSkpKmDFjNpMmvcimTX8FxgL3kJNzJb/5zUkMHXo6WVnR\nZFepqrFjx3LPPfcEXUZa0Wde8/SZ16xVq1ZxzjnnQPh3aazE4rdBBlAbdgkYRwAnO+c2Vzj3DuAf\nFR57Gx9SyqZPlgKNzKxduXUZPQHDLyLdnR0ArVu3Jjc3N8q3IlXRsWNHZs4cyKZNufgBrVy2bm3E\nBx/cSnY2tGsXdIWpLScnR9/jNUyfec3TZx6YmC43iChkmNkfgTnAp0BD/LBKD6BPOGDMxE+FnA5k\nm1nT8FM3Oed+DC/0/KrCNQE+c859AuCcW21m84B/mNklQC3gPqBQd5Ykjm3bSoCJ1Kq1lMzMiTRs\nWMJbb0FuLgwaBOPHK2yIiKS7SO8uaQJMxq/LmI/vldHHObcQOAgfLg4GlgOfA1+E/9llD9es7B7a\nEeVeYzawCLgowloljkpKWtCmjdGrV0cmTDBOOqkFq1fDI4/AO+/4sHHGGbBs2V4vJSIiKSqikOGc\nu8A5d4Rzrq5zrplzrixg4Jz7xDmXWeErI/zPRXu4ZmZZj4xyj21xzp3jnMtxzjV2zl3onPs2urco\nsfbuu7Bly/2MHz+azMxMxowZzbRp95OVBeeeC6tWweTJP4WNQYOguDjoqkVEpKZp7xKJWEEBNGoE\np54K+fn5PzuelQWjRvmw8eij/p95eTBwIBQVBVBwiqnsM5f40mde8/SZp4Zqd/xMFGaWCxQVFRVp\nsVAclZbCkUdCnz7wcBUbvZeUQGEh3HYbfPABDBjg12zk5cW3VhERqZri4mLy/A/lPOdczMaeNZIh\nEVmyBNauhZEjq/6crCz41a/8NMtjj8Hq1dC+vQ8bb74Zt1JFRCRgChkSkVAIDj0UTjwx8udmZcE5\n5/wUNt5/Hzp0UNgQEUlVChlSZT/8AI8/7kcxMqrxnVM+bIRCP4WN00+HN96IXb0iIhIshQypsjlz\nYPPmyKZK9iQz01/r3Xf9YtI1a6BjRzjtNHj99di8hoiIBEchQ6osFIK2beHYY2N73cxMGDHC3/Ja\nUAAffgidOilsiIgkO4UMqZKtW+GZZ/w0R7yUDxtTpsBHH/mwceqp8NqeGsqLiEhCUsiQKpk506/J\nGD48/q+VmQn5+fD22z5srF0LnTsrbIiIJBuFDKmSUAhOOQUOOqjmXrMsbLz1lu+zURY2+veHV1+t\nuTpERCQ6ChmyV//5D7z4YnynSvYkM9OPoLz1FkydCp98Al26KGyIiCQ6hQzZq8JCqF0bzjwz2Doy\nM+Hss38KG59+6sNGv36wdGmwtYmIyM8pZMhehUJ+35F99gm6Eq982Jg2DT77DE44Afr2VdgQEUkk\nChmyR2+9BStXBjdVsicZGTBs2E9hY926n8LGkiVBVyciIgoZskcFBbDffv4Xd6IqCxsrV/qOpOvW\nQdeufhM3hQ0RkeAoZMhulZb6kDFsGNSqFXQ1e5eRAWed5cPG9OnwxRc/hY1XXgm6OhGR9KOQIbu1\naJG/syQRp0r2JCMDhg6FFSt+Chsnngi9eytsiIjUJIUM2a2CAmjRwt/BkYzKh40ZM+DLL33Y6NUL\nXn456OpERFKfQoZUascOPwowciSYBV1N9WRkwJAhsHy5DxsbNkC3bj5sLF4cdHUiIqlLIUMq9eyz\nfr+SWO24mgjKwsayZb5N+oYN0L079OypsCEiEg8KGVKpggJo3x5atQq6ktjLyPCNxZYtgyeegI0b\nfwobixYFXZ2ISOpQyJCf2bTJj2Sk0ihGZTIyYPDgn8LG119Djx5+jxaFDRGR6lPIkJ+ZMQNKSmpm\nx9VEUBY2iovhySd9yCoLGy+9FHR1IiLJSyFDfqagwN/u2axZ0JXUrIwMOOOMn8LG5s1w0klw8sl+\ngzgREYmMQobs4pNP/FRBqk+V7En5sPHUU7Bliw8aJ52ksCEiEgmFDNnFlClQr56fPkh3ZjBokA8b\ns2bBtm0KGyIikVDIkP9xzu+4esYZ0KBB0NUkDjO/C21R0a5ho0cP+Pe//ecmIiI/p5Ah/7NiBbz7\nbnpPlexJ+bDx9NOwfbtfHHrSSQobIiKVUciQ/wmF4IAD/KJP2T0zGDAA3nzTh43//teHjR49YOFC\nhQ0RkTIKGQLAzp1+Pcbw4ZCdHXQ1yaEsbLzxBjzzDHz7rW/o1aMHLFigsCEiopAhgF/I+MUXmiqJ\nhhmcfroPG7Nnw3ff+X1RundX2BCR9KaQIYCfKjnqKOjYMehKkpcZnHYavP66Dxs7dviw0a0bzJ+v\nsCEi6UchQ/juO79h2DnnJP+Oq4mgfNh49ln44Qe/zkVhQ0TSjUKG8Mwz8M03miqJNTM49VR47TUf\nNn780YeNE0+EF15Q2BCR1KeQIYRC0KmTny6R2CsLG6++Cs895/eF6dNHYUNEUp9CRprbuBHmzPFT\nJRJfZtC/vw8bc+b4O3r69IGuXeH55xU2RCT1KGSkuenT/S+3s88OupL0YQb9+sHSpT5slJZC375w\nwgkwb57ChoikDoWMNBcK+V9wBxwQdCXpp3zYmDvXh4t+/RQ2RCR1KGSksY8+giVLNFUSNDMf9MrC\nBviw0aXLT+GjjHOOK6/8A04JRESSgEJGGpsyxW+ENmhQ0JUI/BQ2lizxIxkZGX4NR/mwUVRUxP33\n30dxcXHQ5YqI7JVCRpoq23F18GC/tbskDjO/IPSVV/yC0LKw0bkz3HDDdEpK7uLBB6cHXaaIyF5F\nFDLM7GIzW2FmW8NfS8ysX/hYlplNMLOVZrbdzNaZ2WQza17u+Y3N7P+Z2Woz+9bMPjGze81snwqv\n09jMCsKvsdnMJppZ/di8ZQG/k+h772mqJJGZ+b4avXrdQU5OS1as6M/zz78PXMjTT7/HUUf1o0mT\nltx00x1BlyoiUqlIRzI+A64DcoE8YCEwy8xaA/WAtsAtQDtgMNASmFXu+QcCzYGrgGOBc4F+wMQK\nrzMFaA30BE4DugMPR1ir7EFBATRr5ncPlcQ2btw13HbbtTRqlAU8CRgbNjzJli3ZjBt3LePGXRN0\niSIilYooZDjnnnXOzXXOfeicW+OcuxHYDnR2zm1zzvV1zs10zn3gnHsd+C2QZ2YHh5//jnPuLOfc\nc865j51zLwI3AAPMLAPAzFoBfYHRzrk3nXNLgDHAcDNrFrN3nsZKSqCw0O+4mpUVdDWyN9nZ2YwZ\nM5qcnF17vm/aZLRoMZpsbZsrIgkq6jUZZpZhZsPxIxhLd3NaI8ABW/ZwqUbANudcafjPXYDNzrll\n5c6ZH75Op2jrlZ8sWABffqmpkmRTWlpC3boTad68N3XrTqRevRIGDoQJE3S7q4gkpohDhpkdZ2bf\nAN8DDwCDnXOrKzmvNvAnYIpzbvturrU/cCO7ToU0A74qf55zbiewKXxMqqmgAFq1gtzcoCuRSOTm\ntmDCBGPNmllMmGCcemoLbrgBfv97GDXK7/oqIpJILNL77c0sCzgUyAGGAhcC3csHjfA5T+DXX5xc\nWcgws4b4EYoNwKBwkMDMrgdGOedaVzj/S+Am51ylazPMLBco6t69Ozk5Obscy8/PJz8/P6L3mar+\n+19o2tT/YrrxxqCrkViYOhV+/Wto0waefBKaN9/7c0QkfRUWFlJYWLjLY1u3bmXRokUAec65mN0j\nH3HI+NkFzF4A1jjnLgn/OQuYDhwOnOKc21zJcxoAzwPfAAOccz+UO/Zr4C/Ouf3KPZYJ7ACGOudm\nVbxe+JxcoKioqIhc/RV9t6ZM8butfvghHHFE0NVIrLz5pu93YgazZkFeXtAViUgyKS4uJs//4Ihp\nyIhFn4wMoDbsEjCOAHruJmA0xAeM74CB5QNG2FKgkZm1K/dYT8CA12JQb1oLhfyGXAoYqaV9ex80\nDjoIunWDadOCrkhEJPI+GX80s25mdlh4bcYdQA8gFA4YM/G3t54DZJtZ0/BXdvj5DYEX8ItFL8CH\nibJzMgDC0y7zgH+YWQcz6wrcBxQ659bH5m2np6++8s2dRo4MuhKJh+bN4aWXYMgQf+fQuHF+8zUR\nkaBEegNjE2Ayfq3FVmAl0Mc5t9DMDgNOD5+3PPxPw98VcjKwCB9AOoSPralwTgvg0/BjI4C/4dds\nlAIzgCsirFUqmDbND6cPGxZ0JRIvderAo4/CccfB9dfDO+/4PzdoEHRlIpKOIgoZzrkL9nDsEyBz\nL89/aW/nhM/bgh8NkRgKheDUU2G//fZ+riQvM7juOvjFL2DECD899vTTcNhhQVcmIulGe5ekiQ8+\ngNdf11RJOhkwwO/s+s030KEDvPxy0BWJSLpRyEgTBQXQsKH/xSPp47jjfLj8xS98C/lJk4KuSETS\niUJGGijbcXXoUKhbN+hqpKbtv79f8Hv++XDBBTB2rG8tLyISbwoZaeD1131fDLURT1+1asFDD8H9\n98N998Fpp8GWPTX7FxGJAYWMNBAKwYEHQo8eQVciQbv0Upg3D954Azp1gvffD7oiEUllChkp7scf\nfdvpESMgc6/39Ug66NnTj25lZPig8fzzQVckIqlKISPFvfACbNyoqRLZ1VFHwauvQpcu0L8/3Huv\ndnIVkdhTyEhxoRAceywcf3zQlUiiycmBZ56Bq66CK6+E3/wGfqjY5F9EpBoUMlLYN9/AU0/5UQyz\noKuRRJSZCX/+MzzyiO8M2qsXbNgQdFUikioUMlLYU0/Bd9/59Rgie3LuufDii34haIcOsHJl0BWJ\nSCpQyEhhoRB07w6HHhp0JZIMunTxd53suy+ccIIPqSIi1aGQkaLWr4f587XgUyJzyCGweLFfDDp4\nMNx+uxaEikj0FDJS1NSpkJXlu3yKRKJ+fb9j7803w403+um2774LuioRSUYKGSkqFPJdHRs3DroS\nSUYZGTB+PMyY4Xdw7dYN1q0LuioRSTYKGSlo1SooKtJUiVTfkCHwyivw1VfQvj289lrQFYlIMlHI\nSEEFBdCoEZx6atCVSCpo29YvCG3RwremLygIuiIRSRYKGSnGOf9LYOhQqFMn6GokVTRtCv/+N+Tn\n+xGy66+H0tKgqxKRRJcVdAESW0uWwNq1miqR2KtdG/75TzjuOLj2WnjnHR9oGzYMujIRSVQayUgx\noZC/DbFbt6ArkVRkBldfDbNnw0sv+d4aH30UdFUikqgUMlLIDz/A44/7Ww4z9F9W4qh/f7/B2vff\nQ8eOvluoiEhF+lWUQubOhU2bNFUiNaN1a3+3Sdu20Ls3PPxw0BWJSKJRyEghoRC0aePnzEVqwr77\nwpw5cPHF/uu3v4Uffwy6KhFJFAoZKWLrVr9t98iRQVci6SY7G+67Dx56yI9m9O/vR9RERBQyUsQT\nT/j58fz8oCuRdHXRRX6/nOXLoVMn3xRORNKbQkaKCIXg5JPh4IODrkTSWY8evnFXnTrQuTM891zQ\nFYlIkBQyUsC6db5RkhZ8SiJo0cL3a+nRA04/Hf7yF+3kKpKuFDJSQGEh1KoFZ54ZdCUiXsOG8NRT\n8Pvfw+9+B7/+tZ/OE5H0opCRAkIhGDgQcnKCrkTkJxkZ8Mc/+u/PqVP9dN769UFXJSI1SSEjyb39\nNqxYoakSSVwjR8KiRb7dfYcOsGxZ0BWJSE1RyEhyBQW+V0G/fkFXIrJ7HTv6BaHNmkHXrjB9etAV\niUhNUMhIYqWlPmQMG+bXZIgksoMO8iMagwb579mbb9ZOriKpTruwJrGXX4bPPtNUiSSPunVhyhT4\n5S/hhhv8dN/kyVC/ftCViUg8aCQjiYVCcPjhcMIJQVciUnVm8Ic/wJNP+v12TjwRPv006KpEJB4U\nMpLUjh1+x9WRI/0PbZFkc8YZvp/G5s1+QeiSJUFXJCKxppCRpJ57zu9Xor1KJJkdf7xfENqypb/F\n9ZFHgq5IRGJJISNJFRRAbq7fblskmR1wgN/zZNQo37Trmmtg586gqxKRWFDISEKbN8Ps2VrwKamj\nVi34+9/h3nvhnntgwAA/UiciyU0hIwnNmAElJTB8eNCViMSOGVx+uV8MunSp32Dtgw+CrkpEqkMh\nIwkVFEDPntC8edCViMRe797w2mu+h0anTrBgQdAViUi0FDKSzKefwksvaapEUtsxx/ig0bEj9O0L\nf/ubdnIVSUYRhQwzu9jMVpjZ1vDXEjPrFz6WZWYTzGylmW03s3VmNtnMmle4Rm0zu9/MNprZN2Y2\nw8yaVDinsZkVhF9js5lNNDO168E3MqpbFwYPDroSkfhq1MivPbr8chgzBi6+GH74IeiqRCQSkY5k\nfAZcB+QCecBCYJaZtQbqAW2BW4B2wGCgJTCrwjX+CpwGDAG6AwcCMyucMwVoDfQMn9sdeDjCWlOO\nc74B16BBfittkVSXlQV33w2TJsG//gV9+sDGjUFXJSJVFVHIcM4965yb65z70Dm3xjl3I7Ad6Oyc\n2+ac6+ucm+mc+8A59zrwWyDPzA4GMLN9gPOBsc65l5xzy4BfA13NrGP4nNZAX2C0c+5N59wSYAww\n3MyaxeqNJ6OVK+GddzRVIunn/PNh4UJ4910/hfL220FXJCJVEfWaDDPLMLPh+BGMpbs5rRHggC3h\nP+fh90v531Iu59x7wKdAl/BDnYHN4QBSZn74Op2irTcVhEKw//7+b3Mi6ebEE33jrn32gS5d4Omn\ng65IRPYm4pBhZseZ2TfA98ADwGDn3OpKzqsN/AmY4pzbHn64GfCDc25bhdO/DB8rO+er8gedczuB\nTeXOSTs7d0JhIZx9NmRnB12NSDAOO8xvDNi7t29L/qc/aUGoSCKLZiRjNdAG6Ag8CDxqZq3Kn2Bm\nWcB0/OjDpdUtUvwdJevWaapEpEED3ytm3Di4/nr//8R33wVdlYhUJuKt3p1zJcBH4T8uC6+luAK4\nBHYJGIcAp5QbxQBYD9Qys30qjGY0DR8rO6fi3SaZwL7lztmtsWPHkpOTs8tj+fn55OfnV+0NJqhQ\nCI480vcNEEl3GRlwyy1w7LFw3nm+addTT8GBBwZdmUjiKywspLCwcJfHtsapxa65ao41mtkC4BPn\n3PnlAsYRwMnOuU0Vzt0H2AAMd849GX6sJbAKv3j09fCoyDtA+7J1GWbWB3gOONg5V2nQMLNcoKio\nqIjc3NxqvadE89130KwZjB0LN98cdDUiiaWoyN9x5RzMmgXt2wddkUjyKS4uJi8vDyDPOVccq+tG\n2ifjj2bWzcwOC6/NuAPoAYTCAWMm/vbWc4BsM2sa/soGCI9eTALuNrOTzCwP+CfwSvhuFMLrO+YB\n/zCzDmbWFbgPKNxdwEh1s2fDtm3acVWkMnl5fkHoIYdAt24wdWrQFYlImUinS5oAk4HmwFZgJdDH\nObfQzA4DTg+ftzz8T8OvyzgZWBR+bCywE5gB1AbmApdVeJ0RwN/wd5WUhs+9IsJaU0Yo5G/bO/ro\noCsRSUzNm8OLL8JvfgP5+fDWW3DbbX5aRUSCE1HIcM5dsIdjnwCZVbjG9/i+F2P2cM4W/GhI2vv6\na3juOd+QSER2r04dmDwZjjsOfv9731PmscfUuE4kSMr5CW76dD/XfPbZQVcikvjM4NprfQ+NhQuh\na1dYuzboqkTSl0JGgguFfPOtJk32fq6IeKef7reL/+9/oUMHWLRo788RkdhTyEhgH38Mr7yi3hgi\n0Tj2WHj9dT990qsXTJwYdEUi6UchI4FNmQL16/vb80QkcvvtB88/DxdcABdeCFdcASUlQVclkj4i\nbsYlNaNsx9XBg33QEJHoZGfDAw/AL3/pt4xftQqmTYPGjYOuTCT1aSQjQRUXw+rVmioRiZVLLoEX\nXvDNuzp1gvfeC7oikdSnkJGgCgr8Ys+ePYOuRCR1nHyyX6eRleWDxrx5QVckktoUMhJQSYnfcTU/\n3/8wFJHYOfJIePVVv3X8qafCX/+qnVxF4kUhIwEtXAjr12uqRCRe9tnH73Ny9dV+T6ALLoDvvw+6\nKpHUo5CRgAoK4Jhj/J4MIhIfmZlw552+S2go5Kcmv/oq6KpEUotCRoL573/hiSf8KIZZ0NWIpL5R\no+Cll2DNGt+4a8WKoCsSSR0KGQnm6adh+3btuCpSkzp39ju57r8/nHACPPlk0BWJpAaFjARTUABd\nusARRwRdiUh6OeQQWLzYtyQ/80y/i6tz4Jzjyiv/gNPqUJGIKWQkkA0bYO5cLfgUCUq9ejB1Ktx6\nK9x0EwwfDq+8UsT9999HcXFx0OWJJB2FjAQybZpfhzFsWNCViKQvMxg3DmbMgNmz4YwzplNSchcP\nPjg96NJEko5CRgIpKIB+/fy8sIgE56ab7uCSS1qy33792bz5feBCZs58j6OO6keTJi256aY7gi5R\nJCmo1VOCWLPGNwiaOjXoSkRk3LhrOOCAJtx++1OUlvpVoFu2PMkPPwzg9tuv5bLLRgVcoUhy0EhG\ngigogIYNYcCAoCsRkezsbMaMGU1Ozq73kX/7rTFt2mjWrcsOqDKR5KKQkQDKdlwdMsQvPBORxFBa\nWkLduhNp3rw3detO5OCDS1i/Htq18/1sRGTPFDISwBtv+OkS9cYQSSy5uS2YMMFYs2YWEyYYJ5zQ\ngmXLfHfQIUPgt7+FHTuCrlIkcVmq3PttZrlAUVFREbm5uUGXE5HLL/cr2T/7zLc6FpHE5hw89JDf\n96RVK3j8cb8VgEiyKi4uJs/vZZHnnIvZ/doayQjYjz/6xZ4jRihgiCQLM7jkEnjtNfjuO8jNhcce\nC7oqkcSjkBGw+fN9Ey5NlYgknzZtoKjIT52MGgXnnee3BRARTyEjYKEQ/OIX0LZt0JWISDQaNPA7\nuU6e7Kc9O3SAlSuDrkokMShkBGj7dnjqKe24KpIKRo2CN9+EWrWgY0e/ZiNFlryJRE0hI0BPPQXf\nfuvXY4hI8mvVyjfVO/98v2bj7LNhy5agqxIJjkJGgEIh6NYNDjss6EpEJFbq1oUHHvBTJ88/73tq\nvP560FWJBEMhIyDr18MLL2jHVZFUNWQILFsGTZtC165w111QWhp0VSI1SyEjINOm+VtWhw4NuhIR\niZcWLWDxYt9P45pr/LYBGzcGXZVIzVHICEgoBKedBvvuG3QlIhJP2dlw553w3HN+2qRNG3jppaCr\nEqkZChkBeO89vwpdUyUi6aN/f1i+3HcGPeUUuPVW2Lkz6KpE4kshIwAFBZCT40cyRCR9HHSQb8B3\n001wyy3Quzd8/nnQVYnEj0JGDSvbcXXoUKhTJ+hqRKSmZWbC+PGwYAGsXu0b8c2dG3RVIvGhkFHD\nli6Fjz/WVIlIujvpJFixAtq391Mp113n9zISSSUKGTWsoAAOPhi6dw+6EhEJ2gEHwOzZ8Oc/w913\n+58La9cGXZVI7Chk1KAffvC3ro4YARn65EUE/7Pgmmvg5Zd9/5x27eCJJ4KuSiQ29KuuBs2bB19/\nrakSEfm5Tp18865evXwjr8sugx07gq5KpHoUMmpQKAS//KX/EhGpqFEjePxx35Z80iTo3Bnefz/o\nqkSip5BRQ7Ztg6ef1iiGiOyZmd9c7bXX/EhGbi489ljQVYlERyGjhjzxBHz/PeTnB12JiCSDNm18\n076hQ/028uedB9u3B12VSGQUMmpIKORvWTvkkKArEZFk0aABPPIITJ7sd3Xt0AFWrgy6KpGqiyhk\nmNnFZrbCzLaGv5aYWb9yxweb2Twz22hmpWZ2fCXXaGpmj5nZF2a23cyKzOzMCuc0NrOC8GtsNrOJ\nZlY/+rcZrM8/h4ULYeTIoCsRkWQ0ahQUFUGtWtCxIzz0kG/sJ5LoIh3J+Ay4DsgF8oCFwCwzax0+\nXh9YDFwL7O5/gceAo4HTgeOAJ4DHzaxNuXOmAK2BnsBpQHfg4QhrTRiFhf6Hw5AhQVciIsmqZUu/\nTmP0aL9mY9gw2LIl6KpE9iyikOGce9Y5N9c596Fzbo1z7kZgO9A5fDzknPs/YAFgu7lMF+A+51yR\nc26tc+52YAs+tBAOLH2B0c65N51zS4AxwHAzaxbNmwxaKOS3eG7UKOhKRCSZ1akD99/vp05eeMH3\n1Hj99aCrEtm9qNdkmFmGmQ0H6gFLI3jqK8DZ4SkRC1+jNvBi+HhnYLNzblm558zHj4x0irbeoLzz\njt95UVMlIhIrQ4b4nhpNm0LXrnDXXVBaGnRVIj8Xccgws+PM7Bvge+ABYLBzbnUElzgbqAV8Hb7G\ng+FrfBQoNKtpAAAfKElEQVQ+3gz4qvwTnHM7gU3hY0mloAAaN/Z7E4iIxEqLFrB4MYwd6zuGDhgA\nGzcGXZXIrrKieM5qoA2QAwwFHjWz7hEEjf8LP/cUfNA4A5huZic6596Jop5djB07lpycnF0ey8/P\nJz+Ae0dLS33IGDYMateu8ZcXkRSXnQ133gknn+wXh7ZpA1OmQI8eQVcmiaywsJDCwsJdHtu6dWtc\nXstcNZcom9kLwBrn3CXlHjsM+Bho65xbWe7xI4A1wLHOuVUVrvGBc+5SM/s18Bfn3H7ljmcCO4Ch\nzrlZu6kjFygqKioiNze3Wu8pVhYv9hseLVoE3boFXY2IpLLPP/fTsosWwU03wY03+m3lRaqiuLiY\nvLw8gDznXHGsrhuLPhkZ+DUVFVWWXuqFH99Z4fGd5WpZCjQys3bljvfELyR9rXql1qxQCA47zM+Z\niojE04EHwvz5PmDceiv07u2Dh0iQIu2T8Ucz62Zmh4XXZtwB9ABC4eONw7eiHosPBa3MrI2ZNQ1f\nYjXwIfB3M+tgZkeY2dVAL+BJgPC0yzzgH+FzugL3AYXOufXVf8s14/vv/R4EI0dqx1URqRmZmTB+\nPCxYAO+9B23bwty5QVcl6SzSX39NgMn4sDAff9tpH+fcwvDxgcAy4Bn8iEUhUAxcBOCcKwH6AxuA\np4EVwDnAKOfcvHKvM6Lca8wGFpVdI1nMmePvYdddJSJS0046yd/V1r69X3R+7bXw449BVyXpKKKF\nn865C/ZyfDI+hOzpnA+Bs/ZyzhZ8+EhaoZC/h/0Xvwi6EhFJRwccALNnw913w/XX+7UaU6fC4YcH\nXZmkEw3kx8GWLfDMM9pxVUSClZHhb299+WX48ks/fTJzZtBVSTpRyIiDmTP90OTw4UFXIiICnTr5\n5l29e/tdXS+7zG8jLxJvChlxEApBz55+tbeISCJo1MgvRn/gAZg0CTp39otDReJJISPGPv0UXnxR\nUyUiknjM/OZqr73mRzLy8uCxx4KuSlKZQkaMFRb6TYwGDw66EhGRyrVpA2++6adORo2C886D7duD\nrkpSkUJGjBUUwKBBsM8+QVciIrJ7DRrAI4/A5Ml+V9f27WHlyr0+TSQiChkxtHIlvPWWpkpEJHmM\nGgVFRX5/pY4d4aGHoJq7TYj8j0JGDIVCsN9+0Ldv0JWIiFRdy5Z+ncbo0X7NxrBh/lZ8kepSyIiR\n0lK/++HZZ/udEUVEkkmdOnD//X7q5IUXfDPB118PuipJdgoZMfLSS7BunaZKRCS5DRnie2o0beo3\nd7zrLv+XKJFoKGTESCgERxzh7z0XEUlmLVrA4sUwdqzvGDpgAGzYEHRVkowUMmJgxw4/xDhypL8P\nXUQk2WVnw513wnPP+WmTtm39iK1IJBQyYmD2bNi2TTuuikjq6d8fVqyAY46BU06BW26BnTuDrkqS\nhUJGDIRC0KGDX6EtIpJqDjwQ5s+Hm26CW2+FXr3g88+DrkqSgUJGNW3a5IcTNYohIqksMxPGj4eF\nC+H9933X0Llzg65KEp1CRjVNn+5XXmvHVRFJBz16wPLlvnFX//5w7bV+12mRyihkVFMo5LdPbto0\n6EpERGrGAQfAM8/AX/4C99wD3brB2rVBVyWJSCGjGtauhZdf1lSJiKSfjAy4+mr/M/DLL/3dJzNn\nBl2VJBqFjGqYMgXq1YMzzgi6EhGRYHTq5Jt39e7td3W97DJ/W78IKGREzTk/VTJ4sN/NUEQkXTVq\nBI8/Dg8+CJMm+aaE770XdFWSCBQyorR8OaxapakSERHwjQgvvthvtLZjB+TlwaOPBl2VBE0hI0qh\nkF/81Lt30JWIiCSONm3gzTf91Mm55/qv7duDrkqCopARhZ07/XqM/HzIygq6GhGRxNKgATzyiB/J\nmDkT2reHlSuDrkqCoJARhYULYf16TZWIiOzJr34FRUV+G/mOHf2aDeeCrkpqkkJGFAoK4OijfStx\nERHZvZYt4dVXYfRouPRSOOss2LIl6KqkpihkROjbb/3w3znnaMdVEZGqqFMH7r/f71Y9fz60a+d3\ndpXUp5ARoaef9ouYRowIuhIRkeQyZIi/M69pU+ja1XcMLS0NuiqJJ4WMCBUU+HvAjzoq6EpERJLP\n4YfD4sVw1VXwu9/B6afDhg1BVyXxopARgQ0b/K6D55wTdCUiIskrOxsmTPA7WL/xhm9J/uKLPx13\nznHllX/AaZVo0lPIiMDjj/uV0cOGBV2JiEjy698fVqyAY46Bnj3hllt8i4CioiLuv/8+iouLgy5R\nqkkhIwIFBdCvn2/CJSIi1XfggX4x6PjxcOut0KsX3HXXdEpK7uLBB6cHXZ5Uk0JGFX34ISxdqqkS\nEZFYy8yEkpI72Geflixe3J9p094HLmT27Pc46qh+NGnSkptuuiPoMiUKChlVVFDgu9gNHBh0JSIi\nqWfcuGu49dZr2XffLJx7EjC+/PJJtm/PZty4axk37pqgS5QoKGRUgXM+ZJx5pt/aXUREYis7O5sx\nY0bTuPGuDYh27jR++9vRZGdnB1SZVIdCRhW8+Sa8/76mSkRE4q20tIS6dSfSvHlvsrImsnFjCZdf\n7heESvJRyKiCUAiaNYNTTgm6EhGR1Jab24IJE4w1a2Zx991GXl4LHnjAb0j5/fdBVyeRslS5D9nM\ncoGioqIicnNzY3bdkhI46CC/Gdrdd8fssiIiUkVPPQXDh0OXLv7fc3KCrij1FBcXk5eXB5DnnIvZ\nvcMaydiL+fPhq680VSIiEpQzzvA/i5cvhx494Isvgq5IqkohYy9CIWjd2m/oIyIiwTjxRN+OfONG\nOOEEv05OEp9Cxh5s3w5PPumnSrTjqohIsI47zvcrqlvXb7CmnVwTn0LGHsya5bd2146rIiKJ4ZBD\n4OWXfSvyk0+GOXOCrkj2JKKQYWYXm9kKM9sa/lpiZv3KHR9sZvPMbKOZlZrZ8bu5ThczW2Bm28PX\nedHMapc73tjMCsLHNpvZRDOrH/3bjE4o5IfoWrSo6VcWEZHd2Xdfv0ajVy/fIPHRR4OuSHYn0pGM\nz4DrgFwgD1gIzDKz1uHj9YHFwLVApbetmFkXYA4wF2gf/vobUFrutClAa6AncBrQHXg4wlqr5csv\n4fnn/VSJiIgklrp1YeZMOO88OPdcuPNO3zhREktWJCc7556t8NCNZnYJ0BlY5ZwLAZjZYcDuVjHc\nDfzVOffnco99UPYvZtYK6Iu/jWZZ+LExwLNmdo1zbn0kNUdr2jTfT/+ss2ri1UREJFJZWfD3v0Pz\n5nDddf6uk7vuggwtBEgYUf+nMLMMMxsO1AOWVvE5BwCdgI1m9oqZrQ9PlXQtd1oXYHNZwAibjx8Z\n6RRtvZEKheDUU2G//WrqFUVEJFJmfvfWBx6Ae+/1a+jUtCtxRBwyzOw4M/sG+B54ABjsnFtdxacf\nEf7nePz0R1+gGFhgZkeGjzUDvir/JOfcTmBT+Fjcvf8+vPGGpkpERJLFJZfAjBm+Wddpp8G2bUFX\nJBDhdEnYaqANkAMMBR41s+5VDBploeYh51zZUp2rzKwncD5wQxT17GLs2LHkVGgHl5+fT35+fpWv\nUVAA++wDp59e3WpERKSmnHkmzJsHgwbBSSfBc8/5LSFkV4WFhRQWFu7y2NatW+PyWtVuK25mLwBr\nnHOXlHvsMOBjoK1zbmW5xw8HPgLOcc5NKff4VOBH59yvzOzXwF+cc/uVO54J7ACGOudm7aaOmLQV\ndw6OOsp/g06aFPVlREQkIG+9Bf36Qe3aPnQcfXTQFSW+RG4rngHUruTxn6UX59xa4HOgZYVDxwCf\nhP99KdDIzMr32OyJX0j6WnWL3ZtXX4WPPtJUiYhIsvrlL2HJEqhVyzftevPNoCtKX5H2yfijmXUz\ns8PCazPuAHoAZXeVNDazNsCx+FDQyszamFnTcpf5M3C5mQ0xsyPN7DZ86JgEEJ52mQf8w8w6hBeF\n3gcU1sSdJQUFfkO0Hj3i/UoiIhIvhx0Gr7wCRx7pR6affz7oitJTpCMZTYDJ+HUZ8/G9Mvo45xaG\njw8ElgHP4EcyCvELOy8qu4Bz7l7gDvytrMuBk4FezrmPy73OiHKvMRtYVP4a8fLjjzB1ql+dnJkZ\n71cTEZF42m8/37TrpJP8YtCCgqArSj+R9sm4YC/HJ+NDyN6ucydw5x6ObwFqfN/TefPg6681VSIi\nkirq1/d7UF10kd9Ne/16uPrqoKtKH9HcXZKyCgr8BjzHV9oMXUREklF2tl/I37w5XHMNfP45/PnP\natpVExQywrZt8/dX33yzdlwVEUk1ZnD77f6W1iuu8FtH/POffnGoxI9CRtiTT8KOHdpxVUQklY0Z\nA02bwq9+BRs2+AZeDRsGXVXq0mBRWEGBv6PkkEOCrkREROJp2DCYOxeWLvXbxX/11d6fI9FRyMDP\nzy1Y4BcFiYhI6jv5ZFi0CNatgxNOgA8/DLqi1KSQgb9tNSsLhg4NuhIREakpbdv6pl0ZGT5oFMes\nz6WUUcjA77h6+unQqFHQlYiISE1q0cI37Tr8cD9lPn9+0BWllrQPGe++C8uWaapERCRdHXAALFwI\n3brBqadChb3DpBrSPmQUFPgRjFNPDboSEREJSv36MGuWv8NwxAi4556gK0oNaX0La2mpDxlnneV3\n6xMRkfSVnQ3/+pfvpXHVVfDFF/CnP6lpV3WkdchYsgQ++URTJSIi4pn5YNG8OVx5pW9DPmmSDyAS\nubQOGaEQHHoonHhi0JWIiEgiueIK37Rr1CjfR2PGDGjQIOiqkk/aDgL98AM8/rife9NQmIiIVDR8\nOMyZ40e9TznFdwiVyKTtr9c5c2DzZk2ViIjI7vXsCS+9BJ9+Cl27wscfB11RcknbkBEK+UYsxx4b\ndCUiIpLI2rXzoxnO+aZdy5cHXVHySMuQsXUrPPMMjBwZdCUiIpIMjjjCN+066CDo3t331ZC9S8uQ\nMXOmX5ORnx90JSIikiyaNIEXX4QuXaB/f7+uT/YsLUNGKOQX8Rx0UNCViIhIMmnQwI+En3WWXxh6\n331BV5TY0u4W1v/8xyfRSZOCrkRERJJRrVrw6KO+adfll/umXbff7ntsyK7SLmQUFvrunmeeGXQl\nIiKSrDIy4C9/8U27rrnGB42//11NuypKu5ARCsHAgZCTE3QlIiKS7K6+2o9onHeeb9r1+ON+HxTx\n0mpNxltvwcqVuqtERERiZ+RIePZZ30+jZ0/YuDHoihJHWoWMggLYd1/o1y/oSkREJJX06ePX+330\nkd+qYu3aoCtKDGkTMsp2XD37bL9oR0REJJbat/dNu3780TftWrky6IqClzYhY9Eif2eJ2oiLiEi8\nHHWUDxrNm0O3bn50I52lTcgIhaBFC99ERUREJF6aNvXhomNH6NvX7+CartIiZOzY4f8jjxyp+5hF\nRCT+Gjb0i0GHDIFhw+CBB4KuKBhpcQvrs8/6/Up0V4mIiNSUWrX8KHrTpnDZZfD553Dbben1l920\nCBmhEOTlQatWQVciIiLpJCMD7r4bDjwQrr0W1q+Hhx6CrLT47ZsGIWPTJnjuOZgwIehKREQkHZnB\n737nm3adf75v2jV1KtSrF3Rl8ZfyazJmzICSEr+RjYiISFB+9Su/udqCBdCrF3z9ddAVxV/Kh4xQ\nyP/HbNYs6EpERCTd9esH//43fPCBv8X100+Drii+UjpkfPIJLF6s3hgiIpI4OnaEV16B777zTbve\nfjvoiuInpUPGlCl+zmvw4KArERER+ckxx/imXfvv70c0Fi8OuqL4SNmQ4ZyfKhk0CBo0CLoaERGR\nXTVv7jdVa9cOeveGp54KuqLYS9mQsWIFvPuupkpERCRx5eTAnDn+L8RDhsDDDwddUWyl7C2soRAc\ncIBPhyIiIomqdm0oLPRNuy6+GL74AsaPT42mXSkZMnbu9Osxzj4bsrODrkZERGTPMjLg3nv9FMof\n/uCDxgMPQGZm0JVVT0qGjBdf9P+BNFUiIiLJwgyuv94HjQsu8E27pkyBunWDrix6KbkmIxTy2+12\n7Bh0JSIiIpE57zyYNQvmzfNT/ps2BV1R9CIKGWZ2sZmtMLOt4a8lZtav3PHBZjbPzDaaWamZHb+X\n680JnzewwuONzawg/BqbzWyimdWvSo07dsDMmdpxVUREktdpp8HChbBqlb/F9bPPgq4oOpGOZHwG\nXAfkAnnAQmCWmbUOH68PLAauBdyeLmRmY4GduzlvCtAa6AmcBnQHqrTmdtEi+OYb7bgqIiLJrXNn\n37Rr+3bftOvdd4OuKHIRhQzn3LPOubnOuQ+dc2ucczcC24HO4eMh59z/AQuA3Y4jmFlbYCxwfsXz\nzKwV0BcY7Zx70zm3BBgDDDezvTYHf+456NQJjj46kncmIiKSeFq1gqVLoXFjOPFEHzqSSdRrMsws\nw8yGA/WApRE8ry5QAFzqnPuqklO6AJudc8vKPTYfP+LRaW/XX7JEoxgiIpI6DjzQj9Iff7zfi+vp\np4OuqOoiDhlmdpyZfQN8DzwADHbOrY7gEvcALzvnZu/meDNgl/DhnNsJbAof26PSUsfZZ0dQjYiI\nSIJr1AjmzvVrNQYPhokTg66oaqIZyVgNtAE6Ag8Cj4anOPYqvMDzFPxUSVwcf/xqmjSJ19VFRESC\nUacOTJvmG3ZdeCHcdpvfQiORRdwnwzlXAnwU/uMyM+sIXAFcUoWnnwwcAWy1XW/9eMLMFjnnTgHW\nA7vEBDPLBPYNH9ujtWt/z8CB03Z5LD8/n/z8/CqUJyIikrgyM+Fvf/NTKDfe6HtC3XdfZE27CgsL\nKSws3OWxrVu3xrhSLxbNuDKA2pU8Xlm+ugP4R4XH3saHlLLpk6VAIzNrV25dRk/8AtHX9lpMxlG8\n++4PbNv2MRdffB633np9Vd6DiIhIUjCDG26AZs3gN7+BL7+EggI/0lEVlf3Fu7i4mLy8vJjXGlHI\nMLM/AnOAT4GGwEigB9AnfLwxcChwED4UtDI/ZLHeOfdleKHnVxWuCfCZc+4TAOfcajObB/zDzC4B\nagH3AYXOub2OZGzefBe1ao1n3LhrufjiUZG8PRERkaQxejQ0aQLDhkHfvr6BV6NGQVe1q0jXZDQB\nJuPXZczH98ro45xbGD4+EFgGPIMfySgEioGL9nDNykY8RpR7jdnAor1cYxc5OcaYMaPJ1sYlIiKS\nwgYMgAUL4K23oHt3WLcu6Ip2FdFIhnPugr0cn4wPIZFc82czSc65LUBUO4/Urv0kpaUl0TxVREQk\n6Zxwgu+f0bev//d583x/jUSQcnuXXH65kZvbIugyREREakzr1r5PVMOG0LUrvPpq0BV5KRcyhg8/\ng2nT7g+6DBERkRp18MGweDH84hdwyinw7LNBV5SCIUNERCRdNW4Mzz/vp04GDYJ//jPYehQyRERE\nUkjdujBjBlxwgb8D5fbbg2vaFYs+GSIiIpJAMjPhwQehefOfmnbde29kTbtiQSFDREQkBZnB+PG+\nadell/qmXY89VvWmXbGg6RIREZEUdtFFMHMmPPMM9O8PceogXimFDBERkRR3xhkwfz4sX+6bdn3+\nec28rkKGiIhIGjjxRH+L69df+6Zd770X/9dUyBAREUkTxx3nm3bVq+ebdr0W3nbUxen2E4UMERGR\nNHLoofDyy9CypW/aNWcOrFq1Ki6vpbtLRERE0sy++/o1GsOH+03Wjj/+hbi8jkKGiIhIGrrjjjtY\nuvQR6tU7gmXLdsTlNTRdIiIikobGjbuGceOupV69LOCuuLyGQoaIiEgays7OZsyY0eTkWNxeQyFD\nREQkjZWWllC79pNxubZChoiISBrLzW3B5ZfHZzRDIUNERCSNTZt2P8OHnxGXaytkiIiISFwoZIiI\niEhcKGSIiIhIXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiI\nSFwoZIiIiEhcKGSIiIhIXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiISFwoZIiIiEhcKGSIiIhI\nXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiISFwoZEi1FBYWBl1C2tFnXvP0mdc8feapIaKQYWYX\nm9kKM9sa/lpiZv3KHR9sZvPMbKOZlZrZ8RWe39jM/p+ZrTazb83sEzO718z2qeS8gvBrbDaziWZW\nv3pvVeJBPwhqnj7zmqfPvObpM08NkY5kfAZcB+QCecBCYJaZtQ4frw8sBq4FXCXPPxBoDlwFHAuc\nC/QDJlY4bwrQGugJnAZ0Bx6OsFYREREJUFYkJzvnnq3w0I1mdgnQGVjlnAsBmNlhgFXy/HeAs8o9\n9LGZ3QA8ZmYZzrlSM2sF9AXynHPLwtcbAzxrZtc459ZHUrOIiIgEI+o1GWaWYWbDgXrA0mrU0AjY\n5pwrDf+5C7C5LGCEzcePjHSqxuuIiIhIDYpoJAPAzI7Dh4o6wDfAYOfc6mhe3Mz2B25k16mQZsBX\n5c9zzu00s03hY7tTB2DVqlXRlCJR2rp1K8XFxUGXkVb0mdc8feY1T595zSr3u7NOLK9rzlW2dGIP\nTzDLAg4FcoChwIVA9/JBIzxd8jHQ1jm3cjfXaYgfodgADHLO7Qw/fj0wyjnXusL5XwI3OecqXZth\nZiOAgojejIiIiJQ30jk3JVYXi3gkwzlXAnwU/uMyM+sIXAFcUtVrmFkDYB6wBTizLGCErQeaVDg/\nE9g3fGx35gEjgbXAjqrWIiIiItQBDsf/Lo2ZiENGJTKA2pU8XukQSXgEYx7wHTDQOfdDhVOWAo3M\nrF25dRk98QtJX9tdEc65r/F3pYiIiEjklsT6ghGFDDP7IzAH+BRoiB856AH0CR9vjJ9KOQgfClqZ\nmQHrnXNfhgPGC/jENBIfJsouv8E5V+qcW21m84B/hO9cqQXcBxTqzhIREZHkEdGaDDObCJyC73Wx\nFVgJ/Mk5tzB8/FzgX/x8FOMW59ytZtYD31tjl8uGz2/hnPs0fJ1GwN+AAUApMAO4wjn3bWRvT0RE\nRIIS8cJPERERkarQ3iUiIiISFwoZIiIiEhdJFTLM7DIz+9jMvjOzV82swx7ObRbeZO09M9tpZnfX\nZK2pIsLPfLCZPW9mX5XbQK9PTdabCiL8zLua2cvhTQm/NbNVZnZlTdabCiL5zCs8r6uZ/Whm6hoV\noQi/z3uEN90s/7XTzJrs7jnyc5F+n5tZLTO73czWmtkOM/vIzM6L5DWTJmSY2dnAXcB4oB2wApgX\n7hpamdr4zqG3ActrpMgUE8Vn3h14HuiP30Tv38AzZtamBspNCVF85v/F333VDWiF/37/PzO7oAbK\nTQlRfOZlz8sBJuObCkoEovzMHXA0vvNzM6C5c+6rPZwv5UT5mU8HTgZ+DRwD5APvRfS6ybLw08xe\nBV5zzl0R/rPhd4X9f865O/fy3H8Dy5xzV8W/0tRRnc+83DXeBqY65/4vfpWmjhh95jOB7c65c+NX\naeqI9jM3s0LgffwdcIOcc7k1UW8qiPQzL3dnYmPn3LYaLTZFRPGZ98P3njrCObcl2tdNipEMM8vG\nby2/oOwx59PRfPyGahJjsfjMw9/EDYFN8agx1cToM28XPvfFOJSYcqL9zM3s10AL4JZ415hqqvF9\nbsByM/s8PC17QnwrTR1RfuYDgDeB68zsP+GlB382s4j2NolFx8+asD+QCXxZ4fEvgZY1X05aiMVn\n/jugPvB4DOtKZVF/5mb2GXBA+Pk3O+f+FZcKU0/En7mZHQ38ETjROVdarqGgVE003+dfABfhf+nV\nxu+Z9aKZdXTOaTp876L5zI/AT8PuAM4IX+NB/BYfo6v6wskSMiTJmN+wbhy+dfzGoOtJAycCDYDO\nwAQzW+OcmxZwTSnHzDLwGzGOd859WPZwgCWlBefc+/ipqTKvmtmRwFhA04LxkYGfChzhnNsOYGZX\nAdPN7FLn3PdVuUiyhIyNwE6gaYXHm7LnTdMkelF/5mY2HPg7MNQ59+/4lJeSov7MnXOfhP/1HTNr\nBtwMKGTsXaSfeUOgPdDWzO4PP5aBnx38AejjnHsxTrWmilj9PH8d6BqrolJcNJ/5F8C6soARtgof\nqg8GPqz0WRUkxZoM59yPQBF+ozTgf/P9PYnDhi4S/WduZvnAJGC4c25uvOtMJTH8Ps+k8k0LpYIo\nPvNtwHFAW6BN+OshYHX433e7iaN4Mfw+b4v/RSh7EeVn/gpwoJnVK/dYS/zoxn8iefGk+AKGAd8C\no/C36j0MfA0cED5+BzC5wnPa4L8R3wAeC/+5ddDvJVm+Iv3MgRHAD8DF+IRc9rVP0O8lWb6i+Mwv\nBU4Hjgp/jcbvK3RL0O8lWb6i+dlS4fnjgeKg30cyfUXxfX4FMBA4EjgW+CvwI3BS0O8lWb6i+Mzr\nA5/gR0Rb41sUvAc8FMnrJst0Cc65x8P3896K/8W1HOjrnNsQPqUZcEiFpy3jp83acvG/BD/BL2iR\nvYjiM78Q/7fo+8NfZSYD58e/4uQXxWeegf/hcDhQgh/C/J1z7u81VnSSi/Jni1RDFJ95LXyPhwPx\nvyhXAj2dc4tqrurkFuln7pz7r5n1xvfheQMfSKbh19pVWdL0yRAREZHkkhRrMkRERCT5KGSIiIhI\nXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiISFwoZIiIiEhcKGSIiIhIXChkiIiISFwoZIiIiEhc\nKGSIiIhIXPx/jnhgVOalxLAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdec0718dd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "atTeach = atAll[:100]\n",
    "l=[]\n",
    "# eps = np.linspace(0.001,0.1,10)\n",
    "# for e in eps:\n",
    "#     l.append(likelihood(atTeach,eps=e))\n",
    "    \n",
    "# plt.plot(eps,l,\"-*\")\n",
    "# plt.show()\n",
    "\n",
    "# zeta = [1,2,3,4,5,6]\n",
    "# for z in zeta:\n",
    "#     l.append(likelihood(atTeach,zeta=z,eps=0.01))\n",
    "    \n",
    "# plt.plot(zeta,l,\"-*\")\n",
    "# plt.show()\n",
    "\n",
    "# cutoff = [1.5,2.0,2.5,3.0,3.5]\n",
    "# for c in cutoff:\n",
    "#     l.append(likelihood(atTeach,zeta=1,eps=0.01,cutoff=c))\n",
    "    \n",
    "# plt.plot(cutoff,l,\"-*\")\n",
    "# plt.show()\n",
    "\n",
    "atom_sigma = [0.1,0.2,0.3,0.4,0.5,0.6]\n",
    "for a in atom_sigma:\n",
    "    l.append(likelihood(atTeach,zeta=1,eps=0.01,cutoff=2.5,atom_sigma=a))\n",
    "    \n",
    "plt.plot(atom_sigma,l,\"-*\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}