examples/07. Force Matching.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 07. Force Matching \n",
    "\n",
    "In this tutorial, we show how to use force matching to fit a potential energy function to forces from a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd, hoomd.htf as htf, hoomd.md\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "import math\n",
    "tf.get_logger().setLevel('ERROR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the model graph\n",
    "\n",
    "We build an LJ potential with unknown, trainable parameters (`epsilon`, `sigma`) which start out at 0.9 and 1.1. We then obtain forces from our potential and the simulation. Force matching is used to modify the LJ potential until the forces agree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Build the graph\n",
    "N = 16\n",
    "NN = N - 1\n",
    "true_sigma = 1.0\n",
    "true_epsilon = 1.0\n",
    "graph = htf.graph_builder(NN, output_forces=False)\n",
    "# make trainable variables\n",
    "epsilon = tf.Variable(0.9, name='lj-epsilon')#, trainable=True)\n",
    "sigma = tf.Variable(1.1, name='lj-sigma')#, trainable=True)\n",
    "# get LJ potential using our variables\n",
    "# uses built in nlist_rinv which provides\n",
    "# r^-1 with each neighbor\n",
    "inv_r6 = sigma**6 * graph.nlist_rinv**6\n",
    "# use 2 * epsilon because nlist is double-counted\n",
    "p_energy = 4.0 / 2.0 * epsilon * (inv_r6**2 - inv_r6)\n",
    "# sum over pairs to get total energy\n",
    "energy = tf.reduce_sum(p_energy, axis=1, name='energy')\n",
    "# compute forces\n",
    "computed_forces = graph.compute_forces(energy)\n",
    "# compare hoomd-blue forces (graph.forces) with our \n",
    "# computed forces\n",
    "minimizer, loss = htf.force_matching(graph.forces[:,:3], \n",
    "                                     computed_forces[:,:3])\n",
    "# save loss so we can visualize later\n",
    "graph.save_tensor(loss, 'loss')\n",
    "# Make sure to have minimizer in out_nodes so that the force matching occurs!\n",
    "graph.save(model_directory='CG_tutorial/force_matching',\n",
    "           out_nodes=[minimizer])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the simulation\n",
    "\n",
    "The simulation consists of LJ particles whose sigma and epsilon values are 1.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HOOMD-blue v2.5.1 DOUBLE HPMC_MIXED SSE SSE2 \n",
      "Compiled: 03/04/2020\n",
      "Copyright (c) 2009-2019 The Regents of the University of Michigan.\n",
      "-----\n",
      "You are using HOOMD-blue. Please cite the following:\n",
      "* J A Anderson, C D Lorenz, and A Travesset. \"General purpose molecular dynamics\n",
      "  simulations fully implemented on graphics processing units\", Journal of\n",
      "  Computational Physics 227 (2008) 5342--5359\n",
      "* J Glaser, T D Nguyen, J A Anderson, P Liu, F Spiga, J A Millan, D C Morse, and\n",
      "  S C Glotzer. \"Strong scaling of general-purpose molecular dynamics simulations\n",
      "  on GPUs\", Computer Physics Communications 192 (2015) 97--107\n",
      "-----\n",
      "HOOMD-blue is running on the CPU\n",
      "notice(2): Started TF Session Manager.\n",
      "notice(2): Group \"all\" created containing 16 particles\n",
      "notice(2): Force mode is FORCE_MODE.hoomd2tf \n",
      "notice(2): Starting TensorflowCompute \n",
      "notice(2): completed reallocate\n",
      "notice(2): TF Session Manager has released control. Starting HOOMD updates\n",
      "notice(2): -- Neighborlist exclusion statistics -- :\n",
      "notice(2): Particles with 0 exclusions             : 16\n",
      "notice(2): Neighbors included by diameter          : no\n",
      "notice(2): Neighbors excluded when in the same body: no\n",
      "** starting run **\n",
      "Time 00:00:09 | Step 1000 / 1000 | TPS 127.328 | ETA 00:00:00\n",
      "Average TPS: 127.298\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "74 normal updates / 10 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 10 / n_neigh_max: 14 / n_neigh_avg: 11.375\n",
      "shortest rebuild period: 11\n",
      "-- Cell list stats:\n",
      "Dimension: 2, 2, 1\n",
      "n_min    : 3 / n_max: 6 / n_avg: 4\n",
      "** run complete **\n",
      "notice(2): Sending exit signal.\n",
      "notice(2): Shutting down TF Manually.\n",
      "notice(2): TF Queue is waiting, sending None\n"
     ]
    }
   ],
   "source": [
    "# run the simulation\n",
    "model_dir = 'CG_tutorial/force_matching'\n",
    "hoomd.context.initialize('--mode=cpu')\n",
    "with hoomd.htf.tfcompute(model_dir) as tfcompute:\n",
    "    rcut = 7.5\n",
    "    system = hoomd.init.create_lattice(\n",
    "        unitcell=hoomd.lattice.sq(a=4.0),\n",
    "        n=[int(math.sqrt(N)), int(math.sqrt(N))])\n",
    "    nlist = hoomd.md.nlist.cell(check_period=1)\n",
    "    lj = hoomd.md.pair.lj(r_cut=rcut, nlist=nlist)\n",
    "    lj.pair_coeff.set('A', 'A', epsilon=true_epsilon, sigma=true_sigma)\n",
    "    hoomd.md.integrate.mode_standard(dt=0.005)\n",
    "    hoomd.md.integrate.nve(group=hoomd.group.all(\n",
    "            )).randomize_velocities(kT=2, seed=2)\n",
    "    tfcompute.attach(nlist, r_cut=rcut, save_period=10)\n",
    "    hoomd.run(1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the trained variables and cost\n",
    "\n",
    "Here we load the training progress and fit LJ parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "time = np.arange(0, 1e3, 10)\n",
    "cost = np.empty(len(time))\n",
    "epsilon = np.empty(len(time))\n",
    "sigma = np.empty(len(time))\n",
    "for i, t in enumerate(time):\n",
    "    variables = hoomd.htf.load_variables(model_dir, names = ['loss', 'lj-epsilon', 'lj-sigma'], \n",
    "                                         checkpoint = int(t))\n",
    "    cost[i] = variables['loss']\n",
    "    epsilon[i] = variables['lj-epsilon']\n",
    "    sigma[i] = variables['lj-sigma']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting training progress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3de5zcdX3v8ddnbnvLZjfJbu6EBAxICMglIBetVkQitdg+RIWqpYrl9OHR2lNri0cf1Oo5rdXz0NoWPWK1Vo+F4qWaIhqLgFVukggEkhBIAiFXkmySTbK3uX3OH7/f7E42u5Od2ZlsfjPv5+ORBzu/+WX2O/sj897v9/P9fn/m7oiISOOKTXUDRERkaikIREQanIJARKTBKQhERBqcgkBEpMElproB5erq6vLFixdPdTNERCJl7dq1+929e6znIhcEixcvZs2aNVPdDBGRSDGzbeM9p6EhEZEGpyAQEWlwCgIRkQanIBARaXAKAhGRBqcgEBFpcAoCEZEGpyCIKHfn7jXbGcrmpropIhJxCoKI2rj7CH/+3XU8uGnfVDdFRCJOQRBRhZ5A70BmilsiIlGnIIioXD64s9zRwewUt0REok5BEFHZMAiOKAhEZJIUBBE13CMY0tCQiEyOgiCiMrk8oB6BiEyegiCiCj2CI0MKAhGZHAVBRKlGICLVoiCIqJFZQ6oRiMjkKAgiSj0CEakWBUFEZcNi8VHVCERkkhQEEaUegYhUi4IgokbWEWSHvxYRqYSCIKKyRR/+fWn1CkSkcgqCiMqFNQLQ8JCITI6CIKKKewTaeE5EJkNBEFHFQXBEawlEZBIUBBFVXCDWNhMiMhkKgojK5op7BAoCEamcgiCicvmRYrFqBCIyGQqCiMqoRiAiVaIgiKhc3mlKxIiZtpkQkclJTHUDpDLZnJOKx2hKxFQjEJFJURBEVC6fJx432lIJBYGITIqGhiIqm3cSMaO9OaEagYhMSk2DwMxWmtkmM9tsZreO8fwiM3vAzJ4ws3Vmdm0t21NPcnknHgaBagQiMhk1CwIziwO3A28GlgE3mtmyUad9Arjb3S8EbgC+VKv21JtMzknEYkxr0tCQiExOLXsElwKb3X2ru6eBu4C3jjrHgenh1x3Arhq2p67k8nkScaO9OamhIRGZlFoGwQJge9HjHeGxYp8E3m1mO4B7gQ+N9UJmdouZrTGzNfv27atFWyMnGw4NTdPQkIhM0lQXi28EvuHuC4FrgW+Z2XFtcvc73H2Fu6/o7u4+6Y08FeWKisWHNTQkIpNQyyDYCZxW9HhheKzYzcDdAO7+CNAMdNWwTXUj6BHEaG9KkM7mGcrmprpJIhJRtQyCx4GlZrbEzFIExeBVo855CbgKwMzOIQgCjf1MQDaXD3sESUD7DYlI5WoWBO6eBT4IrAY2EswOWm9mnzKz68LTPgL8oZk9BdwJ/IG76wa8EzBcI2gK1gSqTiAilarpymJ3v5egCFx87LairzcAV9ayDfUql3eS8aBGANqKWkQqN9XFYqlQ8awhUBCISOUUBBEVzBqKMT2sEWgtgYhUSkEQUdlcXjUCEakKBUFEFW86BxoaEpHKKQgiKpd3EvGRGoF6BCJSKQVBRGXDGkFTIk4qEeOwagQiUiEFQUQVtqEGaG9KaEGZiFRMQRBRmXBlMRDenEZBICKVURBEVHGPYJruUiYik6AgiKhs3knEg8vX3pRUsVhEKqYgiKjCNtRQ6BEoCESkMgqCiCosKAPVCERkchQEEVXcI2hvUo1ARCqnIIioTN6Jxws9gqBGoB28RaQSCoKIGl0jyDv0p3WXMhEpn4Iggtx9ePdRYHi/Ic0cEpFKKAgiKJcPhoCGewRNhY3nVCcQkfIpCCIoGwZBoUYwck8C9QhEpHwKggjKju4RaCtqEZkEBUEE5XJhj0A1AhGpAgVBBGXzeQCS8WNrBNqBVEQqoSCIoEKxuLCyOBXuOZTO5aesTSISXQqCCBpdI0iGQZBREIhIBRQEEZQdVSNIJmLHHBcRKYeCIIIKNYKRHkHwXw0NiUglFAQRNLygLAyAZExDQyJSOQVBBI2uEcRiRjxmCgIRqYiCIIJGZg2NXL5k3MioRiAiFVAQRFDhN/9CjwCCmUPqEYhIJRQEETR6HQEEawkUBCJSCQVBBI2uEUBQOM5kNTQkIuVTEETQyKyh4hqBegQiUhkFQQRlxxsayqtHICLlq2kQmNlKM9tkZpvN7NZxznmHmW0ws/Vm9q+1bE+9yOXHKRZn1SMQkfIlavXCZhYHbgeuBnYAj5vZKnffUHTOUuBjwJXuftDMZteqPfUkkzu+R5CIax2BiFSmlj2CS4HN7r7V3dPAXcBbR53zh8Dt7n4QwN331rA9dWP0ymIIegTaYkJEKlHLIFgAbC96vCM8Vuws4Cwze8jMHjWzlWO9kJndYmZrzGzNvn37atTc6BiZNTRy+VLxmDadE5GKTHWxOAEsBV4P3Ah81cw6R5/k7ne4+wp3X9Hd3X2Sm3jqGbNGkNDQkIhUppZBsBM4rejxwvBYsR3AKnfPuPsLwHMEwSAlZMeoEWj6qIhUqpZB8Diw1MyWmFkKuAFYNeqcHxD0BjCzLoKhoq01bFNdyI5RI0jEYqQ1NCQiFahZELh7FvggsBrYCNzt7uvN7FNmdl142mqgx8w2AA8AH3X3nlq1qV6MuY5AQ0MiUqGaTR8FcPd7gXtHHbut6GsH/jT8IxOUCz/wk7FjVxZnFQQiUoGpLhZLBYZ7BPHRNQINDYlI+RQEEZQbY9O5ZNy0jkBEKqIgiKCxagSaNSQilVIQRFBh+mjiuBqBhoZEpHwKgggqLCgr6hBoiwkRqZiCIIKyeScZN8yKt6EOpo8GE7FERCZOQRBBubwfUx+A4CY17iOFZBGRiVIQRFAm58fUByAYGio8JyJSDgVBBOXy+eN6BMlwTUEmrzqBiJRHQRBB2bwfs4YAIJUIewS6S5mIlElBEEG5vB+z4RyMTCXV0JCIlEtBEEFBj2B0jSAcGtIUUhEpk4IggsaaNTQ8NKQgEJEyKQgiKJPLH1cj0KwhEamUgiCCxuoRjASBegQiUh4FQQRl804ifuylKxSPtc2EiJRLQRBBubGmj8Y1fVREKjOhIDCzD5vZdAt8zcx+bWZvqnXjZGzZEkNDWW0xISJlmmiP4H3ufhh4EzADeA/wmZq1SkrKjlks1tCQiFRmokFQ+NS5FviWu68vOiYnWakegYaGRKRcEw2CtWb2U4IgWG1m7YA+cabIWCuLNX1URCqVmOB5NwMXAFvdvd/MZgLvrV2zpJRSK4uz2nRORMo00R7B5cAmdz9kZu8GPgH01q5ZUkouP/6CsrSGhkSkTBMNgi8D/Wb2KuAjwBbgmzVrlZSUzZXaYkJDQyJSnokGQdaDeyC+FfhHd78daK9ds6SU7Ji7j2rTORGpzERrBEfM7GME00Zfa2YxIFm7ZkkpwRYTo2oE2nRORCo00R7BO4EhgvUEe4CFwOdq1iopKZvPkxxvZbGGhkSkTBMKgvDD/9tAh5m9BRh0d9UIpkhujBqBNp0TkUpNdIuJdwC/At4OvAN4zMyur2XDZHyZMWoE8ZhhpiAQkfJNtEbwceASd98LYGbdwH3Ad2vVMBnfWNtQQ9Ar0BYTIlKuidYIYoUQCPWU8XelyoK9ho7/8afiMbKqEYhImSbaI/iJma0G7gwfvxO4tzZNkhMZaxtqCFYXa2hIRMo1oSBw94+a2duAK8NDd7j7v9euWVJKNu/E42MPDSkIRKRcE+0R4O7fA75Xw7bIBI3fI4iRzmpoSETKU3Kc38yOmNnhMf4cMbPDJ3pxM1tpZpvMbLOZ3VrivLeZmZvZikreRCNx93Ab6uMvnYaGRKQSJXsE7l7xNhJmFgduB64GdgCPm9kqd98w6rx24MPAY5V+r0aSC+9ANl6PQLuPiki5ajnz51Jgs7tvdfc0cBfBXkWjfRr4W2Cwhm2pG4VbUY47fVRDQyJSploGwQJge9HjHeGxYWZ2EXCau/+o1AuZ2S1mtsbM1uzbt6/6LY2QQo8gOWaxWENDIlK+KVsLEG5c93mCba1Lcvc73H2Fu6/o7u6ufeNOYSM9grFqBJo1JCLlq2UQ7AROK3q8MDxW0A4sBx40sxeBy4BVKhiXlg0/6MetEWhBmYiUqZZB8Diw1MyWmFkKuAFYVXjS3XvdvcvdF7v7YuBR4Dp3X1PDNkVerlSNIKEtJkSkfDULAnfPAh8EVgMbgbvdfb2ZfcrMrqvV96132RKzhlKqEYhIBSa8oKwS7n4vo7aicPfbxjn39bVsS70Ynj4aPz7DEzHVCESkfNo4LmJK9QiSiZhuTCMiZVMQREyhWDz2OgINDYlI+RQEEVO6RqChIREpn4IgYkrNGkrETUNDIlI2BUHEZIdXFo+zoCyrHoGIlEdBEDG5/Pg1glQ8RkabzolImRQEEVNYOTzeymINDYlIuRQEEXOi3UdzeR+uI4iITISCIGKGZw2Nsfto4ZhmDolIORQEETNSIzj+0qXCArKCoHr601ny6mFJnVMQREzpGoEdc45MTjqb58rP3M/3n9h54pNFIkxBEDG5EkNDyYR6BNV0ZDDDwf4ML/X0TXVTRGpKQRAxmVJ7DYXDRdqKujr607lj/itSrxQEEVOqRpBMFIrFGhqqhkIA9CkIpM4pCCLmROsIgnPUI6iGvnQWCArGIvVMQRAxJWsEcQ0NVVP/UNgjGFKPQOqbgiBiSi0oG5k+qqGhaij0BAYy6hFIfVMQRMzIzevHuEOZFpRV1XCNQD0CqXMKgog50RYTgHYgrRLVCKRRKAgiJldq+mghCLQStipUI5BGoSCImFJ7DaXUI6iqwtDQQEZBIPVNQRAxIz0C1QhqrTAk1DekoSGpbwqCiCkUi8cYGdL00Sor1AiGsnmtzZC6piCImGzeScQMs/GHhrTpXHUUby3Rr+EhqWMKgojJ5X3MGUNQvMWEfnuthv6iIvGAtpmQOqYgiJhCj2AsSd2PoKr6iqaNqk4g9UxBEDG5vJOIj33ZRnYf1dBQNRwzNKQegdQxBUHEZPP58XsEGhqqqv50jrZUHFCPQOqbgiBisrkSNQLtPlpV/eksXe1NwdcqFksdUxBETKkaQeG4hoaqo28oR9e0MAi0uljqmIIgYnJ5Jz7GqmIAMyMZNw0NVUl/Okt3GAR92m9I6piCIGKyeR8uCo8lGY9pi4kqyOedgUyO7sLQkGoEUscUBBGTy+fHrRFAEARZbTo3aYPZHO6MDA2pRiB1rKZBYGYrzWyTmW02s1vHeP5PzWyDma0zs5+Z2em1bE89yJQoFkMQBNpiYvIKO47OaEsSM9UIpL7VLAjMLA7cDrwZWAbcaGbLRp32BLDC3c8Hvgt8tlbtqRfBOoLxgyAVNw0NVUFhw7m2VIK2VEI1AqlrtewRXApsdvet7p4G7gLeWnyCuz/g7v3hw0eBhTVsT13I5p14iRpBIh5TsbgKCgvIWlNxWpvi6hFIXatlECwAthc93hEeG8/NwI/HesLMbjGzNWa2Zt++fVVsYvTk8nmSJYeGTPcsroJCj6C1KUFrKqEagdS1U6JYbGbvBlYAnxvreXe/w91XuPuK7u7uk9u4U0ypBWUQzhpSj2DSCjWCtlSc1lRcs4akriVq+No7gdOKHi8Mjx3DzN4IfBx4nbsP1bA9dSGbd5qT4+d3KqEgqIaRoSHVCKT+1bJH8Diw1MyWmFkKuAFYVXyCmV0IfAW4zt331rAtdeOENYKYhoaqYXhoqFAj0KZzUsdqFgTungU+CKwGNgJ3u/t6M/uUmV0XnvY5YBrwHTN70sxWjfNyEsqV2HQONH20WvoKPYKmcGhIQSB1rJZDQ7j7vcC9o47dVvT1G2v5/evRiWoEqURMO2VWQaEm0JYKi8X6mUodOyWKxTJxubyTLLGOICgWa2hosgo9gJZknLZUfLiHIFKPFAQRc6IagTadq47+dJaWZJxYzGhtSgzXDETqkYIgYkrdmAaCBWWqEUxeXzpHW1NwU5rWZJxMzklrxbbUKQVBxOROVCPQOoKq6B/K0poKSmitTcF/dQN7qVcKgogpdWMaCIaGsqoRTFp/OkdreJvK4dtVanhI6pSCIGJOtOmcVhZXR3EQFHoEqhNIvVIQREzQIyh9YxqNZU9eXzpLWxgArckgELSWQOqVgiBisrkT3ZhGK4uroX+ouEcQDg1pB1KpUwqCiDlxjSBGNq8ewWT1Z0aKxW0pDQ1JfVMQREwuP5HdRx139Qomo7hHUJhGqkVl1fPgpr3sOjQw1c2QkIIgQtw96BHES+8+Cmh4aJKKawQtqcL0UfUIqiGXd2755lr+/mfPT3VTJKQgiJDCPelLLigLn9PMocrl8s5gJn/89FHVCKpi75FB0rk863b0TnVTJKQgiJDCh/uJhoaKz5XyDWRGblMZ/Fc1gmraeTAYEnru5SMM6s5vpwQFQYTkwi5ByWKxhoYmrbDTaCEAUokYybipRlAlO8PaQDbvPLvnyBS3RkBBECnZMAhKbzGhoaHJKnzgF4rEEOxCqi0mqmPXocHhr5/ecWgKWyIFNb0fgVTXRHoEhcVmCoLKFYaAWpIj/zzamhK6z0OV7DzUT0dLkkTMVCc4RSgIIqSwPqDUrKGRoSEFQaX6x+gR6C5l1bPr0CALOluYPb2Jp3cqCE4FGhqKkMJmcqV6BCNDQ6oRVKpvVI0Awh6BisVVsfPgAPM7WzhvQQfP7z2qIbdTgIIgQnITqBFo1tDkjdUjaEmqR1Atuw4NsHBGEAS5vLNh9+GpblLDUxBESKFYfKLdR0FBMBmFD/zWUTUCTR+dvN6BDEeGsszvbOb8hZ0APKPhoSmnIIiQXL6wjmD8y1YIiXRWQ0OVKnzgt46uETT4grKHN+/njZ//OUcGMxW/RmFbiQWdrcyZ3kTXtCYVjE8BCoIIKfQIkie4QxmcvB5BJpcfnhdeLworiNuKawQp1Qh+/MweNu89OqkP7sJisvmdzZgZ5y/s4OmdmkI61RQEEVIoFk+kRlDrHUj3Hhnki/c9z5WfuZ/XffYBth/or+n3O5n601nMoDk58s+jRbOG+PVLBwEmFQS7esMewYwWAM5b0MHmvUc17DbFFAQRUk6NoJZDQz9/bh9XfuZ+vnDfcyye1UY27zy0eX/Nvt/J1p/O0ZqMYzbyc25rCoKgUXd17U9nh1cBT+Y3+J0HB0jFY3S1NQFBEOQdNuxSwXgqKQgiZCI1glSi9iuLv7NmOx0tSR74s9fzb//tMrrbm3hka0/Nvt/J1p/ODt+esqA1lSCXd4Ya9O5v63b0kss7M1qTkxsaOjTAvM5mYmGv9ryFHcOvL1NHQRAhE1lHUOuVxe7OI1t6eM0ruljS1YaZcdkZs3h0a0/d/LbcN5Qb3nG0oLABXaMODxWGhW64dBE7Dg7Qc3SootfZeWiABZ0tw4/nTG9m7vRm1m47WJV2SmUUBBFSzqZz2RotKNv08hF6+tJc8Yqu4WOXnzGLlw8P8cL+vpp8z5OtP50dvgdBQaPfpezX2w6xpKuN31jaDVDxiuBdh4LFZMV+85XdPLhpr3YinUIKggjZEn7QdrU3jXtOsjB9tEY9goc2B0NAVxYFwWVnzASom+Gh/vQYPYKmxu0RuDtPbj/IhYs6Wb5gOgBPVzCUk87m2Xtk6JgeAcA1586lL53j4S31U2eKGgVBhPx0/R6WdLVxRlfbuOfUevrow5v3s3hW6zH/mJd0tTFnehOPbj1Q9us9sqWHR7b0nFK/Dfalc8fVCAo9gkbceG77gQH2H01z0aIZtDcnOaO7jXUV9Aj29A7iznFBcMWZXbQ3JfjJM3uq1WQpkzadi4je/gyPbOnh5tcuOWY2y2hNieA31/0VjuGWks3leeyFA1x3wfxjjhfqBA9tDuoEpdpX7Gu/fIFP37MBCPb8v+C0Tj56zdlcsnhm1dtejv6hLPM7mo851tLANYJCfeCiRTMAOH9BR0WhX1hvUpg6WpBKxHjDObO5b+Nesrl8yU0VpTb0E4+IBzbtJZt3rjl3bsnzWlJxXndWN3f+anvVf3t9akcvR4eyXHlm13HPXX7GLPYfHWLLvonVCe5es51P37OBlefO5Ws3reCmy09n+4F+/uSuJ6e8d9Cfzg1/8BeM1AgaMwhaU3HOmjMNgPMWdrLn8CB7Dw+e4G8eqxAEo2sEEAwPHehL8/iLKhpPBQVBRKxev4fZ7U1cEO7PUsofX7WUA31pvv3Ytqq24eFwrcDlZ8467rnLzgiOTaRO8JNndnPr99bx2qVdfPHGC7jqnDl8/LeW8YV3XsDOQwN89b+2VrXd5epPZ49ZVQzFNYLGGxp64qVDvGph5/Bv6ueHUz7LLRgXtpeYN6q3BfC6s7pJJWKsXq/hoamgIIiAwUyOBzft4+plc4bnX5dy8ekzeO3SLu74r61V3eL3oS37WTZvOjPbUsc9d/qsVuZ1NPPolmB46L4NL/OhO5/ggWf3Dp/j7vzTL7byoTuf4FWndfJ/333x8FAWBGHy5uVz+dKDW9jTW95vm9UU1AjG7hE02g3sB9I5Nu4+zEWnj/wCcu786cSs/Ln/Ow8O0DWtieZk/Ljn2poS/MbSbn66fk/dTEOOEgVBBPzi+f0MZHInHBYq9sdXLWX/0er1CgbSOX697RBXvuL43gCM1Ake3rKfd3/tMd7/zTWsXr+H937jcT7w7bU8u+cw7/+XNfyvH23k9WfP5hvvvZS2puNLVP/z2nPIufPZnzxblXaXK5vLk87mj+sRjNQIJtYj6Dk6xJcf3MIt31zDO77yCG/6ws+54Y5HeGTL+D0md6d3IEP2FNo5dt2OQ2TzPlwfgGBx3dLZ7cM9Andn75HBE36A7+odYEHn8b2BgmvOncOu3kHdrGYM6Wye7Qf66e2vfMO/UmpaLDazlcAXgTjwT+7+mVHPNwHfBC4GeoB3uvuLtWxTFP10/R7amxPDwy8TccnimVxx5iy+8l9befdlp4/5W1g51mw7QDqXP2b9wGiXnzmLf39iJ+t3HeaTv72Md1xyGl//5Qv8w/2buffpPaTiMT7528u46YrF4xaUT5vZyvtfs4QvPbiFS5bMpLMlSTqXp7M1xdLZ05jX0TzhYvREDaRzPP7iAZLxGE3h/kKt4ywoW7vtIIcHn+PZ3YeZPb2Jq86Zw+VnzKIpEWPP4UE27j7Mfzy1mx+t2006l+cVs6cxqy3Fkq421u3o5cavPsrrz+7mpisWM5jOsb8vzc6DA6zf1cszO3s5GP5Db0vFmTktxeVnzOINr5zNa5Z2M21UcA5mchzoSzN3evOYPcVsLs8PntzFlx7YTDbvvHZpF79xVjeXnTGLjpZkyZ9JNpfnR0/v5h/v30zM4MKiIIBgRfCDm/byr4+9xDcefoHnXj7KynPn8unfWU53OL3Z3dnW08+MthQdLUl2HhzglfPax/2ebzxnDvGY8fVfvsBfXbecjtbSbZyoF/f38eyeI1zxillMbz7+NfvTWZ546RBP7TjEWbPbec3Srkn/exmLe7AyfSKvfXgwww+e2MkPn9zFtp7+4ckff/275/F7r15U9bZZrbphZhYHngOuBnYAjwM3uvuGonM+AJzv7n9kZjcAv+vu7yz1uitWrPA1a9bUpM2nomwuzyX/+z5ed1Y3f3fDhWX93Ue39nDDHY9yxZmz+O1Xzec3z57N3DHGZ8fi7vSnczy14xCrntzFvU8HH2xrP3H1mL/JQ/Bby6qndnH1OXOO+Ue8raePf37oRa6/eCHLF3Sc8HsfHcpy9ed/zu4xhofamxKcPbed5Qs6OHf+dM6cPY3OliSdrSmmNSVIxq1kULg7R4eyvHx4iK37jvKjp3fznxtePq4I/Nm3nc87LjntmGPnfXI1RwazxAwWz2pjd+8gA5kcrak4qUSMQ+GH+LSmBG+7aAHvufx0XjF75INvMJPjXx5+kdsf2MzhwZGeRTJunDWnnfMWdHBGdxuDmTy9Axl2HRrgl5v3c2Qw2ARvWipBa1OcpkScg/1pjoSvsWhmK7/36kW8/eKFtKYSbNl3lGd29vLVX2xly74+li+YztzpLTyyZT994ftc0tXGeQs6WNzVxrSmOK2pBO7OviND7Ds6xC+e38+OgwMsnT2Nj7zpLFYun3fMz+Kbj7zIbT9cDwRDRZcumcm3H3uJtlScj17zSvb0DnDPut1sDde+zO9o5uUjQ7z3isV84i3Lxr0+f/Hddfzbmu00JWK85fz5XL1sNvM7W5jX0cKsttS4Q6PuTjqXp7c/wwv7+3ixp48Nuw7z8+f28WJPsCFiWyrO2y5eyPUXL2R37yBrtx3kVy8c4JmdvcP7eBXOe8M5c7jwtE4Wd7WyaGYbM9tStCTjNCVix7Qhl3cyuTxD2TwD6RyHBzMcGczQczTN7t5BdvUOsOPgAFv39fHi/j4GMjkWdLawdM40Tp/ZSixmFD6CzSBmxsG+ND9+Zg8DmRznzp/OeQs6mNvRzPyOFi5ZMpMlJaaPl2Jma919xZjP1TAILgc+6e7XhI8/BuDuf1N0zurwnEfMLAHsAbq9RKMqDYK7H9/OV38xtUXIchR+ANlcnhd7+vnyuy7izefNK/l3xvIPP3ueux7fPjxjo70pQVMyRlMifswupo6TyznZvDOYyXF0KEvh30ZbKs41587l3ZeffswQQS31DmTY1tNHMh4jGTf2H03z/MtHeO7lo2zcfZgNuw+POYMnZsEU2lQiRjxmxMwwg3w+eG9D2RyDmZGhl87WJG9ePo+Vy+eSjBkH+zMMZHKsXD73uN/AN+4+zFA2z9lz2mlJxRnM5Hhkaw/3bwxmdJ0zr51z5k3n3PnTj7nN5XHvrT/D+l29dLam6JqWYkZbanizwNEyuTxrtx3k0a099A5k6B/KMZjN0dmSZPb0ZtpScX78zB4ee+EAiZiRcx/+YFk6exp/evVZrFw+FzMjnc3z65cOsnbbQZ7e0cu6HYfYNSpszWBWWxNLZ0/jfa9ZwlWvnD3mh+/BvjRf/vkWrl42hxWnz8DM2Lz3CH/2nXU8uf0QMQtqPiuXz6VvKMezew7z4v4+/mLlK0v2KiG4Uc2dv3qJHz65i6OjZr4l40ZTIk4iPvIBmss7AyHuHQYAAAkXSURBVJnc8Mr7guZkjMvPmMXrzupm6Zx2vv/rnfzHU7uGF1um4jHOX9jBpUtmcumSmbxqYSfrdvby4/CXg56+9Jjti1nw73MiH52peIz5nc0s6WpjSdc0OlqSbN1/lOdePsqOg/3BCxkYI6+ZiBvXLJvLuy5bNHzznmqYqiC4Hljp7u8PH78HeLW7f7DonGfCc3aEj7eE5+wf9Vq3ALcALFq06OJt28of9/7p+j384Mmdlb6dKWHh/yEdLUlue8uyirur7s7ze4/ywLN7efnw0PCHYX7UtU/EjHjMaErEaG9O0t6c4LSZrfzm2bOPm0451XJ554X9R3npQD+9AxkO9Wc4OpglHf52NpTJkXMn78H7j8eMRCwIle72Jma3NzOvo5kLF80glYh+qWzz3iN8/9c7SSVinDWnnbPmTOOMrmknnFyQDz9EC/damNmamtQ8/lzeeXRrD0vnTGN2+8R6n+PpT2fZsreP3b0D7Dk8yIG+NOlscH0zuTxGUJsyC4buWlMJpjcnWDQrWHQ5v7PluC3b9x8d4oFn97Kkq43lCzrG/Tfl7vT0pdnW089LB/ro7c8wkMkzmMkN/7sxgg0gU4ngT0syTntzgvbmBDPbUifsxZxskQ+CYo02NCQiUg2lgqCWvwrtBIoHWReGx8Y8Jxwa6iAoGouIyElSyyB4HFhqZkvMLAXcAKwadc4q4Kbw6+uB+0vVB0REpPpqNn3U3bNm9kFgNcH00a+7+3oz+xSwxt1XAV8DvmVmm4EDBGEhIiInUU3XEbj7vcC9o47dVvT1IPD2WrZBRERKi/50CRERmRQFgYhIg1MQiIg0OAWBiEiDq9mCsloxs31ApVtqdgGNeGPURnzfjfieoTHfdyO+Zyj/fZ/u7t1jPRG5IJgMM1sz3sq6etaI77sR3zM05vtuxPcM1X3fGhoSEWlwCgIRkQbXaEFwx1Q3YIo04vtuxPcMjfm+G/E9QxXfd0PVCERE5HiN1iMQEZFRFAQiIg2uYYLAzFaa2SYz22xmt051e6rFzE4zswfMbIOZrTezD4fHZ5rZf5rZ8+F/Z4THzcz+Pvw5rDOzi6b2HVTOzOJm9oSZ3RM+XmJmj4Xv7d/C7c8xs6bw8ebw+cVT2e7JMLNOM/uumT1rZhvN7PJ6v9Zm9j/C/7efMbM7zay5Hq+1mX3dzPaGN+wqHCv72prZTeH5z5vZTWN9r9EaIgjMLA7cDrwZWAbcaGbj30E7WrLAR9x9GXAZ8N/D93Yr8DN3Xwr8LHwMwc9gafjnFuDLJ7/JVfNhYGPR478FvuDurwAOAjeHx28GDobHvxCeF1VfBH7i7q8EXkXw/uv2WpvZAuCPgRXuvpxgS/sbqM9r/Q1g5ahjZV1bM5sJ/CXwauBS4C8L4VGSu9f9H+ByYHXR448BH5vqdtXovf4QuBrYBMwLj80DNoVffwW4sej84fOi9Ifgjnc/A94A3ENwC9n9QGL0NSe4J8bl4deJ8Dyb6vdQwXvuAF4Y3fZ6vtbAAmA7MDO8dvcA19TrtQYWA89Uem2BG4GvFB0/5rzx/jREj4CR/5kKdoTH6krYDb4QeAyY4+67w6f2AHPCr+vlZ/F3wJ8D+fDxLOCQu2fDx8Xva/g9h8/3hudHzRJgH/DP4ZDYP5lZG3V8rd19J/B/gJeA3QTXbi31f60Lyr22FV3zRgmCumdm04DvAX/i7oeLn/PgV4O6mSdsZm8B9rr72qluy0mWAC4CvuzuFwJ9jAwVAHV5rWcAbyUIwflAG8cPnzSEWl7bRgmCncBpRY8XhsfqgpklCULg2+7+/fDwy2Y2L3x+HrA3PF4PP4srgevM7EXgLoLhoS8CnWZWuOte8fsafs/h8x1Az8lscJXsAHa4+2Ph4+8SBEM9X+s3Ai+4+z53zwDfJ7j+9X6tC8q9thVd80YJgseBpeFMgxRBsWnVFLepKszMCO79vNHdP1/01CqgMGPgJoLaQeH474ezDi4Deou6npHg7h9z94XuvpjgWt7v7u8CHgCuD08b/Z4LP4vrw/Mj91uzu+8BtpvZ2eGhq4AN1PG1JhgSuszMWsP/1wvvua6vdZFyr+1q4E1mNiPsTb0pPFbaVBdHTmIR5lrgOWAL8PGpbk8V39drCLqL64Anwz/XEoyL/gx4HrgPmBmebwQzqLYATxPMxpjy9zGJ9/964J7w6zOAXwGbge8ATeHx5vDx5vD5M6a63ZN4vxcAa8Lr/QNgRr1fa+CvgGeBZ4BvAU31eK2BOwnqIBmC3t/NlVxb4H3h+98MvHci31tbTIiINLhGGRoSEZFxKAhERBqcgkBEpMEpCEREGpyCQESkwSkIpCGFu3h+IPx6vpl9t4bf6wIzu7ZWry8yWQoCaVSdwAcA3H2Xu19/gvMn4wKCtR0ipyStI5CGZGZ3Eexhs4lgsc457r7czP4A+B2CPW2WEmx4lgLeAwwB17r7ATM7k2BBTzfQD/yhuz9rZm8n2AY4R7Dh2RsJFva0ECz1/xuCHTT/AVgOJIFPuvsPw+/9uwTbIiwA/p+7/1WNfxQiJE58ikhduhVY7u4XhLu23lP03HKCXVybCT7E/8LdLzSzLwC/T7Dz6R3AH7n782b2auBLBHse3QZc4+47zazT3dNmdhvBys8PApjZXxNsffA+M+sEfmVm94Xf+9Lw+/cDj5vZj9x9TS1/ECIKApHjPeDuR4AjZtYL/Ed4/Gng/HCn1yuA7wTb3wDBtgcADwHfMLO7CTZIG8ubCDbN+7PwcTOwKPz6P929B8DMvk+whYiCQGpKQSByvKGir/NFj/ME/2ZiBPvhXzD6L7r7H4U9hN8C1prZxWO8vgFvc/dNxxwM/t7osVqN3UrNqVgsjeoI0F7JX/Tgfg8vhPWAwv1jXxV+faa7P+butxHcROa0Mb7XauBD4W6amNmFRc9dHd6ntoWgVvFQJW0UKYeCQBpSOPzyUHij8M9V8BLvAm42s6eA9QSFZ4DPmdnT4es+DDxFsGXyMjN70szeCXyaoEi8zszWh48LfkVwb4l1wPdUH5CTQbOGRE4R4ayh4aKyyMmiHoGISINTj0BEpMGpRyAi0uAUBCIiDU5BICLS4BQEIiINTkEgItLg/j/skiBHsfE67AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3wc1bn/8c+zq5VkFau7ScaSe6+yMdUGB2MgMZiEhHZpFwgQkpsCgSTcUAIXkssPEhIuCSG0hFACJDFgMGBqCGDLBvfekCzbkm31uuX8/jgjWTayJcu7Wnn2eb+8L2lnZmfO7MjfmTlz5owYY1BKKeVenmgXQCmlVGRp0CullMtp0CullMtp0CullMtp0CullMvFRbsAB8vOzjb5+fnRLoZSSh1Tli5duscYk9PeuB4X9Pn5+RQVFUW7GEopdUwRke2HGqdVN0op5XIa9Eop5XIa9Eop5XI9ro5eKaUA/H4/JSUlNDY2RrsoPUpiYiJ5eXn4fL5Of0aDXinVI5WUlJCamkp+fj4iEu3i9AjGGPbu3UtJSQkFBQWd/lyHVTci8riIlInIqkOMHykiH4tIk4jcdNC4OSKyXkQ2icitnS6VUirmNTY2kpWVpSHfhoiQlZV1xGc5namjfxKYc5jx+4DvAfcfVCAv8DBwFjAauEhERh9R6ZRSMU1D/su68p10GPTGmA+wYX6o8WXGmCWA/6BR04BNxpgtxphm4Dng3CMuYWcFA/DmbVBZHLFFKKXUsSiSrW5ygbapW+IM+xIRuVZEikSkqLy8vGtLq9wOS5+GP8+Duj1dm4dSSoXR/Pnzue+++wC44447uP/++zv4RGT0iOaVxphHjTGFxpjCnJx27+DtWNYQuPg5qCqGZ74BTTXhLaRSSh2huXPncuut0b88Gcmg3wEMbPM+zxkWOYNOhAuegp0r4LlLINAU0cUppdztL3/5C9OmTWPixIl8+9vfJhgMkpKSwg9+8APGjBnDrFmzaKmFeOihhxg9ejTjx4/nwgsvBODJJ5/kxhtv/NJ8P//8c6ZPn8748eOZN28eFRUVAMycOZNbbrmFadOmMXz4cD788MOwrEckm1cuAYaJSAE24C8ELo7g8qwRc+Dch+Ef18EH/wun3xbxRSqlIuvOV1azprQ6rPMcPaA3t39tzCHHr127lueff56PPvoIn8/HDTfcwDPPPENdXR2FhYU8+OCD3HXXXdx555387ne/47777mPr1q0kJCRQWVl52GVfdtll/Pa3v2XGjBn8/Oc/58477+TXv/41AIFAgMWLF7NgwQLuvPNO3n777aNe1w6DXkSeBWYC2SJSAtwO+ACMMb8XkX5AEdAbCInI94HRxphqEbkRWAh4gceNMauPusSdMfEiWPsKFD0Op9wEvsRuWaxSyj0WLVrE0qVLmTp1KgANDQ306dMHj8fDt771LQAuvfRSzj//fADGjx/PJZdcwnnnncd55513yPlWVVVRWVnJjBkzALj88su54IILWse3zG/KlCls27YtLOvSYdAbYy7qYPwubLVMe+MWAAu6VrSjNO0aWP8arH4ZJkb+REIpFTmHO/KOFGMMl19+Offee+8Bw3/xi18c8L6lueNrr73GBx98wCuvvMI999zDypUru7TchIQEALxeL4FAoEvzOFiPuBgbEYNnQvYI+PQPYEy0S6OUOsbMmjWLF198kbKyMgD27dvH9u3bCYVCvPjiiwD89a9/5eSTTyYUClFcXMxpp53GL3/5S6qqqqitrW13vmlpaWRkZLTWv//5z39uPbqPFPd2gSBij+oX3AQlRTBwarRLpJQ6howePZq7776b2bNnEwqF8Pl8PPzwwyQnJ7N48WLuvvtu+vTpw/PPP08wGOTSSy+lqqoKYwzf+973SE9PP+S8n3rqKa677jrq6+sZPHgwTzzxRETXRUwPO9otLCw0YXvwSFMtPDAKhp8JX38sPPNUSnWLtWvXMmrUqGgX40tSUlIOebTeXdr7bkRkqTGmsL3p3Vt1A5CQAhMvgdX/gJrd0S6NUkpFhbuDHmz1TcgPy56KdkmUUi4Q7aP5rnB/0GcNgfxTYOXf9KKsUiomuT/oAcZ+HfZsgN3t9rSslFKuFhtBP2oueOJg5YvRLolSSnW72Aj65CwYfBqselmrb5RSMSc2gh5s9U3VF7ZNvVJKdcHVV1/NmjVrol2MI+beG6YONvJs8CbAqpf05imlVJc89tixeT9O7BzRJ6bBsDNg9d8hFIx2aZRSPVxdXR3nnHMOEyZMYOzYsTz//PPMnDmTlhs6//SnPzF8+HCmTZvGNddc09od8RVXXMH111/P9OnTGTx4MO+99x5XXXUVo0aN4oorrmid//XXX09hYSFjxozh9ttvj+i6xM4RPdjqm3WvwvZ/Q8Ep0S6NUqqzXr8VdnWtk7BD6jcOzrrvkKPfeOMNBgwYwGuvvQbYXicfeeQRAEpLS/nFL37BsmXLSE1N5fTTT2fChAmtn62oqODjjz9m/vz5zJ07l48++ojHHnuMqVOn8vnnnzNx4kTuueceMjMzCQaDzJo1ixUrVjB+/PjwrqMjdo7oAYbPAV+yrb5RSqnDGDduHG+99Ra33HILH374IWlpaa3jFi9ezIwZM8jMzMTn8x3QzTDA1772NUSEcePG0bdvX8aNG4fH42HMmDGtXQ+/8MILTJ48mUmTJrF69eqI1v3H1hF9fBIM+wqsfx3OeQA8sbWfU+qYdZgj70gZPnw4y5YtY8GCBdx2223MmjWr059t6WrY4/G0/t7yPhAIsHXrVu6//36WLFlCRkYGV1xxBY2NjWFfh9blRmzOPdWIc6B2F5R+Fu2SKKV6sNLSUpKSkrj00ku5+eabWbZsWeu4qVOn8v7771NRUUEgEOCll46slqC6uprk5GTS0tLYvXs3r7/+eriLf4DYOqIHe0FWvLB+AeRNiXZplFI91MqVK7n55pvxeDz4fD4eeeQRbrrpJgByc3P56U9/yrRp08jMzGTkyJEHVO10ZMKECUyaNImRI0cycOBATjrppEitBuD2booP5cmvQv1euOHjyC5HKdVlPbWb4ha1tbWkpKQQCASYN28eV111FfPmzeuWZWs3xZ0x4iwoWwP7tka7JEqpY9Qdd9zBxIkTGTt2LAUFBYd9Tmy0xV7VDdigX/hTe1H2hBuiXRql1DHo/vvvj3YROi02j+gzB0POKFtPr5TqsXpa1XJP0JXvJDaDHuxR/fZ/Q/2+aJdEKdWOxMRE9u7dq2HfhjGGvXv3kpiYeESfi82qG4CR58C/HoBNb8P4b0a7NEqpg+Tl5VFSUkJ5eXm0i9KjJCYmkpeXd0Sfid2gHzAZUvrCutc06JXqgXw+HwUFBdEuhivEbtWNxwPDZsPmdyDoj3ZplFIqYmI36MEGfVM1FH8a7ZIopVTExHbQDzkNPD7YsDDaJVFKqYiJ7aBPSIVBJ8LGt6JdEqWUipgOg15EHheRMhFZdYjxIiIPicgmEVkhIpPbjAuKyOfOa344Cx42w2ZD+Vqo/CLaJVFKqYjozBH9k8Ccw4w/CxjmvK4FHmkzrsEYM9F5ze1yKSNp+Jn2p1bfKKVcqsOgN8Z8ABzurqJzgaeN9QmQLiL9w1XAiMsaChkFWn2jlHKtcNTR5wLFbd6XOMMAEkWkSEQ+EZGe2eOPiK2+2foB+BuiU4byDXqHrlIqYiJ9MXaQ023mxcCvRWRIexOJyLXODqEoKnfBDZ8NgQbY+mH3L7uxGv5wKjx8vFYfKaUiIhxBvwMY2OZ9njMMY0zLzy3Ae8Ck9mZgjHnUGFNojCnMyckJQ5GO0KCTwZcEG9/s/mVvfNPuZLzx8Ndvwivfh+a67i+HUsq1whH084HLnNY304EqY8xOEckQkQQAEckGTgIi9/Tbo+FLhMEzYcMb0N0dKK17FZL7wHeL4MTvwtInYdFd3VsGpZSrddjXjYg8C8wEskWkBLgd8AEYY34PLADOBjYB9cCVzkdHAX8QkRB2h3KfMaZnBj3AiLNtt8W7VkD/Cd2zTH+jvQg87hvg6wWz74ady6F4cfcsXykVEzoMemPMRR2MN8B32hn+b2Bc14vWzUacBeKxnZx1V9BveQ+aa2Hk1/YP6zsWip6AUBA83u4ph1LK1WL7zti2krNh4HRY140PI1n3CiT0hoJT9w/rO9bW2e/b0n3lUEq5mgZ9WyPPgd0roWJb5JcVDNhHGQ6bDXHx+4f3G2t/7loZ+TIopWKCBn1bI8+2P7vjqL74E6jfC6O+euDw7BEgXti9OvJlUErFBA36tjIHQ58xtp4+0ta+Ct4EGHrGgcN9iZA9HHa327WQUkodMQ36g408G774N9Ttjdwymmph7Ssw5HRISPny+L5jYJcGvVIqPDToDzbyHDAh26Y+EvZthT/NhppSKLyy/Wn6jYXqEmioiEwZlFIxRYP+YP0nQu9ceyNTuG1+B/54mg3xS/62v+fMg/V1LshqPb1SKgxi9+HghyJib2D66CEo/QwGtNtrQ8eMgcZKqNhuzw7WvmLr3XNGwYXPQFa73f5YbYM+/+SuLV8ppRwa9O055Ufw+bPw6g/h6rf337i08kVY9hQMOxNGnwvpAyHQDPs2Q9kaW6++exWUrYPaXRBsdmYocNwJcOa9MPmy9uvl20rtB70ytYmlUiosNOjbk5gGZ94DL19jg73wKtsS5+VrIbG37dL4zZ9B2nG2rj0UsJ/zxNnmkccdD70HQEpf+yo4FVL6dH75IraeXqtulFJhoEF/KOMugGVPw9t32OD/+3UwYCJc9k+oLYO1822/NJkXQM5I5zUC4hLCs3ztCkEpFSYa9IciAuf8P3jkRHjxKlu3fsmL9oHiCalw8g8iu/y2XSFkD4vsspRSrqatbg4nZwSc/t+2k7P/+DskZXbfsrUrBKVUmGjQd+Tk78O3P4De3fwYXO0KQSkVJhr0PVVLVwi7VkS7JEqpY5wGfU+WNwVKirr/qVdKKVfRoO/J8qZCwz7tm14pdVQ06HuyvKn2Z8mS6JZDKXVM06DvyXJGQnyKBr1S6qhoO/qezOOF3MnhC3pjbM+coSDg/G5CzjUAY3/6ksCrfxZKuYn+j+7p8qbBvx6E5nqIT/ry+KYa+OITKFsL5etsJ2qNVfbVVA1Bv+1zJ+Tv/DLjEqFXBsy5D8acF751UUpFhQZ9T5c3FUwQdn4Og07cP7yxCj59FD7+ne0lE2y/OpmDbWdriWNttU9cAnh94PHZMwTx2Lt+xWNfyP73xoC/AZprYPN78I/rbfVRn5HRWHOlVJho0Pd0eYX2Z/Hi/UG/8kV47Yc27IfPgeOvs3fvhvPO3eqd8IdT4IXL4Jp3Ou5xUynVY+nF2J4uORsyCvbX09eWwSvfh6yhcO37cPHzMOS08HfP0Ls/fP1PsHcjvPr98LTlDzTZ8vckW96DJX/SexUiyd8Au9fodxxFekR/LMibClvft/9R3r3HdnY271HIHhrZ5Q6eAaf9FN6523bBfPIPIWd41+a1YSEsuBkqt8PA42H8N+3ZiMdnq6bik20vod2ptsyesTRWQfGnMPe34et9VFn+RnjmAtj2oe2or/BKGPdN29236jYa9MeCvKmw8gXY+KbtOnnatyMf8i1O/hE0VMKSx2D5szD8LHsdoHoHVJdCsAm88fYVaIT6vVBfAb5ekDvF3t27Y5l9NGP2CJhxq+3i+bUf2VcLjw9m/Tec8F3wdNOJ5pu32Yvcx18Hn/4eKovt078i2Xmdv9FeK/H6jm4+zfW2e4xdK+2rtgzSciH9OPsozKQs6JVud57elus0cfZ6zMGMsRfuGyrsTi9zMKQNbH/aIxEKwt+vtSE//Qb787Ufweu3Qt8xttvvvmMhOadNWePt34I3Dnv9yLmmREtZzP5WYmCHe+Lsd9orY/+OOhS0Z8HrF0BjNWTkQ2YBZA6xT3fz9Tq6dTvGiOlhp1OFhYWmqKgo2sXoWUo/g0dnOke8At/7rHt70gSo22PDfvEfobnOhkrvAbY5ZrDZPmkrLt4JmEx7gXjHUti7yU4z48cw/Tt2GmPsk7iKP7Xr4/HCxrfszmDwaTDv9/YpW5G05X14ei6c+mM4/Wf2usc/brBlaXloTHK2Lbuvl22JRNvgc5qngh3nS7JnJb0HQMYgSO0PNbvsXc17NkLZalt9sW+L/Wximv2uEtMgobft+jou0QmtOLuzE4/t2M7XywnCdBvGWz+w313LE8wS0224V+/Yf2H+aKUOgOOmw7RrDmwE0FnG2DO4JX+E2ffAiTfaYTuWwdp/QunntoFBY1V4yguA2O89/Tj7d1e/x+40ElLs99Z2urSBdrpE57v3Je3fYYj30DvE9polt8yz5bOJafb/Z1KW/RtKzrG/1++zT6Pbt8U+rCghFRLSIDkLUvo5T5bL6PIOVkSWGmMK2x2nQX8MCPrh3jx7xDznlzD9uuiVpeXvpbN/jA0VgNig6mi+y56yR3txCbZaZ8hp9ulcqf2P/uiyrUCTfc5AKAg3fLz/6G7HUljxAtTuhprd9uzE3wD+evuZ/YXdf6RpsNsl2NTekhxij5L7joE+o+3n6vfanWdTtW0i21ht5xEK2hBoud/BhOzy/fX759VvnP1e8k+2v/fO3f/9NFbZM62Giv1H6MFm+zfU8iS0lu+77Xea0NuGTHwylK+H4k/sDqVuj3205sxbO38WUrcXFt1hzz5P/C7Mvrv96YyBmp1OWSud5sBOWYN+Wo/eW3aorV+nc4Qv4nxPAfuqLYfKL+wrtS+MOBuGnWGDt6ESKrbZoN2zyV57qtphW5g11dgzpFDA+c6Dba4nmP3LggPPMA4uR8s2a67t3PfUngGT4dp3u/RRDXo3eGqu/U9x/b+P/rS/JytfD+//Eja/a/v5AXtUlpxjj47ik50mowltmot6OjgCc8JCPPYzteXwxb/h0pdh6KzwlDsYsKFRtcMGTU2pPUrLGmIvpvsSj27+gWZ7tO712UDuDk218MYt8NlfILfQhnbLg3e88fun88TtP+tZOx/eu9d+9oQb4Ct3dV9VXE8R9Nuj9/o9zg693L5PTLc7/MwC+zfcsoOv32P/b9fstmcfky/r0mKPKuhF5HHgq0CZMWZsO+MF+A1wNlAPXGGMWeaMuxy4zZn0bmPMUx0VVoP+EOqd0OvuKptoCYVsHfQXn9j/BHXl9uVvcKqKGu00JmSPwFocfKTa9n4BzP6jv1Fz4Su3d/daHZtWvWxbejV1sppl8Gkw517oMyqy5VIHONqgPxWoBZ4+RNCfDXwXG/THA78xxhwvIplAEVCIPf9ZCkwxxlQcPI+2NOiV6oEaq+yZSnOdPVpvuT4A9q5rf6NtDZY+CAbPDG9Vm+qUwwV9h61ujDEfiEj+YSY5F7sTMMAnIpIuIv2BmcBbxph9TiHeAuYAzx5Z8ZVSUZeYZq8HqGNSOCrPcoHiNu9LnGGHGv4lInKtiBSJSFF5eXkYiqSUUqpFj7hKYox51BhTaIwpzMnJiXZxlFLKVcIR9DuAgW3e5znDDjVcKaVUNwpH0M8HLhNrOlBljNkJLARmi0iGiGQAs51hSimlulGHF2NF5FnshdVsESkBbgd8AMaY3wMLsC1uNmGbV17pjNsnIr8AWp6acVfLhVmllFLdpzOtbi7qYLwBvnOIcY8Dj3etaEoppcKhR1yMVUopFTka9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XIa9Eop5XKdCnoRmSMi60Vkk4jc2s74QSKySERWiMh7IpLXZlxQRD53XvPDWXillFIdi+toAhHxAg8DZwAlwBIRmW+MWdNmsvuBp40xT4nI6cC9wH844xqMMRPDXG6llFKd1Jkj+mnAJmPMFmNMM/AccO5B04wG3nF+f7ed8UoppaKkM0GfCxS3eV/iDGtrOXC+8/s8IFVEspz3iSJSJCKfiMh57S1ARK51pikqLy8/guIrpZTqSLguxt4EzBCRz4AZwA4g6IwbZIwpBC4Gfi0iQw7+sDHmUWNMoTGmMCcnJ0xFUkopBZ2oo8eG9sA27/OcYa2MMaU4R/QikgJ83RhT6Yzb4fzcIiLvAZOAzUddcqWUUp3SmSP6JcAwESkQkXjgQuCA1jMiki0iLfP6CfC4MzxDRBJapgFOAtpexFVKKRVhHQa9MSYA3AgsBNYCLxhjVovIXSIy15lsJrBeRDYAfYF7nOGjgCIRWY69SHvfQa11lFJKRZgYY6JdhgMUFhaaoqKiaBdDKaWOKSKy1Lke+iV6Z6xSSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrmcBr1SSrlcp4JeROaIyHoR2SQit7YzfpCILBKRFSLynojktRl3uYhsdF6Xh7PwSimlOtZh0IuIF3gYOAsYDVwkIqMPmux+4GljzHjgLuBe57OZwO3A8cA04HYRyQhf8ZVSSnWkM0f004BNxpgtxphm4Dng3IOmGQ284/z+bpvxZwJvGWP2GWMqgLeAOUdf7Pa9tLSEyvrmSM1eKaWOSZ0J+lyguM37EmdYW8uB853f5wGpIpLVyc8iIteKSJGIFJWXl3e27AfYuqeOH/1tOSfe9w53vrKaHZUNXZqPUkq5Tbguxt4EzBCRz4AZwA4g2NkPG2MeNcYUGmMKc3JyulSAguxkXv+vU5gzph9//ng7p/7qXV5eVtKleSmllJt0Juh3AAPbvM9zhrUyxpQaY843xkwCfuYMq+zMZ8NpVP/ePPCtibz/49OYkJfGXa+u0aocpVTM60zQLwGGiUiBiMQDFwLz204gItki0jKvnwCPO78vBGaLSIZzEXa2MyyictN7cc+8cVQ3+HnwrQ2RXpxSSvVoHQa9MSYA3IgN6LXAC8aY1SJyl4jMdSabCawXkQ1AX+Ae57P7gF9gdxZLgLucYRE3qn9vLjl+EH/59AvW76rpjkUqpVSPJMaYaJfhAIWFhaaoqCgs86qoa2bm/e8xZkBvnrn6eEQkLPNVSqmeRkSWGmMK2xvn6jtjM5Lj+dHs4fx7814Wrt4d7eIopVRUuDroAS6edhwF2ck89uGWaBdFKaWiwvVBH+f1cNG0gRRtr2BTmdbVK6V6nmDIUFnfTFl1Y0TmHxeRufYw50/O41dvrOf5JcX87JyDe29QSqkjZ4zBHzQ0NAepbQ5Q2xigtslPkz9EUzCEPxDCtJm2tLKRzeW1bC6vpbLeT3MwRJM/RG1TgOpGP8bAlEEZvHT9iWEva0wEfXZKAl8Z1ZeXlu3g5jNHEh/n+hMZpXqMT7fs5Z/LSxmXm8b0wVnkZyWFvWFES6OSlvlu3VPH66t2snDVLspqmkhJiCM1MY7khDh6+bz0ivcS5/EQCIUIBA3+YMh5GZqDIQLBEIGQDfJAMEQwZPCHQjQHnGkCIRr8QYKhI2vM0jsxjiF9UjguM4kEn5d4r4fUxDjSevlIT/KRl5EU1u+lRUwEPcC3pg3kjdW7eHvtbs4e1z/axVEqJny6ZS+XP7GYQNDw10+/ACA9yUdWcjzpSfGkJsbREvkGW4URMoZA0P4MhgxBA6GQ2T8u5ARzIERjIERDc5AGfxCPQEKcF59XqG4MADBhYDonDc2mrilArfMqr2mi0R/EHzTEeYU4jxDn8RAf58HnFXxeD0nxcc44O8zrscPjvXa6+DhP6w4j0eclJcFLSoKP5AT73uf1kBDnoe3+rE9qItkp8VFp/RczQX/qsBz6pyXy/JJiDXql2rGntok1pdUkJ3gZmJlETkrCUYXS58WVXPXkEnLTe/H8t0+gst7Pp1v3srq0mqp6P5UNzeyrO/DOdY/Y4PU44ZsQZ3/3Cng9gogQ77XhG+dtE7ZxHgzQFAjR6A8yKCuZOWP7kZve6yi/FXeImaD3eoQLCgfy23c2sqOyQf8AlGvsq2tmR0UDvriWEPTsr4YIhKhvDtDgD9q65KYAdU0B6pqDVDX4qay3YbumtJrSqgMvBPbyeZk3OZefnDWS1ETfEZVp6fZ9XPnEErJSEnjm6ulkpySQnZLA0D4p4Vx11UkxE/QAF0zJ47fvbORvRcV8/yvDo12cdn3rDx8TDBn++6ujmTAwPdrFUT1MKGRYt6uGjzbtYcm2fawure5yT60JcR4ykuJJT/IxJT+TK3PTGJPbmyZ/iC/21bNqRxXPLf6C99eX88uvj+fkYdkdzrN4Xz2/WrieV5aXkpvei2euPp5+aYldKp8KH1ffGdueCx/9mLKaJhb9cEaPu1O2qsHPhDvfRASMsTumW88aSVZKQrSL5lqP/2srG3bXcNe5Y7v9In1zINR6VF1R72dfnT26rqhvpqKumcoGP5X1fuqbA/aCYShE8b569tTa6o78rCTG5aUzdkBv8rOTCYUMTQF7UdHWN3uc+mZbvZEU7yU53l6QTHLqljuydHsFN/9tOVv21DGyXyrZKQlkJseT6PMQMhByWp7UO/Xfn31RiccD15wymGtPHXzEZwKq6w53Z2xMHdEDnDOuP//9z9VsLKtleN/UaBfnAGtKqwH43UWTWVFSyeMfbaWmMcDv/2NKlEvmTsuLK7n7tTWEDNQ0BvjNhROJ84Y37I0xbNlTx9LtFSzbXsHnxZXsqW2iujFAcyB0yM/18nnJSPKRlhRPcrxzcc8Xx6nDcjhxaDYnDc2if1rkqx+nDMpgwX+dwu/f38yqHVXsrWumuKKe5kAIwbZyiY/zkJzgJSk+jm9NHch3ThuqR/E9TMwF/Zlj+vHz+atZsHJnjwv61aVVAEwryOSc8f2pavCzYOVOQiGDx9Ozzj66YtHa3by1Zjf52V0uOSQAAA8ySURBVMkM75tCRlI8W8rr2FhWS6M/yHUzhnRbQDQHQtzy0gpyUhO4eNogHnx7A4k+L//7jfFH9F0bY6io91Na2cCuqkaqGvzUNNqj85U7qvisuJLKej8Aab18TDounSmDMkhJjCMlPo70JB/pTvVJRlI8WSnxZCTFd+pou7sk+rw9tqpTdU7MBX2f3olMHZTJ6yt39bg/3jWl1fRJTSAn1VbVTM3P5LklxWwoq2Fkv95dmmcwZGjwB2n0BwkETWu74YDTXC3gtA1uaa1Q1xSktslPcyDEWeP6kx2GaiNjDA8t2sSDb28gKd5LffOBz6TxeQVBeHlZCXfMHcO8SbkRr1b7v/c2sW5XDY9dVshXRvcF4MG3N1BZ38yY3DQyk3wkxcfRFAzR5FzIrGkKUN3gp6rBT1lNE+U1TZTVNNLob//IfFifFM4c3a813IfkpLhih62OPTEX9ABnjevHna+sYVNZbY9qBbBmZzVjBuwP9Kn5mQAs2VZx2KA3xrBhdy0fbChneUklu6oa2VnVyJ7aJpoOUz3QkSf/vY1nr5lOn95HdpTdHAhRUd9MIGRvNvnVG+t5beVOzp+cy//MG0eTP8Sm8hr21jYzOCeF/KwkSioauOlvy/nhC8uZv7yUs8f1Z+LAdIbkpOANYzj6gyFW7aji4Xc3MXfCgNaQ/96soYSM4amPt7FoXVm7n030eUhN9JHWy0dOSgKTjksnJyWBAem9GJDei/5piaQn+UhN9JGaGIcvzNVASnVVzF2MBdhZ1cAJ977DTbOHc+PpwyK6rM5q9AcZc/tCrp8xhJvOHAHYAJ9+7yKOL8jioYsmfekzwZDhjx9u4cmPtrHL6SNjYGYvctN70T+tF9kp8STFx7VejIvzeFpvEPE67ZS9HiHBZ2/uSIjzkpoYR0pCHFv31HHN00X0T0vsMOz31DaxsqSKpdsrWLxtH8uLKw/YwXgEfnLWKK4+peCwR+rBkOGJj7by0KKNrTe8JPpsy5CWuxp9HnsTStvwNwYMxvlpW6aEjL3RpjngHJH7g9Q02guGAJnJ8bz1g1PbvdAdDBmqGvzUNQVI8HlI9Hlbvx+leiq9GHuQ/mm9mHxcOq+v2tVjgn7D7hqCIXPAEb2IUJifSdG2Lz+rpaSinh8+v5zF2/YxY3gOPzhjGKcMy2FAmO4PGJDeiyevnMYVTyzmoj9+wq++MZ5xuenEx3loDoT4YEM5r6woZcnWfa3tr70eYcyA3lw6fRAF2cn4vIJHhBH9Uhmf13FTUa9HuPqUwVx1UgFb9tSxvLiSNTurqWrwO/2IBJw7JW378Lb7DEFw/hEf5zng5ppEJ6xTE+NamxOePDT7kK2ZvB4hMzmezOT4cHyVSkVdTAY9wNnj+nP3a2vZvreOQVnJ0S4Oq50WN2MGpB0wfFp+Jq+t2ElJRX1rPxjvrivje89+hgEe+OaEiNVpTyvI5IkrpnLlk0v4+iMfkxDnYWxuWmunTOlJPk4ZlsNVeWmMy01jbG4ayQlH/yfl8QhD+6QwtE8KXw/DeigV62I26OeM7cfdr63l9VW7uG7GkGgXh9WlVaQmxJGXceAReWF+BgBF2yrIy0ii0R/k1pdXMCC9F49dXsjAzMh0gtTi+MFZ/OuW01m8dS9F2yr4rLiSU4flcN6kAZw8NEc7iFPqGBCzQZ+XkcT4vLQeE/RrSqsZNaD3l1pljOzXm9SEOBZv28d5k3J5cWkJu6ubeOCbEyMe8i0yk+OZM7Y/c8ZqH0FKHYti+nBs9ui+LC+upKwmMp39d1YwZFi7s+aA+vkWXo8wJT+Dom378AdDPPLeZiYdl86JQ7KiUFKl1LEopoN+1ijbtO6dte03p+suW/fU0eAPfql+vsXU/Ew27K7lyY+2saOyge+ePrTHdd+glOq5YjroR/ZLJTe9F2+vje6Dw1vuiG3viB72t6f/34XrGd2/N6eN6NNtZVNKHftiOuhFhDNG9+XDjXtoOOhuze60prSaeK/nkDdvjc9LI97roTkY0qN5pdQRi+mgB/jKqL40BUJ8tGlP1MqwZmc1w/ulHPJOykSfl8L8DEb0TeXMMf26uXRKqWNdzLa6aTGtIJPUhDjeXru79Xb47hQKGVbtqGL26MMH+MMXT8aA9pWilDpiMR/08XEeZozI4e21ZVHpJXJ1aTUV9X6mFmQedroMvUtTKdVFMV91A3DG6L7sqW1ieUllty970brdiMBpI3K6fdlKqdigQQ/MHN4Hr0ei0vrmnXVlTBqYrk+RUkpFjAY9kJbkY1p+Jm+u7t6gL6tuZEVJVWt7fqWUioROBb2IzBGR9SKySURubWf8cSLyroh8JiIrRORsZ3i+iDSIyOfO6/fhXoFwOXNMXzaW1bK5vLbblvnuenuj1ukjtV28UipyOgx6EfECDwNnAaOBi0Rk9EGT3Qa8YIyZBFwI/F+bcZuNMROd13VhKnfYnTnWtnp5Y9WublvmorVl5Kb3YmS/nvVIQ6WUu3TmiH4asMkYs8UY0ww8B5x70DQGaLmtMw0oDV8Ru0f/tF5MHJjebUHf6A/y4cY9nD6yj94ApZSKqM4EfS5Q3OZ9iTOsrTuAS0WkBFgAfLfNuAKnSud9ETmlvQWIyLUiUiQiReXl5Z0vfZidNbYfK3dUUbyvPuLL+mTLXhr8QU4fpdU2SqnICtfF2IuAJ40xecDZwJ9FxAPsBI5zqnR+CPxVRL7UoYsx5lFjTKExpjAnJ3rNDOc41TcLV0f+qP6ddWX08nk5YbD2QqmUiqzOBP0OYGCb93nOsLb+E3gBwBjzMZAIZBtjmowxe53hS4HNwPCjLXSkDMpKZlT/3hGvvmkOhFi0toyThmaT6NPnkCqlIqszQb8EGCYiBSISj73YOv+gab4AZgGIyChs0JeLSI5zMRcRGQwMA7aEq/CRcNbYfiz9ooKy6sj0UV/fHODqp4vYUdnAN6YcXAOmlFLh12HQG2MCwI3AQmAttnXNahG5S0TmOpP9CLhGRJYDzwJXGGMMcCqwQkQ+B14ErjPGfPlJ1z3InLH9MAYWrgl/m/rK+mYufexT/rWxnPvOH6dPbFJKdQuxedxzFBYWmqKioqgt3xjDrAfeJyclgee/fUJY5tnoD7Jw9S5++84mvthbz0MXTdSQV0qFlYgsNcYUtjcu5js1O5iIcOHUgfzPgnW8vebLPVoGgiHiDtGdMIA/GGL73jqKKxrYUdHAht01zF9eSmW9n+Myk3jyyqmcODQ70quhlFKtNOjbccWJBby4tISf/3MVJwzJIjkhDmMMt89fzXOLizllWDZzxvZjan4mJRUNbCyrYcPuGlaXVrNuVw3NgVDrvOLjPJwxui8XTzuOEwZnaTfDSqlup0Hfjvg4D/eeP46vP/IxD7y1gdvOGcX/LFjL0x9vZ8bwHNbtqmHRugOfM5ue5GN0/95cfsIgRvXvzaCsJHLTk8hJTcCr4a6UiiIN+kOYMiiTS44/jic+2kplvZ+XlpVw+QmDuGPuGABW7qhidWk1g7KSGNYnleyUeL3DVSnVI2nQH8aP54zkzTW7eWlZCRdMyeP2r41pDfPxeemMz0uPcgmVUqpjGvSHkdbLx/9dMpmPN+/lO6cN1fp1pdQxSYO+A1PzM5maf/jH/CmlVE+mDx5RSimX06BXSimX06BXSimX06BXSimX06BXSimX06BXSimX06BXSimX06BXSimX63H90YtIObD9KGaRDewJU3GOFbG4zhCb6x2L6wyxud5Hus6DjDHtPnS7xwX90RKRokN1vu9WsbjOEJvrHYvrDLG53uFcZ626UUopl9OgV0opl3Nj0D8a7QJEQSyuM8TmesfiOkNsrnfY1tl1dfRKKaUO5MYjeqWUUm1o0CullMu5JuhFZI6IrBeRTSJya7TLEy4iMlBE3hWRNSKyWkT+yxmeKSJvichG52eGM1xE5CHne1ghIpOjuwZHR0S8IvKZiLzqvC8QkU+d9XteROKd4QnO+03O+PxolrurRCRdRF4UkXUislZEToiFbS0iP3D+vleJyLMikujGbS0ij4tImYisajPsiLeviFzuTL9RRC7vaLmuCHoR8QIPA2cBo4GLRGR0dEsVNgHgR8aY0cB04DvOut0KLDLGDAMWOe/BfgfDnNe1wCPdX+Sw+i9gbZv3vwQeNMYMBSqA/3SG/ydQ4Qx/0JnuWPQb4A1jzEhgAnbdXb2tRSQX+B5QaIwZC3iBC3Hntn4SmHPQsCPaviKSCdwOHA9MA25v2TkckjHmmH8BJwAL27z/CfCTaJcrQuv6T+AMYD3Q3xnWH1jv/P4H4KI207dOd6y9gDznD/904FVAsHcKxh283YGFwAnO73HOdBLtdTjC9U0Dth5cbrdvayAXKAYynW33KnCmW7c1kA+s6ur2BS4C/tBm+AHTtfdyxRE9+/9QWpQ4w1zFOUWdBHwK9DXG7HRG7QL6Or+76bv4NfBjIOS8zwIqjTEB533bdWtdb2d8lTP9saQAKAeecKqrHhORZFy+rY0xO4D7gS+AndhttxR3b+u2jnT7HvF2d0vQu56IpAAvAd83xlS3HWfsbt1V7WRF5KtAmTFmabTL0o3igMnAI8aYSUAd+0/jAddu6wzgXOyObgCQzJerN2JCpLavW4J+BzCwzfs8Z5griIgPG/LPGGNedgbvFpH+zvj+QJkz3C3fxUnAXBHZBjyHrb75DZAuInHONG3XrXW9nfFpwN7uLHAYlAAlxphPnfcvYoPf7dv6K8BWY0y5McYPvIzd/m7e1m0d6fY94u3ulqBfAgxzrtLHYy/kzI9ymcJCRAT4E7DWGPNAm1HzgZar7Zdj6+5bhl/mXLGfDlS1OS08ZhhjfmKMyTPG5GO35zvGmEuAd4FvOJMdvN4t38c3nOmPqSNfY8wuoFhERjiDZgFrcPm2xlbZTBeRJOfvvWW9XbutD3Kk23chMFtEMpyzodnOsEOL9oWJMF7gOBvYAGwGfhbt8oRxvU7GnsqtAD53Xmdj6yQXARuBt4FMZ3rBtkDaDKzEtmSI+noc5XcwE3jV+X0wsBjYBPwNSHCGJzrvNznjB0e73F1c14lAkbO9/wFkxMK2Bu4E1gGrgD8DCW7c1sCz2OsQfuwZ3H92ZfsCVznrvwm4sqPlahcISinlcm6pulFKKXUIGvRKKeVyGvRKKeVyGvRKKeVyGvRKKeVyGvRKKeVyGvRKKeVy/x8lr3mc8AkT/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(time, cost) #, label = 'cost')\n",
    "plt.xlabel('timestep')\n",
    "plt.ylabel('loss')\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(time, epsilon, label = 'epsilon')\n",
    "plt.plot(time, sigma, label = 'sigma')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.linspace(0.2, 4, 1000)\n",
    "start_pot = hoomd.htf.compute_pairwise_potential(model_dir, r, 'energy', checkpoint=0)\n",
    "end_pot = hoomd.htf.compute_pairwise_potential(model_dir, r, 'energy', checkpoint=-1)\n",
    "true_pot = 2 * (r**-12 - r**-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFpCAYAAABeVxsLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzdd5ycVdn/8c+Zvr23JJtsGgTSQ0CKKBB6r4oiCogRENDHRx+xwc/6oMaGCghIUUABkR4C4SGAmBBSgDRCetmS7GZ7m90p5/fHbGLA7D272ZlNZvN9v1772inXzH1uZtgr1znnPsdYaxEREZGDm+tAN0BERETiU8IWERFJAUrYIiIiKUAJW0REJAUoYYuIiKQAJWwREZEUMOCEbYwpN8YsMMasMcasNsZ8dR8xxhhzhzFmgzFmhTFmxkCPKyIicijxJOA9wsB/W2uXG2OygGXGmPnW2jV7xZwFjO/5+RhwV89vERER6YMBV9jW2hpr7fKe263A+8Dwj4RdAPzZxrwF5BpjygZ6bBERkUNFQsewjTEVwHRg8UeeGg5s3+t+Jf+Z1EVERKQXiegSB8AYkwk8CXzNWtuyn+8xG5gNkJGRcdSECRMS1TzZD5t3rcFlXIwq6P1z2NXWRU1zkInDsnEZM4itExEZepYtW7bLWlu0r+cSkrCNMV5iyfoRa+0/9hFSBZTvdX9Ez2MfYq29B7gHYObMmXbp0qWJaJ7sp/n3HYfbn8MpV87rNeahhVu47dnVvPL908jP8A1i60REhh5jzNbenkvELHED/Al431r7q17CngU+3zNb/Fig2VpbM9BjS3KdZrI4xfodY7zu2FcoFIkORpNERA5ZiaiwTwCuBFYaY97teew7wEgAa+3dwFzgbGAD0AFcnYDjSpI1u92EIkEKHWK87lg3eHdYCVtEJJkGnLCttW8CjoOXNraH51cGeiwZXLeaRirp4kmHGJ8nVmF3q8IWEUmqhE06k6Hn06582oONjjE+dYmLSAKFQiEqKysJBoMHuilJFQgEGDFiBF6vt8+vUcKWXh3vzYPWJseYPWPYYTsYTRKRIa6yspKsrCwqKiowQ/TKE2st9fX1VFZWMnr06D6/TmuJS6+aXC6qol2OMd49XeKRwWiSiAxxwWCQgoKCIZusAYwxFBQU9LsXQQlbevXbcA2fy+h2jNndJd6tCltEEmQoJ+vd9ucclbClVxf4y/hWm3Pl7PPEvnQawxaRoeo3v/kNHR0d/X7dgw8+SHV1dcLaoYQtvZrmL+TMoHOFreuwRWSo25+EHYlElLBl8DQa2GTCjjHePV3iStgikvra29s555xzmDp1KpMmTeIHP/gB1dXVnHzyyZx88skAXH/99cycOZOJEydy22237XltRUUF3/rWt5gxYwZ//etfWbp0KVdccQXTpk2js7NzwG3TLHHp1SPBbdxTmMl71vY63qLrsEUkWX7w3GrWVO/X1hS9OnJYNredN7HX5+fNm8ewYcN44YUXAGhubuaBBx5gwYIFFBbGlpH6yU9+Qn5+PpFIhFmzZrFixQqmTJkCQEFBAcuXLwfgvvvuY86cOcycOTMhbVeFLb06PX0UP9/V1+uwNelMRFLf5MmTmT9/Pt/61rf45z//SU5Ozn/EPP7448yYMYPp06ezevVq1qxZs+e5T3/600lrmyps6dVhgSIOa2tzjNEYtogki1MlnCyHHXYYy5cvZ+7cuXzve99j1qxZH3p+8+bNzJkzhyVLlpCXl8dVV131ocuzMjIyktY2VdjSqyYirPV5CYd6H3vRWuIiMpRUV1eTnp7O5z73Ob75zW+yfPlysrKyaG1tBaClpYWMjAxycnLYuXMnL774Yq/vtffrEkEVtvRqbsd2/nd4Ga931pPvS99nzO4xbFXYIjIUrFy5km9+85u4XC68Xi933XUXixYt4swzz2TYsGEsWLCA6dOnM2HCBMrLyznhhBN6fa+rrrqK6667jrS0NBYtWkRaWtqA2qaELb06MWsMJaufIcPV+9dkzyxxJWwRGQLOOOMMzjjjjA89NnPmTG666aY99x988MF9vnbLli0fun/JJZdwySWXJKxtStjSq/K0Aso7OnEaOdFa4iIig0Nj2NKrZhthpc9HsLv3MRi3y+B2Ga0lLiKSZErY0qtFHVV8dngplW2VjnE+t0uXdYmIJJkStvTqqJwx/GFHLWW+XMc4r9tolriISJJpDFt6VRQooKgzCMbtGOfzuDTpTEQkyVRhS69ao2GW+f20BJ1XO/O6XYRUYYuIJJUStvRqbXAnVw0rYW3zZsc4n8el67BFZMi44447OOKII8jLy+P222/f7/fJzMxMYKvUJS4ODssazT01Ozn8uBLHOK8mnYnIEHLnnXfyyiuvMGLEiAPdlA9RhS29ygnkclywixyX1zHO63bRpS5xERkCrrvuOjZt2sRZZ53Fr3/9a2688UYgtmrZzTffzPHHH8+YMWP4+9//DkBbWxuzZs1ixowZTJ48mWeeeSZpbVOFLb3qIMKKgJ9xnQ0UOsT53EZd4iKSeC/eAjtWJvY9SyfDWb13c999993MmzePBQsW8Pzzz3/ouZqaGt58803Wrl3L+eefz6WXXkogEOCpp54iOzubXbt2ceyxx3L++ef3uiXxQKjCll7VdDfzpbISljZ94BinMWwRORRceOGFuFwujjzySHbu3AmAtZbvfOc7TJkyhVNPPZWqqqo9zyWaKmzp1bDM4TxQs5MxU0Y6xsXGsJWwRSTBHCrhA8Hv9++5bW1s3s4jjzxCXV0dy5Ytw+v1UlFR8aHtNhNJCVt6lebLZGawC1w+xziv20V7V3iQWiUicvBobm6muLgYr9fLggUL2Lp1a9KOpYQtveoGFqcFGNNZx3CHOK/bRbdmiYvIIeiKK67gvPPOY/LkycycOZMJEyYk7VhK2NKr9miYG0qLuaX5A65wiPN7XHSHtfmHiAwNu7fJvOqqq7jqqquA/9xSs62tDYDCwkIWLVq0z/fZHZMoStjSq6xALg9X72D42NGOcVqaVEQk+ZSwpVcebzpTu7rBOF+H7fe46AopYYuIJJMu65LeuX28mp7Ghk7nSxRUYYuIJJ8StvTO7eXrxYXMbd3gGBYbw1bCFhFJJnWJS++M4W876iksrXAM83m0NKmISLIpYYujCVE38TpifG43kaglHInicavTRkQkGfTXVRy9np7Gys5axxi/N/Y10ji2iKS6pqYm7rzzzgPdjH1SwhZHP8n287dO55V7fD1VtcaxRSTV9Zaww+EDv5qjusTF0T3tLjKynNY5+3eFrXFsEUl1t9xyCxs3bmTatGl4vV4CgQB5eXmsXbuWl19+mXPPPZdVq1YBMGfOHNra2vh//+//sXHjRr7yla9QV1dHeno69957b8JXPVOFLY4qTIAi67zsqCpsEUmWq+ddzdMbngYgFA1x9byreW7jcwB0hju5et7VzNs8D4DW7launnc1r2x9BYDGYCNXz7ua17a/BsCuzl1xj3f77bczduxY3n33XX7xi1+wfPlyfvvb37Ju3TrH182ePZvf/e53LFu2jDlz5nDDDTfs7yn3ShW2OHrTZ3B1N3C8Q4zf6wZUYYvI0HPMMccwerTzao9tbW0sXLiQyy67bM9jXV1dCW+LErY4utvbTVqk1jFh766wu7SeuIgk2ANnPrDnttfl/dD9NE/ah+5n+bI+dD8vkPeh+4Vphf0+fkZGxp7bHo+HaPTfhcnubTSj0Si5ubm8++67/X7//lCXuDiaQxH/G813jPF71CUuIkNDVlYWra2t+3yupKSE2tpa6uvr6erq4vnnnwcgOzub0aNH88QTTwCxvbLfe++9hLdNFbY4KnWnQ7fzjjO7E7a6xEUk1RUUFHDCCScwadIk0tLSKCkp2fOc1+vl1ltv5ZhjjmH48OEfmlT2yCOPcP311/PjH/+YUCjE5ZdfztSpUxPaNiVscbTIFaI12srpDjE+VdgiMoQ8+uijvT538803c/PNN//H46NHj2bevHnJbJYStjh7wjazydPpmLD9Hk06ExFJNiVscfT9wBiijc5jMaqwRUSSTwlbHOV5MyEccozZM+ksolniIiLJkpBZ4saY+40xtcaYVb08f5IxptkY827Pz62JOK4k37JoO3/3Oi/Jt7vC7gqpwhaRgbNxFmsaCvbnHBN1WdeDwJlxYv5prZ3W8/PDBB1Xkmx+uIFfZfocY/Z0iWvzDxEZoEAgQH19/ZBO2tZa6uvrCQQC/XpdQrrErbVvGGMqEvFecnC5MXsi161/2zHGrwpbRBJkxIgRVFZWUldXd6CbklSBQIARI0b06zWDOYZ9nDHmPaAa+Ia1dvVHA4wxs4HZACNHjhzEpklvMn2ZEAqCtWDMPmNUYYtIoni93rhLgR6qBmuls+XAKGvtVOB3wNP7CrLW3mOtnWmtnVlUVDRITRMnq8KtPJSdRSTc+7q4e5YmDWnSmYhIsgxKwrbWtlhr23puzwW8xpj+L+oqg25xdx1zCvLodljtzBiDz+OiSxW2iEjSDErCNsaUGhPrTzXGHNNz3PrBOLYMzOfyZ/DWlu0E2Hd3+G5+t0vXYYuIJFFCxrCNMX8FTgIKjTGVwG2AF8BaezdwKXC9MSYMdAKX26E8BXAI8fvS8VsLkW7nOK9LK52JiCRRomaJfybO878Hfp+IY8ngWh9q4bWcbD4drCc7u6zXOJ8qbBGRpNL2muJobVcDd+Tn0tixyzHO73WrwhYRSSIlbHF0VtF0lm3exki/857YsQpbs8RFRJJFa4mLI483PXYjEmc9cY1hi4gklSpscbQ91MKduTlUt1c7xmkMW0QkuZSwxVFNqI278nKoaq9xjPN5lLBFRJJJCVsczSyYxHubt3F05ijHOL9HXeIiIsmkMWxx5PL27CbjsDQpqMIWEUk2VdjiqD4c5Ld5OXzQus0xzu9x06VZ4iIiSaOELY5abTcP5mSzuWOHY5wqbBGR5FKXuDiqyBnDO1u2w6QxjnF+j0vba4qIJJEqbHHm8cd+h53XEvd5XHSFlLBFRJJFCVscBW2UOfm5LGnb4hin7TVFRJJLCVscRd1eHs/KZF1nrWOc3+OmOxxFm7CJiCSHxrDFUbo/m7e3VsGo+GPYAN2RKH6PezCaJiJySFGFLc6MiY1jx7kOe3fCDmocW0QkKZSwJa45edm80h7nOmxvrKrWtdgiIsmhhC1xvZzm4/3uBseYQE+FrZniIiLJoTFsievlJgv55Y4xgZ4KOxhShS0ikgyqsCU+t09j2CIiB5gStsT1u3TDE13O+2HvqbA1hi0ikhTqEpe4lnpgXKTNMUZd4iIiyaWELXE9FC0CV8AxJuDVpDMRkWRSl7jE14frsNUlLiKSXErYEtf9rg7utvEu69rdJa4KW0QkGdQlLnGtJ0TQhBxjdneJawxbRCQ5lLAlrv/NmAD1bznG+DXpTEQkqdQlLvH1Yy3xrrC6xEVEkkEJW+J6IlTHj9KdY/weF8aowhYRSRYlbIlrJ2E2xNkx0xiD3+NShS0ikiQaw5a4bsyfASvng7Wx7TZ7EfC6VWGLiCSJKmyJz+OP/Y53LbZHCVtEJFmUsCWu+Z3V3FRcSKS73TEu4HXpOmwRkSRRwpa4Womyw+MhHIqXsFVhi4gkixK2xHVx8dE8Ub0Dv3Wunv0eF0FNOhMRSQolbIlv9xh2KOgY5leFLSKSNErYEteSjmquKyliZ1v8PbG7lLBFRJJCCVviChkXLS4X4VCHY1xA12GLiCSNrsOWuI4vns7xNTvBm+kYp0lnIiLJowpb4uvjGLYu6xIRSR4lbIlrfWcd15YWs7p5o2NcwOsmGFaFLSKSDErYEpfx+Ok2EIm30pm6xEVEkkYJW+Ial384f66pZYo/3zHO74l1iVtrB6llIiKHDiVsic8TiP3uQ4UN2hNbRCQZlLAlroZwJ18oK2ZB01rHOL8n9nVSwhYRSTwlbInL7U3HY8EVCTvG7amwNY4tIpJwCUnYxpj7jTG1xphVvTxvjDF3GGM2GGNWGGNmJOK4Mjhy0vL4U20jn/QVOMbtTti6tEtEJPESVWE/CJzp8PxZwPien9nAXQk6rgwWb1ofxrBjXydd2iUikngJSdjW2jeABoeQC4A/25i3gFxjTFkiji2D48rCLB5tXecYE/DsrrCVsEVEEm2wxrCHA9v3ul/Z89iHGGNmG2OWGmOW1tXVDVLTpC9ycZMWdb5cy7+7wlaXuIhIwh1Uk86stfdYa2daa2cWFRUd6ObIXn4X9HORK9sxJs2rCltEJFkGK2FXAeV73R/R85ikCm+gz9dhdyphi4gk3GAl7GeBz/fMFj8WaLbW1gzSsSUBbvIH+VV3pWNMuq8nYXcrYYuIJFpCttc0xvwVOAkoNMZUArcBXgBr7d3AXOBsYAPQAVydiOPK4Ck1PgrijGGn+2Jfpw4lbBGRhEtIwrbWfibO8xb4SiKOJQfGdz3DoKPeMSZNXeIiIklzUE06k4OYNwBh5/2w0/Z0iTuviCYiIv2nhC19clu4mm942x1jvG6D22XUJS4ikgQJ6RKXoa/ck0lneKdjjDGGdK9bXeIiIkmghC19cm3W4bBlRdy4NJ9bs8RFRJJAXeLSN940CHXGDUvzqcIWEUkGJWzpkzva1nF5cQ5EnZcdTfO6NYYtIpIE6hKXPhkZKKA92A3hTvBl9BqnLnERkeRQhS19cmHBNL7d0Bi3WzxdXeIiIkmhhC19402P/e52vrQrzetRl7iISBIoYUufPN60ik+OHE5Hp9O257Euce3WJSKSeErY0iflGcM4tb0j7mpn6V43HVrpTEQk4TTpTPrkuKKpHFffCFHn6jnNp1niIiLJoApb+mb3GHacSWfqEhcRSQ4lbOmThc3rOH7kCNY0bXCMS/e6CUUsoYjz9doiItI/StjSJyWZwzmvrZ3sOHl4z45dqrJFRBJKY9jSJ2PzJ8Suw3b5HOP+vcVmhOyAdzCaJiJySFCFLX3jTYv9jjeG7Y0lbE08ExFJLCVs6ZPtwXqOGTWCuY2rHePS96qwRUQkcZSwpU+y0/L5VFsno0y8LvHYKEtnSNdii4gkksawpU9y/Dl8oz0MJuAYpy5xEZHkUIUtfedNx8ZZS1xd4iIiyaGELX1ireWoAi+/b1/vGBfw6rIuEZFkUMKWPjHGcHXYz1HWeQxbFbaISHJoDFv67EbyIOJ2jNmdsDWGLSKSWKqwpc+sN41QyHkMW13iIiLJoYQtfXYh1dxCvWOM3+PC7TK0d+myLhGRRFKXuPTZFb5ScjuqHGOMMWT43ErYIiIJpoQtffap9Aqo3hg3LtPvoV1j2CIiCaUucemzkCdAR7gjbly636MKW0QkwZSwpc++2baSK/PT4sZl+D20KWGLiCSUusSlz87LGk/L9pUQjYCr98u7Mv0awxYRSTQlbOmzWbmHQ1s7hDrAn9VrXIbPQ31b9yC2TERk6FOXuPRZl9tLo8sVd09sdYmLiCSeErb02V0Ny5k1cjjE2QAkQ13iIiIJpy5x6bOT8idTvOYFbFcbxiEuQ5d1iYgknBK29Nm0giOZ1toWG8N2kOHz0B2OEopE8brViSMikgj6ayp91uX2U+t2E+5qcozL8Mf+HahucRGRxFHClj57seE9Zo0cTk1rtWNcpj92yZcmnomIJI4StvTZtOJp3Lqrnpxo1DFud4WtLTZFRBJHY9jSZxX5E6hobYeIcyLO8MW+VqqwRUQSRxW29FnI46fK46Yj2OAYpzFsEZHEU8KWPtvUVsmZ5cNZ1LLJMS6jZwy7vUtd4iIiiaKELX1WljWMHzZ1cgR+x7jdXeKqsEVEEkdj2NJn2b5sLor4449h7+4S71bCFhFJFFXY0mfWWrYG0qnvanSMy/Rr0pmISKIlJGEbY840xnxgjNlgjLllH89fZYypM8a82/NzbSKOK4PLYjkvPcjfgs7XYQe8LlwGOjSGLSKSMAPuEjfGuIE/AKcBlcASY8yz1to1Hwl9zFp740CPJweOy7i43RQxPk7hbIwhw6cdu0REEikRFfYxwAZr7SZrbTfwN+CCBLyvHITO9pUyvrsrblyG36NJZyIiCZSIhD0c2L7X/cqexz7qEmPMCmPM340x5Qk4rhwAW71utoda4sZl+N2adCYikkCDNensOaDCWjsFmA88tK8gY8xsY8xSY8zSurq6QWqa9Mc3ghv5eSB+Is70e2jTGLaISMIkImFXAXtXzCN6HtvDWltvrd3dj3ofcNS+3shae4+1dqa1dmZRUVECmiaJ9o2cKVzbFL/Czgp4aQ2GBqFFIiKHhkQk7CXAeGPMaGOMD7gceHbvAGNM2V53zwfeT8Bx5QD4WPYYpna0QcS5ys5O89DSqYQtIpIoA54lbq0NG2NuBF4C3MD91trVxpgfAkuttc8CNxtjzgfCQANw1UCPKwdGlYnQ7PNyZHcbpOX2Gpfl99Ia1Bi2iEiiJGSlM2vtXGDuRx67da/b3wa+nYhjyYF1d+MK3iopYn53u3PCDniUsEVEEkgrnUm/fK70BH5aVw/dbY5x2WleOkMRQhHnvbNFRKRvtJa49MvheeMh2AVdzgk7KxD7arUGw+Rn+AajaSIiQ5oqbOmXWhvm7YCfcFezY1xWwAugmeIiIgmihC39Mr9xFV8sK6GtY5djXPZeFbaIiAycErb0y6wRn+T+mp1khJ0T8e4KW5d2iYgkhsawpV9Kc8dSGuyCULtj3O4x7BZV2CIiCaEKW/ql2cDCtADN7bWOcTlpGsMWEUkkJWzpl/WtW/lyaTHvt293jMvSGLaISEIpYUu/TMifwF8aupgYcTvGZfp3d4mrwhYRSQSNYUu/ZPoymebJgu4OxziP20WGz60KW0QkQVRhS7+EoiHeSAuwpTP+9qfasUtEJHGUsKVfojbKV7wtvBKujxubFfDQ0qkKW0QkEdQlLv3id/t51H8Yw+q3xY3NTvPS2qUKW0QkEVRhS79NTh9GQWdL3Djt2CUikjhK2NJvi+himXWedAa7x7CVsEVEEkFd4tJvv2lfR0Gmn6MiIXB7e43LDni0NKmISIKowpZ++8WIs7htVwMEnbvFc9K8NHeGsNYOUstERIYuJWzpt5FZIymJRCDOFpt56T7CUUtbl7rFRUQGSglb+m1VuIUXM9Ih6Jywc9Jj3eVNHeoWFxEZKCVs6bfnm1bzo4L8uF3ieek+QAlbRCQRlLCl37407jKerKqJW2Hn7q6wO7sHo1kiIkOaZolLvxXkjIRIBLqcK+zcni02G1Vhi4gMmCps6bft4VYez8qktd15PfHcni7x5g5V2CIiA6WELf22pm07PyrMZ0f7Dse4HFXYIiIJo4Qt/XbiiE/yak0ToyPO11f7PC4y/R5NOhMRSQCNYUu/pXvTSffnxB3DhliV3aQucRGRAVOFLf3W1t3GI5npbGiviRubl+GlScuTiogMmBK29FtHuIPbfZ0s73aedAaQm+ajURW2iMiAqUtc+q0wrZDX06aStWN13NicdC/VTZ2D0CoRkaFNFbb0m8u4yM8owdvZGDc2L11d4iIiiaCELfvliUgDrxOEaMQxLjfNR1NHN9GoduwSERkIJWzZLw+2reeFzPgbgOSme4laaNWOXSIiA6KELfvl8Qmzub2uHjoaHONy92wAoolnIiIDoYQt+yUjqyz25el0TtgFmbGEvatNCVtEZCCUsGW/vNG+nftzsuJW2IUZfgDq27oGo1kiIkOWErbsl4Ut63kgJ7vPFXZ9uypsEZGBUMKW/fLfR32dN7ZVxa2w9yRsVdgiIgOihVNkv3jTCsC4IM612H6Pm6yAR2PYIiIDpApb9svGls38qqiYuj6sJ16Y6WeXKmwRkQFRwpb9UtNewyPpPnZ21MaNLcjwUa8KW0RkQNQlLvvlhGEnsCxSBl3xK+fCTD+bdrUNQqtERIYuVdiyX4wxkF4Qd9IZxCaeaQxbRGRglLBlv4QiIX5umvhX9664sQWZfho7uglHooPQMhGRoUkJW/aL2+Xmqe6drIt2xN0ApCjTh7XQ2KFdu0RE9pcStuwXl3GxaML1XN3cAh31jrEFmbHVzjRTXERk/ylhy/7LKIr9btvpGFaQsXvxFI1ji4jsr4QkbGPMmcaYD4wxG4wxt+zjeb8x5rGe5xcbYyoScVw5sP7cuIL7crKhzfnSrsIsVdgiIgM14IRtjHEDfwDOAo4EPmOMOfIjYV8EGq2144BfAz8b6HHlwFvZWcN7fh+01znGlWQHANjZEhyMZomIDEmJuA77GGCDtXYTgDHmb8AFwJq9Yi4A/l/P7b8DvzfGGGutTcDx5QD5xYk/g389GrdLPNPvIdPvoaZZCVsOEGtjkyNtBKIRbDQM0QjGRrHRMB3dbXiMwW/c2GiYus5dZLh8ZLh9RCIhtrVXk+/JIMcdIBTpZl3bdkr8uRR6s+iOdLGqdSsjAwUUeLPoCAV5t2UzFf4iCryZtIU6eLd1C+P8JeR70mkNd7K8dTMTAmXku9JpCrezvH0LR/qHUeBKpyHcxtKOLUzzDSfXlU5DuIUlwW0c5R1OriuNulALS7q3c5y3nBwTYEekhSWhSk70jCTH+KmMtrAsVMVJ3gqy8LE90szySA2nuCvIMD62Rpt4J7KD092jScPNpkgT79paznKNxo+bDdFGVtg6zmU0XuNmXbSBFeziAsbgwfC+bWA19VzEGFwWVtkG1poGLrFjwVreo571polLoqOxFt41u9hkWrk4MgqApaaO7aaNiyKjAMti1y6q6eCiSDlgWeSqo9YEuTA8HIB/uutooJsLwsMAywJ3Ha0mxPmhMsDyiqeWThPhvFApWMtL3lrCWM4JlQCWF721YOGsUDEAz/l24LMuzgjFhvSe9u4kw7o4LVQIwJP+HeREvZwaygfgcd8OiqI+TgrlA5bHAjsoi/j5RCgPgNyTbmTMlBMG5WuciIQ9HNi+1/1K4GO9xVhrw8aYZqAA+NA1QcaY2cBsgJEjRyagaZJU/izwpMXtEgcoyfarwh6qrIVICMKdEApCOEh3Vyvh7jbSAcJBtrdWEQl3UOHOhHCQJS0bcYW7OcqbB5FunmnbQJqF072FEOnmgc6t5FrDRa4couFu5oSqGRY1XBbyYyNhvjPHnvcAAB9kSURBVO9tZnwYruyIYqIhrs+Bo7rCzG7pBBvl0tJsZnUEubGxBWOjnFJexkVtbXy1sRmAGRXlXN3cws2NzVjg2NEjuaGxieubWugyhlkV5Xy1oYlrm1voMIbzK8r5Rn0jX2hppdXl4vJRI/j2rgY+29pGo9vNF0YO57Zd9Vza2k6Dx8115cP5cV09F7S1U+318F8jhvHz2l2c1d5Bjc/Ld4aX8euddYzv6KTG5+OHw0v5/Y5aJnYGqfH7+MWwUv5YU8uUYJCagJ87ykq4v2YnY4JdVKcF+GNpMSdWP8O4rm62padxf0kRZ1S9wPjuEJvS03iopIgLN73EuFCI9Rnp/KW4kE9veplR4TBrMzN4tKiAL2yez7BwhFVZmTxWmM91m1+hMBLlvexMnijI56tbXyUnGmVpThb/yM/jm5tfI81aFuZm8VR+Lt/e9AYeDK/nZfFcbhbf27IQgFfzspmfk873tr2FxTA/P5PXs9L5TuXbWODFgkwWZQT4n6rlADxbkMmyNC/frH4HgCcLM1jl9/JfNe9hMfytMINNPjc371iBxVBdlE61x035ztUAVBWl0+Q2DN+5FoDK4jS6jGF47QcAbCsJYIBhdbvvp5FuLWV16wDYWhqgIGIp27UeTOz+sFCUsvp1WGBLWQBXOEpZQyx+U1mAtFCUsob1AOxs2JHE/7k+zAy0yDXGXAqcaa29tuf+lcDHrLU37hWzqiemsuf+xp6YXi/inTlzpl26dOmA2ibJ9XbN2zz+/Be5Lf9osi653zH2ivveoqM7wlM3DM6/RKWPIiHobILORmxHA6arGTob2daylV2ddcxwZUJXG6+1b6Mq1MwVkQB0tfGQbaTSdvPdxhbobufW/Gx2ut38cWdseOTLJUW0uVw8UhPrfbmmtJgo8OCO2D/uriwrwW8t9+2oJYyHK4YVkRux3F4bJIybG0rTKAnDTbWGkHXz4zJLWbebS3elEcLDPaXtDOvy8YnGHEJ4+EdRPaVdaUxuLiCKiwWFOykKpjOurYAILt4urKakO5uKziJwuVmWs42SUD7Du4vAuHkvcwul4QJKo8VEjeGDtK2URIspooioMWz0bqWYYgpcRVg3bDXbKHEVk+fOJ2IiVNkqitxFZLlziRChxlZT6C4i25NL2ETYGamhwFtMlieLMGF2heso8BWR7skgYiM0RuvJ9RWS4UknbMM0RZvJ9xUQ8AQI2zBt0XbyfLl4PX6iNkKH7STbm4Xb7Sdiw3RFu8jwZuF2uYnYKGHbTcATwG3cRGyEiAkTcAdwGTdRGyZsLX63F7fbQ9RGsUTxuX0YY7BYwOJ2eXC7XIDFABhwu1wYYgsnGcBlDBjo+YXLmJ7bPb/3vt3zlTMmduvf9z/8+KHOGLPMWjtzX88losKuAsr3uj+i57F9xVQaYzxADuB8LZAc9Jq7m/nA66a1bQdZcWJLsgO8tVEfedJFI9C+C9prsa07MO110LaTzU2b+aBtG2d2RaCzkWciTbzmDvHrHbHq4Jd5uTyblcHr22L/6z5QkMeC9HQWVNVifZnMz81ksd9wVms+QRNggz+Dre4ArwVm0uTx0RStJRiN8CN7Pi1hN02Nu+g2Lj7fPZag9dFUW083Xj7ZNYqg9RHZ0Y7Lk87H/cMI+L2kdUfA5+W/RwVI87mp8LrJ8Ht4eqybdK+bU31uAl431R4XPo+Lz3vc+L0u/G4Xfq+LW9yx+76e+xe4Xfi9bvweFx6XUTKQISERCXsJMN4YM5pYYr4c+OxHYp4FvgAsAi4FXtX4deo7bdRpnOY/Ehq3xI0tzQ5Q29pFNGpxufTHc78FW6B5OzRtp6NxE/7matwtVaxu3sT8UB3X19Xgj0Z5JDuTOfl5LNpaScBa5ucX8rucdE4Ol+LPHkG7u4CacD0NIz9HE5lkd+5kWncj9xw2i21BP2s7G+lq6GJcsJxIh4WmMOBixl7zVHPSvPwww0d+ho+8DB9FaV7saC/DAh4mpHnJDnjITvOSFfCQHfCS03M70+/B49YVpSL9NeCE3TMmfSPwEuAG7rfWrjbG/BBYaq19FvgT8BdjzAaggVhSl6Egswgq344bVpoTIBy17GrvojgrMAgNS2GRcOwfQfXrYdd6qN9Adf1aXm7fxnkNOymIRnkuI53vFBfyfHUdozLK2JiTy0NeD5cc+2XKc8cywXbwhY4q2k+7kk2hArLrG7m8oY7bWnLZ3hhka30H1U2dzNjzz+bDAHjV76EkJ0BZTiknFAconRKgONtPfk9iLsiI3c5N9+JV0hUZVAnZrctaOxeY+5HHbt3rdhC4LBHHkoNHR6iD73Wu52zbzqmRMLh7/zrtvrRrR3NQCXtvwRbYuRp2rKSlejnuuvfJqF3LBybC94oK+E59A9PdWVQVjuSXWV4mHPZFCoYdzxSvj6+2biDj0isgo4hTuoKMrG3n7dpOHq5tZcPONtbVtvLbl7dh7bY9hyvICFGen86MkXlcNH045fnpjMhNoyQnQGl2gAy/NvATOVjp/07Zb363n83RDppdBtp2QM6IXmNL90rYU3oPG9qshfoNsH0x0W1vsb56MRn1GxkRjrDN4+Gc8mH80F/IRUdfS37eSAp2voY583oY+QmmRrp5M9xJjj+HYChC445W/A2T+PmL1ayofJ/1tW1EorFy2es2jCnMZOqIXC6dUc74kkwqCjIoz08jK+A9wP8RRGR/KWHLfnO73Dx11Hdh7aXQXOmYsMtyDsHFU6yFhk2waQF2w6usrV6M7WrmyO4QoUAOl5fl8PnS0/mviVczomQiX9/2MpNHfALyxlEE3M2XAWjuDLFsayOLNzewZPMqVlY1E4rEknN+ho9Jw3M49YgSJg7LZnxJFhUF6RojFhmClLBlYHYn6eZKx7CCTD9ul2HHUE/YoU7YuADWzaNp0wK2d+xgcnc3JqecrxfnMy5jMr878Wf4C8bz2+o3OSzvMMgoxQVcPfkaALrCEZZsbuS1D2r518Z61u5owdpY5TxlRC7XfHw008tzmTQ8h+G5aZoBLXKIUMKWAblr+0vUFuRxW8tHr+T7MLfLUJodoKqxc5BaNoiCLbD+ZXj/OVo2vEJ2dxv4s7lleDlVJbk8d/qDkD+G23etoCyjDNJjKy59YsQn9rxFTXMnr6zZyWsf1LFwYz2doQg+t4uZFXl8bdZhHDM6n2nluaT53AfoJEXkQFPClgHpAjo8fmh2TtgAI/LS2D5UEnY0ApsWwLt/hbXPQzjIg4Vl3Dm8kNeOvZf0sadyQ+Na3MYNBWMBmFo09UNvUd3UydyVNcxdWcPybU0AlOencelRIzjp8CKOG1tAuk//i4pIjP4ayIB87aivwVt/gzgVNkB5fjpvrHPeKOSg17QNlt4P7/6V7cE67i4o4ropF1M+9fMcnZ7JNdX/IjL6RPD4mFI05T9e3tYV5oUV1Ty+tJJlWxsBOKIsm2+cfhhnTipjbFGGurhFZJ+UsGXgcobHFvOIY2R+OrWtXQRDEQLeFOratRa2LoTFd9G99gXaXC7yx56K94izeX3tvcyadhnlI49jIjCxaPI+Xm5ZsqWRx5du54UVNXSGIowrzuSbZxzO2ZPLGF2YMfjnJCIpRwlbBuSd2neYww5+0raTijix5flpAFQ2djKuODPpbRuwaBTWPgdvzIEdK4im5XHJ2AkcWTaTn836HaXAgqmfxeva96VSwVCEZ9+r5v43N7N2RyuZfg8XTh/GZTPLmV6eq0paRPpFCVsGxO/2k+7NoDvYFNupydv7oijleekAbG/sOLgTdjQCa56BN35BW937vF48mnPO+y2uyZ/i6q3zYhPHeuwrWTe0d/PQwi08sngru9q6mVCaxc8vmcK5U8s0Ji0i+01/PWRAjiw4knsPvwreXxwbx+6ZYLUv5fmxhF3Z0DFIrdsPGxfAy9+Dnaug8HCeOOFqflX1ChPHnUSFL52Lx1/c60sb27u595+beGjhFtq7I8yaUMwXPz6a48YWqJoWkQFTwpaB230tdtM2x4RdlOnH73Gx7WBM2HUfwMvfx65/iQVFIyk88wdMOeYmPhUJckzztVTkVPT60ubOEPe+sYkH/rWZjlCEc6cM4+ZTxjG+JN4eZiIifaeELQN207q/MDovl683bgZO7jXO5TKxS7saDqJLu7o74LX/hUV/AF8G3bNu46c75jKju5qfu9xkuDKYWDhxny8NR6I8+vY2fj1/HU2dIc6eXMZXZ43nMCVqEUkCJWwZsNKcURRgoH5j3NiR+elsPVgq7I0L4PmvEW7cwtyJp3PuWX/An1nMfc0XMyLLecHz19fV8ePn17C+to3jxhTw/XOP5Mhh2YPUcBE5FClhy4B999jvwbLnY+tmxzGmKJNFm+oP7L7YwRaY921492HIH8ur5/yI7665l9ymtXwis9ix+7u2NcgPnlvDCytqqChI554rj+K0I0s0Ri0iSaeELYlRMDa2E1Uc44ozCYaiVDV17pmENqi2L4Env4ht3s7OY6+jdNYPOM3j508jjuWYsmN6fZm1lieWVvLjF9YQDEf5xumHMfsTY/F5tMmGiAwO/bWRAXt+0/PMCq2jpXFL7NplB7sv59pQ2zYILdtLNAKv/wLuPwOw/PLEa7i8eTHN0S6MMY7JekdzkCv/9Db/8+QKJpRl8+JXT+TGU8YrWYvIoFKFLQNWkl7CcdljCUc3xC7tyi3vNXZc0b8T9skTigengZ2N8OS1sOEVmPwpOGcO53XsoGTH22T5nCeIvbiyhlv+sZLucJSfXDSJzxw98sB15YvIIU0JWwbs6NKjOXry9fDOi9Cw0TFh52X4KMjwDV6FvWMVPHYFNFcx78Tr2Vo4mi8Hcjg8kMPh+Yf3+rL2rjA/fG4Njy3dzpQROfz28ulaQlREDij16UliFIzFQp9mio8tzmRD3SAk7DXPwp9Og3AXXD2Xt9ICLKxeSCgacnzZpro2LrrzXzy+bDtfOXksT15/vJK1iBxwqrAlIS567SaOKSrm23Vr48aOK87khRU1WGuTN7t60Z3w0ncIjTiK9gvvIrfwML4zbBoG0+va3wDz1+zk64+9i8dt+Ms1H+Pj4wuT0z4RkX5ShS0JMWvkLKaklcLO1XFjDy/JorkzxI6WYOIbEo3GLtl66dsw4Rz+Z8xEZi/6HqFoCJ/bh9e972QdjVp+9fIHfOnPS6kozOC5mz6uZC0iBxVV2JIQN06/Eao2wsonY9tROlTOk4bnALCyspmynLTENSLcBf/4Umzjjo9dB2f8lIurF7Krc5djVR0MRfjvx9/jhZU1XHbUCH504aTU2v5TRA4JStiSMKHiCbi7mnE1VzpOPDuyLBuXgVVVzZw+sTRBB++Exz4HG16h4ZTvsm7ciRzrcnPiiBMdX9bQ3s2X/ryUZVsb+e7ZR3DtiaO1CIqIHJTUJS4JMW/LPGa+/wcqPZ643eJpPjfji7NYWdWcmIN3t8Mjl8GG/4Pzf8cc08TXX/9vWrpbHF+2eVc7F9/5L1ZVNXPnFTP40ifGKFmLyEFLFbYkxGF5hzF74lUEtv0Ydq6Ew890jJ80PIfX19UNfOJZsAUe/RRsXwwX/RGmfppvdTXzqcM/Rbav97W9V1Q28YX738YYw19nH8uMkXn73wYRkUGgClsSYkzOGL4y8+sUZ5fDjpVx4ycPz2ZXWxc7W7r2/6BdrfDwxVC5hNYL7+SPNBGOhsnx5zCteFqvL1u2tZEr7l1MZsDDUzccr2QtIilBCVsSJhwNs2vYNNj+dmzimYOp5bkALN/WuH8HC3XCXz8DVcvhsgd5JT3A3e/dzdoG58vKFm+q5/N/Wkxhlp/HZh/HqAJdXy0iqUFd4pIw171yHUFbw8OtNdC0FfIqeo2dNDyHDJ+bRRvrOXtyWf8OFO6Gx78AW96Ei++BI87jImB68XTHnbbeXL+La/+8hBF56Tx67ccozg7077giIgeQKmxJmM9M+AyfO+xTsTvb3nKM9bpdzKzI561N9f07SDQCT82G9S9hz/klvw5Vsa1lG4Bjsl7wQS3XPLSEioIM/jb7WCVrEUk5StiSMLNGzuLMKdeAPwe2LYobf+yYAtbXtlHX2sdx7GgUnrsZVj8Fp/2ImiPO4qn1T7Fg+wLHl728egdf/vMyDivJ5K9fOpbCTH/fjicichBRl7gkjLWWmo4d+EbMoHDrwrjxx47JB+CtTfWcN3VYvDePrV72zsPwyW/BCTczDPjHBf+gIFDQ68teWFHDV//2DhOH5/Dna44hJ633BVRERA5mqrAlYTrDnZz55Jk8kZcHu9ZB4xbH+MnDc8hJ87JgbW38N1/wE1h8Nxx7A6snX8AT654AoDCtsNfLwp5+p4qb/rqcaeW5PPxFJWsRSW1K2JIw6d50fnriTzlj8tWxB9bOdYz3uF2cekQJr7y/k+5wtPfAN38Db/wCZnwezvgpj33wGPetuI+OUEevL3l8yXb+6/F3+djoAh665hiyAkrWIpLalLAloc4dcy5jRn0Cio6AD5wTNsBZk0ppCYZZ1Nvks7fvhVdug4kXw7m/AWO47bjbeOish0j3pu/zJQ+/tZX/eXIFHx9XyP1XHU2GXyM/IpL6lLAloYLhIAurF7Jr/CmwdSG073KM//j4QjJ8buauqPnPJ995BOZ+Aw4/m+hFd3Pf6gdo7W7F7XJTmrHvNcjvf3Mz33t6FbMmFHPv52eS5tMmHiIyNChhS0JVt1Xz5flf5p8Fw8FG4N1HHOMDXjfnTCnj2feqae4M/fuJVf+AZ2+EMSfDpQ+wpmk9f3jnD8zfOn+f72Ot5fevrueHz6/hzIml3PW5o7TjlogMKUrYklCjc0bzx9P+yOlTroJRH4clf4pdO+3g88dV0BmK8PdllbEH3n8utk1m+cfg8kfAG2BS4SSevOBJLhp30X+83lrL7fPWMufldVw8fTi//+x0fB59tUVkaNFfNUkoYwzHDzueDG8GHP3F2IpnH7zo+JpJw3M4alQeDy7cTHj5o7FVzIZNh88+zrtN61m6YykQW6/8ozPCI1HL959ZxR9f38QVHxvJnMum4nHray0iQ4/+sknCNXc18/Cah9k+YhoUjIf/+wFEQo6vuenksZzS/DSeZ6+Hio/DlU9j/Vn8Zvlv+MninxDZR5XeGgzxxYeW8PBb2/jyJ8fw4wsn4XJpe0wRGZqUsCXhguEgP1vyMxbuWAKn/zh2TfZr/9v7C8JdnLTuJ/zA+xCv2plUnfMQ+DMxxnDHKXdwxyl34HZ9eDx6a307l9y1kH+u38VPLprEt886QntZi8iQputdJOFKMkp4+ZKXKcvs2dRj+pXwz19CVhkcfS3snVi3LoTnvga7PqDpqJv52tLjGfW3FVx2UjVXHBnb03rvfa2jUcvfl1Xyg+dW43G7+PM1x3DCuMJBPkMRkcGnhC1JsSdZA5zzS2irjV2i9f5zMP602PaYG/4Ptr8FOSPhiifJHX8qvxq7k688cz+/WPYoHlvIZyefBkA4EuWN9XX89pX1vFfZzLFj8vnVp6YxLDftAJ2hiMjgMjbOvsUHysyZM+3SpUsPdDNkP1lr+fmSn1OUXsQ1k66JzRRffDe8dRc0b48FFR0BM66Eo64C37/3pV7wQS3/9dQzNDUNY3RhBlkBD5t3tdMaDFOWE+Abpx/ORdOHa7xaRIYcY8wya+3MfT2nCluSwhhDVVsVXlfPkqAuNxz3FTj2BuhsBLcP/Jkfes3TG55mZslMTj58BK9/9SqeWLqdJVsa6OiOcP7UYZw4vpBTJpToki0ROSSpwpakCUfDeFx9+zdhc1cz5z51LqeNOo1bj7s1yS0TETk4Ja3CNsbkA48BFcAW4FPW2sZ9xEWAlT13t1lrzx/IcSU17E7WG5s2kuZJY1jmf26hufsfjDn+HP5y1l/2GSMiIgO/rOsW4P+steOB/+u5vy+d1tppPT9K1oeQYDjINS9dw8/e/tl/PBeKhLj51Zv3bJVZkVOBz+0b7CaKiKSEgY5hXwCc1HP7IeA14FsDfE8ZQgKeAL/85C8ZlT0KgA2NG6jtrOX4YcfjdXuJ2AhR67C1poiIAAMcwzbGNFlrc3tuG6Bx9/2PxIWBd4EwcLu19ule3m82MBtg5MiRR23dunW/2yYHp9vfvp1/rP8Hb3z6DQKeAJFo5D8WRREROVQ5jWHHTdjGmFeAfe1l+F3gob0TtDGm0Vqbt4/3GG6trTLGjAFeBWZZazc6HVeTzoam+s56mrqaGJk98t8zyEVEBBjgpDNr7akOb7zTGFNmra0xxpQBtb28R1XP703GmNeA6YBjwpahqSCtgIK0ggPdDBGRlDPQSWfPAl/ouf0F4JmPBhhj8owx/p7bhcAJwJoBHldEROSQMtCEfTtwmjFmPXBqz32MMTONMff1xBwBLDXGvAcsIDaGrYQtIiLSDwOaJW6trQdm7ePxpcC1PbcXApMHchwREZFDndZ4FBERSQFK2CIiIilACVtERCQFKGGLiIikACVsERGRFKCELSIikgKUsEVERFKAEraIiEgKUMIWERFJAUrYIiIiKUAJW0REJAUoYYuIiKQAJWwREZEUoIQtIiKSApSwRUREUoAStoiISApQwhYREUkBStgiIiIpQAlbREQkBShhi4iIpAAlbBERkRSghC0iIpIClLBFRERSgBK2iIhIClDCFhERSQFK2CIiIilACVtERCQFKGGLiIikACVsERGRFKCELSIikgKUsEVERFKAEraIiEgKUMIWERFJAUrYIiIiKUAJW0REJAUoYYuIiKQAJWwREZEUoIQtIiKSApSwRUREUoAStoiISApQwhYREUkBStgiIiIpQAlbREQkBShhi4iIpAAlbBERkRQwoIRtjLnMGLPaGBM1xsx0iDvTGPOBMWaDMeaWgRxTRETkUDTQCnsVcDHwRm8Bxhg38AfgLOBI4DPGmCMHeFwREZFDimcgL7bWvg9gjHEKOwbYYK3d1BP7N+ACYM1Aji0iInIoGYwx7OHA9r3uV/Y8JiIiIn0Ut8I2xrwClO7jqe9aa59JZGOMMbOB2T1324wxHyTy/Q8ChcCuA92IBBuK5wRD87x0TqljKJ7XUDwnSPx5jertibgJ21p76gAPXgWU73V/RM9j+zrWPcA9AzzeQcsYs9Ra2+vkvFQ0FM8JhuZ56ZxSx1A8r6F4TjC45zUYXeJLgPHGmNHGGB9wOfDsIBxXRERkyBjoZV0XGWMqgeOAF4wxL/U8PswYMxfAWhsGbgReAt4HHrfWrh5Ys0VERA4tA50l/hTw1D4erwbO3uv+XGDuQI41RAzF7v6heE4wNM9L55Q6huJ5DcVzgkE8L2OtHaxjiYiIyH7S0qQiIiIpQAk7weItw2qMucoYU2eMebfn59oD0c7+MMbcb4ypNcas6uV5Y4y5o+ecVxhjZgx2G/dHH87rJGNM816f1a2D3cb+MsaUG2MWGGPW9Cwb/NV9xKTU59XHc0rFzypgjHnbGPNez3n9YB8xfmPMYz2f1WJjTMXgt7Tv+nhOKfc3EGKrdhpj3jHGPL+P5wbnc7LW6idBP4Ab2AiMAXzAe8CRH4m5Cvj9gW5rP8/rE8AMYFUvz58NvAgY4Fhg8YFuc4LO6yTg+QPdzn6eUxkwo+d2FrBuH9/BlPq8+nhOqfhZGSCz57YXWAwc+5GYG4C7e25fDjx2oNudgHNKub+BPe3+OvDovr5ng/U5qcJOrD3LsFpru4Hdy7CmNGvtG0CDQ8gFwJ9tzFtArjGmbHBat//6cF4px1pbY61d3nO7ldiVGR9dWTClPq8+nlPK6fnv39Zz19vz89FJRRcAD/Xc/jswy8RZC/pA6uM5pRxjzAjgHOC+XkIG5XNSwk6svi7DeklPV+TfjTHl+3g+1Qzl5WeP6+nee9EYM/FAN6Y/errlphOrcvaWsp+XwzlBCn5WPd2s7wK1wHxrba+flY1dItsMFAxuK/unD+cEqfc38DfA/wDRXp4flM9JCXvwPQdUWGunAPP597/K5OCzHBhlrZ0K/A54+gC3p8+MMZnAk8DXrLUtB7o9iRDnnFLys7LWRqy104itAHmMMWbSgW7TQPXhnFLqb6Ax5lyg1lq77EC3RQk7seIuw2qtrbfWdvXcvQ84apDalkx9Xn42lVhrW3Z379nYWgJeY0zhAW5WXMYYL7HE9oi19h/7CEm5zyveOaXqZ7WbtbYJWACc+ZGn9nxWxhgPkAPUD27r9k9v55SCfwNPAM43xmwhNsx5ijHm4Y/EDMrnpISdWHGXYf3IWOH5xMbjUt2zwOd7Zh8fCzRba2sOdKMGyhhTunscyhhzDLH/Xw7qP5Y97f0T8L619le9hKXU59WXc0rRz6rIGJPbczsNOA1Y+5GwZ4Ev9Ny+FHjV9sxsOhj15ZxS7W+gtfbb1toR1tqK/9/OHaNEEARRGP4fewUTEy/iDTYwMthAQUMTD2AieAoTBTEQDMXYE5gtZnsONxkog14DF4MRwbHh/+IJuiioN8w0RZvpL1V1tPXYn/TpV5vO9FVVDUk+17DOgNuqektyBbxW1RNwnuQAGGgXnk4mO/BISR5ot3B30lbRXtIuk1BV17QtdnNgBbwDp9Oc9GdG1HUInCUZgDWw+M/DcmMfOAaWm/+IABfAHnTbrzE19dirXeAuyYz2gvFYVc9b8+IGuE+yos2LxXTHHWVMTd3NwO9M0Sc3nUmS1AE/iUuS1AEDW5KkDhjYkiR1wMCWJKkDBrYkSR0wsCVJ6oCBLUlSBwxsSZI68AFZMulYI/iD/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(r, start_pot[0], label='start')\n",
    "plt.plot(r, end_pot[0], label='final')\n",
    "plt.plot(r, true_pot, label='true', linestyle=':')\n",
    "plt.legend()\n",
    "plt.ylim(-1,2)\n",
    "plt.show()"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}