{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "88a4b6e6",
   "metadata": {},
   "source": [
    "# Tutorial 4 – Predicting Energies and Nonadiabatic Couplings\n",
    "\n",
    "In Tutorial 3, we showed how to train a NN model for predicting energies, forces and nonadiabatic couplings.\n",
    "Here we demonstrate, how to use the trained model, to predict the energies and nonadiabatic couplings for a geometry of the sample molecule, which was not included in the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "30b9b5fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import ase\n",
    "from ase.db import connect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c207ca6b",
   "metadata": {},
   "source": [
    "## Molecule out of Dataset\n",
    "\n",
    "In the following we create an `Atoms` object from the symbols and positions of our target molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1d7d384f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# molecular symbols of target structure (same as database)\n",
    "symbols = 'CNHHHH'\n",
    "# positions of atoms in target structure (not included in database)\n",
    "positions = np.array(\n",
    "    [[ 0.0000,  0.0000,  0.0000 ],\n",
    "     [ 2.4321,  0.0000,  0.0000 ],\n",
    "     [-1.0111,  1.7951,  0.0000 ],\n",
    "     [ 3.4373,  1.6202, -0.2566 ],\n",
    "     [ 3.4373, -1.6202,  0.2566 ],\n",
    "     [-1.0111, -1.7951,  0.0000 ]]\n",
    ")\n",
    "# create ase Atoms object\n",
    "target_mol = ase.Atoms(symbols=symbols, positions=positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67f20190",
   "metadata": {},
   "source": [
    "## Define the calculator\n",
    "\n",
    "In this step, we define a calculator for predicting for the target geometry the properties. \n",
    "Therefore, we have to specify the path to the trained model (see Tutorial 3) as well as the neighbor list for building the representation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ff3065ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "from schnetpack.transform import MatScipyNeighborList\n",
    "from spainn.interface import NacCalculator\n",
    "\n",
    "# set up calculator for predictions\n",
    "calc = NacCalculator(\n",
    "    model_file=os.path.join(os.getcwd(), 'train', 'best_model_E_F_C'), \n",
    "    neighbor_list=MatScipyNeighborList(cutoff=10.0)\n",
    ")\n",
    "# apply calculator to target molecule\n",
    "target_mol.calc = calc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c59c277",
   "metadata": {},
   "source": [
    "In the next step, we get the predictions from our model, by calling the `get_properties` function for the set of properties, we are interested in (here energies and nonadiabatic couplings).\n",
    "Note, the `NacCalculator` always returns nonadiabatic couplings $\\mathbf{C}_{ij}$ rather than the smoothed nonadiabatic couplings $\\mathbf{C}_{ij}^s$ employed in the training. \n",
    "Nevertheless, the predictions, need to be requested with the same property key as used in the training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4531486c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'energy': array([-94.67488565, -94.42807255, -94.33222212]),\n",
       " 'smooth_nacs': array([[[-5.55834436e-02, -1.67016947e-17, -7.63236560e-19],\n",
       "         [-1.71284160e-01,  4.90413892e-17, -3.03844639e-19],\n",
       "         [-9.52767548e-01,  1.31373348e-16,  2.18561979e-17]],\n",
       " \n",
       "        [[-9.73515189e-02, -4.93364229e-17,  2.94756393e-18],\n",
       "         [-2.03955834e-01,  1.57202331e-18, -3.02399188e-18],\n",
       "         [ 7.18859946e-01, -4.18493612e-17,  5.87670367e-18]],\n",
       " \n",
       "        [[ 3.81299868e-02, -4.82265658e-02, -1.96809142e-04],\n",
       "         [-1.94629845e-02,  3.90448386e-02,  1.13046205e-03],\n",
       "         [ 1.16888676e-01, -9.81919929e-02,  4.69116636e-03]],\n",
       " \n",
       "        [[-3.24237656e-02, -1.53159303e-02,  2.94397172e-03],\n",
       "         [ 1.06732541e-01,  4.91115399e-02, -8.72029081e-03],\n",
       "         [ 2.22506380e-02, -3.17968002e-02,  3.09087171e-04]],\n",
       " \n",
       "        [[-3.24237656e-02,  1.53159303e-02, -2.94397172e-03],\n",
       "         [ 1.06732541e-01, -4.91115399e-02,  8.72029081e-03],\n",
       "         [ 2.22506380e-02,  3.17968002e-02, -3.09087171e-04]],\n",
       " \n",
       "        [[ 3.81299868e-02,  4.82265658e-02,  1.96809142e-04],\n",
       "         [-1.94629845e-02, -3.90448386e-02, -1.13046205e-03],\n",
       "         [ 1.16888676e-01,  9.81919929e-02, -4.69116636e-03]]])}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# predict energies and nacs from NN model\n",
    "pred = target_mol.get_properties(['energy', 'smooth_nacs'])\n",
    "pred"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ae3f05c",
   "metadata": {},
   "source": [
    "## Mean absolute error for energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8e336c54",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE: 0.03294638125980498\n"
     ]
    }
   ],
   "source": [
    "# reference energies from (6,4)-CASSCF/cc-pVDZ calculations\n",
    "ref = {'energy': np.array([-94.69603647, -94.43388056, -94.40410243])}\n",
    "\n",
    "mae = np.mean(\n",
    "    np.abs(np.asarray(pred['energy']).flatten() \n",
    "           - np.asarray(ref['energy']).flatten()\n",
    "          ))\n",
    "\n",
    "print('MAE: '+str(mae))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ec5aed3",
   "metadata": {},
   "source": [
    "## Visualization of prediction quality for energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4a46cd29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAHeCAYAAAC7RCO/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWiUlEQVR4nO3dd1iUV6I/8O8w9DaDdATpINGgGNdCikbEqFGzm6ybaHTFAtnkejUkMSvqWtibawzXkkTWcF0xkJiYTUwxZiU/ozHRqICFrEFDsUtTurQBmff3B5dZR9pUpvD9PM88j/O+55w574zw5bzvmfOKBEEQQEREREbLwtAdICIiot4xrImIiIwcw5qIiMjIMayJiIiMHMOaiIjIyDGsiYiIjBzDmoiIyMgxrImIiIwcw5qIiMjIMayJiIiMHMOaSI8mTpwIkUgEkUiEsLCwPstPmDBBUX748OF9lt+5c6ei/J/+9Kc+y8fFxSnK9/W4evWqKoeolpMnT+KZZ56Bp6cnbG1tERgYiJdeegklJSVatVtSUoIXX3wRAQEBsLW1haenJ5555hmcOnWqxzrHjx/Hq6++isceewz+/v5wcHCAnZ0dgoODsWDBAuTk5PRYt6CgAFu2bMETTzwBb29vWFtbQyKR4De/+Q2Sk5NRU1PTa3/r6uqQnJyMkSNHwtHREY6OjoiMjMRf//pXNDU1afw+kBkTiEhvJkyYIABQPI4dO9Zj2eLiYkEkEinKDhs2rM/2x40bpygvkUiEpqamXssvWLBAACB4eHgIDz/8cK+PsrIytY+3Nzt37hQsLCwEAIKbm5vw0EMPCc7OzgIAwcXFRcjLy9Oo3bNnzwpSqVQAINjb2wujRo0SfHx8BACCWCwW0tPTu623ePFiAYBgYWEheHt7Cw899JAQHh4u2NjYCAAEkUgkbNq0qUu94uJipc/Ux8dHGD16tDB48GDFNm9vb+Ff//pXt69bXFwsBAYGKl7jgQceEKKiogRra2sBgBARESHcunVLo/eCzBfDmkiPOsM6IiJCACAsXry4x7KrV69WKttXWF+8eFHxC18ikQgAhA8//LDXOp1hvWDBAk0OR2P/+te/BLFYLAAQ/vznPwutra2CIAhCY2Oj8PzzzwsAhMDAQKGlpUWtdpubm4UhQ4YIAITJkycLVVVVgiAIglwuF9555x0BgGBpaSnk5+d3qfvFF18I+/fvF+rq6pS2V1dXC3/6058U721OTo7S/qKiIsHDw0NITk4WLl26pLTv+PHjgr+/vwBACAoK6nI87e3tQmRkpABACAsLEy5evKjYV1FRITz++OOKYyG6F8OaSI86w3rt2rWCvb294OTkJDQ2NnYp197eLvj5+QlisVj47//+b5XCesWKFQIAYcKECcLy5csFAMKkSZN6rWOosH7mmWcEAEJ0dHSXfS0tLYqR5o4dO9Rq99133xUACE5OTkJlZWWX/XPmzBEACM8++6xa7crlcuGBBx5Q/HFxr+bmZqGhoaHHuj/99JNihP3VV18p7Ttw4ECvZ1lKSkoEBwcHAYBw+PBhtfpM5o3XrIn6gZOTE37/+9/jzp072LdvX5f9hw4dwo0bNzB16lR4e3v32d7du3fxwQcfAAD++Mc/YsGCBQCA77//Xi/XmrXR2NiIb775BgDw4osvdtlvY2ODuLg4AMDevXvVavsf//gHAOAPf/gDXF1du+zvvI6/f/9+ta4Fi0QiDB06VNH/e9na2sLBwaHHutHR0ZBIJACAixcvKu07duwYAGDw4MF45JFHutT18fHBY489BgD46KOPVO4vmT+GNVE/WbhwIQDg/fff77Jv9+7dSmX6cvDgQZSXl8Pe3h6zZ89GVFQUIiMjIQhCt+3rwvr16yESiRAQEKBWvXPnzqGlpQUAFEF0vwkTJgAAsrOzIZfLVWq3vb1dMQmsp3bHjRsHa2trNDc3Iy8vT+U+NzU1ITc3FwDwm9/8RuV6nf1qa2sDANjb2yvtq6qqAtAR1j3x8/MD0DEBjqgTw5qon0yYMAFBQUH4/vvvce3aNcX2mpoafPnll3B1dcXMmTNVais9PR0A8Nvf/hZOTk4AoBhdv//++xAEQce911xhYSEAwNraWhFE9wsODgYAtLS0KL03vbl27RpkMhkAICQkpNsy975mQUFBn23W1NTghx9+wPTp03Hjxg089thjeP7551XqT6cvvvhCMYp/9NFHlfZJpVIA6HX2+40bNwAAxcXFuHv3rlqvTeaLYU3UT0QiERYsWABBEJCRkaHY/tFHH0Emk+H555+HtbV1n+3cvn1bcVq5M6ABYN68ebC0tMS1a9dw+PDhXtvIyMjo9WtbEydO7FJHKpXC398fvr6+Kh5xh+rqagCAi4sLRCJRt2UGDRqk+HdfX3u6v9376/fUdk/t5uXlKY570KBBmDhxIi5evIhNmzbh0KFDEIvFKvWns0+vvvoqAGD69OkYOXKk0v5x48YB6AjrEydOdKlfXl6uOFXe3t6O+vp6lV+bzBvDmqgfxcXFwcLCAhkZGYrRr7qnwD/44AO0tbXBx8cHkydPVmz38PDA1KlTldrsiYeHBx5++OEeHw8++GCXOi+//DKuXr2q9unZ5uZmAOj1DxFbW1vFv1W9ttzZrqpt99Suo6Oj4riDg4NhZWWFW7du4ZNPPsHp06dV6gsAtLW1Yfbs2bh+/Trc3d2xc+fOLmVmzpypuIwQFxeHn3/+WbGvpKQEf/jDH9DQ0KDYxu9cUydLQ3eAaCAZMmQIJk2ahO+++w4//vgjBg0ahDNnziAqKqrLKKwnnUE8b948WFgo/70dFxeHAwcO4PPPP0dtba3itOv9pk2bprdr2/ezs7MDALS2tvZYpvOaNtD1Om9f7aradk/thoSEKP0BUltbi02bNmHTpk14/PHHcezYMYwZM6bXvrS3t2PevHk4cuQInJyc8PXXX8PHx6dLOWtra+zbtw9PPPEEioqKEBUVBX9/f9jZ2aGwsBByuRxxcXGKz8bZ2bnX16WBgyNron5270SzzmvPqo6qc3Jy8MsvvwDomAV+v5kzZ2LQoEFoaWnBxx9/rKMea8fFxQVAx2nonq6l33tKu7O8qu0C/5641VvbqrYrlUqxceNGxMfHo7W1FWvWrOm1vFwux6JFi/CPf/wDDg4O+OabbzB27Ngey48aNQo///wzEhMTERoairKyMpSWliImJgaHDh3CjBkzAAAODg6K+QhEDGuifva73/0OEokEn376KT788ENYW1tj7ty5KtW99/T28OHDu1xrtrGxUYRTX6fC+0t4eDiAjtHv9evXuy1z6dIlAB2nrP39/VVqNyAgADY2NgA6JmN1p7W1VTFhq7Mfquqc7NfbqXBBEPDCCy8gMzMT9vb2OHDgQJdJZd3x8fHBli1bUFBQgJaWFtTW1uLbb79FTEyMYhb66NGje7zGTwMPw5qon9nZ2eG5555DY2MjKisrMWvWrG6/I3y/lpYWxfeQXVxc4Onp2eMDAHJzcxWjcEOKiopSXDf+8ccfuy3zww8/AADGjBnT5dR+T8RiseJrVT21e+rUKbS2tsLOzk7lywydOmdit7e391jmP/7jP/D3v/8ddnZ22L9/f7cT89Qhl8vx6aefAuiY6U/UiWFNZAAvvPACYmJiEBMTg//4j/9Qqc6+fftQW1sLKysrFBYWory8vMdH5wSxztPshuTg4IDp06cDANLS0rrsl8lkimu0zz77rFpt/+EPfwAAfPrpp92eCn/vvfcAADNmzFD5Wninzz77DEDHHxvdWbZsGXbs2AFbW1t89dVXiImJUav97vztb3/D5cuX4erqqlgohggAb+RBpE+dy42mpKSoXGf37t3dLjcaExMjABCefvrpPtvYvHmzAEBwd3dXrMMtCNotN7p161bB399fePjhh9Wum5eX1+Pa4PPmzRMACP7+/t2uDf7qq68K/v7+3S4Z2tTUJPj5+fW5Nvj58+eV6t28eVP4z//8T+Hs2bNd2rx165bw8ssvK5YF3b9/f5cynUu92traCt9++61a78WhQ4eEw4cPC3K5XLGtublZSElJUbxHe/fuVatNMn+cDU5kAq5du4YjR44AUG0y2rx58/DnP/8Zt2/fxoEDB/C73/1Oaf/Bgwe7Xe7yXhs3blS6/lpbW6vygiX3GzFiBFJTU/HSSy9h06ZN2LVrF/z9/VFUVIT6+npIpVJ88cUXimvQ96qsrMS1a9e6XTnNzs4On3/+OSZPnozvvvsOfn5+GDp0KMrLy1FaWgoLCwu89957XW432tbWhnfffRfvvvsunJ2dERgYCDs7O1RXV+PSpUtob2+HlZUVNm3a1GWhmpMnTyIlJQVAx2zt5ORkJCcnd3vcixYtwqJFi5S2HT9+HBs2bICjoyMCAgJgaWmJwsJCNDU1wdraGn/729/UPsNA5o9hTWQCOlcl8/T0VHyXujceHh548skn8dVXXyE9Pb1LWN+6dQu3bt3qtY3eZlhr4oUXXsCDDz6IlJQU/PTTTzh//jy8vLwwd+5crF69Wu3FVjqNHj0a58+fxxtvvIGDBw/il19+gUQiwW9/+1u8/vrrGD9+fJc6Xl5e2LlzJ44ePYpz587hxo0bqKurg6OjI0aMGIGJEyciISGh20lpnaumAX2/j/d+D77T1KlTceXKFZw6dQpXr15Fe3s7fH19MWXKFCxfvhyhoaEavQ9k3kSCYETrEhIREVEXnGBGRERk5BjWRERERo5hTUREZOQY1kREREaOYU1ERGTkGNZERERGjt+zNgC5XI7S0lI4OTlxoX4iogFMEATcuXMHPj4+va6Lz7A2gNLSUvj5+Rm6G0REZCRu3LjR68JADGsD6LxH7Y0bN3hzeSIiM1VVVdXjbWE7NTY24sknn+zz3uUMawPoPPXt7OzMsCYiMkO3b99GdXU1HB0dVSrf1yVRTjAjIiLSodu3b/c5olYXw5qIiEhH9BHUAE+DExER6cStW7dw48YNvbTNkTUREZGW9BnUAMOaiIhIK/oOaoCnwYmIiDRWUVGBmzdvalzf399fpXIMayIi6pYgCKhpakOj7C4cbCzhYm/FVRfvoU1Qi0QiBAYGQiwWq1SeYU1ERErqmtuw78xNZJy4imvVTYrt/oPssSA6AM885AuJnZUBe2h4ughqFxcX1NfXq1ZHEARBo1cjjdXX10MikaCuro6LohCRUfmh8DZe/PAMmlvbAQD3BkTnmNrOWowd8x7ChDD3fu+fMSgvL0dJSYlGde8NakD1POAEMyIiAtAR1At356C5rR0ClIMa//dcANDc1o6Fu3PwQ+Ht/u+kgekyqNXBsCYiItQ1t+HFD890BHIf51sFoSO0X/zwDOqa2/qje0ZB26AOCgrSKKgBhjUREQHYd+Ymmlvb+wzqToIANLe24/Ozms+ENiVlZWVaB7VUKtX49RnWREQDnCAIyDhxVaO67/90FeY+9amsrAylpaUa1dVFUAMMayKiAa+mqQ3Xqpu6XKPuiwDgWnUTapvM91S4tkEdHBysdVADDGsiogGvUXZXq/oNWtY3VqWlpVoHtUQi0UlfGNZERAOcg412S244alnfGJWWlqKsrEyjuroOaoBhTUQ04LnYW8F/kD3UXZtMhI6FUqT25rVAirEFNcCwJiIa8EQiERZEB2hUN+7hALNaglTboA4JCdF5UAMMayIiAvDMQ76wsxZD1dy1EHWsZPb0KF/9dqwflZSUaB3U+lqVkmFNRESQ2Flhx7yHIAL6DOzO/e/Ne8hs1ggvKSlBeXm5RnUtLCz0GtQAw5qIiP7PhDB37F44BnZW4o7Qvm9/5zY7KzHeXzgGj5nJ2uDaBnVwcLDe7/NgflP4iIhIYxPC3HEyKQafn72J939SvuvWkEH2iHu4465bzrbmMaK+efMmKioqNKrbOaJ2cnLSca+64l23DIB33SIiUyAIAmqb2tAguwtHG0tIzex+1sYQ1KrmAUfWRETULZFIBBcHa7g4WBu6KzpnDEGtDoY1ERENKDdu3MCtW7c0qmuIoAYY1kRENIBoG9ShoaFwdHTUca9UeO1+f0UiIiIDMNWgBkw0rM+ePYvY2FhIpVK4uroiISEBDQ0N3ZatqqqCr68vRCIRamtre233hRdeQHBwMOzs7ODu7o6nnnoKv/76q1KZgICO1XrufaxcuVJXh0ZERHpw/fp1kw1qwATDurS0FJMnT0ZISAiys7ORlZWF/Px8xMXFdVt+8eLFiIyMVKnthx56CLt378bFixfx7bffQhAETJkyBe3t7UrlkpOTUVZWpnisWbNG28MiIiI9uX79Om7fvq1RXWMIasAEr1kfOHAAVlZWSE1NhYVFx98aqampiIqKQnFxMUJCQhRld+zYgdraWqxduxYHDx7ss+2EhATFvwMCAvBf//VfGDFiBK5evYrg4GDFPicnJ3h5eenwqIiISB+0CWqxWIyQkBCDBzVggiNrmUwGa2trRVADgJ2dHQDg+PHjim0XLlxAcnIyMjMzlcqqqrGxEbt370ZgYCD8/PyU9m3atAmurq4YOXIk3njjDbS2tvbZ5/r6eqUHERHpl7ZBbQwj6k4mF9aTJk1CeXk5UlJS0NraipqaGqxatQoAFAuwy2QyzJkzBykpKRgyZIha7f/tb3+Do6MjHB0dkZWVhUOHDsHa+t/fMVy+fDn27t2L77//HkuXLsW2bdvw0ksv9drmxo0bIZFIFI/7w5+IiHTr2rVrWge1g4ODjnulOaMJ6/Xr13eZuHX/4/Tp0xg2bBgyMjKwefNm2Nvbw8vLC0FBQfD09IRYLAYAJCUlISIiAvPmzVO7H88//zzOnTuHH374AaGhofjDH/6AlpYWxf7ExERMmDABkZGRWLJkCd577z3s2rULVVVVPbaZlJSEuro6xePGjRvqv0FERKSSa9euobKyUqO6xhjUgBEtN1pZWdnnmxsQEABbW1vF84qKCjg4OEAkEsHZ2Rl79+7F7NmzMXLkSJw/f16xLJ4gCJDL5RCLxVi9ejU2bNigUp9aW1vh4uKCv//975gzZ063ZUpKSuDr64tTp05h7NixKrXL5UaJiPTD1ILa5JYbdXNzg5ubm1p1PD09AQDp6emwtbVFbGwsAGDfvn1obm5WlMvNzcWiRYtw7NgxpYliqhAEATKZrMf9586dAwB4e3ur1S4REemOIAi4du1ar2c5eyMWixEWFgZ7e3sd90w3jCas1bF9+3ZER0fD0dERhw4dwooVK/Dmm29CKpUCQJdA7vwrKyIiQlGmpKQEMTExyMzMxJgxY3D58mV88sknmDJlCtzd3VFSUoJNmzbBzs4O06dPBwCcPHkSp06dwuOPPw6JRILc3FwkJiZi1qxZal8bJyIi3TD3oAZMNKxzcnKwbt06NDQ0YOjQoUhLS8P8+fPVaqOtrQ0FBQVoauq4/ZutrS2OHTuGbdu2oaamBp6ennjsscdw4sQJeHh4AABsbGzwySefYMOGDZDJZPD390d8fDxef/11nR8jERH1bSAENWBE16wHEl6zJiLSniAIuHr1KqqrqzWqb2lpidDQUIMGtap5YDSzwYmIiFRlDkGtDpM8DU5ERAOXLoI6LCxMsaCWKeDImoiITMZADGqAI2siIjIRgiDgypUrqKmp0ai+qQY1wJE1ERGZgIEc1ADDmoiIjNxAD2qAp8GJiMiI6SKow8PDlZaqNkUMayIiMkqCIODy5cuora3VqL6VlRXCwsJMPqgBngYnIiIjxKBWxrAmIiKjwqDuiqfBiYjIaAiCgEuXLqGurk6j+uYY1ABH1kREZCQY1D3jyJqIiAxOF0EdHh4OGxsbHffMODCsiYjIoORyOS5fvqxxUFtbWyMsLMxsgxpgWBMRkQHJ5XJcunQJ9fX1GtUfCEEN8Jo1EREZCINadQxrIiLqdwxq9TCsiYioXzGo1cdr1kRE1G/kcjmKi4tx584djepbW1sjPDwc1tbWOu6ZcePImoiI+gWDWnMcWRMRkd5pG9Q2NjYICwsbkEENcGRNRER6xqDWHkfWRESkN3K5HEVFRWhoaNCoPoO6A0fWRESkFwxq3eHImoiIdK69vR3FxcVaBXV4eDisrKx03DPTxLAmIiKd0jaobW1tERYWxqC+B0+DExGRzjCo9YMjayIi0on29nYUFRWhsbFRo/oM6p5xZE1ERFpjUOsXw5qIiLTCoNY/hjUREWmMQd0/eM2aiIg00t7ejsLCQjQ1NWlU387ODmFhYbC0ZBT1hSNrIiJSG4O6f/FdIiIitdy9exdFRUUM6n7EkTUREamMQW0YDGsiIlIJg9pw+I4REVGf7t69i8LCQjQ3N2tU397eHqGhoQxqDZnkyPrs2bOIjY2FVCqFq6srEhISelzarqqqCr6+vhCJRKitrVWpfUEQMG3aNIhEInz55ZdK+2pqajB//nxIJBJIJBLMnz9f5XaJiEwRg9rwTC6sS0tLMXnyZISEhCA7OxtZWVnIz89HXFxct+UXL16MyMhItV5j27ZtEIlE3e6bO3cu8vLykJWVhaysLOTl5WH+/PnqHgYRkUlgUBsHk3v3Dhw4ACsrK6SmpsLCouNvjdTUVERFRaG4uBghISGKsjt27EBtbS3Wrl2LgwcPqtT+zz//jC1btiA3Nxfe3t5K+y5evIisrCycOnUKY8eOBQDs3LkT48ePR0FBAcLDw3V0lEREhqeLoA4LC4NYLNZxzwYekxtZy2QyWFtbK4Ia6Ji0AADHjx9XbLtw4QKSk5ORmZmpVLY3TU1NmDNnDrZv3w4vL68u+0+ePAmJRKIIagAYN24cJBIJTpw4oekhEREZHW2D2sHBgUGtQyYX1pMmTUJ5eTlSUlLQ2tqKmpoarFq1CgBQVlYGoCPQ58yZg5SUFAwZMkTlthMTExEdHY2nnnqq2/3l5eXw8PDost3DwwPl5eU9tiuTyVBfX6/0ICIyVm1tbSgoKNAqqENDQxnUOmQ0Yb1+/XqIRKJeH6dPn8awYcOQkZGBzZs3w97eHl5eXggKCoKnp6fiP0ZSUhIiIiIwb948lV9///79OHLkCLZt29Zrue6uZQuC0OM1bgDYuHGjYkKaRCKBn5+fyv0iIupPbW1tKCwsREtLi0b1GdT6IRIEQTB0JwCgsrISlZWVvZYJCAiAra2t4nlFRQUcHBwgEong7OyMvXv3Yvbs2Rg5ciTOnz+vCFBBECCXyyEWi7F69Wps2LChS9svv/wy3nnnHaVT5u3t7bCwsMCjjz6Ko0ePIj09Ha+88kqX2d9SqRRbt27FwoULu+23TCaDTCZTPK+vr4efnx/q6urg7Ozc53tDRNQfGNT9r76+HhKJpM88MJoJZm5ubnBzc1OrjqenJwAgPT0dtra2iI2NBQDs27dP6fRNbm4uFi1ahGPHjiE4OLjbtlauXIklS5YobXvwwQexdetWzJw5EwAwfvx41NXVIScnB2PGjAEAZGdno66uDtHR0T3208bGBjY2NmodGxFRf2JQGzejCWt1bN++HdHR0XB0dMShQ4ewYsUKvPnmm5BKpQDQJZA7R+wRERGKMiUlJYiJiUFmZibGjBkDLy+vbieVDRkyBIGBgYr6U6dORXx8PNLS0gAACQkJmDFjBmeCE5HJ0jaoHR0dERISwqDWI5MM65ycHKxbtw4NDQ0YOnQo0tLS1P6uc+cECnWXzduzZw+WLVuGKVOmAABmzZqF7du3q9UGEZGx6PxdeO+lOnUwqPuH0VyzHkhUvUZBRKRPDGrDM7lr1kRE1H9aW1tRWFioVVCHhoaqvI4FaYdhTUQ0wGgb1E5OTggJCWFQ9yO+00REAwiD2jTx3SYiGiAY1KaLp8GJiAaA1tZWFBQUoLW1VaP6DGrDYlgTEZk5bYPa2dkZwcHBDGoDYlgTEZkxmUyGwsJCBrWJ47tPRGSmGNTmg58AEZEZYlCbF34KRERmhkFtfnjNmojIjMhkMhQUFKCtrU2j+hKJBEFBQQxqI8OwJiIyE7oI6uDgYIhEIh33jLTFsCYiMgMtLS0oLCxkUJsphjURkQkQBAE1TW1olN2Fg40lXOytFMGqbVBLpVIEBQUxqI0Yw5qIyIjVNbdh35mbyDhxFdeqmxTb/QfZY0F0AJ4c5oaKG1cY1GaO97M2AN7PmohU8UPhbbz44Rk0t7YDAO79ZS0CIL/bCvGdcqx8IhQP+buo3T6D2vBUzQNO9yMiMkI/FN7Gwt05aG5rhwDloAY6grqtugQtra3Y8HU+zlyrUat9BrVpYVgTERmZuuY2vPjhmY6Q7ubcp7xNhrbqmxDk7RCEjiDfePAiGlruqtS+i4sLg9rEMKyJiIzMvjM30dza3mNQ360pgSCXK7YJAiBrk+PIrxV9tu3i4oLAwEAGtYlhWBMRGRFBEJBx4mq3+7oL6nt9/XMZepuGxKA2XQxrIiIjUtPUhmvVTV2vUfcR1AKAsvoW3OnhVDiD2rQxrImIjEijrGvY9hXU92r6v5nj9xo0aBCD2sTxe9ZEREbEwUb517K8rQV3a0pVCmoAsLcWKz0fNGgQAgICGNQmjiNrIiIj4mJvBf9B9h3fo1YjqEUAvJ1t4WT777BnUJsPhjURkRERiURYEB3QEdTVqo+oAWDmCG9FMLu6ujKozQjDmojIyEyLGATxnXIAqgW1hQiwsbLApKGeADqC2t/fn0FtRhjWRERGpKmpCRU3rmDl1HCIAPSVt537V02LgKOtJYPaTDGsiYiMRGNjIwoLC9He3o6H/F2wbuYw2FhadIT2fWU7t9lYWmD9zGEY5e/CoDZjnA1ORGQEGhsbUVRUhPb2f3/16iF/F7wfNwZHfq3A1z+Xoay+RbHPy9kWM0d4Y1KEJxxtLOHm5gZ/f39DdJ36AcOaiMjAugvqTo62lpg1cjBmjvDBnZa7aGpth721GE62looRNIPa/DGsiYgMqLegvpdIJIKznRWc7ayUtjOoBwaGNRGRgaga1D1xd3fHkCFDdNwrMkYMayIiA2hoaEBRURHkanyP+l4M6oGFs8GJiPoZg5rUxZE1EVE/0jaoPTw84Ofnp+NekbHjyJqIqJ8wqElTHFkTEfWDO3fuoLi4mEFNGuHImohIzxjUpC2TDOuzZ88iNjYWUqkUrq6uSEhIQENDQ7dlq6qq4OvrC5FIhNraWpXaFwQB06ZNg0gkwpdffqm0r/MuNvc+Vq5cqeUREZG50jaoPT09GdRkemFdWlqKyZMnIyQkBNnZ2cjKykJ+fj7i4uK6Lb948WJERkaq9Rrbtm3rdW3d5ORklJWVKR5r1qxRq30iGhh0EdS+vr467hWZIpO7Zn3gwAFYWVkhNTUVFhYdf2ukpqYiKioKxcXFCAkJUZTdsWMHamtrsXbtWhw8eFCl9n/++Wds2bIFubm58Pb27raMk5MTvLy8tD8YIjJb9fX1uHTpEoOadMLkRtYymQzW1taKoAYAOzs7AMDx48cV2y5cuIDk5GRkZmYqle1NU1MT5syZg+3bt/caxps2bYKrqytGjhyJN954A62trRoeDRGZI22D2svLi0FNSkwurCdNmoTy8nKkpKSgtbUVNTU1WLVqFQCgrKwMQEegz5kzBykpKWotHJCYmIjo6Gg89dRTPZZZvnw59u7di++//x5Lly7Ftm3b8NJLL/XarkwmQ319vdKDiMyTLoJ68ODBOu4VmTqjCev169d3mbh1/+P06dMYNmwYMjIysHnzZtjb28PLywtBQUHw9PSEWCwGACQlJSEiIgLz5s1T+fX379+PI0eOYNu2bb2WS0xMxIQJExAZGYklS5bgvffew65du1BVVdVjnY0bN0IikSgenCxCZJ7q6+u1ukbNoKaeiARBEAzdCQCorKxEZWVlr2UCAgJga2ureF5RUQEHB4eOu9E4O2Pv3r2YPXs2Ro4cifPnzysmiQmCALlcDrFYjNWrV2PDhg1d2n755ZfxzjvvKJ0yb29vh4WFBR599FEcPXq02z6VlJTA19cXp06dwtixY7stI5PJIJPJFM/r6+vh5+eHuro6ODs793rMRGQaOoNa01+p3t7e8PHx0XGvyNjV19dDIpH0mQdGM8HMzc0Nbm5uatXx9PQEAKSnp8PW1haxsbEAgH379qG5uVlRLjc3F4sWLcKxY8cQHBzcbVsrV67EkiVLlLY9+OCD2Lp1K2bOnNljH86dOwcAPU5GAwAbGxvY2NiodlBEZHLq6upw6dIlBjXpjdGEtTq2b9+O6OhoODo64tChQ1ixYgXefPNNSKVSAOgSyJ0j9oiICEWZkpISxMTEIDMzE2PGjIGXl1e3k8qGDBmCwMBAAMDJkydx6tQpPP7445BIJMjNzUViYiJmzZrFRfWJBigGNfUHkwzrnJwcrFu3Dg0NDRg6dCjS0tIwf/58tdpoa2tDQUEBmpqaVK5jY2ODTz75BBs2bIBMJoO/vz/i4+Px+uuvq3sIRGQGtA1qHx+fXs/KEXUymmvWA4mq1yiIyHgxqEkXTO6aNRGRqaitrcXly5cZ1NRvGNZERGrQNqgHDx7MFRBJbUbzPWsiImPHoCZD4ciaiEgFNTU1uHLlCoOaDIJhTUTUB22D2tfXV7EuBJEmeBqciKgXDGoyBgxrIqIeMKjJWPA0OBFRN6qrq3HlyhWN6zOoSZcY1kRE99E2qP38/ODh4aHDHtFAx9PgRET3YFCTMWJYExH9HwY1GSueBiciAlBVVYWrV69qXH/IkCFwd3fXXYeI7sGRNRENeAxqMnYMayIa0BjUZAp4GpyIBqzKykpcu3ZN4/r+/v5wc3PTYY+IuseRNRENSAxqMiUMayIacBjUZGoY1kQ0oDCoyRRpdc36woULuH79OlpbW5W2z5o1S6tOERHpw+3bt3H9+nWN6wcEBMDV1VWHPSJSjUZhffnyZfzud7/D+fPnIRKJFIvci0QiAEB7e7vuekhEpAMMajJlGp0GX758OQIDA1FRUQF7e3vk5+fjxx9/xOjRo3H06FEdd5GISDsMajJ1Go2sT548iSNHjsDd3R0WFhawsLDAI488go0bN2LZsmU4d+6crvtJRKSRW7du4caNGxrXDwwMxKBBg3TYIyL1aTSybm9vh6OjIwDAzc0NpaWlADomXhQUFOiud0REWmBQk7nQaGQ9fPhw/Otf/0JQUBDGjh2Lt956C9bW1vjf//1fBAUF6bqPRERqY1CTOdEorNesWYPGxkYAwH/9139hxowZePTRR+Hq6opPPvlEpx0kIlJXRUUFbt68qXF9BjUZG5HQOZVbS9XV1XBxcVHMCKee1dfXQyKRoK6uDs7OzobuDpFZ0SaoRSIRAgMD4eLiouNeEXVP1TxQa2RdX1/fc0OWlrhz5w4AMICIyCAY1GSu1AprqVTa68hZEASIRCJ+z5qI+l15eTlKSko0qsugJmOnVlh///33in8LgoDp06fj73//OwYPHqzzjhERqYpBTeZOrbCeMGGC0nOxWIxx48ZxBjgRGYy2QR0UFASpVKrbThHpGO9nTUQmq6ysTLHOg7oY1GRKeNctIjJJDGoaSLQeWfOrWkTU37QN6uDgYEgkEh33ikh/1Arrp59+Wul5S0sL/vSnP8HBwUFp++eff659z4iIulFaWoqysjKN6jKoyVSpFdb3/wefN2+eTjtDRNQbBjUNVGqF9e7du/XVDyKiXjGoaSDjbHAiMnraBnVISAhXViSTxrAmIqNWUlKC8vJyjeoyqMlcmORXt86ePYvY2FhIpVK4uroiISEBDQ0N3ZatqqqCr68vRCIRamtre2134sSJEIlESo/nnntOqUxNTQ3mz58PiUQCiUSC+fPn99kuEWlGm6C2sLBgUJPZMLmwLi0txeTJkxESEoLs7GxkZWUhPz8fcXFx3ZZfvHgxIiMjVW4/Pj4eZWVlikdaWprS/rlz5yIvLw9ZWVnIyspCXl4e5s+fr80hEVE3tA3q4OBgBjWZDbVOg69atQq//e1vMWbMGH31p08HDhyAlZUVUlNTYWHR8bdGamoqoqKiUFxcjJCQEEXZHTt2oLa2FmvXrsXBgwdVat/e3h5eXl7d7rt48SKysrJw6tQpjB07FgCwc+dOjB8/HgUFBQgPD9fy6IgIAG7evImKigqN6naOqJ2cnHTcKyLDUWtkXVZWhhkzZsDb2xsJCQn45ptvIJPJ9NW3bslkMlhbWyuCGgDs7OwAAMePH1dsu3DhApKTk5GZmalUti979uyBm5sbhg0bhtdee01x208AOHnyJCQSiSKoAWDcuHGQSCQ4ceJEr32ur69XehBR9xjURF2pFda7d+9GRUUF/vGPf0AqleLVV1+Fm5sbnn76abz//vuorKzUVz8VJk2ahPLycqSkpKC1tRU1NTVYtWoVAChmi8pkMsyZMwcpKSkYMmSIym0///zz+Pjjj3H06FH85S9/wb59+5QWgikvL4eHh0eXeh4eHr2ertu4caPiGrdEIoGfn5/KfSIaSBjURN1T+5q1SCTCo48+irfeegu//vorcnJyMG7cOOzcuRODBw/GY489hv/5n/9R+y4469ev7zK56/7H6dOnMWzYMGRkZGDz5s2KU9ZBQUHw9PSEWCwGACQlJSEiIkLtRVvi4+MxefJkDB8+HM899xw+++wzfPfddzh79qzS8d+v8z7ePUlKSkJdXZ3icePGDbX6RTQQ3Lhxg0FN1AORIAiCrhq7ffs29u/fj/379+PRRx/Fa6+9pnLdysrKPkfmAQEBsLW1VTyvqKiAg4MDRCIRnJ2dsXfvXsyePRsjR47E+fPnFQEqCALkcjnEYjFWr16NDRs2qNQnQRBgY2ODDz74AM8++yzS09PxyiuvdJn9LZVKsXXrVixcuFClduvr6yGRSFBXV8cJMEToCOpbt25pVNfCwgKhoaFwdHTUca+I9E/VPNDp96zd3d2xePFiLF68WO26bm5ucHNzU6uOp6cnACA9PR22traIjY0FAOzbtw/Nzc2Kcrm5uVi0aBGOHTuG4OBgldvPz89HW1sbvL29AQDjx49HXV0dcnJyFJPssrOzUVdXh+joaLX6TkQdGNREfTPJRVG2b9+O6OhoODo64tChQ1ixYgXefPNNxe3u7g/kzhF7RESEokxJSQliYmKQmZmJMWPG4NKlS9izZw+mT58ONzc3XLhwAa+++iqioqLw8MMPK+pPnToV8fHxiq90JSQkYMaMGZwJTqSB69ev4/bt2xrVZVDTQGKSYZ2Tk4N169ahoaEBQ4cORVpamtrfdW5ra0NBQQGampoAANbW1jh8+DDefvttNDQ0wM/PD08++STWrVunuBYOdMwWX7ZsGaZMmQIAmDVrFrZv3667gyMaIBjURKrT6TVrUg2vWdNAp01Qi8VihISEMKjJLBjkmjURUV+0DerQ0FA4ODjouFdExk3jsG5ra0N5eTmamprg7u6OQYMG6bJfRGSGrl27pvF6DAxqGsjU+p51Q0MD0tLSMHHiREgkEgQEBOCBBx6Au7s7/P39ER8fj9zcXH31lYhMGIOaSHMqh/XWrVsREBCAnTt3YtKkSfj888+Rl5eHgoICnDx5EuvWrcPdu3cRGxuLqVOnoqioSJ/9JiITwqAm0o7KE8xmz56NtWvX4sEHH+y1nEwmw65du2BtbY0lS5bopJPmhhPMaKAQBAHXrl1DVVWVRvXFYjHCwsJgb2+v454RGQdV84CzwQ2AYU0DAYOaqG+q5oHJ3c+aiIwfg5pIt/jVLSLSKUEQcPXqVVRXV2tU39LSEqGhoQxqonuoFdYWFha93l0K6Lgr1d27d7XqFBGZJgY1kX6oFdZffPFFj/tOnDiBd999F7wETjQw6SKow8LCYGdnp+OeEZk+tcL6qaee6rLt119/RVJSEr7++ms8//zz+Otf/6qzzhGRaWBQE+mXxhPMSktLER8fj8jISNy9exd5eXnIyMjAkCFDdNk/IjJygiDgypUrDGoiPVI7rOvq6vDnP/8ZISEhyM/Px+HDh/H1119j+PDh+ugfERmxzqCuqanRqD6Dmkg1ap0Gf+utt7Bp0yZ4eXnh448/7va0OBENDAxqov6j1qIoFhYWsLOzw+TJk5Xu8Xy/zz//XCedM1dcFIVMnS6COjw8HLa2tjruGZFp0cstMv/4xz/2+dUtIjJvgiDg8uXLqK2t1ai+lZUVwsLCGNREalArrN9//309dYOITAGDmsgwuNwoEamEQU1kOFzBjIj6JAgCLl26hLq6Oo3qM6iJtMMVzIioVwxqIsPjCmZE1CNdBHV4eDhsbGx03DOigYUrmBFRt+RyuVZBbW1tzaAm0hGuYEZEXegiqMPCwhjURDrCFcyISElnUNfX12tUn0FNpHtcwcwAuIIZGSsGNVH/4gpmRKQWBjWR8eIKZkQEuVyO4uJi3LlzR6P6nZPJrK2tddwzIgLUmGB2/fp1tRouKSlRuzNE1P8Y1ETGT+Ww/s1vfoP4+Hjk5OT0WKaurg47d+7E8OHDed2ayARoG9Q2NjYMaqJ+oPJp8IsXL+K///u/MXXqVFhZWWH06NHw8fGBra0tampqcOHCBeTn52P06NFISUnBtGnT9NlvItKSLoI6LCyMQU3UD9SaDQ4ALS0t+Oc//4ljx47h6tWraG5uhpubG6KiovDEE0/w+9Yq4GxwMjS5XI6ioiI0NDRoVJ9BTaQbquaB2mFN2mNYkyExqImMh16+ukVEpq29vR3FxcVaBXV4eDisrKx03DMi6o1W97N+++23AQAFBQWQy+U66RAR6Ye2QW1ra8ugJjIQrUbWndenExMTUVxcDEdHRwwbNgzDhw/H8OHD8eSTT+qkk0SkHV0EdVhYGIOayEB0es26vr4ev/zyC3755RdcuHAB27Zt01XTZoXXrKk/tbe3o6ioCI2NjRrVZ1AT6Y+qeaDWafDp06cr3YXnjTfeQG1treJ5W1sblixZgoSEBAY1kRFgUBOZB7XC+ttvv4VMJlM837RpE6qrqxXP7969i4KCAt31rgdnz55FbGwspFIpXF1dkZCQ0OPpvaqqKvj6+kIkEin9YdGdiRMnQiQSKT2ee+45pTIBAQFdyqxcuVJXh0akMwxqIvOhVljff8bcEN/6Ki0txeTJkxESEoLs7GxkZWUhPz8fcXFx3ZZfvHgxIiMjVW4/Pj4eZWVlikdaWlqXMsnJyUpl1qxZo+nhEOkFg5rIvJjcV7cOHDgAKysrpKamwsKi42+N1NRUREVFobi4GCEhIYqyO3bsQG1tLdauXYuDBw+q1L69vT28vLx6LePk5NRnGSJDaW9vR2FhIZqamjSqb2dnh7CwMFhamtyvByKzpdbIuvO07/3b+pNMJoO1tbUiqIGOXy4AcPz4ccW2CxcuIDk5GZmZmUpl+7Jnzx64ublh2LBheO2117pdinHTpk1wdXXFyJEj8cYbb6C1tVWLIyLSHQY1kXlS6ydSEATExcUp7lfb0tKCP/3pT3BwcAAApevZ+jJp0iS88sorSElJwfLly9HY2IhVq1YBAMrKyhT9mDNnDlJSUjBkyBBcvnxZpbaff/55BAYGwsvLC7/88guSkpLw888/49ChQ4oyy5cvx6hRo+Di4oKcnBwkJSXhypUr+Pvf/95juzKZTOm90fR+wUS9uXv3LoqKihjURGZIrZ/KBQsWKD2fN29elzJ//OMfNerI+vXrsWHDhl7L5ObmYvTo0cjIyMArr7yCpKQkiMViLFu2DJ6enhCLxQCApKQkREREdNu/3sTHxyv+PXz4cISGhmL06NE4e/YsRo0aBaDjO+WdIiMj4eLigt///veK0XZ3Nm7c2OexEWmDQU1k3oxmbfDKykpUVlb2WiYgIAC2traK5xUVFXBwcIBIJIKzszP27t2L2bNnY+TIkTh//rziFL0gCJDL5RCLxVi9erXKwSkIAmxsbPDBBx/g2Wef7bZMSUkJfH19cerUKYwdO7bbMt2NrP38/Pg9a9IJBjWR6TK5tcHd3Nzg5uamVh1PT08AQHp6OmxtbREbGwsA2LdvH5qbmxXlcnNzsWjRIhw7dgzBwcEqt5+fn4+2tjZ4e3v3WObcuXMA0GsZGxsbxaUDIl26e/cuCgsLlf6/q8Pe3h6hoaEMaiIjZ5I/odu3b0d0dDQcHR1x6NAhrFixAm+++SakUikAdAnkzhF7RESEokxJSQliYmKQmZmJMWPG4NKlS9izZw+mT58ONzc3XLhwAa+++iqioqLw8MMPAwBOnjyJU6dO4fHHH4dEIkFubi4SExMxa9YsDBkypN+OnwhgUBMNJCb5U5qTk4N169ahoaEBQ4cORVpaGubPn69WG21tbSgoKFCcOrS2tsbhw4fx9ttvo6GhAX5+fnjyySexbt06xbVwGxsbfPLJJ9iwYQNkMhn8/f0RHx+P119/XefHSNQbBjXRwGI016wHEq4NTtrQRVCHhYUp/gglIsMxuWvWRNQ3bYPawcEBoaGhDGoiE8OwJjIRbW1tKCwsREtLi0b1GdREpkutFcyIyDAY1EQDG8OayMgxqImIYU1kxBjURATwmjWR0dI2qB0dHRESEsKgJjIDDGsiI9S5DoCmN8dhUBOZF4Y1kZFpbW1FYWEhg5qIFBjWREZEF0EdGhqq1j3cicj4MayJjIS2Qe3k5ISQkBAGNZEZ4k81kRFgUBNRb/iTTWRgDGoi6gtPgxMZUGtrKwoKCtDa2qpRfQY10cDAsCYyEG2D2tnZGcHBwQxqogGAYU1kADKZDIWFhQxqIlIJf9KJ+hmDmojUxZ92on7EoCYiTfAnnqifMKiJSFO8Zk3UD2QyGQoKCtDW1qZRfYlEgqCgIAY10QDFsCbSM10EdXBwMEQikY57RkSmgmFNpEctLS0oLCxkUBORVhjWRHqibVBLpVIEBQUxqImIE8yI9IFBTUS6xLAm0jEGNRHpGsOaSIcY1ESkD7xmTaQjzc3NKCwsxN27dzWq7+LigsDAQAY1EXXBkTWRDjCoiUifGNZEWmJQE5G+MayJtMCgJqL+wGvWRBrSNqgHDRqEgIAABjUR9YlhTaSBpqYmFBUVMaiJqF/wNDiRmhjURNTfOLImUkNTUxMKCwvR3t6uUX1XV1f4+/szqIlILRxZE6mIQU1EhsKwJlIBg5qIDImnwYn60NjYiKKiIgY1ERkMw5qoF9oGtZubG/z9/XXcKyIaaHganKgHDGoiMhYmGdZnz55FbGwspFIpXF1dkZCQgIaGhm7LVlVVwdfXFyKRCLW1tX22ffLkSUyaNAkODg6QSqWYOHEimpubFftramowf/58SCQSSCQSzJ8/X6V2ybQwqInImJhcWJeWlmLy5MkICQlBdnY2srKykJ+fj7i4uG7LL168GJGRkSq1ffLkSUydOhVTpkxBTk4OcnNzsXTpUlhY/Pttmjt3LvLy8pCVlYWsrCzk5eVh/vz5ujg0MhLaBrW7uzuDmoh0yuSuWR84cABWVlZITU1VhGhqaiqioqJQXFyMkJAQRdkdO3agtrYWa9euxcGDB/tsOzExEcuWLcPKlSsV20JDQxX/vnjxIrKysnDq1CmMHTsWALBz506MHz8eBQUFCA8P19VhkoE0NDSgqKgIcrlco/ru7u4YMmSIjntFRAOdyY2sZTIZrK2tlUa7dnZ2AIDjx48rtl24cAHJycnIzMxUKtuTW7duITs7Gx4eHoiOjoanpycmTJig1ObJkychkUgUQQ0A48aNg0QiwYkTJ3RxeGRADGoiMlYmF9aTJk1CeXk5UlJS0NraipqaGqxatQoAUFZWBqAj0OfMmYOUlBSVf3levnwZALB+/XrEx8cjKysLo0aNQkxMDIqKigAA5eXl8PDw6FLXw8MD5eXlPbYtk8lQX1+v9CDjom1Qe3h4MKiJSG+MJqzXr18PkUjU6+P06dMYNmwYMjIysHnzZtjb28PLywtBQUHw9PSEWCwGACQlJSEiIgLz5s1T+fU7f0m/8MILWLhwIaKiorB161aEh4cjPT1dUa6778oKgtDrd2g3btyomJAmkUjg5+encr9I/3QR1PxMiUifjOaa9dKlS/Hcc8/1WiYgIABAxySvuXPnoqKiAg4ODhCJRNiyZQsCAwMBAEeOHMH58+fx2WefAegIU6Bjhu7q1auxYcOGLm17e3sDAB544AGl7REREbh+/ToAwMvLCxUVFV3q3r59G56enj32OykpCa+88orieX19PX+5G4k7d+6guLiYQU1ERs1owtrNzQ1ubm5q1ekMyPT0dNja2iI2NhYAsG/fPqWvW+Xm5mLRokU4duwYgoODu20rICAAPj4+KCgoUNpeWFiIadOmAQDGjx+Puro65OTkYMyYMQCA7Oxs1NXVITo6usd+2tjYwMbGRq1jI/1jUBORqTCasFbH9u3bER0dDUdHRxw6dAgrVqzAm2++CalUCgBdArmyshJAxyi5s0xJSQliYmKQmZmJMWPGQCQSYcWKFVi3bh1GjBiBkSNHIiMjA7/++qtihB4REYGpU6ciPj4eaWlpAICEhATMmDGDM8FNjLZB7enpCV9fXx33ioioeyYZ1jk5OVi3bh0aGhowdOhQpKWlqf1d57a2NhQUFKCpqUmx7eWXX0ZLSwsSExNRXV2NESNG4NChQ0rhv2fPHixbtgxTpkwBAMyaNQvbt2/XzYFRv2BQE5GpEQmdF3Sp39TX10MikaCurg7Ozs6G7s6AUl9fj0uXLjGoicgoqJoHJjmyJtKEtkHt5eWFwYMH67hXRER9M5qvbhHpE4OaiEwZR9Zk9urr61FcXAxNr/gwqInI0BjWZNa0DWpvb2/4+PjouFdEROphWJPZqqurw6VLlxjURGTyeM2azBKDmojMCUfWZHa0DWofHx/F8rNERMaAI2syKwxqIjJHHFmT2aitrcXly5cZ1ERkdhjWZBa0DerBgwfDy8tLx70iItINngYnk8egJiJzx5E1mbSamhpcuXKFQU1EZo1hTSZL26D29fVV3BOdiMiY8TQ4mSQGNRENJAxrMjkMaiIaaHganExKdXU1rly5onF9BjURmSKGNZkMbYPaz88PHh4eOuwREVH/4GlwMgkMaiIayBjWZPQY1EQ00PE0OBm1qqoqXL16VeP6Q4YMgbu7u+46RERkABxZk9FiUBMRdWBYk1FiUBMR/RtPg5PRqaysxLVr1zSu7+/vDzc3Nx32iIjIsDiyJqPCoCYi6ophTUaDQU1E1D2GNRkFBjURUc94zZoM7vbt27h+/brG9QMCAuDq6qrDHhERGReOrMmgGNRERH1jWJPBMKiJiFTD0+BkELdu3cKNGzc0rh8YGIhBgwbpsEdERMaLI2vqdwxqIiL1MKypXzGoiYjUx9Pg1G8qKipw8+ZNjeszqIlooGJYU7/QJqhFIhECAwPh4uKi414REZkGngYnvWNQExFphyNr0qvy8nKUlJRoVJdBTUTUgSNr0hsGNRGRbnBkTXqhbVAHBQVBKpXqtlNERCbKJEfWZ8+eRWxsLKRSKVxdXZGQkICGhoZuy1ZVVcHX1xcikQi1tbV9tn3y5ElMmjQJDg4OkEqlmDhxIpqbmxX7AwICIBKJlB4rV67U1aGZhbKyMgY1EZEOmVxYl5aWYvLkyQgJCUF2djaysrKQn5+PuLi4bssvXrwYkZGRKrV98uRJTJ06FVOmTEFOTg5yc3OxdOlSWFgov03JyckoKytTPNasWaPtYZmNsrIylJaWalSXQU1E1D2TOw1+4MABWFlZITU1VRGiqampiIqKQnFxMUJCQhRld+zYgdraWqxduxYHDx7ss+3ExEQsW7ZMaaQcGhrapZyTkxO8vLx0cDTmRdugDg4OhkQi0XGviIhMn8mNrGUyGaytrZVGu3Z2dgCA48ePK7ZduHABycnJyMzM7DIy7s6tW7eQnZ0NDw8PREdHw9PTExMmTFBqs9OmTZvg6uqKkSNH4o033kBra2uffa6vr1d6mJvS0lIGNRGRnphcWE+aNAnl5eVISUlBa2srampqsGrVKgAdIzugIxznzJmDlJQUDBkyRKV2L1++DABYv3494uPjkZWVhVGjRiEmJgZFRUWKcsuXL8fevXvx/fffY+nSpdi2bRteeumlXtveuHEjJBKJ4uHn56fJoRut0tJSxXuvLgY1EVHfjCas169f32Xi1v2P06dPY9iwYcjIyMDmzZthb28PLy8vBAUFwdPTE2KxGACQlJSEiIgIzJs3T+XXl8vlAIAXXngBCxcuRFRUFLZu3Yrw8HCkp6cryiUmJmLChAmIjIzEkiVL8N5772HXrl2oqqrqse2kpCTU1dUpHtqsjW1sGNRERPpnNNesly5diueee67XMgEBAQCAuXPnYu7cuaioqICDgwNEIhG2bNmCwMBAAMCRI0dw/vx5fPbZZwAAQRAAAG5ubli9ejU2bNjQpW1vb28AwAMPPKC0PSIiotd7Lo8bNw4AUFxc3OO9lW1sbGBjY9PrsZkibYM6JCQEzs7OOu4VEZH5MZqwdnNzg5ubm1p1PD09AQDp6emwtbVFbGwsAGDfvn1KX7fKzc3FokWLcOzYMQQHB3fbVkBAAHx8fFBQUKC0vbCwENOmTeuxD+fOnQPw77AfKEpKSlBeXq5RXQY1EZF6jCas1bF9+3ZER0fD0dERhw4dwooVK/Dmm28qvvJzfyBXVlYC6Bgld5YpKSlBTEwMMjMzMWbMGIhEIqxYsQLr1q3DiBEjMHLkSGRkZODXX39VjNBPnjyJU6dO4fHHH4dEIkFubi4SExMxa9Ysla+NmwNtgtrCwgLBwcEMaiIiNZhkWOfk5GDdunVoaGjA0KFDkZaWhvnz56vVRltbGwoKCtDU1KTY9vLLL6OlpQWJiYmorq7GiBEjcOjQIUX429jY4JNPPsGGDRsgk8ng7++P+Ph4vP766zo9PmPGoCYi6n8iofOCLvWb+vp6SCQS1NXVmVRw3bx5ExUVFRrVtbCwQEhICJycnHTcKyIi06VqHhjNbHAybgxqIiLDYVhTnxjURESGZZLXrKn/3LhxA7du3dKoLoOaiEg3GNbUI22DOjQ0FI6OjjruFRHRwMPT4NQtBjURkfHgyJq6uH79Om7fvq1RXQY1EZHucWRNShjURETGhyNrUtAmqMViMUJCQhjURER6wLAmAMC1a9cUy7KqSywWIzQ0FA4ODjruFRERATwNTmBQExEZO46sTYQgCKhpakOj7C4cbCzhYm8FkUikdbsMaiIi48ewNnJ1zW3Yd+YmMk5cxbXqf990xH+QPRZEB+CZh3whsbPSqG0GNRGRaeCNPAxA1YXbfyi8jRc/PIPm1nYAwL0fVOeY2s5ajB3zHsKEMHeVX18QBFy7dg1VVVUa9L4jqMPCwmBvb69RfSIi6sAbeZi4HwpvY+HuHDS3tUOAclDj/54LAJrb2rFwdw5+KFRtFjeDmojI9DCsjVBdcxte/PBMRyD3cd5DEDpC+8UPz6Cuua2PsgxqIiJTxLA2QvvO3ERza3ufQd1JEIDm1nZ8fvZmL2UEXL16VeOgtrS0ZFATERkIw9rICIKAjBNXNar7/k9X0d0UhM6grq6u1qhdS0tLhIaGMqiJiAyEYW1kapracK26qcs16r4IAK5VN6G2SflUuC6CmiNqIiLDYlgbmUbZXa3qN9xTX1dBbWdnp1WfiIhIO/yetZFxsNHuI3H8v/qCIODKlSuoqanRqB0GNRGR8eDI2si42FvBf5A91F2bTISOhVKk9lYMaiIiM8OwNjIikQgLogM0qhv3cEc9BjURkXlhWBuhZx7yhZ21GKou/W0h6ljJ7HdRg7UO6vDwcAY1EZGRYVgbIYmdFXbMewgioM/A7ty/4/lRqCq7oXFQW1lZITw8HLa2thrVJyIi/WFYG6kJYe7YvXAM7KzEHaF93/7ObXZWYuyO+w18Le+gtrZWo9eysrJCWFgYg5qIyEhxNrgRmxDmjpNJMfj87E28/5PyXbeGDLJH3MMBeHrUYFSW3mBQExGZMd51ywBUvcvKvQRBQG1TGxpkd+FoYwmpfcdtMS9duoS6ujqN+sGgJiIyLFXzgCNrEyESieDiYA0XB2sAHeHNoCYiGhgY1iZIF0EdHh4OGxsbHfeMiIj0gWFtYuRyOS5fvqxxUFtbWyMsLIxBTURkQhjWJkQul+PSpUuor6/XqD6DmojINPGrWyaCQU1ENHAxrE0Ag5qIaGBjWBs5BjUREfGatRGTy+UoLi7GnTt3NKpvbW2N8PBwWFtb67hnRETUnziyNlIMaiIi6sSRtRHSNqhtbGwQFhbGoCYiMhMmObI+e/YsYmNjIZVK4erqioSEBDQ0NHRbtqqqCr6+vhCJRL2un3316lWIRKJuH59++qmiXE1NDebPnw+JRAKJRIL58+drvC53dxjURER0P5ML69LSUkyePBkhISHIzs5GVlYW8vPzERcX1235xYsXIzIyss92/fz8UFZWpvTYsGEDHBwcMG3aNEW5uXPnIi8vD1lZWcjKykJeXh7mz5+vk2OTy+UoKipiUBMRkRKTOw1+4MABWFlZITU1FRYWHX9rpKamIioqCsXFxQgJCVGU3bFjB2pra7F27VocPHiw13bFYjG8vLyUtn3xxRd49tln4ejoCAC4ePEisrKycOrUKYwdOxYAsHPnTowfPx4FBQUIDw/X+Lg6g7qnMwR9YVATEZkvkxtZy2QyWFtbK4IaAOzs7AAAx48fV2y7cOECkpOTkZmZqVRWVWfOnEFeXh4WL16s2Hby5ElIJBJFUAPAuHHjIJFIcOLECU0OBwDQ3t6udVBzMhkRkfkyubCeNGkSysvLkZKSgtbWVtTU1GDVqlUAgLKyMgAdgT5nzhykpKRgyJAhGr3Orl27EBERgejoaMW28vJyeHh4dCnr4eGB8vLyHtuSyWSor69XenRqb29HcXGxxkFta2uL8PBwWFlZaVSfiIiMn9GE9fr163uc4NX5OH36NIYNG4aMjAxs3rwZ9vb28PLyQlBQEDw9PSEWiwEASUlJiIiIwLx58zTqS3NzMz766COlUXUnkUjUZZsgCN1u77Rx40bFhDSJRAI/Pz8AugnqsLAwBjURkZkTCYIgGLoTAFBZWYnKyspeywQEBCjdf7miogIODg4QiURwdnbG3r17MXv2bIwcORLnz59XBKggCJDL5RCLxVi9ejU2bNjQ6+t88MEHWLx4MUpKSuDu7q7Ynp6ejldeeaXL7G+pVIqtW7di4cKF3bYnk8kgk8kUz+vr6+Hn54fc3NxeQ743DGoiItNXX18PiUSCuro6ODs791jOaCaYubm5wc3NTa06np6eADpC1NbWFrGxsQCAffv2obm5WVEuNzcXixYtwrFjxxAcHNxnu7t27cKsWbOUghoAxo8fj7q6OuTk5GDMmDEAgOzsbNTV1SmdLr+fjY1Nt8t9NjY2KiavqYNBTUQ0sBhNWKtj+/btiI6OhqOjIw4dOoQVK1bgzTffhFQqBYAugdw5Yo+IiFCUKSkpQUxMDDIzMxXBCwDFxcX48ccf8c9//rPL60ZERGDq1KmIj49HWloaACAhIQEzZszQaia4OhjUREQDj0mGdU5ODtatW4eGhgYMHToUaWlpan/Xua2tDQUFBWhqalLanp6ejsGDB2PKlCnd1tuzZw+WLVum2D9r1ixs375dswNRE4OaiGhgMppr1gNJ5zWKo0ePqnwa3M7ODmFhYbC0NMm/r4iIqBuqXrM2mtng1DMGNRHRwMawNnIMaiIiYlgbMQY1EREBJjrBzNR1ThNobGzssYydnR28vLy6TIAjIiLz0bmiZV/TxzjBzABu3rypWMWMiIjoxo0b8PX17XE/w9oA5HI5SktL4eTk1OcKZp2rnd24caPXmYKkX/wcjAM/B8PjZ6BbgiDgzp078PHx6fWmUzwNbgAWFha9/gXVHWdnZ/5gGAF+DsaBn4Ph8TPQHYlE0mcZTjAjIiIycgxrIiIiI8ewNnI2NjZYt25dtzcCof7Dz8E48HMwPH4GhsEJZkREREaOI2siIiIjx7AmIiIycgxrIiIiI8ewJiIiMnIMawM7e/YsYmNjIZVK4erqioSEBDQ0NHRbtqqqCr6+vhCJRKitre2xzatXr0IkEnX7+PTTT/V0JKZNH59Dp5MnT2LSpElwcHCAVCrFxIkT0dzcrOMjMA/6+hwmTpzY5Wfhueee08MRmAd9/jwAHat2TZs2DSKRCF9++aXuOm7GGNYGVFpaismTJyMkJATZ2dnIyspCfn4+4uLiui2/ePFiREZG9tmun58fysrKlB4bNmyAg4MDpk2bpuOjMH36+hyAjqCeOnUqpkyZgpycHOTm5mLp0qW9Lis4UOnzcwCA+Ph4pZ+JtLQ0HfXcvOj7cwCAbdu29bnUMt1HIINJS0sTPDw8hPb2dsW2c+fOCQCEoqIipbJ/+9vfhAkTJgiHDx8WAAg1NTVqvdbIkSOFRYsW6aLbZkefn8PYsWOFNWvW6KPbZkefn8OECROE5cuX66HX5kffv5fy8vIEX19foaysTAAgfPHFFzo+AvPEP+8NSCaTwdraWmmUZWdnBwA4fvy4YtuFCxeQnJyMzMxMjUZkZ86cQV5eHhYvXqx9p82Qvj6HW7duITs7Gx4eHoiOjoanpycmTJig1Cb9m75/Hvbs2QM3NzcMGzYMr732Gu7cuaO7zpsRfX4OTU1NmDNnDrZv3w4vLy/ddtzMMawNaNKkSSgvL0dKSgpaW1tRU1ODVatWAQDKysoAdPzgzJkzBykpKRgyZIhGr7Nr1y5EREQgOjpaZ303J/r6HC5fvgwAWL9+PeLj45GVlYVRo0YhJiYGRUVF+jkYE6bPn4fnn38eH3/8MY4ePYq//OUv2LdvH55++mm9HIep0+fnkJiYiOjoaDz11FN66bs5Y1jrwfr163uc4NX5OH36NIYNG4aMjAxs3rwZ9vb28PLyQlBQEDw9PSEWiwEASUlJiIiIwLx58zTqS3NzMz766KMBOao29Ocgl8sBAC+88AIWLlyIqKgobN26FeHh4UhPT9fLMRsjQ38OQMf16smTJ2P48OF47rnn8Nlnn+G7777D2bNn9XHIRsnQn8P+/ftx5MgRbNu2TU9HaN643KgeVFZWorKystcyAQEBsLW1VTyvqKiAg4MDRCIRnJ2dsXfvXsyePRsjR47E+fPnFZMxBEGAXC6HWCzG6tWrsWHDhl5f54MPPsDixYtRUlICd3d37Q/OhBj6c7hy5QqCgoLwwQcfKP1Se/bZZ2FpaYk9e/bo6EiNm6E/h+4IggAbGxt88MEHePbZZzU/OBNi6M/h5ZdfxjvvvKN0yry9vR0WFhZ49NFHcfToUd0cqJni/az1wM3NDW5ubmrV8fT0BACkp6fD1tYWsbGxAIB9+/Ypfc0nNzcXixYtwrFjxxAcHNxnu7t27cKsWbMGXFADhv8cAgIC4OPjg4KCAqXthYWFA2pWvqE/h+7k5+ejra0N3t7eavXLlBn6c1i5ciWWLFmitO3BBx/E1q1bMXPmTLX6NSAZcHIbCYLw7rvvCmfOnBEKCgqE7du3C3Z2dsLbb7/dY/nvv/++y6zLmzdvCuHh4UJ2drZS2aKiIkEkEgkHDx7UV/fNhr4+h61btwrOzs7Cp59+KhQVFQlr1qwRbG1theLiYn0ejsnSx+dQXFwsbNiwQcjNzRWuXLkifPPNN8LQoUOFqKgo4e7du/o+JJOkz99L9wJng6uMI2sDy8nJwbp169DQ0IChQ4ciLS0N8+fPV6uNtrY2FBQUoKmpSWl7eno6Bg8ejClTpuiyy2ZJX5/Dyy+/jJaWFiQmJqK6uhojRozAoUOH1BoFDiT6+Bysra1x+PBhvP3222hoaICfnx+efPJJrFu3TnENlpTp8/cSaYbXrImIiIwcZ4MTEREZOYY1ERGRkWNYExERGTmGNRERkZFjWBMRERk5hjUREZGRY1gTEREZOYY1ERGRkWNYExERGTmGNRH1qqqqCh4eHrh69apW7QiCgISEBAwaNAgikQh5eXn4/e9/jy1btuimo0RmjGFNNIDExcV1ex/jqVOn9lhn48aNmDlzJgICArR67aysLLz//vs4cOAAysrKMHz4cKxduxZvvPEG6uvrVe7/ypUre9z329/+tsv2o0ePQiQSoba2VoveExkWb+RBNMBMnToVu3fvVtpmY2PTbdnm5mbs2rUL//znP3tsr7W1FdbW1n2+7qVLl+Dt7Y3o6GjFtsjISAQEBGDPnj148cUXe60vl8vxzTffYP/+/X2+FpG54ciaaICxsbGBl5eX0sPFxaXbsgcPHoSlpSXGjx+v2DZx4kQsXboUr7zyCtzc3BT3OBYEAW+99RaCgoJgZ2eHESNG4LPPPgPQMer9z//8T1y/fh0ikUhplD5r1ix8/PHHffb7p59+goWFBcaOHavF0XeM8B955BFIpVK4urpixowZuHTpklZtEukbw5qIevTjjz9i9OjRXbZnZGTA0tISP/30E9LS0gAAa9aswe7du7Fjxw7k5+cjMTER8+bNww8//IC3334bycnJ8PX1RVlZGXJzcxVtjRkzBjk5OZDJZL32Zf/+/Zg5cyYsLLT7tdXY2IhXXnkFubm5OHz4MCwsLPC73/0Ocrlcq3aJ9ImnwYkGmAMHDsDR0VFp25///Gf85S9/wXfffYfz588jMTERAHD16lX4+Ph0aSMkJARvvfWW4nljYyO2bNmCI0eOKEbhQUFBOH78ONLS0vDRRx/ByckJYrEYXl5eSm0NHjwYMpkM5eXl8Pf377Hf+/fvx//8z/+ofWzt7e1Kz5955hml57t27YKHhwcuXLiA4cOH99o+kaEwrIkGmMcffxw7duxQ2jZo0CAAwOTJkzF58mTF9ubmZtja2nZp4/7R9oULF9DS0qI4Jd6ptbUVUVFRvfbHzs4OANDU1NRjmYsXL+LmzZtKfetOd8eWnZ2NefPmKZ5funQJf/nLX3Dq1ClUVlYqRtTXr19nWJPRYlgTDTAODg4ICQnpdt+0adOwZcsWREREAADc3NxQU1PTbRv36gy8b775BoMHD1ba19PktU7V1dUAAHd39x7L7N+/H7GxsYpg70l3x3bz5k2l5zNnzoSfnx927twJHx8fyOVyDB8+HK2trb22TWRIDGsiUigqKkJoaKjieVRUFD788MM+6z3wwAOwsbHB9evXMWHCBLVe85dffoGvry/c3Nx6LPPVV19hyZIlarXbnaqqKly8eBFpaWl49NFHAQDHjx/Xul0ifWNYEw0wndeH72VpaQkrKys4OjrC0vLfvxaeeOIJJCUloaampscZ4wDg5OSE1157DYmJiZDL5XjkkUdQX1+PEydOwNHREQsWLOix7rFjxzBlypQe99+6dQu5ubn48ssvVT/IHri4uMDV1RX/+7//C29vb1y/fr3H720TGROGNdEAk5WVBW9vb6Vt4eHh2LVrF4YNG6a0/cEHH8To0aPxj3/8Ay+88EKv7f71r3+Fh4cHNm7ciMuXL0MqlWLUqFFYtWpVj3VaWlrwxRdf4Ntvv+2xzNdff42xY8fCw8NDhaPrnYWFBfbu3Ytly5Zh+PDhCA8PxzvvvIOJEydq3TaRPokEQRAM3QkiMry0tDRUV1cjKSlJafs///lPvPbaa/jll1+0/trU/VJTU/HVV1/h//2//9djmVmzZuGRRx7B66+/rtPXJjIlHFkTEQAgPz+/y2xuAJg+fTqKiopQUlICPz8/nb6mlZUV3n333V7LPPLII5gzZ45OX5fI1HBkTUREZOS4ghkREZGRY1gTEREZOYY1ERGRkWNYExERGTmGNRERkZFjWBMRERk5hjUREZGRY1gTEREZOYY1ERGRkWNYExERGTmGNRERkZFjWBMRERm5/w8TkafIZgaaXgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# define plot\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.scatter(ref['energy'], pred['energy'], s=100.0)\n",
    "plt.rcParams[\"figure.autolayout\"] = True\n",
    "plt.rcParams[\"font.size\"] = 14\n",
    "\n",
    "# get limits of axis and aplly same for x and y-axis (square plot)\n",
    "limits = [*plt.xlim(), *plt.ylim()]\n",
    "lower_lim = min(limits)\n",
    "upper_lim = max(limits)\n",
    "plt.xlim(lower_lim, upper_lim)\n",
    "plt.ylim(lower_lim, upper_lim)\n",
    "\n",
    "# plot y=x line to guide eye\n",
    "plt.plot([lower_lim-2.5, upper_lim+2.5], [lower_lim-2.5, upper_lim+2.5], c='k', alpha=0.2, linewidth=8, linestyle='-')\n",
    "\n",
    "plt.xlabel('E$_i$(ref) / Ha')\n",
    "plt.ylabel('E$_i$(NN) / Ha')\n",
    "plt.title('MAE: '+str(round(mae, 4)))\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}