{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 04. EDS biasing using HTF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the collective variable (CV) being biased is the average distance to center of mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd\n",
    "import hoomd.md\n",
    "import hoomd.dump\n",
    "import hoomd.group\n",
    "import hoomd.htf as htf\n",
    "import tensorflow as tf\n",
    "tf.get_logger().setLevel('ERROR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the computational graph\n",
    "\n",
    "We first create a CV, the distance from the center of mass for each atom. Then we use EDS to bias the CV to match our target value. Note we add the forces ourselves to the compute graph and they depend on the atom positions, not the distance between atoms. Hoomd-TF requires you to be extra sure you want position dependent forces, because often you can accidentally implicit create position dependent forces. So we add `positions=True` to our `compute_forces` call to say we're sure about computing position dependent forces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Build training graph ####\n",
    "\n",
    "def make_eds_graph(NN, set_pt):\n",
    "    # N= Number of atoms in the system \n",
    "    # NN=Number of nearest neighbors \n",
    "    # set_pt=set point in EDS method\n",
    "    graph =htf.graph_builder(NN,output_forces=True)\n",
    "    #calculate center of mass\n",
    "    com = tf.reduce_mean(graph.positions[:, :2], 0) \n",
    "    #calculate distance of each atom from center of mass\n",
    "    rs = graph.safe_norm(tf.math.subtract(graph.positions[:, :2], com), axis=1) \n",
    "    #calculate the average distance from center of mass. This is the collective variable (CV)\n",
    "    real_cv = tf.reduce_mean(rs) \n",
    "    #calculates the running mean of the CV and value\n",
    "    graph.running_mean(tensor=real_cv,name='cv_run')\n",
    "    graph.save_tensor(tensor=real_cv,name='cv')\n",
    "    #calculate the EDS alpha value every 300 steps. \n",
    "    eds_alpha = htf.eds_bias(real_cv, set_point=set_pt, period=300,learning_rate=5.0)\n",
    "    eds_energy = eds_alpha * real_cv #computes EDS energy\n",
    "    #compute EDS forces\n",
    "    eds_forces = graph.compute_forces(eds_energy, positions=True)\n",
    "    graph.save('eds-model',force_tensor=eds_forces,virial=None)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the simulation\n",
    "\n",
    "This simulation is 64 LJ particles in an NVT ensemble. We save data every 10 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: You did not provide a virial for eds-model, so per particle virials will not be correct\n",
      "Note: Backed-up eds-model previous model to eds-model/previous_model_2\n",
      "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 64 particles\n",
      "notice(2): -- Neighborlist exclusion statistics -- :\n",
      "notice(2): Particles with 0 exclusions             : 64\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:00 | Step 3000 / 3000 | TPS 28167.7 | ETA 00:00:00\n",
      "Average TPS: 28059.1\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "103 normal updates / 30 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 0 / n_neigh_max: 30 / n_neigh_avg: 15.3906\n",
      "shortest rebuild period: 5\n",
      "-- Cell list stats:\n",
      "Dimension: 2, 2, 1\n",
      "n_min    : 10 / n_max: 22 / n_avg: 16\n",
      "** run complete **\n",
      "notice(2): Force mode is FORCE_MODE.tf2hoomd \n",
      "notice(2): Starting TensorflowCompute \n",
      "notice(2): completed reallocate\n",
      "notice(2): Setting flag indicating virial modification will occur\n",
      "notice(2): TF Session Manager has released control. Starting HOOMD updates\n",
      "** starting run **\n",
      "Time 00:00:11 | Step 8566 / 18000 | TPS 556.543 | ETA 00:00:16\n",
      "Time 00:00:22 | Step 14500 / 18000 | TPS 550.859 | ETA 00:00:06\n",
      "Time 00:00:28 | Step 18000 / 18000 | TPS 595.259 | ETA 00:00:00\n",
      "Average TPS: 562.765\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "973 normal updates / 150 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 22 / n_neigh_max: 57 / n_neigh_avg: 42.2188\n",
      "shortest rebuild period: 9\n",
      "-- Cell list stats:\n",
      "Dimension: 2, 2, 1\n",
      "n_min    : 12 / n_max: 22 / n_avg: 16\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": [
    "#### Hoomd-Sim code ####\n",
    "\n",
    "make_eds_graph(32, 4.0)\n",
    "\n",
    "hoomd.context.initialize(\"--mode=cpu\")\n",
    "with htf.tfcompute('eds-model', device='CPU:0') as tfcompute:\n",
    "    #cut off radius: must be less than the box size\n",
    "    rcut = 6.0 \n",
    "    #initialize the lattice\n",
    "    system = hoomd.init.create_lattice(unitcell=hoomd.lattice.sq(a=2.0),n=[8, 8])\n",
    "    nlist = hoomd.md.nlist.cell(check_period=1)\n",
    "    #enable lj pair potential\n",
    "    lj = hoomd.md.pair.lj(rcut, nlist) \n",
    "    #set lj coefficients\n",
    "    lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0) \n",
    "    hoomd.md.integrate.mode_standard(dt=0.005)\n",
    "    # set up NVT simulation\n",
    "    hoomd.md.integrate.nvt(kT=1.0, tau=0.5,group=hoomd.group.all()) \n",
    "    #equilibrate\n",
    "    hoomd.run(3000)\n",
    "    #simulation\n",
    "    tfcompute.attach(nlist, r_cut=rcut,save_period=250)\n",
    "    hoomd.run(15000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "Now we plot the CV value and its running average to assess if EDS converged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xUVfr48c+ZSSa9kAQSQksIEFpCIBEEBIIFQWxYKJaVtbu6ll3dlZ8rqLuurvpd66orNqygoOKqYEVBqgFCDz2UQEgCgfRMO78/ZjIESEJIZjLJ8Lxfr3llZu6de565SZ45c+65z1Vaa4QQQvgeg7cDEEII4RmS4IUQwkdJghdCCB8lCV4IIXyUJHghhPBRft4OoLaYmBidkJDg7TCEEKLNWL16dZHWun1dy1pVgk9ISCArK8vbYQghRJuhlNpT3zIZohFCCB/l0QSvlIpUSs1VSuUopbYopYZ6sj0hhBDHeXqI5kVgodb6GqWUCQj2cHtCCCGcPJbglVIRwEhgKoDW2gyYPdWeEKJuFouF/fv3U1VV5e1QRDMEBgbSuXNn/P39G/0aT/bgE4FC4B2l1ABgNXCf1rq89kpKqduB2wG6du3qwXCEODvt37+fsLAwEhISUEp5OxzRBFprDh8+zP79+0lMTGz06zw5Bu8HDAJe01oPBMqBh09eSWv9htY6Q2ud0b59nTN9hBDNUFVVRXR0tCT3NkwpRXR09Bl/C/Nkgt8P7Ndar3Q+nosj4QshWpgk97avKb9Djw3RaK3zlVL7lFLJWuutwAXAZk+1d7JlO4rI2lNMkL+RQJORQD8DQSYjvePC6NEhrKXCEEIIr/H0LJo/Ah86Z9DsAn7v4fYA+C33CL97exVW+6m17v0Mir+N78NNw2Q8UoiWkp+fz/33389vv/1GZGQksbGxvPDCC4wdO5YFCxaQnJzsWvf++++nY8eO/PWvf3U9l5ubS58+fUhOTsZsNpORkcFbb711RgccG2P69OmMHDmSCy+80C3by87OZuDAgSxYsICxY8e6ZZtnQrWmC35kZGTo5p7Jeqikiktf/pUQk5HP/jAck5+BSrONKouN95fn8unq/RRXWOjcLohhSdFEBpuYNq43SimsNjt+Rjn3S/iWLVu20KdPH6+1r7Vm2LBh3HTTTdx5550ArFu3jpKSEhYsWEBAQAAzZswAwG6307VrV5YuXUq3bt1c28jNzeXSSy9l48aN2Gw2LrroIm655Rauv/56r7ynxvrrX//KsmXL6N69O7NmzWr29ur6XSqlVmutM+pa36eymdlq5w8frqG82sp/b8wgIsifZTuKCDIZ6RIVTFKHULpGBdM+NID9xZXMW5PHRyv3unryf5m3nuFP/8Qd72fx+i87OVxW7eV3JETbt2jRIvz9/V3JHWDAgAGMGDGCKVOmMGfOHNfzixcvplu3bick95MZjUYGDx5MXl4e4ChxUlRUBEBWVhaZmZkAPPbYY9x8881kZmbSvXt3XnrpJeD4t4HbbruNfv36MWbMGCorKwGYOnUqc+fOdW13xowZDBo0iJSUFHJycgAoLCzkoosuol+/ftx6661069bN1X5tWms+/fRT3n33Xb7//nuqqqrIyclh8ODBrnVyc3NJSUkB4JtvvqF3796kp6dz7733cumll57Zjq5Dq6pF01x//2ozq/cU88p1A8nJL+Hmd38j72glT1+VwuTBXZl0juMG8PPWAu6bnY1da37KOcT5vWMZlhSD2WpnY94xvt10iBd/2M6do5K478KeXn5nQrjPpP8uP+W5S1M7cuPQBCrNNqa+s+qU5dekd+bajC4cKTdz1werT1g2546GT1DfuHEj6enpdS5LSUnBYDCwbt06BgwYwOzZs5kyZUqD26uqqmLlypW8+OKLDa4HkJOTw6JFiygtLSU5OZm77roLgO3bt/Pxxx8zc+ZMJk6cyLx587jhhhtOeX1MTAxr1qzh1Vdf5bnnnuPNN9/k8ccf5/zzz2fatGksXLiQt956q862ly1bRmJiIklJSWRmZvL1119z9dVXYzab2b17N4mJicyZM4dJkyZRVVXFHXfcweLFi0lMTDztPmgsn+nBf5q1j/dX7OH2kd25NDWeNxbvAuD1GwZxTXrnU9bPTO7AV388j65Rwdz23mqW7SjimvTOvHLdIH5+aDQ//GkkY/vHUTNMb7drCkrkRBEh3G3KlCnMnj0bq9XKF198wbXXXlvnejt37iQtLY3Y2Fg6duxIamrqabc9fvx4AgICiImJoUOHDhw6dAiAxMRE0tLSAEhPTyc3N7fO11911VWnrPPrr78yefJkAMaOHUu7du3qfO3HH3/sWm/y5Ml8/PHHAEycONH1raUmwefk5NC9e3fXHHd3JXif6MFv2H+MR77YyNDu0fzlYsfBmmqrnbQukYzt37He13WJCmbOHUOZ8J+l3PPxWr68Zzid2zmqKfToEMbzk9KoOUbx3eZD3Dt7LVPO6cJdmT2Iiwj0/BsTwgMa6nEHmYwNLo8KMZ22x36yfv36uYY96jJ58mTGjBnDqFGjSE1NJTY2ts71kpKSyM7OpqioiOHDh/Pll19y+eWX4+fnh91uBzhlnnhAQIDrvtFoxGq11vl8zRDNyWrWq/3axrDZbMybN4/58+fz5JNPuk5UKi0tZdKkSVx77bVcddVVKKXo2bMn2dnZjd72mWjzPfhjFRbu/GA1MSEmXrluoOsgqdlqx994+lkyoQF+vPG7DCw2O3e8v5pKs+2E5TXj8/3iw7lqYCc+XLmXkc8sYsb8jeQfkx69EKdz/vnnU11dzRtvvOF6bv369SxZsgRwJO6YmBgefvjhRvVcY2JiePrpp3nqqacAx1j56tWOYaN58+Z54B2caPjw4XzyyScAfPfddxQXF5+yzo8//khqair79u0jNzeXPXv2cPXVV/P555+TlJSE0Wjk73//O5MmTQIgOTmZXbt2ub4l1D4u0RxtPsGHBfox6ZwuvHZDOtGhxz+V3556Dg9enNzAK49LjAnhpckD2XywhGmfraeumUVdooJ5+upUFj2YydXpjkR/3Zsr6lxXCHGcUorPP/+cH374gaSkJPr168e0adOIi4tzrTNlyhRycnJcQyKnc+WVV1JRUcGSJUuYMWMG9913HxkZGRiNRk+9DZcZM2bw3Xff0b9/fz799FPi4uIICzvx3JqPP/6YCRMmnPDc1Vdf7RqmmTRpEh988AETJ04EICgoiFdffZWxY8eSnp5OWFgYERERzY7V56ZJNscrP23nue+28bfxfbh1RPcG1913pIKDx6oYnBhFlcXGMwu3ctOwbnSLDmmhaIVoHG9Pk/Q11dXVGI1G/Pz8WL58OXfddZdbhljKysoIDQ1Fa83dd99Nz549eeCBB05Y50ynSfrEGHxdPlixhz4dw0jvFtXo19w9ugcb80p4akEOfTuGM6xHTL3rdokKpkuUY7w+e99RPli5h3eX7eaKtE7cPTpJzpYVwkft3buXiRMnYrfbMZlMzJw50y3bnTlzJrNmzcJsNjNw4EDuuOOOZm/TZ3vwvf62gN8PT2DauDPruZRVW5nwn6UUlVWz4L6RjT6YWlBSxcwlu/hgxV6qrDbG9I3l+UlpBJt89jNUtBHSg/cdZ/WJTjW01pitdgKacFZqaIAf/70xnUqLjUc+39DoMfYO4YE8Mr4vSx8+n7sze1BlsRPk7xgP3HTgGLY6yiYIIYQn+WSCt9gcydTk17S31719KA+OSebHnAK+yM47o9dGhZh48OJkZt08GKUUxeVmrn5tGaOf+5mPVu6VRC+EaDE+meDNNsec2KYmeIDfD09kUNdIHvtyMwWlTZ8OGR7kzwuT0ogONfH/Pt/Alf9ZSva+o03enhBCNJZvJnirM8E3o3CY0aB45poBVFpszJi/qVnbGdu/I5/dNYyXpwykoLSKq19bRt7Ruk+sEEIId/HJBB8R5M+Sv4xmwqBTSxSciR4dQrn/wp4s2JjP1+sPNmtbSikuGxDPj3/O5OUpA+kUGQRAVu4RmUsvfF5+fj6TJ08mKSmJ9PR0LrnkEgwGA1u3bj1hvfvvv59//etfXorS9/hkgjcaFF2igokIan6t6NtHdCelUwTT52/kSHnzrxkeGuDHJSmO8glbDpZwzevLufOD1RytkOuRC9+ktWbChAlkZmayc+dOVq9ezVNPPcWoUaOYPXu2az273c7cuXNd9VtE8/lkgj9SbuY/i3awo6Cs2dvyMxp49tpUSqosPP6/pg/V1CU5Noy/je/DTzkFjHtxCSt3HXbr9oVoDeorF/zSSy+dcalgcWZ8cpJ2/rEqnv12K0ntQ+jRIbTZ2+sdF849o3vy/A/bGJ4Uw8RzurghSjAYFLeO6M6QxGhHIbOZK/jzmGTuHt3DLdsX4mSP/28Tmw+UuHWbfePDmXFZv3qX11cuuCmlgsWZ8ckevDtm0ZzsrswkzusRw8OfrWf+GU6dPJ2UzhF89cfzmDCwM8Emz9fSEKK1aGypYNE0PtmDt9QkeDcWHjL5GZj5uwymvrOKP32yDpPRwLiU+ksRn6mQAD/+b+IA1wHX7zcfwt+oyEzu4LY2hGiop+0pDZULbmypYNE0vtmDt7q/Bw+OWtlvTz2HtC6R/PHjtfyw+ZBbtw+O2TZaa2Yu3sXUd37juW+3YnV+YAnRFjVULvhMSwWLM+PTCb4x9eDPVEiAH+/8/hz6xYfzhw/X8PPWAre3oZRi1s2DmZTRhVcW7eCGt1bK1aREm3W6csFnWipYNJ5PFhuz2OyUVlkJC/TDvxknOzXkWIWFKTNXsLOwjH9c2Z+rBnXGaHD/B8rc1fv52xcbCA3w57sHRhIVYnJ7G8K3SbEx3yHFxgB/o4GoEJPHkjtARLA/H9w6hD4dw3lo7noueXEJP2w+5PaTlq5J78yX95zHdYO7uJJ7tdV2mlcJIYSPJvj1+4/yzMIcjlVYPNpOVIiJz+4axivXDcRss3Pre1lc8/pyVu0+4tZ2esWG8acxjqtTbTtUyvCnF/Fp1j45A1YI0SCfTPAb80p49eedVFo839M1GBSXpsbz3QMj+eeEFPYXVzDxv8v55zdbPNKev9FAYkwwD81dz41vrWJnYfNP5hJC+CaPJnilVK5SaoNSKlsp1WLX4rN4YB786fgbDVw3pCs/PziaiRmdeWPxLpbtKHJ7O4kxIcy5fSj/uLI/6/YdZczzi3nif5vd3o4Qou1riQw4WmudVt9BAE/w5Cya0wkyGXn88v4kRAfzl3nrKa+2ur0Ng0Fxw7ndWPRQJjcM6eo6OUprTaVZxueFEA4+OUTjiTNZz0SQycgz1wwg72glzyzM8Vg7MaEBPH5Ffx682DE+v3h7Eef96yfeW57r+hYjhDh7eToDauA7pdRqpdTtda2glLpdKZWllMoqLCx0S6PuqAffXIMTo5g6LIFZy/ewooWKiMWEmujRIZTp8zdx0b9/4av1B+RArGgVjEYjaWlprtvTTz8NQGZmJsnJyaSmptK7d2/uuecejh49fkGcJ598kn79+pGamkpaWhorV648ZdsrVqxgyJAhpKWl0adPHx577LEGY8nOzuabb76pc9nPP/9MRESEa1uPP/54g9uaPn06P/zwQ4Pr/PzzzyxbtqzBdTxGa+2xG9DJ+bMDsA4Y2dD66enp2h1sNruuNFvdsq3mKK+26JHP/KRH/OsnXV5taZE27Xa7/mnLIT3m37/obn/9St/8zqoWaVe0Xps3b/Z2CDokJKTO50eNGqV/++03rbXW1dXV+k9/+pMeOXKk1lrrZcuW6XPPPVdXVVVprbUuLCzUeXl5p2yjV69eOjs7W2uttdVq1Zs2bWowlnfeeUfffffddS5btGiRHj9+vNZa67KyMt2jRw+9evXqRrzD+s2YMUM/++yzzdpGjbp+l0CWrienerSLq7XOc/4sAD4HBnuyvRoGgyLQ3/tFu4JNfjxzdSp7j1TwzMKtp3+BGyilGN27A9/cN4Lnrh3AZQPiAce3mqxc907fFMKdTCYTzzzzDHv37mXdunUcPHiQmJgYAgICAIiJiSE+Pv6U1xUUFNCxo6MulNFopG/fvgCUl5dz8803M3jwYAYOHMj8+fMxm81Mnz6dOXPmkJaWdkK54pOFhISQnp7Ojh07yM7O5txzzyU1NZUJEyZQXFwMwNSpU111dhISEpgxYwaDBg0iJSWFnJwccnNzef3113n++edJS0tjyZIlbt1np+OxYmNKqRDAoLUudd4fAzzhqfZqm5+dx9b8Uv4ytndLNNegId2jmTosgXeX5TKufxxDuke3SLtGg+Ka9ONXtJqfncdDc9czJDGKP57fk+E9olGq5Q9CCy9b8DDkb3DvNuNSYNzTDa5SWVlJWlqa6/G0adOYNGnSKesZjUYGDBhATk4O48eP54knnqBXr15ceOGFTJo0iVGjRp3ymgceeIDk5GQyMzMZO3YsN910E4GBgTz55JOcf/75vP322xw9epTBgwdz4YUX8sQTT5CVlcUrr7zSYMyHDx9mxYoVPProo0yZMoWXX36ZUaNGMX36dB5//HFeeOGFU14TExPDmjVrePXVV3nuued48803ufPOOwkNDeXBBx9ssD1P8GQPPhb4VSm1DlgFfK21XujB9lyW7TjMvDX7W6KpRvnL2GS6RgUz48tNXhsTvzQ1numX9iX3cDk3vLWSK19dxsKN+TJGL1pEUFAQ2dnZrltdyb1Gzd9kaGgoq1ev5o033qB9+/ZMmjSJd99995T1p0+fTlZWFmPGjOGjjz5i7NixAHz33Xc8/fTTpKWlkZmZSVVVFXv37j1trEuWLGHgwIGMGTOGhx9+mM6dO3P06FHXh8tNN93E4sWL63xtTT2d9PR0cnNzT9uWp3msB6+13gUM8NT2G2K22b02g6YuwSY/7h6dxF/nbSBrTzHnJES1eAxBJiM3n5fI9ed2Ze7q/bz2805eWbSdsf0dBZ/yj1URFxHY4nGJFnaanra32Ww2NmzY4Kq3YjQayczMJDMzk5SUFGbNmsXUqVNPeV1SUhJ33XUXt912G+3bt+fw4cNorZk3bx7JycknrFvXgdraRowYwVdffeV6fOzYsUbHXzOcZDQasVrdP0X6TLWeLOhGZqvdqzNo6nLZgHjCAvz4cMUer8YR4Gfk+iHdWPRgJq/f4LjKTkmVhcznFjH+pSW8/etu8o9J5UrR8iwWC9OmTaNLly6kpqaydetWtm/f7lqenZ1d5+X8vv76a1evf/v27RiNRiIjI7n44ot5+eWXXcvWrl0LQFhYGKWlpY2OKyIignbt2rnGz99///06h4rqc6btuVPryoJuUm21Y/Lz/kHW2oJNfkwY1IlvNua75eLdzeVvNNC5XTAARqV4eGxvDErxxFebOfepH7n6tWWs3lPs5SiFr6gZg6+5Pfzww65l119/PampqfTv35/y8nLmz58PQFlZGTfddBN9+/YlNTWVzZs31zkF8v333yc5OZm0tDRuvPFGPvzwQ4xGI48++igWi4XU1FT69evHo48+CsDo0aPZvHnzaQ+y1jZr1iweeughUlNTyc7OZvr06Y1+75dddhmff/65Vw6y+mS54Nvey+JwWTWf/WG4G6Jyn5z8Esa+sIRHLunDbSO7ezucOu0oKGPhxoN8syGff08aQO+4cFbuOkzWnmLG9I2lR4dQOTjbxki5YN9xpuWCffKSfTN/12JVEc5I77hwMrq146NVe7l1RGKrTJQ9OoRyz/k9uef8nq7nlu08zIs/bufZb7eSEB3MRX1juahvHOcktGuV70EI4eCTQzSt2XVDurK7qJzlO1vm7FZ3eOCiXiyfdj5/v7I/XaNDeHdZLg/MyXYtX7K9kF2FZTIjR4hWxid78P/33VYig03ccl6it0M5xSUpHXniq818uHIvw3rEeDucRusYEcSN53bjxnO7UVplIbeownX92AfmrKOorJpOkUEM7xHN8B4xDO8RQ0xogLfDFk5aa/m21cY1pQPlkz347zYdYtXu1tlDDvQ3cs2gzny7KZ/C0mpvh9MkYYH+pHSOABxnzs67y1G+OKVTBAs35nPf7Gxe+WkHAFabnUU5BZRUefbiK6J+gYGBrmmDom3SWnP48GECA89sKrNP9uAd8+Bb1yya2qYM6cqbv+7mk6x93D26h7fDabZu0SF0iw7hhnO7YbNrNuYdIyzQ8ae1Ie8Yv3/3N5SCPnHhDE6MYnBiFMN7xBAR5O/lyM8OnTt3Zv/+/birmJ/wjsDAQDp37nz6FWvxzQRvtXulFnxjJbUPZWj3aD5etZe7RiVh8MDFur3FaFAM6BLpetynYzgf3TqElbuP8FvuEWb/tpd3l+Xy0a1DGNYjhq35pWzMO8bgxCg6twuSYQQP8Pf3JzGx9Q1XCs/zzQRvsxPQis5krcv153blno/Wsnh7IZnJHbwdjscE+hsZ1iPGdbzBbLWzIe8Y/eLDAfhmw0Fe/NFxMktseAD94yPoFx/OXZk9CDK13m9hQrQFPpngwwL9CG/lX//H9I0jJtTEhyv3+nSCP5nJz0B6t3aux/dd0JNxKXGs2n2ENXuK2XSghFW5R7j/wl4A/PObLazZU0yPDqGuW8/YMDpFBnnrLQjRZvhkgv/pz5neDuG0TH4GJmZ04fVfdrLvSAVdooK9HZJXGAyK3nHh9I4L53dDEwBHL79m2CoqxIRBKb7bfIjZv+0DICE6mJ8fGg3Au0t3A5AcF06v2FCiZeaOEC4+meDbihuHduONxbt4e+luZlzWz9vhtBq1C8XdOSqJO0clAXC4rJodBWVUWY9fjnD2b/vIyT9e5yMqxMSVaZ2YfpmjJnhW7hFiwwOJjwzC6EPHOoRoDJ9L8Ha75o4PVjNhYCcuSeno7XAa1DEiiMsGxPPJb/u4/8JeMqvkNKJDA07poS+4bwQFpdXk5Jey/VApOwrKiI90TCWz2uxMmbkCi01jMhroGh1MQnQIEwZ2YnxqR2x2zcrdh4kLDyQ2PJCQAJ/7dxBnOZ/7izbb7Hy/+RADu0aefuVW4JbzEvl8bR6zV+3lDmdPVTSeUopYZ4Ie1av9Kcvfu3kIuYfLyS0qd/6sIO9oBQBFZdVcN/N46diwAD9iIwK5Z3QPrhzYiWOVFr5cd4D4iEDiIgKJjwgiMtjfrTN9NuYdY9OBY1yT3kW+YQi387kEX90KLrh9Jvp3imBo92jeXZbLzecl4t9G4m4L/IwGhiZFMzSp7qtoRQT589FtQygoqSa/pIr8Y45beJDj32J3UTmPfrHxhNcE+ht4fmIa41I6sruonDm/7SMuPICYsACiQkxEhwTQLTr4tJeM1Frz0aq9PP7lZsw2O59m7effE9PoGn12HosRnuFzCd7sTPCtfZpkbbeNTOTmd7P4ev1BrhzYydvhnDUC/Y0MS6q/XERKpwhWTLuAg8cqOXisynE7Wkli+xAAcovKeevXXVhsJ54h+tFtQxiWFMPX6w8yff5GAv2NBPobCDIZ8TMYePqqFN76dTefrt5PWIAf7cMCWLO3mMznFhEfGcQnd5xLfGQw87PzmL1qHyY/A/5GAwF+BkICjDwyvi8RQf6s3VvM1vxSQgP9CAnwI9R56xUbhtGgMFvtGA1KvhmcxXwvwducPfg2lOAze3UgqX0IM5fs4oq0eDnZp5UwGhRxzuGZgXUsH927A1v/Po4jFWYOl5k5XF7NkXIzveMcc/zjIwMZlxJHpdlOlcVGpcVGaZWFe2evZduhMi5L7UhJlQWlFF2jgthysJT9xZU8PG8D/56UhtZgs2uOVlqwWO1UW22UV9t4ZLyj/YWb8vnvL7tOiWvrP8ZiNBj5x9ebeW/5Hkx+BoKcHzIhAX6uWWYv/rCdpTuLCDYZCTH5EWQyEh1qYto4RznaxdsKKSqrdn1whAQ4ph8nxjg+4Gx2LR8erZzPJXitNV2jgtvUAUuDQXHriO5M+2wDy3cdbrBXKVoXg0ERExrgLKwWdsKygV3bMbDr8Tn/y3YW8YcP12Cza978XQYX9o09YX27XfP20t088+1Wxr6whP/9cXiD3+juu6AnvxuaQHm1lbJqK2VVVirMVgKcZTpG9+5AVIiJSouNSrONaoud2n0Hfz+FAo6Um9l3pIJKs43QQD9Xgn9n6W4WbT2xvEFiTAiLHswE4Po3V5CVW+z48DAZCfI3ktIpgv9cPwiARz7fwIGjlZj8DPgZDRiVIjkuzFWe49/fbeVopQWjQeFvNGA0KHrFhjJhoON0/PnZedi1JsjfSLDJ8QETGx7gulCN2dq6Ls3ZGvnkBT/aoiqLjeFP/0Ral0jemnqOt8MRbra/uIJxLy4hLjyQmb/LIMHZC67L5gMlXP3aMjIS2vHezYO99o3uWIWFo5VmSquslFZZqbRYMRoMroPZn2TtY8/hcirNdiotNqosNuIiAvnr2N4A/GlONjsKyzBb7Vhsduwa0rpE8vykNACueOVX9hypwGbTWO0aq93OBb1jef1Gx6UkBz7xHcUVJxapu2pgJ/7tfH3y3xZg19r17SI0wI8JAztxx6gk7HbNQ3PXExJgJMhkJNDPSKC/kXMS2pGREIXZaueXbYUEm4zOmx/BJiNRIaY2N5vqrLvgR1sU6G/khnO78eKP29lRUEaPDqHeDkm4ic2ueWBONlrDWzedc9oDqX3jw/l/l/Tm0fmb+GjVXq4fcup1SFtCRLA/EcH1fxOemNGlwdfXJOL6zL/nvFOeq93hXHDfSCotNirMVirNNsrNNqJDTK717r2gp+ubS3m1lZIqK6HOIndVVhsrdh2m3PnamskX917Qk4yEKI5WmrntvVM7kw+P682do5LYe7iCcS8uJsjkR4CfgQB/AwF+Ru49v4frAPv0+RvxMyj8jAb8nMc6pg5LICMhit1F5byzdDf+RsfxE5NRYfIzMD41nsSYEPYXV7Bke5FzuWJ07w6EB7p/1MHnEvy2Q6U89uUmHh7Xm9TObWOqZI0bh3bjtV928tavu3nqqhRvhyPc5D+LdvBbbjHPTxrQ6Fky1w/pxrebDvHk11sY0aN9m55dMz87z3WBm9oDBkkdQvjd0IQTZhzV/rYSF1F/aVylVIOVWINNfix9+HzXY7tdU209PkQVGWTif/ecRyNalTQAACAASURBVIXZSoXZRoXZRrnZSqqzDHaQycjkwV2pMNuottowW+1UW+2u3r3Nbqes2oq15tuHzY7Nrl1lsQtLq5mffQCrzY7Fpl3HBvvFR5AYE8KmAyVM+2yDK74f/zzKIwne54Zolu88zJSZK1zVCtuaaZ+t57M1efzwp1FnbfkCX7J6TzET/7ucy1I78sLkug7V1i/vaCVjn19Mn47hzL793DZZdfTnrQVMfec3IoP9XTPbFAqN5lBJNV2jgpl+ad9Tjkf4Gq0dHwQG5ejpV1lsHK2wYLE5hq86twtu8vGEs2qIxtIGZ9HU9ofMHny17iD3fLSGT+8c1mbfh4CSKgv3zV5Lx4hAnriy/xm/vlNkENMv68tDc9fz9tLd3Dqi6Rdqr7LYOFxuRmuN1mDXGruGDmEBHhtzPniskgfmZNM7Lowv7h5+yrkBS3cUMePLTdz6Xhbn9+7A9Ev7Nnhsoi1TSp1QwjzQ30hchOerpXo8wSuljEAWkKe1vtTT7dXMg2+ribFLVDDPXJPKXR+u4Z/fbOGxy6VGTVs1/YuNHDxWxSd3DG3y1+9r0h1X/3rm261kJrenR4ew07/IKbeonEVbC1i0tZAVuw67/jdqiwox8cEtQ+jrLN/sLhabnT9+tBaz1c5/rh9U54lfw3vEsOC+Eby7NJcXftjGmOcXc/foHtx7QY82MVVYa83+4krW7jtKzsESRvVqz5DudZ9U5y0t0YO/D9gCuPcvqB5tcR78ycaldOT3wxN4Z2kuQxKjGFdHTR2rzc7XGw4ytHs0HcLP7DJewvM+X7ufL7IP8MCFvU4oj3ymlFL886oUxjy/mD9/so73bh7S4IHP3KJyPl61l+82H2J3UTkA3WNCuGFIN3rFhmIwOKZGGpTCrjX//n4b1725gg9vHUK/+Igmx3my//tuG1l7inlxchpJ7eufMOBvNHDbyO5cnhbP37/azPM/bMPP2PD4ujdVWWy8v3wPq3KPsHbvUYrKjl9289Wfd3Jhnw78dWxvesY2/oPYkzya4JVSnYHxwJPAnzzZVo1gk5HecWGEmNr26NO0cX1Ys/cof5m7nn7xESccZNtysISH5q5jY14JnSKDmHXzYJl104pszDvGI59v5JyEdtw9uvn1hTqEBfLklSnc/dEaznnyBy7s24EJAzuTmdwef6MBm13zy7YCZi3bwy/bCvEzKIb3iGHqsAQyk9vTLbr+YY/BiVFMeWMF181cyYe3DqF/p+Yn+Z9yDvH6Lzu5bkhXrkhr3JnZseGBvDxlIP5GA89+u5Vu0cFcmhrf7Fjc6VilhdtmZbEq9wiJMSGM7BXjONehSyQJMSG8v3wPry7awcUvLGbSOV24/8JexHq58+XRg6xKqbnAUzjOAHnwdEM0GQkROmvGqVOnzlZVVhsb8o4R6Gd0XQEp72glB45W4mdQxEcGceBoFRpNclwYYQFt5+QuX1VttbHxQAkKR50hd9ZEKjdbKSytpqisGqtd42dQtAs2UVJlodpqx99oIDYsgA7hgWfUbpXVxuYDJdi0pk9cOKENjMlrNGXVVgpKqimuMBPgZyAs0N9xkZ1Af+xasyHvGCY/A/3jIzCc4VCLXWu2HCyhzGylb8fwVvM3bbbZyTlYQqXFRlKHUGJC6r7ugMVuJ6+4kkMlVSiliAzyJ8jkmIsf5O+41Xx7stkdB15tdo1d6yYP46mbv2n5g6xKqUuBAq31aqVUZgPr3Q7cDpAaL1fpqS3Qz0hS+1C2HSplR2EZVWYbFRYbMSEmusWE4G8w0C7ExJaDJWw5WErP2FDaBZm8HXabodGUV9soKqum2monJtREuxATBpo2/mu129maX4rdrukXH+72gnchJj9Cov3oGh3MsQoLhWXVHC43ExJgpGtUcJNjD/Qz0jc+3PF3lF9y/BuwwrU1q01TWFZNYWk1lRYbBqVoF+yPxWbnUGkV+SWO9QzKMQTUq0PYGSf3mtf3ig1j44FjbMsvpV+nCAL9vHvpxkqLjZz8Eiw2R0cqsoH/MX+DgYToEOIiAskrrqS0ysqRCvMJ6xgU2PWpr2vOUF59PNaDV0o9BdwIWIFAHGPwn2mtb6jvNe6YJjk/O4+3l+by4a1DGuyJtCVPfr2ZmUt2ExsewJNXppwypaywtJqp76wiJ7+UZ65O5er0M7vyOjjOtPzvL7v4KaeAkAAj4YH+RAQ5bjFhAVzcL45BXSPbxMGv09lRUMaX2XnMX3eAPYcrMPkZiAzyp6C0mvZhAUw+pwuTB3c9o8sCWm12bpmVxa87inhn6jmMrKN0cWu3v7iCKTNXsO9IZb3rDOwayaSMLlw6IN71/1Vznd3fco+wdm8xEzO6cEGf5k173FFQxoRXl9IxIpC5dw1z6xzxvKOV+BtUo45dbdh/jKnvrEID7/7+nCadW1NlsZF7uJzth8rYUVBGpcVGRJA/4YF+hAX6Ex7kR0SQqckJvqFpki0yD97Zgz/9EI0bEvzrv+zk6QU5bHlirM9ctNlis7NgYz6jerWvt8ZOaZWFO95fzbKdh7nvgp7ce0HPRhWCyi0q57WfdzJvzX6Uggt6O/4xj1VaXLfCsmrMVjvJsWFMHtyFqwZ2rvNAn9a6VX8AFJRW8fC8DfyUU4BBwbCkGC5Pi2ds/zhCTH78sq2AD1fs5aetBSjg/N6xPDK+j6u4Vn201jzyxUY+WrmXp69KYfLgri3zhjygoKSK+dkHMNvs2O0am3M6pcmoGNMvjl4tePBw6Y4ibnp7FRkJ7bh8QCcig/2JDHKcXdsu2ESHsAD8GvktqazayjcbDjJv9X5W7j6Cyc/AAxf24rYRiXVuw2bXzF29jyf+t5nIYBPv3zKY7g0cLPamsyrBv/Tjdv79/TZ2PDmu0b98X1FttTHtsw18tiaPcxLa8fykNFdhppNtO1TKaz/vZH52Hn5GA1PO6cIdo5KIr6PXWlZt5X/rDvDxqr2s33+MAD8DF/eLI9hkpLC0moLSatfYcJeoYEb0jGFkz/acmxTdpG9RFWYruUUV5B4uZ39xBWXVNiqqrZSbbVSarZhtdoZ2j2Z8ajxRIY0bkvp+8yH+Om895dVW7r2gJ9emd663B7e/uII5v+3j3WW5mK12/jymF7ec173OD0y7XfP64p08s3Ard2UmueqwCPf4JGsfj3y+4ZSSzOCo9hkfGUiXdsGOW1QQIQF+J5QOUMCynYdZuDGfSouNxJgQrhrYic0HS1iwMZ+UThE8e22qqwIoOKpo/vObLeTkl5LerR2vXj/I6wdLG+L1BN9Y7kjwz327lVd/3sGup8a7Kaq25/O1+3n0i00oBU9OSOHyAY7ZCFprft1RxMwlu1m8rZAgfyPXD+nK7SO7N3qq5aYDx5i9ah9frT+A0WCgQ1gA7cMC6BAWQFSoiW35pazYdYRKiw0/g2JQt3aM7RfHNRmd6/2aXV5t5ZOsfXy7KZ/dReUcKqk+ZZ3aRaFsdk2e80BzZnJ7rkjrxEV9Y+uca11ptvGPrzfz4cq99O0YzktT0ho9l/xQSRV/+2Ij328+xIAukTx7TaqrB7u/uIK5q/fzadZ+8o5WMj6lIy9PGdgmzzZt7WrO+jxaaXb8rLBQXGEmr7iSvUcq2Fdcwb4jlSdMWawtLNCPS1PjuSa9E4O6tnN9y/xmw0Ee/WIjJVUW7h7dgwv7xPLMt1tZvK2QLlFBPDy2D5ekxLXqb6VwliX4p77ZwqzlueT8fZx7gmqj9h6u4P45a1mz9yhXDerEuYnRvL10Nzn5pcSEBnDT0G5cf263RveAz0S11cbq3GIWby/il22FbDlYQrDJyDXpnfnd0ATXlM78Y1W8uyyXj1buoaTKMWuiT8dwEmOCSYgJISE6hK7RwYSa/E5InFprthwsZX52HvOzD5BfUkVogB/JcWGu66vGRQQQGWzi9V92squwnDtGdudPY3q5Suk2ltaar9YfZMaXmyitsnDjuQlsLyjl1x1FAJzXI4Zr0jszPqXjWfeNsbWpcla0tNgcM1MszvowcRGB9V5h60i5mSf+t4kvsg8AEB7ox70X9OTGod3O+G/FW86qBD9rWS4/bDnE+7cMcVNUbZfVZueln3bwyk/bsWvoHRfGLeclcnlafIv+8a7ff5R3l+Xy1bqDmG12RvSMISY0gP+tO4Bda8b178itIxJPqJ3eWDa7ZuWuw3y14SC7C8s5VFJFfkkVFWYbAHHhgfzfxAEMb2ZdosNl1Tz+v818ue4AnSKDuDajM9ekd653CEy0LT/lHGLdvmP8fngCkcFtaybaWZXgxam2HCyhpNLC4MQor37dLCqrZvaqvby/Yg+lVVYmZnThlvMS3V5UTWtNabWVgpIq4iODCHbjSW8FJVXEhAbIUIxoNSTBi1bF5jy5oy2XkxCitWgowfvcf9jfv9rM1HdWeTsM0QCjQUlyF6IF+Nx/WV6x41R+IYQ42/lcgjfb5EK8QggBPpjgLTZH0SUhhDjb+VwmrLba3V7kSQgh2iLfqMZVS1qXSPxkCpsQQpw+wSulegGvAbFa6/5KqVTgcq31PzweXRP8v0v6eDsEIYRoFRozljETmAZYALTW64HJngxKCCFE8zUmwQdrrU+eWG71RDDucOV/lvLYl5u8HYYQQnhdYxJ8kVIqCdAASqlrgIMejaoZ8o9VUV7daj9/hBCixTTmIOvdwBtAb6VUHrAbqPeqTN5mkXnwQggBNCLBa613ARcqpUIAg9a61PNhNZ3ZKgleCCGgcbNopp/0GACt9RMeiqlZqqUHL4QQQOOGaMpr3Q8ELgW2eCac5hvXP45+8RHeDkMIIbyuMUM0/1f7sVLqOeBbj0XUTC9OHujtEIQQolVoylhGMNDZ3YEIIYRwr9MmeKXUBqXUeudtE7AVeMHzoZ25oxVm+k5fyIcr93g7FCGE8LrGjMFfWuu+FTiktW6VE82rrXbXtTiFEOJsV2+CV0pFOe+ePC0yXCmF1vqI58JqGrPVDiDVJIUQgoZ78KtxnL1aV2lGDXT3SETNUF2T4GWapBBC1J/gtdaJLRmIO0gPXgghjmtUPXilVDugJ4558ABorRef5jWBwGIgwNnOXK31jKaHenrhQX5MGdyFrtHBnmxGCCHahMacyXorcB+OqZHZwLnAcuD807y0Gjhfa12mlPIHflVKLdBar2hmzPXq3C6Yp65K9dTmhRCiTWnMWMZ9wDnAHq31aGAgcPR0L9IOZc6H/s6bbmqgjWG3a7T2aBNCCNFmNCbBV2mtqwCUUgFa6xwguTEbV0oZlVLZQAHwvdZ6ZR3r3K6UylJKZRUWFp5J7Kf4ZXshidO+Ye3e4mZtRwghfEFjEvx+pVQk8AXwvVJqPtCoM4m01jatdRqO4Z3BSqn+dazzhtY6Q2ud0b59+zOJ/RQ1B1n95SCrEEI0OA/+IeBjrfUE51OPKaUWARHAwjNpRGt91PnascDGpgZ7OjUJPkCmSQohRIM9+HhguVJqiVLqD0qp9lrrX7TWX2qtzafbsFKqvbPnj1IqCLgIyHFP2HUzyzx4IYRwqTcTaq0fALoCfwNSgPVKqYVKqZuUUmGN2HZHYJFSaj3wG44x+K/cEXR9zDZJ8EIIUaPBaZLaMSXlF+AXpdQ9wIXA08BrOKpKNvTa9Thm3LSYXrFh3HJeImGB/i3ZrBBCtEqNPdEpBZgMTAKKgGmeDKqp0ru1I71bO2+HIYQQrUJDB1l74kjqkwEbMBsY47xGa6tUZbFhs2uCTUbXpQWFEOJs1dBg9UIcZQYmaa1Ttdb/bM3JHeC1n3fSb8a3yLlOQgjRcLGxpJYMxB3MNjv+RoXBIL13IYTwqekmZqtdTnISQggnn8qGFptdpkgKIYRTvdlQKfWQUqpNXVzbbLVLLXghhHBqaJpkzZmsucDHwKda6+ZVA/OwC/rE0qNDqLfDEEKIVsGTZ7K2uIv6xnLriFZ3JUEhhPCKBscznDXdf9Fa34WjIuTzwP3AoZYI7kwdLqvmSPlpy+QIIcRZwafOZP3zp+soLjcz/57zvB2KEEJ43enOZJ2CI6m3iTNZzVaZRSOEEDUa6sEvxHFwdZLW2mM13N1JErwQQhzXUIIfC8SenNyVUsOBfK31To9G1gRmm53QwEaNOgkhhM9rqLv7PHCsjudLgBc8E07zyDx4IYQ4rqHubqzWesPJT2qtNyilEjwWUTPcNqI74UFSC14IIaDhBB/ZwLIgdwfiDlent6kTb4UQwqMaGs/IUkrddvKTSqlbgdWeC6npdhaWUVRW7e0whBCiVWioB38/8LlS6nqOJ/QMwARM8HRgTXH1a8u4fEA8T1zR39uhCCGE1zVUD/4QMEwpNRqoyZhfa61/apHImkAOsgohxHGnnVOotV4ELGqBWJpN5sELIcRxPpMN7XaN1a4lwQshhJPPZEOzzQ4gCV4IIZx85rRPo0Hxr6tTSOnU0OxOIYQ4e/hMgvc3Gph0TldvhyGEEK2Gx8YzlFJdlFKLlFKblVKblFL3eaotgCqLjbV7iymWevBCCAF4dgzeCvxZa90XOBe4WynV11ON5R+rYsKry1i0tcBTTQghRJvisQSvtT6otV7jvF8KbAE6eao9OcgqhBAnapFs6CxONhBYWcey25VSWUqprMLCpl/T22x1JHh/OdFJCCGAFkjwSqlQYB5wv9a65OTlWus3tNYZWuuM9u3bN7kd6cELIcSJPJoNlVL+OJL7h1rrzzzZVk0PPkB68EIIAXh2Fo0C3gK2aK3/7al2avToEMrLUwbSKy7M000JIUSb4Mnu7nDgRuB8pVS283aJpxqLCQ3gsgHxxIQGeKoJIYRoUzx2opPW+ldAeWr7JysorWL7oTIGdo0k2OQz528JIUST+cyA9fKdh7n+zZUcPFbl7VCEEKJV8JkEX+08yCr14IUQwsFnsqFrFo1MkxRCCMAHE7zMgxdCCAefyYZyopMQQpzIZ6abjE/pSHJsGAF+Rm+HIoQQrYLPJPguUcF0iQr2dhhCCNFq+Mx4xuYDJfy45ZC3wxBCiFbDZxL8J1n7uH9OtrfDEEKIVsNnErzFZpc58EIIUYvPZESz1S4zaIQQohafyYhmmyR4IYSozWcyotkqQzRCCFGbz0yT/POYZCrNNm+HIYQQrYbPJPgeHUK9HYIQQrQqPjOmsWhrAUu2N/2i3UII4Wt8pgf/0o/bCQ3wY0TPpl+4WwghfInP9ODlIKsQQpzIZzKizIMXQogT+UxGNNvs+EsPXgghXHwmI1qkBy+EECfwmYOs790yhCCT1IIXQogaPpPgZR68EEKcyGfGNN5bnsuavcXeDkMIIVoNjyV4pdTbSqkCpdRGT7VR22NfbmJRTkFLNCWEEG2CJ3vw7wJjPbh9F6vNjl0j8+CFEKIWj2VErfVi4Iintl+b2WYHkFk0QghRi9czolLqdqVUllIqq7CwabVkzFZHgpd58EIIcZzXM6LW+g2tdYbWOqN9+6bVkalJ8NKDF0KI43ximmRUiIklfxlNRLC/t0MRQohWwycSvJ/RQJeoYG+HIYQQrYonp0l+DCwHkpVS+5VSt3iqrYLSKv6zaAe7Css81YQQQrQ5HuvBa62neGrbJ8s/VsWz326ld1wY3dvLGa1CCAGt4CCrO8gsGiGEOJVPZESZRSOEEKfyiYwoJzoJIcSpfCIjunrwMkQjhBAuPjFNcnTvDqx99CJCA33i7QghhFv4REb0NxpoF2LydhhCCNGq+MSYxtq9xfxrYQ6lVRZvhyKEEK2GTyT4DXnHeO3nnVQ7x+KFEEL4SIKXaZJCCHEqn8iI1TKLRgghTuETGVGmSQohxKl8IiNabHb8DAqDQXk7FCGEaDV8Yprkg2OSufeCnt4OQwghWhWfSPAGgyLQYPR2GEII0ar4xBDNZ2v289y3W70dhhBCtCo+keB/3VHE52vzvB2GEEK0Kj6R4M1WOwEyB14IIU7gE1nRbLXLSU5CCHESn8iKZpskeCGEOJlPZEWDUgT5yywaIYSozSemSb499RxvhyCEEK2OT/TghRBCnMonEvy/Fubw7tLd3g5DCCFaFZ9I8N9uyue3PcXeDkMIIVoVn0jwZqudAKkkKYQQJ/BoVlRKjVVKbVVK7VBKPeypdsxWO/6S4IUQ4gQey4pKKSPwH2Ac0BeYopTq64m2LDIPXgghTuHJrDgY2KG13qW1NgOzgSs80VB4kD/hQT4x41MIIdzGk1mxE7Cv1uP9wJCTV1JK3Q7cDtC1a9cmNfTLQ6Ob9DohhPBlXh/X0Fq/obXO0FpntG/f3tvhCCGEz/Bkgs8DutR63Nn5nBBCiBbgyQT/G9BTKZWolDIBk4EvPdieEEKIWjw2Bq+1tiql7gG+BYzA21rrTZ5qTwghxIk8OvVEa/0N8I0n2xBCCFE3rx9kFUII4RmS4IUQwkdJghdCCB8lCV4IIXyU0lp7OwYXpVQhsKeJL48BitwYjidIjO7TFuKUGN2jLcQI3ouzm9a6zrNEW1WCbw6lVJbWOsPbcTREYnSfthCnxOgebSFGaJ1xyhCNEEL4KEnwQgjho3wpwb/h7QAaQWJ0n7YQp8ToHm0hRmiFcfrMGLwQQogT+VIPXgghRC2S4IUQwke1+QTfUhf2rqftLkqpRUqpzUqpTUqp+5zPRymlvldKbXf+bOd8XimlXnLGul4pNajWtm5yrr9dKXWTB2I1KqXWKqW+cj5OVEqtdMYyx1nSGaVUgPPxDufyhFrbmOZ8fqtS6mIPxBiplJqrlMpRSm1RSg1tbftSKfWA83e9USn1sVIqsDXsS6XU20qpAqXUxlrPuW3fKaXSlVIbnK95SSml3BTjs87f93ql1OdKqchay+rcR/X9z9f3e2hujLWW/VkppZVSMc7HXtmPZ0Rr3WZvOMoQ7wS6AyZgHdC3BdvvCAxy3g8DtuG4wPgzwMPO5x8G/uW8fwmwAFDAucBK5/NRwC7nz3bO++3cHOufgI+Ar5yPPwEmO++/DtzlvP8H4HXn/cnAHOf9vs79GwAkOve70c0xzgJudd43AZGtaV/iuAzlbiCo1j6c2hr2JTASGARsrPWc2/YdsMq5rnK+dpybYhwD+Dnv/6tWjHXuIxr4n6/v99DcGJ3Pd8FR+nwPEOPN/XhG78eTG/f0DRgKfFvr8TRgmhfjmQ9cBGwFOjqf6whsdd7/LzCl1vpbncunAP+t9fwJ67khrs7Aj8D5wFfOP66iWv9Yrv3o/CMe6rzv51xPnbxva6/nphgjcCRPddLzrWZfcvw6w1HOffMVcHFr2ZdAAicmT7fsO+eynFrPn7Bec2I8adkE4EPn/Tr3EfX8zzf0N+2OGIG5wAAgl+MJ3mv7sbG3tj5EU9eFvTt5IxDn1++BwEogVmt90LkoH4h13q8vXk+/jxeAvwB25+No4KjW2lpHe65YnMuPOdf3dIyJQCHwjnIMJb2plAqhFe1LrXUe8BywFziIY9+spvXtyxru2nednPc9He/NOHq1TYmxob/pZlFKXQHkaa3XnbSote5Hl7ae4FsFpVQoMA+4X2tdUnuZdnxUe20uqlLqUqBAa73aWzE0kh+Or8avaa0HAuU4hhVcWsG+bAdcgePDKB4IAcZ6K54z4e19dzpKqUcAK/Cht2OpTSkVDPw/YLq3Y2mKtp7gvX5hb6WUP47k/qHW+jPn04eUUh2dyzsCBc7n64vXk+9jOHC5UioXmI1jmOZFIFIpVXNFr9rtuWJxLo8ADns4RnD0ZvZrrVc6H8/FkfBb0768ENittS7UWluAz3Ds39a2L2u4a9/lOe97JF6l1FTgUuB65wdRU2I8TP2/h+ZIwvGBvs75P9QZWKOUimtCjB7dj3Xy5PiPp284en27cPwCag649GvB9hXwHvDCSc8/y4kHt55x3h/PiQdlVjmfj8Ix/tzOedsNRHkg3kyOH2T9lBMPSP3Bef9uTjww+Inzfj9OPOi1C/cfZF0CJDvvP+bcj61mXwJDgE1AsLPdWcAfW8u+5NQxeLftO049OHiJm2IcC2wG2p+0Xp37iAb+5+v7PTQ3xpOW5XJ8DN5r+7HR78WTG2+JG44j2dtwHFl/pIXbPg/H1971QLbzdgmO8cAfge3AD7V+uQr4jzPWDUBGrW3dDOxw3n7voXgzOZ7guzv/2HY4/zECnM8HOh/vcC7vXuv1jzhj34oHjv4DaUCWc39+4fznaFX7EngcyAE2Au87E5DX9yXwMY7jAhYc34Zucee+AzKc73kn8AonHQxvRow7cIxX1/z/vH66fUQ9//P1/R6aG+NJy3M5nuC9sh/P5CalCoQQwke19TF4IYQQ9ZAEL4QQPkoSvBBC+ChJ8EII4aMkwQshhI+SBC/aPKVUtFIq23nLV0rlOe+XKaVe9UB7yUqpn51tbFFKveF8Pk0pdYm72xOiqfxOv4oQrZvW+jCOOfQopR4DyrTWz3mwyZeA57XW851tpjifT8Mxz/kbD7YtRKNJD174LKVUpjpe//4xpdQspdQSpdQepdRVSqlnnLW5FzpLTtTU6/5FKbVaKfVtzan+J+lIraJRWusNztrjTwCTnD37SUqpEGd98VXOAmpXONuYqpSa7/wWsF0pNcPze0OcjSTBi7NJEo5aPJcDHwCLtNYpQCUw3pnkXwau0VqnA28DT9axneeBn5RSC5TjAiCRWmszjoJUc7TWaVrrOTjOxPxJaz0YGA0866yQCTAYuBpIBa5VSmV46k2Ls5cM0YizyQKttUUptQFHXZOFzuc34Kg/kgz0B753XmjHiOO09RNord9RSn2Lo47KFcAdSqkBdbQ3BkehtwedjwOBrs773zuHllBKfYaj7EVWs9+hELVIghdnk2oArbVdKWXRx+t02HH8Lyhgk9Z66Ok2pLU+gKOH/7bz8m7961hNAVdrrbee8KRSQzi1dK/UDBFuJ0M0Qhy3FWivlBoKjlLQSql+J6/kvCZozZh9HI6iXnlAKY5LN9b4FvhjzXU3lVIDay27SDmuyrBwrQAAAJFJREFUmRoEXAks9cQbEmc3SfBCODnH0a8B/qWUWoejuuGwOlYdA2x0rvMt8JDWOh9YBPStOcgK/B3wB9YrpTY5H9dYheM6AuuBeVprGZ4RbifVJIVoYc4LXGRore/xdizCt0kPXgghfJT04IUQwkdJD14IIXyUJHghhPBRkuCFEMJHSYIXQggfJQleCCF81P8HioSUSE4Fsf4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np \n",
    "\n",
    "cv_value = []\n",
    "cv_avg = []\n",
    "# we saved every 10 steps\n",
    "for i in range(0, 15000, 250):\n",
    "    variables = htf.load_variables('eds-model', ['cv_run', 'cv'], i)\n",
    "    # sum energy across particles\n",
    "    cv_avg.append(variables['cv_run'])\n",
    "    cv_value.append(variables['cv'])\n",
    "plt.plot(range(0,15000, 250), cv_avg, label='CV Running Avg', linestyle='--')\n",
    "plt.plot(range(0,15000, 250), cv_value, label='CV', color='C0')\n",
    "plt.axhline(4.0, label='EDS Set Point', color='C1')\n",
    "plt.ylabel('CV Value')\n",
    "plt.xlabel('Time Step')\n",
    "plt.legend()\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
}