{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PhyGNN Test with Pythagorean's Theorem\n",
    "\n",
    "In this example test we create a dataset with pythagorean's theorem `a^2 + b^2 = c^2` where `a` and `b` are input features and we are trying to predict `c`. We train on a noisy and biased `c` dataset resulting in a predictably biased neural network. We then train with an augmented loss function that includes 80% weight on the predicted vs. physical calculation of `c`. The neural network loses its bias and is able to predict much more accurately. \n",
    "\n",
    "Obviously this is a contrived situation where we know the solution for `c`, but the physical loss function can be created using whatever benchmark might be applicable to a given physical domain. This example is simply intended to show how to build and train a physics guided neural network. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "\n",
    "from phygnn import PhysicsGuidedNeuralNetwork"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2.1.0'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pythag Inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we set up the training features and known output values based on the pythagorean theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((10000, 1), (10000, 2), (10000, 2))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 100\n",
    "\n",
    "a = np.linspace(-1, 1, N)\n",
    "b = np.linspace(-1, 1, N)\n",
    "a, b = np.meshgrid(a, b)\n",
    "\n",
    "a = np.expand_dims(a.flatten(), axis=1)\n",
    "b = np.expand_dims(b.flatten(), axis=1)\n",
    "\n",
    "y = np.sqrt(a ** 2 + b ** 2)\n",
    "x = np.hstack((a, b))\n",
    "p = x.copy()\n",
    "\n",
    "y.shape, x.shape, p.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we make a y_noise dataset which is noisy and systematically biased."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y_noise')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3zU9Z3v8dcnIUBACFbQKopYl3pBwEsUFdtqt9YroogIaLvb0z082nN6WrWlta2Pimhbdz27257txVLX7XarXIURtyq6tV1bFQQaIGJFkVplaEWU3CCS2+f88ZuEScgkk2R+M7+ZeT8fDx5k5vf7ZT5cZt75/r43c3dERES6U5LrAkREJLoUEiIikpJCQkREUlJIiIhISgoJERFJaVCuC8ik0aNH+/jx43NdhohIXtm0adNedx/T3bGCConx48ezcePGXJchIpJXzOxPqY7pdpOIiKSkkBARkZQUEiIikpJCQkREUlJIiIhISgU1uklEpBjEquLct3Y7u2saOW5UOQsuO4VrzxobymspJERE8kisKs7XV1XT2NwKQLymka+vqgYIJSh0u0lEJI/ct3Z7R0C0a2xu5b6120N5PYWEiEge2V3T2KfnB0q3m0REIuSOWDVL1r9FqzslBkMGlfB+c1tH38Nxo8qJdxMIx40qD6UehYSISETcEavmF+ve7Hjc5tDY3AYEfQ+3LNvc7XXlZaUsuOyUUGpSSIiIRMTD69/s/aQuxoY8ukl9EiIiEdHmua7gcAoJEZE81j4ENlYVD+X7h3q7ycweBK4G9rj7Gd0cXwDclFTLacAYd3/PzN4A6oFWoMXdK8OsVUQkW7pOhrvk1DH8+pV3+v392ofAhnHLKew+iZ8BPwB+3t1Bd78PuA/AzKYDt7r7e0mnXOLue0OuUUQka7qbDJfcWd1fYQ2BDfV2k7s/C7zX64mBucCSEMsREcmpWFWcLy/fcthkuEwIawhsJPokzGwYcDnwSNLTDjxlZpvMbH4P1843s41mtvGdd/rfXBMRCVOsKs6ClVto9cz3ThfDENjpwHNdbjVNc/fdZnY08LSZvZJomXTi7ouBxQCVlZURHBsgIgJ3PbaN5tbMf0QVyxDYOXS51eTuuxO/7wFWA+floC4RkYzYd6A549+zPSDuW7udk27/JdPufSbjo5xyHhJmVgF8DHg06bnhZjai/Wvgk8BLualQRCR6ystKueTUMXx9VTXxmkaccIbDhhoSZrYEeAE4xcx2mdlnzexzZva5pNOuA55y9/1Jzx0D/M7MtgAvAr909yfDrFVEJF+MHVXOd2dO4tevvBP6irCh9km4+9w0zvkZwVDZ5Od2AlPCqUpEJDNiVXEWrtlGTWPnW0lHDivjzukTgWBp7+4W5OuvYWUlPHf7xwG4NcVaTpkcDhuVjmsRkbwRq4rzzdXV7G/qfijrvgPN3LY8+ADP9FIbgweVdnydjRVhc94nISKST9qHsqYKiHZtHs5aTLVJrZYFl51CeVlpp+OZHg6rloSISBpiVXHuemxbKKOU+qKivKzj6/Zhr2Hud62QEBHpQbCMxtaOfR1yzazz42vPGhvaHAnQ7SYRkZRiVXFuW7Y5MgEBQX9HGPMhUlFIiIiksHDNNqITD4eEvTx4Mt1uEhFJ0r6MdyaHrYYhzOXBkykkREQSui7jHXVhLQ+eTCEhIkXvjlg1S9a/FcoKrWEKa3nwZAoJESlqd8SqM7LpT7aFuTx4MoWEiBSlWFWcb6zayoEIjVxKV9jLgydTSIhIUYnKpLj+MuhYuykbFBIiUjTyrWO6O9noh0imeRIiUjTuemxbXgdEtvohkqklISIFL2pLa/RHNvshkikkRKSgxariLFixheYwlmTNgvKyUr47cxIQLOR367LNoSzkl4pCQkQKUr53UMOh1gPQqS+lfVkOQDOuRUR60zUQykogj+8sAVBq1jGKadq9z6TcpjTskAh7j+sHzWyPmb2U4vjFZlZrZpsTv76VdOxyM9tuZjvM7PYw6xSR/HVHrJpblm3u1GLI94AAaHUnVhVn2r3PpFxHqhCW5fgZ8APg5z2c81t3vzr5CTMrBX4IXArsAjaY2Rp3fzmsQkUk/9z00xd47vX3cl1GKMzgtmWbe1yFNhvDYUNtSbj7s0B//gXPA3a4+053bwKWAjMyWpyI5K07YtWMv/2XBRsQAO70GBBG0DcR9t4SUZgncYGZbTGzJ8xsYuK5scBbSefsSjx3GDObb2YbzWzjO++8E3atIpJj+brWUiYZ0D5WK+y9JXIdEr8HTnT3KcC/ALHE89bNud2OX3P3xe5e6e6VY8aMCalMEcm1WFWcM+96qugDAg7/MGzvxA5DTkPC3evcvSHx9eNAmZmNJmg5nJB06vHA7hyUKCIR0L6NaE1j/g5nDVtYndg5DQkz+6BZsK23mZ2XqOddYAMwwcxOMrPBwBxgTe4qFZFc+saqrZHcRjQXhpV1/7EdVid2qKObzGwJcDEw2sx2AXcCZQDufj8wC/i8mbUAjcAcd3egxcy+AKwFSoEH3X1bmLWKSHS0byG6u6aRUcPK8nI570wrNWPu1BOoPPEDhy1SGOaaTqGGhLvP7eX4DwiGyHZ37HHg8TDqEpFoilXF+ebqavY3HfoAzOcZ05kydlT5YcuDt4do2Et0aMa1iERCrCrOl1dsoTVP11gKS3ethGvPGpu1hf4UEiKSU3fEqnlo3ZvdD18scqVmXH/O4YGQfDsu7JZErofAikgRa5/zoIDoXqs7j2yKd5oD0b5xUrymESeYJ3HLv6xgwT/8JJQa1JIQkaxr/0k41ZpEckjXhfzuW7u9o9P64O7t1Dz3MO/v3MQPRh/P33/lf1JSktmf/RUSIpJVhbzeUliS50Dsrmnk4O7t1D63hMadGykpH8moj/0tI8++KuMBAQoJEcmSYl1OI3kJjf5qnwPx4osvUhO7m5rt6xPh8DeMOPtqSgaXMzYf50mIiBRrOEDnuQ0LVm6hubVzXJSVWK875pWXlXLt8Y1UfuQTbPrdrygZOiIIh7OuomTIsI5z8nKehIgUt2IOCAg6nh9K/PnvmzWFb6za2mli4OBBJTQ3taa6nJa/7KBl2yq+es+zlJSPYNRHPx20HIYM62ihhL33tUJCREJx6T/9htf27M91GTnn0BEU3mXt0v0pAuLgn1+j9rmHaXx9AyXlIzj2rz/DoElXdLQc2r9vNigkRCTjJt/5JHUHU/+EXGwceHj9m/Q2T/DgX3YE4bDjRUqGHtGp5ZBKvKaRBSu2AOHsd62QEJGMuCNWzZL1b9HqmvXQnZ4C4rBw+MinGHHO9B7DIVlzm7NwzTaFhIhET6wqzq3LNmtCXD8E4bCExh3rew2HslI7rOM7WVjLqCskRKTfYlVxblm2Oddl5J2mt1+n5rklNL62jpIhw6n4yM2MPGc6JUOGp7ymp4AIk0JCRPpMM6b7p+ntndQ89/ChcLjoJkZWXtNjOKTryGFlGajwcAoJEekTLcjXd2GGAwS3ou6cPjEj36srhYSIpCVWFWfhmm3aQrQPmt7eSc3zS2h89QWsPRzOmU7J0CMy9hqaJyEiOde+8mjybmiSWtOenUGfQ3s4TJsXtBwyGA7tum5GlGkKCRFJSX0PfdO054/UPreEA68+36dwGMj6TtPufSZ/WxJm9iBwNbDH3c/o5vhNwNcSDxuAz7v7lsSxN4B6oBVocffKMGsVkc60Wmv6OoXD4GFUTJvLyMoZabccBtK/E69p5OurqoFwJtOZhzjxxcw+SvDh//MUIXEh8Ad332dmVwAL3X1q4tgbQKW770339SorK33jxo2ZKV6kSCkc0tc1HEZWzmDEuTMozcBtpfKyEhqT1nlKR3/7J8xsU6ofxENtSbj7s2Y2vofjzyc9XAccH2Y9IpJarCrOl5dvJkfD8fNK0ztvBOGw/bmg5XDh3IyFQ7vG5jbKy0o79QN1fdxVGK2KKPVJfBZ4IumxA0+ZmQM/cffF3V1kZvOB+QDjxo0LvUiRQlTsq7WmKwiHpRzY/jtscDkVF85hROUMSstHhPJ63505qWMv64ryMszodfBA153sBioSIWFmlxCExEVJT09z991mdjTwtJm94u7Pdr02ER6LIbjdlJWCRQqEwiE9h4XDBTcy4txrQwuHdvet3d6xT0RfRpftzuBAg5yHhJlNBh4ArnD3d9ufd/fdid/3mNlq4DzgsJAQkb5TOKSn6Z0/Ufv8Ug688jts8NCshUO79hVejxg6qE/Dj4/L4C51OQ0JMxsHrAI+5e6vJj0/HChx9/rE158EFuWoTJGCoo7p3nUNh5EXzGbkuTMoLR+Z9Vqa25x9B/o2gfGSU8dk7PXDHgK7BLgYGG1mu4A7gTIAd78f+BZwFPAjM4NDQ12PAVYnnhsEPOzuT4ZZq0gxuCNWrYDoQdPeN4MO6TTCocR6Xv47l379yjsZ+15hj26a28vxvwP+rpvndwJTwqpLpJhoQlzvmva+GbQc/vDbRDjcwMhzr+2x5RDVgAAy+m+d8z4JEQmPltPoWfPet6h5fikH/vBsEA7nzwrCYVhFrksbkFKz3k9Kk0JCpACp9dCzTuFQNqRgwqFdJncHVEiIFJip336at+ubcl1GJDW/mwiHl9vD4XpGnntdwYRDu7GFMrpJRDJr8p1PUndQt5a6an53FzXPLzkUDlOvZ+R5hREOXRcHLC8r7ZhbkQkKCZE8p1tLqTW/u4va55ey/w/PYoMGF0w4JK/R1P7vv7umkeNC2FtCISGSp2JVce56bFufx9AXg+Z3d1H7wjL2v/zf2KAyRp53HSPPm5n34VBeVsp3Z07qFALXnjU2tGXCQSEhkpc0Y7p7ze/Fg5ZDezice20QDsNH5bq0AQt7B7pUFBIieSZWFVdAdFHI4dBd6yGbFBIieUJ7TB+u+b14cFtp22+w0sIKB8hd6yGZQkIkD8Sq4ixYsYXmKE/zzaLmfbuDlkN7OFTOYOTUmZQOPzLXpWVU2PtXp0MhIZIHFqzYTB83KStIQTgsY/+2Xxd0OEBm5zoMhEJCJII0rLWzIByWs3/bM1hpGSMqr6HivOspPaLwwgEyP9dhIBQSIhGiYa2dNe/7c6Ll8AxWOogR50ynYuqsgguHUYld52oONIcy12EgFBIiEaFhrYc07/tz0CH9UmGHAwS3laLQ95CKQkIkx2JVcb68fDOt6pOmueYvQcvhpV91hMPIqdcz6IgP5Lq0UETptlIqaYWEmX0Y+DFwjLufkdhy9Bp3vyfU6kQKnFoPgeRwoKSUEWdfzcjzZxVkOAwrK6GxuS1yt5VSSbcl8VNgAfATAHffamYPAwoJkX7SpLggHOpeWE7DS78CKwnCYer1DBpxVK5Ly7hR5WUsvGZi5EOhq3RDYpi7v2idN7JoCaEekYIWbAK0lcYiH896WDicdSUjp84qyHA4clgZd07Pv3Bol25I7DWzk0msSGtms4A/93aRmT0IXA3scfczujluwPeBK4EDwN+6++8Txy5PHCsFHnD3e9OsVSSSdGsJWmrfpvb5ZQUbDvnaWuhJuiHxv4HFwKlmFgf+CNycxnU/A34A/DzF8SuACYlfUwn6PaaaWSnwQ+BSYBewwczWuPvLadYrEinFvs9DS+3b1L6wnIbq/0oKh+sZNGJ0rkvLmJvPH8c9107KdRkZl1ZIuPtO4BNmNhwocff6NK971szG93DKDODn7u7AOjMbZWbHAuOBHYnXxcyWJs5VSEjeKeaAaKndQ+0LyxLhYIw484qgQ7oAwqHUjFb3SKyvFKZ0Rzd9Cfg3oB74qZmdDdzu7k8N8PXHAm8lPd6VeK6756emqG0+MB9g3LhxAyxHJHPuiFXz0Po3yeB2w3kjCIf2lgNBOEydxaCR+R8OuV6VNdvSvd30P9z9+2Z2GXA08BmC0BhoSFg3z3kPzx/+pPtiglthVFZWFuHbUaKk2JfTaKndQ+265TRsbQ+Hyxg59Ya8D4fhg0s50NSaN8NWMyndkGj/0L4S+Dd332Jdhjr10y7ghKTHxwO7gcEpnheJrGLumG6pS7QcEuFwxJTLqDh/FoNGjsl1af1SLLeS0pFuSGwys6eAk4Cvm9kIIBNj+NYAX0j0OUwFat39z2b2DjDBzE4C4sAcYF4GXk8kFDf99AWee/29XJeRdUE4rKBh69MFEQ6F2vk8EOmGxGeBM4Gd7n7AzI4iuOXUIzNbAlwMjDazXcCdQBmAu98PPE7QOtlBMAT2M4ljLWb2BWAtwRDYB919Wx/+XCJZU4wB0VL3DrXrVtCwJbjjfMSUT1Jx/g15Gw4G3KSA6FaPIWFmp7r7KwQBAfChvtxlcve5vRx3guG13R17nCBERCJr6ref5u36plyXkTWHhcPkS6m44AYGjTw6x5Wlb1hZCUPKSiO54moU9daSuI1g5NA/dnPMgeguXSgSgmJdyrulbm8QDlvXgudnOLxx71W5LiEv9RgS7j4/8fsl2SlHJLqKsWO6pW4vdetXUL+lPRw+QcX5sxlUkT/hAMFMaOmfdOdJlAGfBz6aeOo3wE/cvbh+nJKiVWyL8XUOB+eISZdScUH+hQNAWYmx8JqJuS4jb6Xbcf1jgg7nHyUefyrx3N+FUZRIlFz6T7/htT37c11GVrTU76Vu3UrqtzyZFA43MKjimFyX1i8awjpw6YbEue4+JenxM2a2JYyCRKLk1G8+zvtFsBtQS/1e6tY/Qv3mJ8HbOGLSJxIth/wIh2knf4A33m1kd02jOqMzLN2QaDWzk939dQAz+xBQnIvRSEGLVcVZuGYbNY3FcSe1pf5d6tavPBQOZ/w1FRfemDfhUFZi3HfDFAVCiNINiQXAr81sJ8GQ4hNJY56ESD4ppo7p7sJh5AWzKRv1wVyXlrZCXJY7itJdBfZXZjYBOIUgJF5x94OhViaSRcUSEC0N71G3biUNW57EW1s4YtIn8iYcNBs6N9JtSQCcQ7CE9yBgipnh7qn2iRDJG7GqOA8VeEB0DYfhidtK+RAO6nzOrXSHwP4HcDKwmUN9EU7qzYREIi9WFeerK7fQVMAd0y0N71G3/hEaNj+RN+GgUIiWdFsSlcDpiWU0RPJarCrOghWbKeRtplsb9lG7fmVSOHycigtupOzIY3NdWkq6nRRN6YbES8AHSWNfa5GoilXFuW3Z5owsXxxV+RgO6oCOtnRDYjTwspm9CHR0WLv7NaFUJZJhhd4x3dqwj9oXH6Gh6gm8tZnhEz9OxYWzKTvyuFyX1q3yshK+O3OygiEPpBsSC8MsQiQshT7voXX/PmrXJ4fDJUGfQ0TD4chhZdw5Xa2GfJLuENj/7um4mb3g7hdkpiSRzCjkfR5a9++jbv0q6qsePxQOF8ym7APR/PBVZ3T+6ssQ2J4MzdD3EcmIQl1vqXV/DXUvrqL+979MhMPFQZ9DhMLhyGFl2quhgGQqJDTqSXKukPd66AiHql/iLc0MP/1jVFw4J1LhoH6GwpSpkBDJqYJvOUQ4HAaXGv8wS+snFap0J9N9AXjI3felOiVzJYn0zeQ7n6TuYGGtN9l6oDZxW+k/D4XDBTdSdtTxuS6tQ1kJ3HfDmQqHApduS+KDwAYz+z3wILC2y8S6T6W60MwuB74PlAIPuPu9XY4vAG5Kquc0YIy7v2dmbwD1BLO8W9y9Ms16pcAV6q2lw8LhtI8Go5WOOiHXpTHh6OE8fdvFuS5DsszSnURtZgZ8kmD110pgOfCv7cuHp7imFHgVuBTYBWwA5rr7yynOnw7c6u4fTzx+A6h0973p1FhZWekbN25M688j+StWFefLK7bQ2lY4XWFBOKwOwqH5YOK2UjTCQS2Gwmdmm1L9EJ52n4S7u5n9BfgL0AIcCaw0s6fd/aspLjsP2OHuOxOFLAVmAN2GBDAXWJJuTVKc7npsW8EEROuBWuo2rKZ+UxAOw07/KKMunBOJcNCcBoH0+yS+CPwNsBd4AFjg7s1mVgK8BqQKibHAW0mPdwFTU7zGMOBy4AtJTzvwlJk5wZ7ai7u5bj4wH2DcuHHp/HEkDxXanIfDwuG0RDiMVjhItPRlWY6Z7v6n5Cfdvc3Mru7huu46tFP9CDgdeM7dkz8Jprn7bjM7GnjazF5x92e71LAYWAzB7abe/iCSfwpp5FJrY92h20pN70cmHDTZTVJJd8b1t3o49oceLt0FJP/vPx7YneLcOXS51eTuuxO/7zGz1QS3r57t5lopUDf99IWCCIjWxjrqNsSo3/RYIhw+QsWFcxg8OretX628Kr0Je57EBmCCmZ0ExAmCYF7Xk8ysAvgYcHPSc8OBEnevT3z9SWBRyPVKhBRCC+KwcDj1Iiqmzc16OKjzWfor1JBw95bEHIu1BENgH3T3bWb2ucTx+xOnXgc85e7JnwjHAKuDQVUMAh529yfDrFeiIVYV55Zlm3NdxoC0NtYnwmHNoXC4cA6Dx5yY1To0C1oGKu0hsPlAQ2DzV6wqztdXbaUxz3cCam2sp35DjLpO4XAjg8eMz1oNuoUkfZWRIbAiYSiUSXGdw6GRYadcRMW0OVkNB1BASOYpJCRnpn77ad6ub8p1GQPS+n5DEA4b1+BNB3ISDuVlpXx35iTdUpJQKCQk6wphl7ggHB6lbuOjiXCYFnRIZykchg8u5UBTq5bjltApJCRrYlVxbl22Oa/Xlc91OJQYzJuqW0qSPQoJyYp8bz20vd9A3cZHg9tKB/cz7MMXBuFw9Emhv/a0kz/AQ/9TGz9KbigkJFSxqjjfWLWVA3k6aqn7cJjD4KM/lJXXV0BIrikkJDT53HoIwmFNcFspS+GgOQ0SRQoJCUWsKp6XAdF2cD91G9dQvyFG28H9lH/4AkZdOJfBx4QTDsPKSviOgkEiTCEhGZePK7YeFg4TzmfUtLkMPubkUF5PG/hIvlBISEbcEavmoXVv5t3IpWyHw/DBpXz7Os1pkPyhkJABiVXF+fLyzbTmWTq0HTxA3aZEOLzfEGo4aNiq5DOFhPRbPs57OCwc/moqFdPmMuSDf5Xx19IeDVIIFBLSZ/nZ53CA+k2PUbdhdajhYMBNWj9JCohCQtIWq4qzYMVm8mnKQ9vBA9T//j+pe3E1be/XU/5X51ExbV7Gw0EL60mhUkhIWmJVcW5bvpm2PLm3dFg4nHwuFRfdlNFw0PBVKQYKCUkpVhVn4Zpt1DTmzzLe3YbDtHkMOXZCxl5Ds6ClmCgkpFvBraUtNOdJ06Ht4AHqq34ZhENjXSIc5jLk2A9n5PuPKi9j4TUT1WqQoqOQkG7lS9/DYeHwoUoqLpqXkXDQ6CSRLISEmV0OfJ9gj+sH3P3eLscvBh4F/ph4apW7L0rnWsm8fOmcbmtqpP73v6TuxVWHwmHaXIYcd8qAvq/mNIh0FmpImFkp8EPgUmAXsMHM1rj7y11O/a27X93PayUD8mVYa1tTY9ByWB+Ew9APncOoafMGFA5aWE8ktbBbEucBO9x9J4CZLQVmAOl80A/kWklTXrUcMhQO2u5TJH1hh8RY4K2kx7uAqd2cd4GZbQF2A19x923pXmtm84H5AOPGjctQ2cUhH1oPbU3vJ8LhkSAcTjqHUdPmMmTsqf36fuqAFumbsEPCunmu63CZ3wMnunuDmV0JxIAJaV6Luy8GFgNUVlbmx1CcCLj0n37Da3v257qMlIJweJy6Fx+h7UAtQ086OxEOp/Xr+6kTWqR/wg6JXcAJSY+PJ2gtdHD3uqSvHzezH5nZ6HSulb6L+kZAAw2HYWUlDCkrpeZAM8cpGEQGLOyQ2ABMMLOTgDgwB5iXfIKZfRB4293dzM4DSoB3gZrerpX0Rb3voa3pfRo2P07t+lW0Hahh6PizqJg2j6HH9x4O2ptBJDyhhoS7t5jZF4C1BMNYH3T3bWb2ucTx+4FZwOfNrAVoBOa4uwPdXhtmvYUqVhXntmWbiWI+tDW/T0NV13CYy9DjT+/1Wt1CEgmfBZ/HhaGystI3btyY6zIiJVYV55Zlm3NdxmGCcHiC2vWP9DkcAL5345kKB5EMMbNN7l7Z3THNuC5QUd0prq35fRo2P0nt+pW07a9h6IlnUnHRXIYePzGt60sN/nG2AkIkWxQSBSRWFeeux7ax70D0FuRraz5Iw+YnksJhChUzbmfoCWekdb2GrorkhkKiQMSq4ixYuYXmiO0jGoTDk9StX0nr/n19Dgf1O4jklkIij8Wq4ty3dju7axojeFvpIA1bnqRuXXs4TGb0jK+lFQ7awEckOhQSeSqqS3l3DYch49ILh+GDS/n2dVoqQyRqFBJ5KFYV59ZlmyPVegjCYW1wW6nhvSAcrvkqQ8f13CI4clgZd05XX4NIVCkk8kzUltPwlibqt6ylbt2KRDhMYvT0rzB03ORuzzfgj/deld0iRaTfFBJ55KafvhCZgDgsHE44o8dw0JafIvlJIZFHorBia7rhMGRQCU0tbVo/SSTPKSQi7I5YNQ+tf5MoTIr3liYatj5F7QsraG14lyHHT2T01V9hyLhJmB1asFfzGUQKi0Iiok795uO8H4E5D4eHw+kcdfVtDB03uSMcSs2YO/UEDVsVKUAKiQhpn/cQr2nMdSl4S3MiHJanDAf1M4gUPoVERERlIT5vaaah+ukgHOr3MmTs6Rx11a0MPXFKRzgYcJMmvIkUBYVEjkWl9dBtOFx5C0NPnMLwwaV8Z+Zk9TOIFCGFRA5FYZe4Q+Gwgtb6dxgy9rSOcCgx45+1JLdIUVNI5ECsKs5XV26hKYcd097aTMPWpHA47lSOuuKLDB1/JmammdAiAigksi7XrQdvbaah+r+C20p1QTiMv+5W7vvSp7ju7ONzVpeIRJNCIktyvQlQEA6/ovaFZYlwOIU7vvs97vz83E7zHEREkoUeEmZ2OfB9gn2qH3D3e7scvwn4WuJhA/B5d9+SOPYGUA+0Ai2ptteLsqiFw+BjT+HEGbfwf2/5tFoOItKrUEPCzEqBHwKXAruADWa2xt1fTjrtj8DH3H2fmV0BLAamJh2/xN33hllnWHK5GN+hcFhOa90ehhx7Ct/8zj+z8H/NU8tBRNIWdkviPGCHu+8EMLOlwAygIyTc/fmk89cBBfHjba4W4/PWFhpeSoRD7dsMPvbDjJv+Rf7xtr9Ry0FE+izskBgLvKZuupkAAApfSURBVJX0eBedWwldfRZ4IumxA0+ZmQM/cffFXS8ws/nAfIBx48YNuOCByOUe013DYcIZZ/L9Jf/G5ZdfrpaDiPRb2CHR3adTt7fnzewSgpC4KOnpae6+28yOBp42s1fc/dlO3ywIjsUAlZWVObn1n9vbSi00vPQMdS8so6X2bc4991zuukvhICKZEXZI7AJOSHp8PLC760lmNhl4ALjC3d9tf97ddyd+32NmqwluXz3b9fpcmvrtp3m7vinrr+utLezf9gz165bTtO8vnHvuuSxc+CBXXHGFwkFEMibskNgATDCzk4A4MAeYl3yCmY0DVgGfcvdXk54fDpS4e33i608Ci0Kut09iVfGsB4S3tuCvPUvTxpW8G3+TyspKFi58gCuvvFLhICIZF2pIuHuLmX0BWEswBPZBd99mZp9LHL8f+BZwFPCjxIdc+1DXY4DViecGAQ+7+5Nh1puubN9eMuC+mROp3/Zr7rnnHnbu3EllZSX/+pMfKRxEJFTmUdjRJkMqKyt948aNob5GtgOizNq4YugOHv+PH/L6669zzjnnsHDhQq666iqFg4hkhJltSjUPTTOu05St1VqPGTGY9d+8lJaWFn7xi19wzz338C+vv87ZZ5/NmjVruPrqqxUOIpI1Cok0ZGO9pfbd3RZefRr//u//zt13383rCgcRyTGFRC9iVfFQAmJwqfEPs6Z0rLLa0tLCww8/zGmnzWTHjh2cddZZPProo0yfPl3hICI5o5DoRpi3lm7usqNbezjcfffd7NixgzPPPJNYLMY111yjcBCRnFNIdBHmNqLJAdHS0sKSJUu4++67ee211xQOIhJJCokuMhUQY0eVs+CyUw7btKelpYWlS5dy99138+qrrzJlyhRWr17NjBkzFA4iEjkKCTI3rLW98zn5dlK71tbWjpZDcjhcc801lJSUDPi1RUTCUPQhkamA+F6KvaBbW1tZunQpixYt4tVXX2Xy5MmsWrWKGTNmKBxEJPKK/lNqIAFhBLeVuguI1tZWHnroISZOnMjNN9/M0KFDeeSRR6iqquK6665TQIhIXij6lkRfDSsr4TszJ3fbaoAgHJYtW8aiRYvYvn07kyZN4pFHHuHaa69VMIhI3lFIpKmn/gYIwmH58uUsWrSIV155hTPOOIOVK1eq1SAiea3oQ2LC0cO7veU04ejhPH3bxb1e3104rFixgpkzZyocRCTvFf2n2NO3XcyEo4d3ei6dgGjvkJ40aRLz5s1j0KBBrFixgi1btjBr1iwFhIgUhKJvSQBptRjatba2snLlShYtWsTLL7/MxIkTWb58Oddff72CQUQKjj7V0tTW1sayZcuYPHkyc+bMAWDZsmVs3bqVG264QQEhIgVJn2y9aGtrY/ny5YeFQ3V1NbNnz1Y4iEhB0ydcCm1tbaxYsYLJkydz44030tbWxtKlS9m6davCQUSKhj7pukgOh9mzZ9PW1saSJUuorq7mxhtvpLS0NNcliohkTeghYWaXm9l2M9thZrd3c9zM7P8ljm81s7PTvTaT2traWLlyJVOmTGH27Nkday1VV1czZ84chYOIFKVQRzeZWSnwQ+BSYBewwczWuPvLSaddAUxI/JoK/BiYmua1GfHWW29x1VVXUV1dzamnnsrDDz/M7NmzFQwiUvTCHgJ7HrDD3XcCmNlSYAaQ/EE/A/i5uzuwzsxGmdmxwPg0rs2I4447jpNOOonbb79dt5RERJKEHRJjgbeSHu8iaC30ds7YNK/FzOYD8wHGjRvXryJLS0t59NFH+3WtiEghC7tPortddDzNc9K5Fndf7O6V7l45ZsyYfpQoIiKphN2S2AWckPT4eGB3mucMTuNaEREJUdgtiQ3ABDM7ycwGA3OANV3OWQN8OjHK6Xyg1t3/nOa1IiISolBbEu7eYmZfANYCpcCD7r7NzD6XOH4/8DhwJbADOAB8pqdrw6xXREQ6s2BQUWGorKz0jRs35roMEZG8Ymab3L2yu2OacS0iIikpJEREJCWFhIiIpFRQfRJm9g7wpz5eNhrYG0I5YcinWiG/6s2nWiG/6lWt4clUvSe6e7cTzQoqJPrDzDam6rCJmnyqFfKr3nyqFfKrXtUanmzUq9tNIiKSkkJCRERSUkjA4lwX0Af5VCvkV735VCvkV72qNTyh11v0fRIiIpKaWhIiIpKSQkJERFIqipAYyD7buZBGvTcl6txqZs+b2ZRc1JmoJa19yM3sXDNrNbNZ2ayvmzp6rdfMLjazzWa2zcz+O9s1JtXR2/+DCjN7zMy2JGr9TC7qTNTyoJntMbOXUhyP2nust3qj9B7rsdak88J5j7l7Qf8iWEH2deBDBHtUbAFO73LOlcATBBsdnQ+sj3i9FwJHJr6+Ilf1plNr0nnPEKz4Oyvif7ejCLbIHZd4fHSEa/0G8PeJr8cA7wGDc1TvR4GzgZdSHI/MeyzNeiPxHkun1qT/L6G8x4qhJdGxz7a7NwHte2Un69hn293XAe37bOdCr/W6+/Puvi/xcB3Bhky5kM7fLcD/AR4B9mSzuG6kU+88YJW7vwng7rmqOZ1aHRhhZgYcQRASLdktM1GI+7OJ108lSu+xXuuN0Hssnb9bCPE9VgwhkWoP7b6eky19reWzBD+h5UKvtZrZWOA64P4s1pVKOn+3HwaONLPfmNkmM/t01qrrLJ1afwCcRrBjYzXwJXdvy055fRal91hf5fI91quw32Nhb18aBQPZZzsX0q7FzC4h+A98UagVpZZOrd8DvuburcEPvDmVTr2DgHOAvwbKgRfMbJ27vxp2cV2kU+tlwGbg48DJwNNm9lt3rwu7uH6I0nssbRF4j6Uj1PdYMYTEQPbZzoW0ajGzycADwBXu/m6WausqnVorgaWJ/7yjgSvNrMXdY9kpsZN0/y/sdff9wH4zexaYAmQ7JNKp9TPAvR7clN5hZn8ETgVezE6JfRKl91haIvIeS0e477FcdcZksdNnELATOIlDHYATu5xzFZ071V6MeL3jCLZ7vTDqf7ddzv8Zue24Tufv9jTgV4lzhwEvAWdEtNYfAwsTXx8DxIHROfz7HU/qjuDIvMfSrDcS77F0au1yXsbfYwXfkvAB7LMd4Xq/BRwF/Cjx00OL52DlyjRrjYx06nX3P5jZk8BWoA14wN17HHqYq1qBu4GfmVk1wYfv19w9J8tcm9kS4GJgtJntAu4EypJqjcx7DNKqNxLvsTRrDff1E+kjIiJymGIY3SQiIv2kkBARkZQUEiIikpJCQkREUlJIiIhISgoJERFJSSEhIiIpKSREQmRmd5vZl5Ief9vMvpjLmkT6QpPpREJkZuMJlh4/28xKgNeA8zzaawGJdCj4ZTlEcsnd3zCzd83sLIL1laoUEJJPFBIi4XsA+Fvgg8CDuS1FpG90u0kkZGY2mGBToDJggru35rgkkbSpJSESMndvMrNfAzUKCMk3CgmRkCU6rM8Hbsh1LSJ9pSGwIiEys9MJ9lD4lbu/lut6RPpKfRIiIpKSWhIiIpKSQkJERFJSSIiISEoKCRERSUkhISIiKf1/cvW4d9BBFBkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_noise = y * (1 + (np.random.random(y.shape) - 0.5) * 0.5) + 0.1\n",
    "plt.scatter(y, y_noise)\n",
    "plt.plot((y.min(), y.max()), (y.min(), y.max()), 'k-')\n",
    "plt.xlabel('y')\n",
    "plt.ylabel('y_noise')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example P_Fun\n",
    "\n",
    "This is an example physics loss function that supplements the normal y_predicted vs. y_true neural network loss function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_fun_pythag(y_predicted, y_true, p):\n",
    "    \"\"\"Example function for loss calculation using physical relationships.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y_predicted : tf.Tensor\n",
    "        Predicted y values in a 2D tensor based on x values in this batch.\n",
    "    y_true : np.ndarray\n",
    "        Known y values that were given to the PhyGNN fit method.\n",
    "    p : np.ndarray\n",
    "        Supplemental physical feature data that can be used to calculate a \n",
    "        y_physical value to compare against y_predicted. The rows in this \n",
    "        array have been carried through the batching process alongside y_true \n",
    "        and the features used to create y_predicted and so can be used 1-to-1 \n",
    "        with the rows in y_predicted and y_true.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    p_loss : tf.Tensor\n",
    "        A 0D tensor physical loss value.\n",
    "    \"\"\"\n",
    "    \n",
    "    p = tf.convert_to_tensor(p, dtype=tf.float32)\n",
    "    y_physical = tf.sqrt(p[:, 0]**2 + p[:, 1]**2)\n",
    "    y_physical = tf.expand_dims(y_physical, 1)\n",
    "    \n",
    "    p_loss = tf.math.reduce_mean(tf.math.abs(y_predicted - y_physical))\n",
    "    \n",
    "    return p_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NN Model Structure\n",
    "\n",
    "Here we define the model layers using a simple list of kwargs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hidden_layers = [{'units': 64},\n",
    "                 {'activation': 'relu'},\n",
    "                 {'units': 64}, \n",
    "                 {'activation': 'relu'},\n",
    "                ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train the Model without P_Fun\n",
    "\n",
    "Here we train the model with loss weights (1.0, 0.0) which fully weights the mean absolute error of y_predicted vs. y_noise and does not weight the p_fun calculation at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From C:\\Users\\GBUSTER\\AppData\\Local\\Continuum\\anaconda3\\envs\\mlclouds\\lib\\site-packages\\tensorflow_core\\python\\ops\\array_grad.py:563: _EagerTensorBase.cpu (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use tf.identity instead.\n"
     ]
    }
   ],
   "source": [
    "PhysicsGuidedNeuralNetwork.seed(0)\n",
    "model = PhysicsGuidedNeuralNetwork(p_fun=p_fun_pythag, \n",
    "                                   hidden_layers=hidden_layers, \n",
    "                                   loss_weights=(1.0, 0.0), \n",
    "                                   n_features=2, n_labels=1)\n",
    "model.fit(x, y_noise, p, n_batch=4, n_epoch=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXyU1b348c83kz2ZhJCEJCRAWEWQsIhoRXHBDUVRqxZr3ar1avW69LaVe29r1eqt16r1en9Wqq3W9tq6W5fiUqkK1I2lEBZBFlkCCEmArGSb+f7+eJ6EIUzCTDKTAfJ9v17zmuc5z3meOfMQ8s055znniKpijDHGhCou1gUwxhhzeLHAYYwxJiwWOIwxxoTFAocxxpiwWOAwxhgTlvhYF6An5OTkaHFxcayLYYwxh5XFixdXqGpu+/ReETiKi4tZtGhRrIthjDGHFRHZFCzdmqqMMcaExQKHMcaYsFjgMMYYE5Ze0cdhjOk5zc3NlJWV0dDQEOuimBAlJydTVFREQkJCSPktcBhjIqqsrAyv10txcTEiEuvimINQVSorKykrK2Pw4MEhnWNNVcaYiGpoaCA7O9uCxmFCRMjOzg6rhmiBwxgTcRY0Di/h/ntZ4OjMmndg/iOxLoUxxhxSLHB0ZsOHFjiMMaYdCxyd8eZBUw001cW6JMaYEO3Zs4df//rXYZ937rnnsmfPnk7z3HXXXbz//vtdLVpQ6enpEb1eT7DA0Zn0fOe95uvYlsMYE7KOAofP5+v0vDlz5tCnT59O89x7772cccYZ3SrfkcAex+2MN895r/kasofGtizGHIbueXMlq7ZVR/Sao/pn8LPzR3d4fNasWaxfv55x48aRkJBAeno6BQUFLF26lFWrVnHhhReyZcsWGhoauO2227jhhhuAfXPa1dbWMm3aNE466SQ+/vhjCgsLef3110lJSeGaa65h+vTpXHLJJRQXF3P11Vfz5ptv0tzczEsvvcTIkSMpLy/n29/+NpWVlRx33HG88847LF68mJycnE6/l6ry4x//mLfffhsR4Sc/+Qnf+ta32L59O9/61reorq6mpaWFJ554ghNPPJHrrruORYsWISJ897vf5Y477ojofe6M1Tg601rjqLUahzGHiwceeIChQ4eydOlSfvnLX/L5559z//33s2rVKgCefvppFi9ezKJFi3jssceorKw84Bpr167l5ptvZuXKlfTp04dXXnkl6Gfl5OSwZMkSbrrpJh566CEA7rnnHk4//XSWLFnCRRddxObNm0Mq96uvvsrSpUtZtmwZ77//Pj/60Y/Yvn07f/rTnzj77LPbjo0bN46lS5eydetWVqxYwfLly7n22mu7eLe6xmocnfG2NlXtiG05jDlMdVYz6CmTJk3ab2DbY489xmuvvQbAli1bWLt2LdnZ2fudM3jwYMaNGwfAsccey8aNG4Ne++KLL27L8+qrrwKwYMGCtuufc845ZGVlhVTOBQsWcPnll+PxeMjLy+OUU05h4cKFHHfccXz3u9+lubmZCy+8kHHjxjFkyBA2bNjAv/7rv3Leeedx1llnhX5DIsBqHJ1JyQJPotU4jDmMpaWltW1/+OGHvP/++3zyyScsW7aM8ePHBx34lpSU1Lbt8XhoaWkJeu3WfIF5VLVL5ezovClTpjBv3jwKCwu58sor+cMf/kBWVhbLli3j1FNP5fHHH+f666/v0md2lQWOzog4zVVW4zDmsOH1eqmpqQl6rKqqiqysLFJTU1m9ejWffvppxD//pJNO4sUXXwTgvffeY/fu3SGdN2XKFF544QV8Ph/l5eXMmzePSZMmsWnTJvr168f3vvc9rrvuOpYsWUJFRQV+v59vfvOb/PznP2fJkiUR/x6dsaaqg/HmWY3DmMNIdnY2kydP5phjjiElJYW8vLy2Y+eccw6zZ8+mpKSEo446ihNOOCHin/+zn/2Myy+/nBdeeIFTTjmFgoICvF7vQc+76KKL+OSTTxg7diwiwoMPPkh+fj7PPvssv/zlL9s6+v/whz+wdetWrr32Wvx+PwC/+MUvIv49OiNdrVYdTiZOnKhdXgHw+Sugch3c/FlkC2XMEeqLL77g6KOPjnUxYqaxsRGPx0N8fDyffPIJN910E0uXLo11sQ4q2L+biCxW1Ynt81qN42C8+bBxQaxLYYw5TGzevJnLLrsMv99PYmIiTz31VKyLFHEWOA4mPR8a9kBzAyQkx7o0xphD3PDhw/nnP/+5X1plZSVTp049IO/cuXMPeKLrcBDVwCEi5wD/A3iA36rqA+2OXwHc6e7WAjep6jL32EagBvABLa3VJRHpC7wAFAMbgctUNbTep65ofSS3dgdkDYraxxhjjlzZ2dmHRXNVqKL2VJWIeIDHgWnAKOByERnVLttXwCmqWgL8HHiy3fHTVHVcuza2WcBcVR0OzHX3oycwcBhjjInq47iTgHWqukFVm4DngRmBGVT144DawqdAUQjXnQE8624/C1wYofIGl9467cj2qH6MMcYcLqIZOAqBLQH7ZW5aR64D3g7YV+A9EVksIjcEpOep6nYA971fsIuJyA0iskhEFpWXl3fpCwA2etwYY9qJZh9HsCWlgj77KyKn4QSOkwKSJ6vqNhHpB/xNRFar6rxQP1xVn8Rt+po4cWLXnzlOzQHx2FgOY4xxRbPGUQYMCNgvAra1zyQiJcBvgRmq2jbbmKpuc993Aq/hNH0B7BCRAvfcAmBnVErfKi7Oaa6yGocxR6TW9TC2bdvGJZdcEjTPqaeeysHGgj366KPU19e37Yeyvkc4rrnmGl5++eWIXa87ohk4FgLDRWSwiCQCM4E3AjOIyEDgVeBKVf0yID1NRLyt28BZwAr38BvA1e721cDrUfwODhs9bswRr3///t36xdw+cISyvsfhKmpNVaraIiK3AO/iPI77tKquFJEb3eOzgbuAbODX7mLprY/d5gGvuWnxwJ9U9R330g8AL4rIdcBm4NJofYc26flQVRb1jzHmiPP2LPh6eWSvmT8Gpj3Q4eE777yTQYMG8f3vfx+Au+++GxFh3rx57N69m+bmZu677z5mzNjvWR02btzI9OnTWbFiBXv37uXaa69l1apVHH300ezdu7ct30033cTChQvZu3cvl1xyCffccw+PPfYY27Zt47TTTiMnJ4cPPvigbX2PnJwcHnnkEZ5++mkArr/+em6//XY2btzY4bofBzN37lx++MMf0tLSwnHHHccTTzxBUlISs2bN4o033iA+Pp6zzjqLhx56iJdeeol77rkHj8dDZmYm8+aF3OLfoaiO41DVOcCcdmmzA7avBw6Y1lFVNwBjO7hmJXDgSJpo8uZB2cIe/UhjTNfMnDmT22+/vS1wvPjii7zzzjvccccdZGRkUFFRwQknnMAFF1yA+8fpAZ544glSU1MpLS2ltLSUCRMmtB27//776du3Lz6fj6lTp1JaWsqtt97KI488wgcffHDAgk2LFy/mmWee4bPPPkNVOf744znllFPIyspi7dq1/PnPf+app57isssu45VXXuE73/lOp9+voaGBa665hrlz5zJixAiuuuoqnnjiCa666ipee+01Vq9ejYi0NZPde++9vPvuuxQWFkas6cxGjociPR/qK8DXDJ6EWJfGmMNHJzWDaBk/fjw7d+5k27ZtlJeXk5WVRUFBAXfccQfz5s0jLi6OrVu3smPHDvLz84NeY968edx6660AlJSUUFJS0nbsxRdf5Mknn6SlpYXt27ezatWq/Y63t2DBAi666KK26d0vvvhi5s+fzwUXXBDyuh+B1qxZw+DBgxkxYgQAV199NY8//ji33HILycnJXH/99Zx33nlMnz4dgMmTJ3PNNddw2WWXta0f0l02rXooWpeQrY1uP7wxJjIuueQSXn75ZV544QVmzpzJc889R3l5OYsXL2bp0qXk5eUFXYcjULDayFdffcVDDz3E3LlzKS0t5bzzzjvodTqbSDbUdT9CuV58fDyff/453/zmN/nLX/7COeecA8Ds2bO577772LJlC+PGjQu64mG4LHCEwlvgvFsHuTGHhZkzZ/L888/z8ssvc8kll1BVVUW/fv1ISEjggw8+YNOmTZ2eP2XKFJ577jkAVqxYQWlpKQDV1dWkpaWRmZnJjh07ePvtfUPPOloHZMqUKfzlL3+hvr6euro6XnvtNU4++eQuf7eRI0eyceNG1q1bB8Af//hHTjnlFGpra6mqquLcc8/l0UcfbZviZP369Rx//PHce++95OTksGXLls4uHxJrqgpF2+hxeyTXmMPB6NGjqampobCwkIKCAq644grOP/98Jk6cyLhx4xg5cmSn5990001ce+21lJSUMG7cOCZNckYDjB07lvHjxzN69GiGDBnC5MmT28654YYbmDZtGgUFBXzwwQdt6RMmTOCaa65pu8b111/P+PHjQ2qWCiY5OZlnnnmGSy+9tK1z/MYbb2TXrl3MmDGDhoYGVJVf/epXAPzoRz9i7dq1qCpTp05l7Nig3cdhsfU4QlG9DR45Gs57BI67LnIFM+YI1NvX4zhchbMehzVVhSKtHyA20aExxmBNVaHxxENaDtRYH4cxJrpuvvlm/vGPf+yXdtttt3HttdfGqEQHssARKm++1TiMCZGqdjhGwnTu8ccf7/HPDLfLwpqqQpWebzUOY0KQnJxMZWVl2L+MTGyoKpWVlSQnh77CqdU4QuXNgx0rDp7PmF6uqKiIsrIyurWcgelRycnJFBWFshySwwJHqNLdpiq/D+I8sS6NMYeshIQEBg8eHOtimCiypqpQefNB/VBXEeuSGGNMTFngCFXrIEAbPW6M6eUscISqddoRGz1ujOnlLHCEyms1DmOMAQscobP5qowxBrDAEbr4JEjJgprtsS6JMcbElAWOcKTb6HFjjLHAEQ6vjR43xhgLHOGw+aqMMSa6gUNEzhGRNSKyTkRmBTl+hYiUuq+PRWSsmz5ARD4QkS9EZKWI3BZwzt0islVElrqvc6P5HfaTnufUOGwOHmNMLxa1KUdExAM8DpwJlAELReQNVV0VkO0r4BRV3S0i04AngeOBFuDfVHWJiHiBxSLyt4Bzf6WqD0Wr7B3y5oO/Gep3QVp2j3+8McYcCqJZ45gErFPVDaraBDwPzAjMoKofq+pud/dToMhN366qS9ztGuALoDCKZQ2NjR43xpioBo5CIHBV9DI6/+V/HfB2+0QRKQbGA58FJN/iNm89LSJZ3S9qiLz5zrt1kBtjerFoBo5gq7gE7RwQkdNwAsed7dLTgVeA21W12k1+AhgKjAO2Aw93cM0bRGSRiCyK2PTOrYHDOsiNMb1YNANHGTAgYL8I2NY+k4iUAL8FZqhqZUB6Ak7QeE5VX21NV9UdqupTVT/wFE6T2AFU9UlVnaiqE3NzcyPyhUi3GocxxkQzcCwEhovIYBFJBGYCbwRmEJGBwKvAlar6ZUC6AL8DvlDVR9qdUxCwexHQc6srJaZCUobVOIwxvVrUnqpS1RYRuQV4F/AAT6vqShG50T0+G7gLyAZ+7a5P3KKqE4HJwJXAchFZ6l7yP1R1DvCgiIzDafbaCPxLtL5DUOl5Nu2IMaZXi+oKgO4v+jnt0mYHbF8PXB/kvAUE7yNBVa+McDHD4823iQ6NMb2ajRwPlzffHsc1xvRqFjjClZ7n1Dhs9LgxppeywBEubz607IXG6oPnNcaYI5AFjnC1PZJr/RzGmN7JAke4WpeQtSerjDG9lAWOcKXb6HFjTO9mgSNcNl+VMaaXs8ARriQvJKRajcMY02tZ4AiXyL4FnYwxpheywNEVtva4MaYXs8DRFel5NnrcGNNrWeDoCm+BjeMwxvRaFji6wpsHTTXQVBfrkhhjTI+zwNEVtqCTMaYXs8DRFa2jx+2RXGNML2SBoyvaahw27YgxpvexwNEVXpvo0BjTe1ng6IqULPAk2iO5xpheyQJHV4g4zVVW4zDG9EIWOLrKa4MAjTG9kwWOrrL5qowxvVRUA4eInCMia0RknYjMCnL8ChEpdV8fi8jYg50rIn1F5G8istZ9z4rmd+iQzVdljOmlohY4RMQDPA5MA0YBl4vIqHbZvgJOUdUS4OfAkyGcOwuYq6rDgbnufs/z5kPDHmhuiMnHG2NMrESzxjEJWKeqG1S1CXgemBGYQVU/VtXd7u6nQFEI584AnnW3nwUujOJ36JitBGiM6aWiGTgKgS0B+2VuWkeuA94O4dw8Vd0O4L73C3YxEblBRBaJyKLy8vIuFP8gvBY4jDG9UzQDhwRJ06AZRU7DCRx3hntuR1T1SVWdqKoTc3Nzwzk1NOnutCPWz2GM6WWiGTjKgAEB+0XAtvaZRKQE+C0wQ1UrQzh3h4gUuOcWADsjXO7Q2NrjxpheKpqBYyEwXEQGi0giMBN4IzCDiAwEXgWuVNUvQzz3DeBqd/tq4PUofoeOpeaAeGwshzGm14mP1oVVtUVEbgHeBTzA06q6UkRudI/PBu4CsoFfiwhAi9u8FPRc99IPAC+KyHXAZuDSaH2HTsXFuWM5rI/DGNO7RC1wAKjqHGBOu7TZAdvXA9eHeq6bXglMjWxJu8hGjxtjeiEbOd4dNl+VMaYXssDRHVbjMMb0QhY4uiM9H+rKwdcc65IYY0yPscDRHW2DAGPzRLAxxsSCBY7uaAsc1lxljOk9LHB0R9vocesgN8b0HhY4usNqHMaYXsgCR3ek9QPEph0xxvQqFji6wxMPaTkWOIwxvYoFju7y5tvU6saYXsUCR3el2xKyxpjeJaTAISJpIhLnbo8QkQtEJCG6RTtMePOsxmGM6VVCrXHMA5JFpBBnne9rgd9Hq1CHlfR8ZwCg3xfrkhhjTI8INXCIqtYDFwP/q6oXAaOiV6xDR0PzQQKCNx/UB3UVPVMgY4yJsZADh4h8A7gC+KubFtUp2Q8Fd7+xkqkPf9R5JhvLYYzpZUINHLcD/w685i7GNAT4IHrFOjQU9klh6569VNQ2dpwpvXUJWevnMMb0DiHVGlT1I+AjALeTvEJVb41mwQ4FY4oyAVi+tYrTjuoXPJPXnXbEahzGmF4i1Keq/iQiGSKSBqwC1ojIj6JbtNg7pjATEVheVtVxJpuvyhjTy4TaVDVKVauBC3GWcx0IXBm1Uh0i0pPiGZqbTmnZno4zxSdBShbUbO+5ghljTAyFGjgS3HEbFwKvq2ozoNEr1qGjpDCT0s5qHOA+kms1DmNM7xBq4PgNsBFIA+aJyCCgOlqFOpSMKcpkZ00jO6obOs7ktdHjxpjeI6TAoaqPqWqhqp6rjk3AaQc7T0TOEZE1IrJORGYFOT5SRD4RkUYR+WFA+lEisjTgVS0it7vH7haRrQHHzg3j+4atxO0g77TWYfNVGWN6kVA7xzNF5BERWeS+HsapfXR2jgd4HJiGM1jwchFpP2hwF3Ar8FBgoqquUdVxqjoOOBaoB14LyPKr1uOqOieU79BVowoy8cQJyzvr50h3px3RXtF6Z4zp5UJtqnoaqAEuc1/VwDMHOWcSsE5VN6hqE/A8MCMwg6ruVNWFQHMn15kKrHdrOT0uJdHD8H7pLDtYjcPXBHt391zBjDEmRkINHENV9WduENigqvcAQw5yTiGwJWC/zE0L10zgz+3SbhGRUhF5WkSygp0kIje01pDKy8u78LH7lBRlsnxrFdpRjaLtkVx7ssoYc+QLNXDsFZGTWndEZDKw9yDnSJC0sNpyRCQRuAB4KSD5CWAoMA7YDjwc7FxVfVJVJ6rqxNzc3HA+9gBjivqwq66JrXs6+MreAufdOsiNMb1AqPNN3Qj8QUQy3f3dwNUHOacMGBCwXwRsC694TAOWqGpbz3Pgtog8BbwV5jXDVlLojiAvq6IoK/XADG2jx62D3Bhz5Av1qaplqjoWKAFKVHU8cPpBTlsIDBeRwW7NYSbwRpjlu5x2zVQiUhCwexGwIsxrhm1kgZcEj3Tcz9E2X5XVOIwxR76wZrh1R4+3+gHwaCd5W0TkFuBdwAM87U6QeKN7fLaI5AOLgAzA7z5yO0pVq0UkFTgT+Jd2l35QRMbhNHttDHI84pLiPYzMz2D51g6erEpMhaQMq3EYY3qF7kyNHqwPYz/uo7Jz2qXNDtj+GqcJK9i59UB2kPSYTHUypiiTN5dtQ1URCfLV0/OsxmGM6RW6s+Z4rxq0UFKYSU1DC5sq64NnsNHjxpheotMah4jUEDxACJASlRIdokqK+gBQurWK4pwgYx+9+VC2sIdLZYwxPa/TGoeqelU1I8jLq6pH/AqAgYbnpZMUH0fplg76OdLznKnVbfS4MeYI152mql4lwRPHqP4ZlG7t4Mkqbz607IXGyM39WFXfzKptvWIuSWPMYcQCRxhKCjNZubUKnz9IrSIKS8g++O5qLvz1P6ja29mMLMYY07MscIRhTFEf6pp8fFVRe+BBb2SnHVFVPlxTTlOLn/dWWqe7MebQYYEjDGPdKdaXbQnSXNU67UiExnJsrKxvm+LkzVKbA8sYc+iwwBGGIbnppCZ6WB6sn6NtosPI1A4WrHUmZjx3TD7/WFfBrrqmiFzXGGO6ywJHGDxxwjH9M4OvQZ7khYTUiNU45q2tYEDfFG4+bRg+v/LOCmuuMsYcGixwhGlMUSYrt1XT4vPvf0AkYqPHm31+Pl1fyUnDchlVkMGQnDTeXBbu/JDGGBMdFjjCVFKUSWOLn7U7g3WQR2YJ2WVb9lDT2MKU4TmICNPH9uezryrZWdPJuufGGNNDLHCEqW0EebDmqvS8iDxVNW9tBXECJw7NAeD8kgL8Cm8vt+YqY0zsWeAI06C+qXiT4ykNNsW6tyAi4zgWrC2npKgPmakJAAzP83JUnteaq4wxhwQLHGGKixPGFGYGf7LKmwdNNdBU1+XrV+1tZumWPUwZnrNf+vljC1i0aTfbOlqF0BhjeogFji4YU5TJF9uraWzx7X8gAgs6fbK+Ar/CySP2X+52ekl/AOYstzEdxpjYssDRBWOL+tDsU778ul0HeQSWkJ2/toL0pHjGDeizX3pxThrHFGZYc5UxJuYscHTBGHcN8mXtO8gjUOOYv7aCE4Zkk+A58J/m/JL+LCurYnNHa4IYY0wPsMDRBUVZKWSlJrC8fQe5t3uBY1NlHZt31XNyu/6NVueVONOavFlqtQ5jTOxY4OgCEWFMUZ8Dp1hPyQJPEtR2LXDMX1sB0GHgKMpKZcLAPrxlc1cZY2LIAkcXlRRm8uWOGhqaAzrI20aPd62PY/7acgr7pDA42AqDrukl/fliezXrgg1ANMaYHhDVwCEi54jIGhFZJyKzghwfKSKfiEijiPyw3bGNIrJcRJaKyKKA9L4i8jcRWeu+Z0XzO3SkpCgTn19Z2X6hJW9el2ocLT4/H6+v5GR3tHhHzispQATesuYqY0yMRC1wiIgHeByYBowCLheRUe2y7QJuBR7q4DKnqeo4VZ0YkDYLmKuqw4G57n6Pax1BvvyADvKu1TiWlVVR09DCycNzO82Xl5HMpOK+vFW6HbVlao0xMRDNGsckYJ2qblDVJuB5YEZgBlXdqaoLgXCWuJsBPOtuPwtcGInChisvI4lcb9KB/Rze/C5NOzJ/bTkiMHlY9kHzTh/bn3U7a1mzoybszzHGmO6KZuAoBLYE7Je5aaFS4D0RWSwiNwSk56nqdgD3vV+wk0XkBhFZJCKLysvLwyz6wYkIJYWZwZ+satgDzeFNSLhgbQUlhZn0SU08aN5px+QTJ9iYDmNMTEQzcARrqA+nbWWyqk7Aaeq6WUSmhPPhqvqkqk5U1Ym5uZ03/3TVmKJM1pXXUtfYsi+xdSxHGIMAqxua+eeWPQdtpmqVk57E5GE51lxljImJaAaOMmBAwH4REPKfyKq6zX3fCbyG0/QFsENECgDc950RKW0XjC3qgyqsCGyu8oYfOD5dX4nPr5zUwWO4wUwvKWBTZT0rtlYfPLMxxkRQNAPHQmC4iAwWkURgJvBGKCeKSJqIeFu3gbOAFe7hN4Cr3e2rgdcjWuowHOOOIN9vwsMuLCE7f20FqYkeJgwM/QGxs0fnEx8nNhjQGNPjohY4VLUFuAV4F/gCeFFVV4rIjSJyI4CI5ItIGfAD4CciUiYiGUAesEBElgGfA39V1XfcSz8AnCkia4Ez3f2YyPUm0T8zef8p1rtQ45i/tpxvDMkmMT70f44+qYlMGZHLX0u34/dbc5UxpufER/PiqjoHmNMubXbA9tc4TVjtVQNjO7hmJTA1gsXsljFF7aZYT80B8YT8ZNWWXfVsrKzn6hOLw/7s6SUF/H31Tv65ZTfHDuob9vnGGNMVNnK8m0qK+vBVRR1Ve90niuPiwhrLsW+akfA78M8clUdifBxvLrMpSIwxPccCRzeVFDn9HPt3kIc+enz+2nL6ZyYzNLfjaUY64k1O4LSjcpmzfDs+a64yxvQQCxzd1DrF+n79HOn5IdU4fH7lH+sqOOkg04x0ZnpJf3bWNPL5V7u6dL4xxoTLAkc39UlNZGDfVJZvDZh6JMQaR2nZHqpDmGakM1OP7kdKgsfmrjLG9BgLHBEwpijzwBpHXTn4Op9JZcHaCneakdDHb7SXmhjP1KP78faKr2nx+bt8HWOMCZUFjggYW5RJ2e69VNY2Ogltj+R2PjZx/toKjumfSd+0g08z0pnpJf3ZVdfEx+sru3UdY4wJhQWOCBhT6M6U29pB3hY4Om6uqm1sYcnm3R0u2hSOU4/KJT0p3pqrjDE9wgJHBBxTmAGwb8LDttHjHXeQf7q+kpYwpxnpSHKCh7NG5fHOiq9parHmKmNMdFngiABvcgJDctP2TbEeQo1j/tpyUhI8HDuok2lGGmtg0yewZzP4fR3nA6aPLaC6oYX5ayM/E7AxxgSK6sjx3qSkMJNPNrh9DGn9AOm0xjF/XQUnDOlLUrwneIY9W+CPF0HlWmffkwh9BkLWYOg7BPoOdrcHQ59BnDQsl8yUBN4q3c7Uo/Mi++WMMSaABY4IKSnqw1+WbmNHdQN5GcmQltvhtCNlu+vZUF7HFccPCn6xnavh/y52ahwXPgEtjbD7K9j1lfO++VNoClzESUjM6M+rKbksXdWX5g8nk5DjBpecEZAY/uBCY4zpiAWOCGkdQb68rIq8UcnuWI7gNVtlkToAABrzSURBVI4F7jQjU4L1b2xZCH+6FOIS4Jq/QkHJgXlUob7SCSS7NrQFlZxtXzKlbjEJH/59X96MQrjpY0jp0+3vaIwxYIEjYkb1zyBOoHRrFWeMynNHjwfv45i/roK8jCSG9Uvf/8C69+GFKyG9H1z5mtMkFYwIpOU4rwHHtSWn+fyc/l9zOWVwKo+ckQnbl8Hr34ePH4Opd0XqqxpjejnrHI+Q1MR4hvfzsrzMHUHeQY2jdZqRk4fn7j/NyPKX4U/fgr5D4bvvdRw0OhHviWPamHzmrKmmLmskjL8CjrkEPvl1WOuDGGNMZyxwRFCJO4JcVZ0aR+3OA56GWrG1ij31zfuP3/h0NrxyHQw4Hq79qxN0umh6SX8amv3MXe0OPjz9P8HfDB892OVrGmNMIAscEVRSlEllXRPbqhqcR3LV5/RFBFiwzunfmDwsx+mr+Pt98M6dMHI6fOdVSM7sVhmOK+5LXkYSby1zBwP2HQLHXgNLnoXK9d26tjHGgAWOiBpT5I4gL9uzbyxHuyer5n1Zzuj+GeSkxsNbt8O8X8L4K+HSZyEhudtl8MQJ544p4MM15VQ3uHNlTfmx8zjv3+/r9vWNMcYCRwSNzPcSHyfOhIfprYFjXz9HnTvNyKlDM+Clq2Hx7+GkH8AF/wueyD2ncP7Y/jT5/PxtpfvZ3jw44fuw8lXYtjRin2OM6Z0scERQcoKHo/K9TuBo7acIGD3+2VeVJPnq+N7mH8MXb8LZ/wVn/Mx5SiqCxg/oQ2GflP3nrpp8K6Rkwdx7IvpZxpjexwJHhJUU9aG0bA+a1s9JCKhxLF75JS8k3Udm+SK46En4xs1RKYOIML2kgPlrK9hd1+QkJmfCyf8G6/8OGz6KyucaY3oHCxwRVlKUSXVDC5urfc5f+K01jt0buXzl9xgWtx2Z+WcY+62oluP8sf1p8Svvrgx4DPe47zkDAufe43TMG2NMF0Q1cIjIOSKyRkTWicisIMdHisgnItIoIj8MSB8gIh+IyBcislJEbgs4dreIbBWRpe7r3Gh+h3Dtt5Sst8AZP7FjJb7fnkW6r5p3j/0NjDgr6uUY3T+D4uxU3ioN6JxPSIZT/x22LnaayowxpguiFjhExAM8DkwDRgGXi8iodtl2AbcCD7VLbwH+TVWPBk4Abm537q9UdZz7mhOdb9A1I/K8JMbHOWtzpOfB1iXwzDQafXBp0884auIZPVIOEeHC8YUsWFfBoo0B65GPvRxyjoK594KvpUfKYow5skSzxjEJWKeqG1S1CXgemBGYQVV3qupCoLld+nZVXeJu1wBfAIVRLGvEJMbHcXRBBsu2uI/k1myDtFwe6P8/VKUPZURe+sEvEiHfO3kIhX1SmPXqchpb3IGInniY+lNn1t1lf+qxshhjjhzRDByFwJaA/TK68MtfRIqB8cBnAcm3iEipiDwtIkEXtBCRG0RkkYgsKi/v2TUqxhZlsmJrFf5hZ8Hws/Ff8w5vbYrnpOE5+08zEmVpSfHcd9ExrNtZy68/CBj8N3I6FE6ED34BzXt7rDzGmCNDNANHsN+QYfXIikg68Apwu6pWu8lPAEOBccB24OFg56rqk6o6UVUn5ubmhvOx3TamMJO6Jh8b8s6CK15kVXUiu+qamDK8Z8sBcNpR/bhwXH9+/eE6vtzhTsUuAmfc7dSGPn+qx8tkjDm8RTNwlAEDAvaLgJAXxRaRBJyg8Zyqvtqarqo7VNWnqn7gKZwmsUNKSesI8q3OhIfz3FX5Jg/r/jKxXfHT6aNIT4pn1iul+P1u7B58Mgw7A+Y/DHv3xKRcxpjDUzQDx0JguIgMFpFEYCbwRignitOe8zvgC1V9pN2xgoDdi4AVESpvxAzNTSMlwcOyLc5SsvO/rODoggxyvUkxKU92ehJ3nT+KJZv38H+fbdp3YOrPoGGPM+26McaEKGqBQ1VbgFuAd3E6t19U1ZUicqOI3AggIvkiUgb8APiJiJSJSAYwGbgSOD3IY7cPishyESkFTgPuiNZ36Kp4TxzHFGawfGsV9U0tLN60O/iiTT3ownGFTBmRy3+/vZpte9x+jYISm3bdGBO2qI7jUNU5qjpCVYeq6v1u2mxVne1uf62qRaqaoap93O1qVV2gqqKqJe0fu1XVK1V1jHvsAlUNvj5rjI0p7MPKbVV8vK6SJp+fk2IcOESE+y88Br/CT/+ywpn6HWzadWNM2GzkeJSUFGXS0OznmY+/Iik+juOK+8a6SAzom8q/nTWCuat37hsYaNOuG2PCZIEjSsa4a5D/Y10lkwb3JTnBE+MSOa6dPJixRZnc8+ZK9tS781i1Trv+wf2xLZwx5rBggSNKBmen4U1ypko/OcbNVIE8ccIvLi5hT30z9//1Cyexddr1Fa/YtOvGmIOywBElcXHCMe68VSfHYPxGZ0b1z+BfThnCS4vLWLDWWZFw37Tr98a2cMaYQ54Fjig6bWQuR+V5GZnvjXVRDvCvpw9ncE4a//HacvY2+QKmXZ8LX82LdfGMMYcwCxxRdMOUobx7x5QenWYkVMkJHn5x8Rg276rn0fe/dBJbp11//26bdt0Y0yELHL3YCUOyuXzSAJ6av4EVW6ts2nVjTEgscPRys6YdTXZ6Ene+UkqLzx/xadfVai7GHHEscPRymSkJ/HzGaFZuq+Z3C76K2LTrK7ZWceXvPmPkT9/h0tkf88Dbq5n7xY59jwAbYw5b8bEugIm9c44p4OzReTzyty85e3Q+xYHTro+5FBJSQr7W5sp6Hv7bGl5fuo0+qQlcPKGQ1V/X8Nv5G5j9kVP7GJGXzsTivhxXnMXEQX0pyko5JPuBjDHBSW9oSpg4caIuWrQo1sU4pO2obuCMhz9iTFEmz11/PLJxATw7Hc78ufOo7kFU1jbyv39fx3OfbcITJ3x38mD+5ZShZKYkALC3yceysj0s2riLhRt3s2TTbmoanaaw/IxkJhZncVxxXyYWZzEyPwNPnAUSY2JNRBar6sT26VbjMADkZSQz69yR/OdrK3hpcRmXTXSnXf/7fbDzCxj/HRh0orOWR4D6phZ+N/8rfjNvA/VNLXzruAHcNnUE+ZnJ++VLSfRwwpBsThiSDYDPr6z5uoZFm5xAsmjjrrZpUNKT4hk/sA/HFfdlyohcxg3o0zM3wRgTEqtxmDZ+vzLzyU9Zs6OGv/1gCv10N3z0ACx/BZpqnHmtxn8Hxn6b5rQ8Xli4hUffX0tFbSNnj87jR2cfxbB+XR+zsnXPXrdGsotFG3ezZkcNqvCvpw/jjjNGEGe1EGN6VEc1DgscZj/ry2uZ9uh8zhydx+PfnuAkNtXBqjfgn3+ETf9AieNTzwR+v/dkqgecxg/PLeHYQUFX8O06v5+aslU889EqHlmZxtSR/fjVzHFkJCdE9nOMMR2ywGGBI2T/7+9reei9L/ntVRM5Y1ReW/qnGyr5/ZtzGVP+JjMTFpCtu9DUHGTsTBh/JfQb2fUPrfnaGT9Stgi2LnLmzGp0VgveknsKV269mLjsYp66aiJDc9O7+xWNMSGwwGGBI2RNLX4u+H8LqNrbzHt3TKFs914efGc1H6wppyAzmTvOHME3x+Xj2fB3pxay5m3wt0DRcU5T1uiLITmjkw+ocwLD1kVusFgM1WXOsbh4yBvtPNVVNBFqd8JHD+Lz+3jCfxFP+6fz0OXHcfrIvI6vb4yJCAscFjjC8s/Nu7n4iY8ZmpvO+vJavEnx3HzaMK4+sfjAKeJry6H0BSeIlK+GhFQYdaETRAYcDxVr3JrEYue1cxWo3zm3zyAnQBQe6wSLgpIDH/+t2grv/jusep0tcUXMariKE8/8Jt8/dag9xmtMFFngsMARtv+a8wW//3gj104u5vunDCMz9SD9C6pOYPjnH/d1qMfFO7URcCZSbA0QRROh/wRID2Pm4LXv45/zQ+J2f8XrvhP5ZNgd/HTm6aQl2cOBxkSDBQ4LHGFTVRqa/aQkdmERqtYO9R0rIH+MEyyyhx7wOG/YmhvQBb/CN/8R9vo8/F/qd5h+3V0MyOmkaSwGquqb2bSrjpqGFkqKMvFap745DFngsMBxZKlcz66Xb6Pv9vmsYRCNZz9EyTfO6rGPV1XKaxrZtKuejRV1bN5Vz6bKejZV1rFpVz176pvb8nrihLFFmZw4NIcTh2UzYWDWIbMipDGdiUngEJFzgP8BPMBvVfWBdsdHAs8AE4D/VNWHDnauiPQFXgCKgY3AZaq6u7NyWOA4QqlS/vlL8M4scrWSNf0vYsQVDyNp2RG5vN+vbN2zl42VdfuCQmV9W5DY2+xry+uJEwr7pDAoO5WBfVMpzk5jYHYqKQkePv9qF/9YX0FpWRU+v5IUH8fE4ixOHJrD5GE5HNM/g3iPTRtnDj09HjhExAN8CZwJlAELgctVdVVAnn7AIOBCYHdr4OjsXBF5ENilqg+IyCwgS1Xv7KwsFjiObLXVu/n4dz/m9D0v0xDvJemce0k49iqIC++XcWOLjxVbq9pGsi/atHu/mkNSfBwD+6YyKDuNQdmp7iuNQX1TKcxKIeEgv/xrGpqdILKuko/XV7D66xoAvEnxHD8km8nDsjlxaA4j8tKt098cEmIROL4B3K2qZ7v7/w6gqr8IkvduoDYgcHR4roisAU5V1e0iUgB8qKpHdVYWCxxHPr9f+dNb7zBi0c+YFLeGpoKJJM541Olf6UDV3maWbNrdNlJ9adkemlqcp72G5KQxsTiL8QOzGJzjBIo8b3JER69X1DbyyfpKPl7vBJJNlfUA5KQnceLQbE4c6kzRMrBvqo2aNzERi7mqCoEtAftlwPERODdPVbcDuMGjX7ALiMgNwA0AAwcODKPY5nAUFyd854Jp/G3oeP7jxf/hh9v/j6zfTEHGXu488puUzm5fMuuqlFWVyrKdPr7YpdRoCg2SQnH/PK46YRAT3YkWc9KTwiuAqvP0WEsj+Jrc90bwNTvpfh+oz31XUB85fh/nZ/o4f5wPxirl1crqbXtYvW0PX66r4qPljczFw+b4YlLzhjKyIIOj8rwclZ/BUfle+qYlRudm9nINzT4+Xl/Beyt38PlXu4j3COlJ8aQlxe/33radHE96koe0RDc9eV8eb3I8qYlH3lN/0fxGwf5ECrV6051zncyqTwJPglPjCOdcc/g6c3Q+g2/+D6559kRmVv+eS0pfJdG/F4As4Dj3BUDg790KYHcSrEqHxHRIyoCkdIhPdn75+xrbBYUg7+H9iB4g132dHKR8tRXprCwvZmnLIP7sH8wKLaY2bRBHFWS6wcTLyPwMhuelR6Xjvdnnp7ahhdrGFmrc97rGFmoaW9z0ZmobnP26xn359jb5iPcICZ449+VsJ7r7rccS451j8XH7thM8caQlxjO60AmY0ewHqm5o5oPVO3lv5Q4+XLOTuiYf6UnxfGNoNnECdY0+qhta2F7V0Pb96hpb8IfwT56W6KFfRjL9vEnkBb5nJNHPm0xeRhL9MpJJP4weK49mScuAAQH7RcC2CJy7Q0QKApqqdna7pOaIMqyflz/ecg63P5/Hf6wpp783nskDkzm+fyLj8jwM9iqe5lpoqoXGGmisdcactG23ptc4255EJ5ikZjvb8UngSYL4xA7ekwLyJYLEQZwHxBPwHtduv4P05nr4ejnp25cxafsyJu14H/E1AtDgT2Hd9iEs3jSQT1oG8ZQO5iv6MyDbqZGMyPNSnJNKi09paPHT2OyjodlHQ7OfvQHbDS0+91hgunus2UdtYwuNbhNeZ0QgPbH1L3DnPSXBQ4tfqW1sodnnp8WnNPn8NPv8NLcozT4/TW56s89PSwe/iVMTPYwt6sOEQX2YMNBpQuxujWtndQPvrdrBe6t28Mn6Cpp9Sk56EheMK+Ts0Xl8Y2g2SfEdB2FVZa97f+oafW2Bta6xhbomZ7t6bwvlNY3srGlgZ3Ujy8r2sKO6gYbmA+9nRwFmQN8UBvZ1HrY4VIJLNPs44nE6uKcCW3E6uL+tqiuD5L2b/fs4OjxXRH4JVAZ0jvdV1R93Vhbr4+idVJU99c30SU04cjqbfc3O6Pzty9pe+vVypNnpH2mOS2JLwlBK/YP4pL6Itf5CABLwES8tJOAjQXykeZRUj5+UeD/JcX5SPM57cpyP5Dg/SXF+kuJ8JMb5SY5TZ98Die6xRFES4pQE8REvfhJE8eBD/D6naa61Wc7vcwKniBsM49xXwHbcvjSVOPwIfpz3Jr9QWd/CzlofO2qb2Vnno1nj8BNHekoyeZmp5GWlU5CVTm5GKnGeeGfQqcQ573EBv/hVKa9tZNW2KlZuq2bzrnoEJTstkVH9Mxld4GVA31TiAmuOIvv+GIhPcmqgB3v3JHX6YIaqUtPYws7qRnZWN7DDDSo7qvcFmB01DUEDTHZaIgP7pjAkK4HirAQGZSYwsE88hd44cpJB2mrHTU4t2NcE/cdDWk6Xftxi9TjuucCjOI/UPq2q94vIjQCqOltE8oFFQAbgB2qBUapaHexc95rZwIvAQGAzcKmq7uqsHBY4zBHN74PKdfuCybal8HVp2ySR3RKX4P4Cdn8JB3sXT7s8AfsS19ang/qdlz9gO2ha67a26x9qQf0+/L4W8Lcg6oSXQ5InoBYqcbS1vouEvK0oquBraULd5tA4fzPx2kw4ai55Hu8x07r0NWwAoAUO05v4/bD7K6hc7/z1G5cAnoR9723b8QFpic4v+9ZjcZ7uj/SPMvX72FxRw9LNlZRuqqR0yy427KxC/D7i8JOVLFQ3+JxBmAMymTKiH1NG5FKQGTgfWsB3bPu+0voB+/qxWhr2PfTQuh3Ke9vvWO1km+DpsO/fpq051Hm1SAJVTUJFA1TsVXbWKV/X+dle52drtZ/aljiaiOcH357OyWOGden+WuCwwGFMr1Df1MKyLVUs2bybjRV1TBrclzOOziOrFz2Fpm6z3JZd9Qzr521bwjlctnSsMaZXSE10nob6xtDIzCBwOBIR+nmT6edNPnjmLrB5DowxxoTFAocxxpiwWOAwxhgTFgscxhhjwmKBwxhjTFgscBhjjAmLBQ5jjDFhscBhjDEmLL1i5LiIlAObunh6Ds6k24cqK1/3WPm6x8rXfYdyGQepam77xF4ROLpDRBYFG3J/qLDydY+Vr3usfN13OJSxPWuqMsYYExYLHMYYY8JigePgnox1AQ7Cytc9Vr7usfJ13+FQxv1YH4cxxpiwWI3DGGNMWCxwGGOMCYsFDpeInCMia0RknYjMCnJcROQx93ipiEzowbINEJEPROQLEVkpIrcFyXOqiFSJyFL3dVdPlc/9/I0istz97AOWW4zx/Tsq4L4sFZFqEbm9XZ4evX8i8rSI7BSRFQFpfUXkbyKy1n3P6uDcTn9Wo1i+X4rIavff7zUR6dPBuZ3+LESxfHeLyNaAf8NzOzg3VvfvhYCybRSRpR2cG/X7122q2utfgAdYDwwBEoFlwKh2ec4F3sZZjPgE4LMeLF8BMMHd9gJfBinfqcBbMbyHG4GcTo7H7P4F+bf+GmdgU8zuHzAFmACsCEh7EJjlbs8C/ruD8nf6sxrF8p0FxLvb/x2sfKH8LESxfHcDPwzh3z8m96/d8YeBu2J1/7r7shqHYxKwTlU3qGoT8Dwwo12eGcAf1PEp0EdECnqicKq6XVWXuNs1wBdAYU98dgTF7P61MxVYr6pdnUkgIlR1HrCrXfIM4Fl3+1ngwiCnhvKzGpXyqep7qtri7n4KFEX6c0PVwf0LRczuXysREeAy4M+R/tyeYoHDUQhsCdgv48BfzKHkiToRKQbGA58FOfwNEVkmIm+LyOgeLRgo8J6ILBaRG4IcPyTuHzCTjv/DxvL+AeSp6nZw/lgA+gXJc6jcx+/i1CCDOdjPQjTd4jalPd1BU9+hcP9OBnao6toOjsfy/oXEAodDgqS1f045lDxRJSLpwCvA7apa3e7wEpzml7HA/wJ/6cmyAZNVdQIwDbhZRKa0O34o3L9E4ALgpSCHY33/QnUo3Mf/BFqA5zrIcrCfhWh5AhgKjAO24zQHtRfz+wdcTue1jVjdv5BZ4HCUAQMC9ouAbV3IEzUikoATNJ5T1VfbH1fValWtdbfnAAkiktNT5VPVbe77TuA1nCaBQDG9f65pwBJV3dH+QKzvn2tHa/Od+74zSJ5Y/xxeDUwHrlC3Qb69EH4WokJVd6iqT1X9wFMdfG6s7188cDHwQkd5YnX/wmGBw7EQGC4ig92/SmcCb7TL8wZwlft00AlAVWuzQrS5baK/A75Q1Uc6yJPv5kNEJuH821b2UPnSRMTbuo3TibqiXbaY3b8AHf6lF8v7F+AN4Gp3+2rg9SB5QvlZjQoROQe4E7hAVes7yBPKz0K0yhfYZ3ZRB58bs/vnOgNYraplwQ7G8v6FJda984fKC+epny9xnrj4TzftRuBGd1uAx93jy4GJPVi2k3Cq06XAUvd1brvy3QKsxHlK5FPgxB4s3xD3c5e5ZTik7p/7+ak4gSAzIC1m9w8ngG0HmnH+Cr4OyAbmAmvd975u3v7AnM5+VnuofOtw+gdafwZnty9fRz8LPVS+P7o/W6U4waDgULp/bvrvW3/mAvL2+P3r7sumHDHGGBMWa6oyxhgTFgscxhhjwmKBwxhjTFgscBhjjAmLBQ5jjDFhscBhzCHOnbn3rViXw5hWFjiMMcaExQKHMREiIt8Rkc/ddRR+IyIeEakVkYdFZImIzBWRXDfvOBH5NGBtiyw3fZiIvO9OtrhERIa6l08XkZfd9TCeax3lbkwsWOAwJgJE5GjgWzgT1I0DfMAVQBrO/FgTgI+An7mn/AG4U1VLcEY7t6Y/BzyuzmSLJ+KMPgZnRuTbgVE4o4snR/1LGdOB+FgXwJgjxFTgWGChWxlIwZmk0M++Ce3+D3hVRDKBPqr6kZv+LPCSO0dRoaq+BqCqDQDu9T5Xd34jd+W4YmBB9L+WMQeywGFMZAjwrKr++36JIj9tl6+zOX46a35qDNj2Yf93TQxZU5UxkTEXuERE+kHb+uGDcP6PXeLm+TawQFWrgN0icrKbfiXwkTprrJSJyIXuNZJEJLVHv4UxIbC/WoyJAFVdJSI/wVm5LQ5nVtSbgTpgtIgsBqpw+kHAmTZ9thsYNgDXuulXAr8RkXvda1zag1/DmJDY7LjGRJGI1KpqeqzLYUwkWVOVMcaYsFiNwxhjTFisxmGMMSYsFjiMMcaExQKHMcaYsFjgMMYYExYLHMYYY8Ly/wF1ffWVZjjxxgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU9dn/8fcN4lNEKCrqT5EIrah1wS2tVVtFrXV7qtWCG+1zPT5WBNSEVVZZRZFNQQWKllqrBXFptAVBXMC60IKyREAUNyD0qoKgrYCScP/+mAkdYzJLPOfMTObzui4vMjknM7e5Mrlzvt/P+X7N3RERkcLVKNsFiIhIdqkRiIgUODUCEZECp0YgIlLg1AhERArcXtkuIFOtWrXytm3bZrsMEZG88vrrr2929wNrO5Z3jaBt27YsXbo022WIiOQVM/uwrmMaGhIRKXBqBCIiBS60RmBmM8zsIzN7M8k5Hc1suZmtMrNFYdUiIiJ1C/OK4EHggroOmllLYApwibsfC3QOsRYREalDaI3A3V8CPklyyjXAk+6+Pn7+R2HVIiIidctmauhIoImZLQSaA5Pc/aHaTjSzrkBXgKKiosgKFBHJBWXLKhg3fy2btu3g0JZN6Xf+Ufz8pNaBPX82G8FewCnAuUBT4DUzW+zub9c80d2nA9MBiouLtVyqiBSMsmUVDHyynB27qgCo2LaDgU+WAwTWDLLZCDYCm939c+BzM3sJOAH4WiMQESk0Q8rKmfm3DVTVslXAjl1VjJu/NrBGkM346FPAj81sLzPbBzgVWJPFekREcsKQsnIeXry+1iZQbdO2HYG9XmhXBGY2E+gItDKzjcAwoAmAu09z9zVmNg9YCewGHnD3OqOmIiINXZf7X+OVd5NlbP7j0JZNA3vd0BqBu1+dxjnjgHFh1SAiki8yaQJNGhn9zj8qsNfOu7WGREQaiuohoIxZsHVoiQkRkSyodxMAdlU54+avDawWXRGIiEQokyGgZPJislhERP6jbFkF/R5bzq7dwTxf0ybBDeioEYiIhOzU0Qv457++DPQ5d1QG1FFQIxARCc03mQdIJcktBhlTIxARCcHRg+eysyq8FXEaW3DRITUCEZGAhHkFUNPVp7YJ7LnUCEREAhBVE2hsxtWntuG2nx8f2HOqEYiIfEPnTVzIOx99HvrrnPHd/Xnk+tMCf141AhGReqjeI6AiwDx/MmE1AVAjEBHJWFA3haUjzAZQTY1ARCRNsU1iVrIjqLvCUoiiCYAagYhIWqKaBwC4+8oTA92KMhU1AhGRJMqWVdDz0eWRvFbTJo244/IOkTYBUCMQEalVlA0Aor8KSBTaMtRmNsPMPjKzpLuOmdn3zazKzDqFVYuISCbOm7gw0quAVE3A3fnLX/7CunXrQqkhzP0IHgQuSHaCmTUG7gTmh1iHiEhahpSV03bAnMjuCfhgzMWsGXVh0iawdu1aLrroIn72s58xceLEUGoJc6vKl8ysbYrTbgaeAL4fVh0iIqkEvUR0Mo0Mrjm1KOWdwZ9++imjRo1i0qRJ7LPPPkyYMIGbbroplJqyNkdgZq2By4BzSNEIzKwr0BWgqKgo/OJEpGBEuT5QOvMAu3fv5sEHH2TgwIF8/PHHXHvttdx+++0cfPDBodWVzcniu4H+7l5lKVbRc/fpwHSA4uLi8JbzE5GCULasgkFPrmR7RPcDpDsRvHjxYkpKSliyZAmnnXYac+bMobi4OPT6stkIioFZ8SbQCrjIzCrdvSyLNYlIAxflXcGQXhPYtGkTAwYM4A9/+AOHHHIIDz30EF26dKFRo2i2lc9aI3D3dtUfm9mDwF/UBEQkDFGvCwTwyx+mngf44osvuPvuu7ntttv48ssvGTBgAIMGDaJ58+YRVRkTWiMws5lAR6CVmW0EhgFNANx9WlivKyKSaEhZOY8sXk+UY8qpmoC7M2fOHHr16sW6dev2JIKOOOKICKv8jzBTQ1dncO7/hlWHiBSuqIeB0rkKWLt2LT179mTevHkcddRRzJs3j/PPPz+iCmunO4tFpMEpW1ZB70eXE81UcEyqJlBXHHTvvfeOsMraqRGISINQtqyC4U+vYtuOXZG+bqrJ4Jpx0P/7v/9j9OjRocZBM6VGICJ5L8p7ARKlagLZioNmSo1ARPJa1PMATRrBuM7JG0DNOOgf/vAHunTpQqp7prJFjUBE8lI2rgJSbRSTK3HQTKkRiEheifoKAFLvE1C9Omjv3r1Zt24dl1xyCRMmTMhaHDRTagQikheycQXQyGDiFcmHgd566y169erFvHnzOProo3MiDpopNQIRyWllyyoY/KdyPv+yKtLXbdm0CcMvObbOJvDpp58ycuRIJk+ezD777MPEiRO56aabaNKkSaR1BkGNQERyUtmyCkb8eRVbt0cbB011P0BtcdDbb7+dgw46KMIqg6VGICI5JVsNwIC7UsRBX3vtNUpKSli6dCmnnXYac+fO5ZRTTomuyJCoEYhITshWA4DUVwGJcdBDDz2Uhx9+mGuuuSZn46CZUiMQkazL1g1hjQ0mJJkM/uKLL7jrrru47bbb2LVrFwMHDmTQoEHsu+++EVcaLjUCEcmabE0EpxsH7dWrF++++y6XXHIJEydO5Lvf/W6kdUZFjUBEIpetBgCpl4V466236NmzJ/Pnz+foo49m/vz5/PSnP42wwuipEYhIpLI1DNT+oGYs6N2xzuM146B33XUXN954Y17GQTOlRiAikcjVq4Ddu3fzu9/9jkGDBvHxxx9z3XXXMXr06LyOg2YqzB3KZgD/DXzk7sfVcrwL0D/+8N9Ad3dfEVY9IpI92VgWItVEMMCrr75KSUkJr7/+OqeffnqDiYNmKsydkR8ELkhy/H3gLHfvAIwCpodYi4hkSYdh8yJvAr/8YRHv3nFxnU1g06ZN/OpXv+KMM87gH//4Bw8//DAvv/xyQTYBCHerypfMrG2S468mPFwMHBZWLSISvWzMBaRKAxVKHDRTuTJHcB3wTLaLEJFgnDp6Af/815eRvV6qxeHcnT//+c/07t27IOKgmQpzaCgtZnY2sUbQP8k5Xc1sqZkt/fjjj6MrTkQyUrasgrYD5kTaBFo2bZK0CaxZs4YLLriASy+9lCZNmjB//nyeeuopNYEEWb0iMLMOwAPAhe6+pa7z3H068TmE4uJij6g8EUlTLm4S8+mnnzJixAjuueeegouDZiprjcDMioAngV+5+9vZqkNE6q9sWQU9H10e+eumEwcdOHAgmzdvLsg4aKbCjI/OBDoCrcxsIzAMaALg7tOAocABwJT4wk2V7p57uzqLSK2ingeA1HcF14yDPvPMMwWbBMpEmKmhq1Mc/zXw67BeX0TCkY2rgIOb783fBp9X5/GKigr69+/PI4880iBXBw1brqSGRCQPRH0V0KQRjOtc91XAzp07ueuuuxg9ejS7du1i0KBBDBw4sODjoJlSIxCRlLJxFZBsj4CacdBLL72UCRMmKAlUT2oEIpLUeRMX8s5Hn0f2eqnSQGvWrKFnz548++yzfO973yuI1UHDpkYgInVqO2BOZK+VapewxDhos2bNFAcNkBqBiHxNlENBqRaHqxkH/fWvf81tt92mOGiA1AhEZI8obwxLZ7N4xUGjoUYgIpQtq6DXo8uJ6rb9VJvE1IyDPvLII1x99dWKg4ZEjUCkwEW5V4DioLlJjUCkgB0xcA6VEV0GJLsruDoO2qtXL9577z1+/vOfM2HCBL7zne9EU1yBUyMQKUBRRkIzjYM+++yznHde3XcRS/DUCEQKTLsBcyKZC0g1Gbxt2zZGjBjBvffeS7Nmzbj77rvp0aOH4qBZoEYgUiCinAtINgxUVVW1Z7P46jjo6NGjOfDAAyOpTb5OjUCkgYvynoBUi8O98sorlJSU8MYbb3DGGWcwb948Tj755Ehqk7qpEYg0UGXLKuj32HJ27Q7/tVo2bcLwS46t8yogMQ7aunVrxUFzjBqBSAMU1Y1hqSaCd+7cycSJE7n99tuprKxk8ODBDBgwQHHQHKNGINLAdBg2j8++qAr9dVKtDvr000/Tu3dvxUHzgBqBSAORK1cBa9asobS0lAULFnDMMcewYMECfvKTn4Rel9RfmFtVzgD+G/jI3Y+r5bgBk4CLgO3A/7r7G2HVI9JQRTUZnOquYMVB81eYVwQPAvcCD9Vx/EKgffy/U4Gp8X9FJE1R7BjW4r8as3LEBXUer6qqYsaMGQwePJjNmzdz/fXXc9tttykOmkfC3LP4JTNrm+SUS4GH3N2BxWbW0swOcfd/hFWTSEMR1VWA4qCFIZtzBK2BDQmPN8Y/97VGYGZdga4ARUVFkRQnkquiuApo2qQRd1zeIWkc9JZbbuGPf/wjrVu35o9//CNXXXWV4qB5KpuNoLafmFrvfHf36cB0gOLi4qhWyhXJKVFcBaTaJKa2OOjAgQNp1qxZqHVJuLLZCDYCbRIeHwZsylItIjkrigaQaiLY3Xnqqafo06cP7733Hpdddhnjx49XHLSByGYjeBq4ycxmEZsk/lTzAyJfFcUqoan2Cl69ejU9e/ZUHLQBCzM+OhPoCLQys43AMKAJgLtPA+YSi46uIxYfvTasWkTyTdmyCvrMXk5ViAOhqYaBtm3bxvDhw7n33nvZd999FQdtwJI2AjPbP9lxd69zKUN3vzrF1zpwY9LqRApQ2KuENgImplgddMaMGQwaNIgtW7YoDloAUl0RvE5sAteAImBr/OOWwHqgXajViRSQKOYCki0PDbE46M0338yyZcs444wzmDx5suKgBSBpI3D3dgBmNg142t3nxh9fCGiQUCQgYV8FpFoWYuPGjfTv319x0AKV7hzB9929W/UDd3/GzEaFVJNIwYhifaBkVwE7d+5kwoQJ3H777VRVVSkOWqDSbQSbzWwI8DCxoaJfAltCq0qkAIR9Y1iqzeKfeuopevfuzfvvv684aIFLtxFcTSz18ydijeCl+OdEJENhXwW0P6gZC3p3rPP46tWrKS0t5bnnnlMcVIA0G0E8HVRqZvu6+79DrkmkQQp7MjjV4nCJcdDmzZszadIkunfvrjio0Cidk8zsdDNbDayOPz7BzKaEWplIAzKkrDzUJtD+oGZ1NoGqqiruv/9+2rdvz+TJk7nuuut4++23KSkpURMQIP2hobuA84ndDYy7rzCzM0OrSqSBCHsYKJM46I9+9CMmT57MSSedFFo9kp/SvrPY3TfUiJKFvxeeSB4Lc8vIVMtC1IyDzpw5kyuvvFJxUKlVuo1gg5mdDriZ7Q2UAGvCK0skf4U5F5BqIrhmHHTIkCEMGDBAcVBJKt1G0I3YtpKtia0a+izQI6yiRPJRmMNAqTaIqS0OOmHCBNq1083/klq6jeAod++S+AkzOwN4JfiSRPLP0YPnsjOkFeJSzQMkxkGPPfZYnnvuOc4999xQapGGKa3UEHBPmp8TKThHDJwTShM4uPnefDDm4qSrg/bs2ZMOHTqwdOlSJk+ezPLly9UEJGOpVh89DTgdONDMeiccagE0DrMwkVwX5vpAya4Cqqqq+O1vf8vgwYPZsmULXbt2ZdSoUVodVOot1dDQ3sC+8fOaJ3z+M6BTWEWJ5LKwJoMNeH/MxUnPefnllykpKVEcVAKVavXRRcAiM3vQ3T+MqCaRnFS2rIKBT65kx67dgT93OquD3nLLLcycOZPDDjtMcVAJVLqTxQ+YWWd33wZgZvsBs9z9/GRfZGYXEEsbNQYecPcxNY5/m9hCdkXxWsa7++8y/H8QCV1Yw0CpGkDNOOitt95K//79FQeVQKXbCFpVNwEAd99qZgcl+wIzawzcB5xHLHK6xMyedvfVCafdCKx295+Z2YHAWjN7xN3DW5JRJENh3BiW6n4Ad6esrIw+ffrw/vvvc/nllzN+/HjFQSUU6TaC3WZW5O7rAczscGKrkCbzA2Cdu78X/5pZwKXE1yuKc6C5xa5v9wU+ASozqF8kFGXLKhj+9Cq27dgV+HOns1m84qASpXQbwWDgZTNbFH98JtA1xde0BjYkPN4InFrjnHuJrV+0idhk9JXu/rUBWDPrWv16RUVFaZYsUj9hDQOlugrYunUrw4cP57777qN58+ZMnjyZ7t27s9deaa8EI1Ivad1H4O7zgJOBR4HZwCnuPj/Fl9U2i1XzKuJ8YDlwKHAicK+Ztajl9ae7e7G7FysiJ2EKqwn88odFdTaBqqoqpk+fzpFHHsk999zD9ddfzzvvvMPNN9+sJiCRSHUfwdHu/paZVe9evSn+b1F8qOiNJF++EWiT8PiwhK+vdi0wxt0dWGdm7wNHA39P+/9AJABhLQ+R6q7gxDjoj3/8YyZPnsyJJ54YeB0iyaT6c6MPcD0woZZjDpyT5GuXAO3NrB1QAVwFXFPjnPXAucBfzexg4CjgvTTqFglM2wFzAn/ObzU23hp9UZ3Ha8ZBZ82axRVXXKE4qGRFqvsIro//e3amT+zulWZ2EzCfWHx0hruvMrNu8ePTgFHAg2ZWTmwoqb+7b870tUQyla3J4J07dzJ+/HjuuOMOxUElZ6QaGro82XF3fzLF8bnA3Bqfm5bw8Sbgp6nLFAlOWMNAyRpAdRy0d+/efPDBB/ziF79g/PjxtG3bNvA6RDKVamjoZ/F/DyK25tAL8cdnAwuBpI1AJNeEcU9AqjjoqlWrKC0t5fnnn+fYY4/l+eef55xzko2qikQr1dDQtQBm9hfgGHf/R/zxIcRuFhPJC2GsD6Q4qDQU6f5Etq1uAnH/BI4MoR6RQGWjAVRVVfHAAw8wePBgPvnkE2644QZGjRpFq1atAq1DJCjpNoKFZjYfmEksLXQV8GJoVYkEIIy5gFRx0L/+9a+UlJSwfPlyxUElb6TVCNz9JjO7jNgdxQDT3f1P4ZUlUn9lyyro/ehyglwjNNVWkRs2bOCWW25h1qxZioNK3slksPIN4F/u/pyZ7WNmzd39X2EVJlIfp45ewD//FdyahXs1MsZ3PqHOq4AdO3YwYcIExUElr6XVCMzsemJr/ewPfJfYOkLTiN0MJpJ1YcwFpIqD/ulPf6JPnz6Kg0reS/eK4EZiq4n+DcDd30m1DLVIFMJoAKnmAd58801KS0t54YUXFAeVBiHdRvCFu39ZPd5pZnuRehlqkVAdMXAOlQH+FKaaB9i6dSvDhg1jypQpNG/enHvuuYdu3bopDip5L92f4EVmNghoambnAT2AP4dXlkjdor4KSIyDbt26dc9m8YqDSkORbiPoD/waKAduILZsxANhFSVSl6DvDE61VWRiHPTMM89k0qRJioNKg5OyEZhZI2Clux8H3B9+SSK1C3KV0FQNIDEO2qZNG8VBpUFL2QjcfbeZrUjcqlIkKkEPAxlwV5JhoMQ46O7duxk6dCj9+/dnn332CawGkVyT7tDQIcAqM/s78Hn1J939klCqEiH4fQKSzQMoDiqFLN1GMCLUKkQSBL00RCZx0OOOO05xUCk4qfYj+BbQDTiC2ETxb929MorCpPAEPQyUSRy0RYsWioNKwUr1E/97YBfwV+BC4BigNN0nN7MLgEnEdih7wN3H1HJOR+BuoAmw2d3PSvf5peEIemmID8ZcXOexmnHQG264gZEjRyoOKgUrVSM4xt2PBzCz35LBpvJm1pjYngXnEdvIfomZPe3uqxPOaQlMAS5w9/W6W7nwBD0MlOoqoGYcdPLkyZxwwgmBvb5IPkrVCPZs6BrfgziT5/4BsM7d3wMws1nApcDqhHOuAZ6sTiO5+0eZvIDkr6CHgVLtErZhwwb69evHo48+Sps2bXj00Ufp3Lmz4qAipG4EJ5jZZ/GPjdidxZ/FP3Z3b5Hka1sDGxIebwROrXHOkUATM1sINAcmuftD6RYv+anL/a/xyrufBPJcexmsu6PuYaAdO3bs2Sze3RUHFalFqq0qG3+D567tT62aK8PsBZxCbBXTpsBrZrbY3d/+yhOZdSW2+ilFRUXfoCTJpvMmLuSdjz5PfWKakt0UVjMO2qlTJ8aNG6c4qEgtwoxHbATaJDw+DNhUyzmb3f1z4HMzewk4AfhKI3D36cB0gOLiYi12l4eCvCcg1TBQzTjoCy+8wNlnnx3Y64s0NGE2giVAezNrB1QQ297ymhrnPAXcG1/NdG9iQ0d3hViTRCzoq4BkaaBPPvmEYcOGMXXqVFq0aMG9997LDTfcoDioSAqhvUPik8s3AfOJxUdnuPsqM+sWPz7N3deY2TxgJbCbWMT0zbBqkugEOQ8AyYeBqqqquP/++xkyZMieOOioUaM44IADAnt9kYbM3PNrpKW4uNiXLl2a7TIkiSCHgZo0gnGd674z+KWXXqKkpIQVK1Zw1llnMWnSJMVBRWphZq+7e3Ftx3TNLIEoW1bBuPlrqdi2I5Dna/FfjVk54oI6jysOKhIcNQL5xoKeB0i2NtCOHTsYN24cY8aMwd0ZNmwYt9xyi+KgIt+AGoF8I1HtEeDuPPnkk/Tp04cPP/yQTp06MX78eA4//PDAXl+kUKkRSL0EvUR0skhoeXk5paWlvPjiixx//PGKg4oETI1AMhLlHgHVcdApU6bw7W9/W3FQkZDoHSVpCXp10PYHNWNB7461HquqqmL69OnceuutioOKRECNQJIKeiI41eqgioOKRE+NQGoV9PLQkHweYP369fTr14/Zs2crDioSMTUC+ZoOw+bx2RdVgT2f4qAiuU2NQIDYDWG9Zy9nd4A3mqfaLP6JJ56gb9++ioOKZJkagQQ+EdzIYOIVdTcBxUFFcosaQQELeh4g1fLQn3zyCUOHDmXq1KmKg4rkEL0DC9QRA+dQGdAwUCNgYpJhoOo46JAhQ9i2bRvdunVj5MiRioOK5Ag1ggIS9D7Bye4FqFYzDjp58mQ6dOgQWA0i8s2pERSIdgPmfG2f0G8i2UQwfDUOWlRUxOzZs+nUqZPioCI5SI2ggQt6IjjV8tCKg4rkHzWCBiroNYGSrQwKX4+Ddu7cmXHjxikOKpIHGoX55GZ2gZmtNbN1ZjYgyXnfN7MqM+sUZj2FYEhZeaBNoJHFhoGSNYHy8nLOPfdcOnfuTIsWLXjxxReZPXu2moBIngjtisDMGgP3AecBG4ElZva0u6+u5bw7ie1tLPUU9BBQqisA+GoctGXLltx333107dpVcVCRPBPmFcEPgHXu/p67fwnMAi6t5bybgSeAj0KspcEqW1ZB2wFzAm0Cv/xhUdImUFVVxdSpU2nfvj1Tp06lW7duvP322/To0UNNQCQPhfmubQ1sSHi8ETg18QQzaw1cBpwDfL+uJzKzrkBXgKKiosALzUdd7n+NV979JNDnNOCuFGmgRYsWUVJSwsqVK+nYsSOTJk1SHFQkz4XZCGrLCdZMMN4N9Hf3qmSxQnefDkwHKC4uDjIFmZeOHjyXnVXBfRtSRUHh63HQxx57jF/84heKg4o0AGE2go1Am4THhwGbapxTDMyK/zJpBVxkZpXuXhZiXXmpbFkFvR5dHui9AAa8P+bipOfs2LGDsWPHcuedd+LuDB8+nH79+ikOKtKAhNkIlgDtzawdUAFcBVyTeIK7t6v+2MweBP6iJvBVQd8NDKnvBYD/xEH79OnD+vXrFQcVacBCawTuXmlmNxFLAzUGZrj7KjPrFj8+LazXbijCmAdIJw1UXl5OSUkJCxcupEOHDvz+97+nY8eOgdYhIrkj1IiHu88F5tb4XK0NwN3/N8xa8kUYv/wh9cqg8PU46JQpU7j++uuVBBJp4PQOzyFB3wsA6V0BVFZW7tksftu2bXTv3p2RI0ey//77B1qLiOQmNYIcEEYDgPTSQAsXLqS0tFRxUJECpkaQRWFsEJ/OL3+ADz/8kH79+vHYY48pDipS4NQIsiCseYB0msD27dv3rA5qZowYMYJ+/frRtGnTwOsRkfygRhChMKKgrVs2pd/5R6VsAO7O448/Tt++fVm/fj1XXHEF48aN053aIqJGEIUwrgC+1dh4a/RFaZ27cuVKSktL98RBH3roIc4666xA6xGR/KVGEKKwhoDSSQIBbNmyhaFDhzJt2jTFQUWkTvqNELCyZRUMenIl23ftDvy507kXABQHFZHMqBEE6LyJC3nno88Df950GwDE4qAlJSWUl5crDioiaVEj+AbKllUwbv5aKrbtCOX5M2kAH374IX379uXxxx9XHFREMqJGUE9hJICqtT+oGQt6d0zr3O3bt+9ZHVRxUBGpDzWCDIR5BWBAlwyuABQHFZGgqBGkKawEEGQ2BASxOGhJSQmLFi1SHFREvjE1gjqEmf6B9G8ES7RlyxZuvfVWfvOb3ygOKiKB0W+QWuTK+H+1yspKfvOb33Drrbfy6aef0qNHD0aMGKE4qIgEQo0gLowF4BLVpwHAV+OgZ599NpMmTeL449MfRhIRSSXURmBmFwCTiO1Q9oC7j6lxvAvQP/7w30B3d18RZk01DSkr55HF6wPdCzhRfRtAYhz08MMP5/HHH+fyyy9XHFREAhdaIzCzxsB9wHnENrJfYmZPu/vqhNPeB85y961mdiEwHTg1rJqqhbX+f6L99mnCsJ8dm9EcACgOKiLRC/OK4AfAOnd/D8DMZgGXAnsagbu/mnD+YuCwMAoJM/GTyID3x1xcr6+tGQe98sorGTt2rOKgIhK6MBtBa2BDwuONJP9r/zrgmdoOmFlXoCuQ8S/GKJpAuovA1SUxDnrCCScoDioikWoU4nPXNphd61C8mZ1NrBH0r+24u09392J3Lz7wwAMzKiLMJrDfPk24+8oT690EtmzZQo8ePTjppJMoLy9n6tSpvP7662oCIhKpMK8INgJtEh4fBmyqeZKZdQAeAC509y0h1hOIxmZcfWqbjG4AqykxDvrZZ58pDioiWRVmI1gCtDezdkAFcBVwTeIJZlYEPAn8yt3fDrGWb6Q+N3/V5cUXX6S0tJTy8nLOOeccJk2axHHHHRdAlSIi9RNaI3D3SjO7CZhPLD46w91XmVm3+PFpwFDgAGBKPBZZ6e7FQdaxd2Pjy6r6hUMzXfohGcVBRSRXmXtYCfpwFBcX+9KlS9M+v2xZBb0eXZ70PoEmjaByNxwa4F/+1bZv386dd97J2LFjMTMGDhxI3759FQcVkUiZ2et1/aHd4O8srv6lPnE8kB8AAAe7SURBVG7+WjZt2xHKL/vauDuPPfYYffv2ZcOGDYqDikjOavCNAGLNIOxf/IlWrFhBaWnpnjjoww8/zJlnnhnZ64uIZCLM+GjBqY6Dnnzyybz55pt74qBqAiKSywriiiBslZWVTJs2jaFDh/LZZ59x4403Mnz4cMVBRSQvqBF8Qy+++CIlJSW8+eabioOKSF7S0FA9ffDBB3Tq1IlzzjmHf//73zzxxBM899xzagIiknd0RZChmnHQkSNHKg4qInlNjSBN7s7s2bPp168fGzZs4KqrrmLs2LG0adMm9ReLiOQwDQ2lYcWKFXTs2JGrrrqK/fffn0WLFjFz5kw1ARFpENQIkti8eTPdu3fn5JNPZtWqVUybNk1xUBFpcDQ0VIva4qAjRoxgv/32y3ZpIiKBUyOo4YUXXqC0tFRxUBEpGBoaiquOg5577rmKg4pIQSn4K4Lt27czZswYxo0bpzioiBSkgm0EioOKiMQU5NDQ8uXLFQcVEYkrqEZQHQc95ZRTFAcVEYkLtRGY2QVmttbM1pnZgFqOm5lNjh9faWYnh1XLnDlzaN++Pffffz833ngj77zzDjfccAONGzcO6yVFRPJCaHMEZtYYuA84D9gILDGzp919dcJpFwLt4/+dCkyN/xu4o446itNOO42xY8cqCSQikiDMyeIfAOvc/T0AM5sFXAokNoJLgYc8tnHyYjNraWaHuPs/gi7miCOOYO7cuUE/rYhI3gtzaKg1sCHh8cb45zI9BzPramZLzWzpxx9/HHihIiKFLMxGYLV8zutxDu4+3d2L3b34wAMPDKQ4ERGJCbMRbAQS85iHAZvqcY6IiIQozEawBGhvZu3MbG/gKuDpGuc8DfxPPD30Q+DTMOYHRESkbqFNFrt7pZndBMwHGgMz3H2VmXWLH58GzAUuAtYB24Frw6pHRERqF+oSE+4+l9gv+8TPTUv42IEbw6xBRESSK6g7i0VE5OvUCERECpzFRmfyh5l9DHyYwZe0AjaHVE4Y8qnefKoV8qvefKoV8qvefKoVgqv3cHevNX+fd40gU2a21N2Ls11HuvKp3nyqFfKr3nyqFfKr3nyqFaKpV0NDIiIFTo1ARKTAFUIjmJ7tAjKUT/XmU62QX/XmU62QX/XmU60QQb0Nfo5ARESSK4QrAhERSUKNQESkwDWYRpBL22KmkkatXeI1rjSzV83shGzUmVBP0noTzvu+mVWZWaco66tRQ8pazayjmS03s1VmtijqGmvUkupn4dtm9mczWxGvN2vrcZnZDDP7yMzerON4Lr3HUtWaa++xpPUmnBfOe8zd8/4/YovavQt8B9gbWAEcU+Oci4BniO2B8EPgbzlc6+nAfvGPL8xWrenWm3DeC8TWluqUq7UCLYntklcUf3xQLn9vgUHAnfGPDwQ+AfbOUr1nAicDb9ZxPCfeY2nWmjPvsXTqTfh5CeU91lCuCPZsi+nuXwLV22Im2rMtprsvBlqa2SFRF0oatbr7q+6+Nf5wMbF9GrIlne8twM3AE8BHURZXQzq1XgM86e7rAdw91+t1oLmZGbAvsUZQGW2Z8ULcX4q/fl1y5T2WstYce4+l872FEN9jDaURBLYtZgQyreM6Yn9lZUvKes2sNXAZMI3sSud7eySwn5ktNLPXzex/Iqvu69Kp917ge8Q2bCoHSt19dzTlZSxX3mOZyvZ7LKWw32OhLkMdocC2xYxA2nWY2dnEfkh/FGpFyaVT791Af3eviv3hmjXp1LoXcApwLtAUeM3MFrv722EXV4t06j0fWA6cA3wXWGBmf3X3z8Iurh5y5T2Wthx5j6Uj1PdYQ2kE+bQtZlp1mFkH4AHgQnffElFttUmn3mJgVvwHtBVwkZlVuntZNCXuke7PwWZ3/xz43MxeAk4AstEI0qn3WmCMxwaJ15nZ+8DRwN+jKTEjufIeS0sOvcfSEe57LJsTJAFOtOwFvAe04z+TbsfWOOdivjqR9fccrrWI2K5tp+fD97bG+Q+SvcnidL633wOej5+7D/AmcFwO1zsVGB7/+GCgAmiVxZ+HttQ9AZsT77E0a82Z91g69dY4L/D3WIO4IvA82hYzzVqHAgcAU+J/AVR6llZLTLPenJBOre6+xszmASuB3cAD7p40spfNeoFRwINmVk7sF2x/d8/KEspmNhPoCLQys43AMKBJQq058R6DtGrNmfcYpFVvuK8f7zAiIlKgGkpqSERE6kmNQESkwKkRiIgUODUCEZECp0YgIlLgGkR8VCRMZnYAsXsPAP4fUAV8HH/8A4+tEySStxQfFcmAmQ0H/u3u4xM+t5e7Z2UhOJEg6IpApB7M7EFiq0WeBLxhZv8ioUHE15X/b3f/wMx+CZQQu3v4b0APd6/KTuUiX6c5ApH6OxL4ibv3qesEM/secCVwhrufSGxYqUtE9YmkRVcEIvX3WBp/2Z9LbLXTJfGlDJqS3T0bRL5GjUCk/j5P+LiSr15hfyv+rwG/d/eBkVUlkiENDYkE4wNiWw0S36u3XfzzzwOdzOyg+LH9zezwrFQoUgc1ApFgPAHsb2bLge7E9zdw99XAEOBZM1sJLACysn2jSF0UHxURKXC6IhARKXBqBCIiBU6NQESkwKkRiIgUODUCEZECp0YgIlLg1AhERArc/wdg6GnPIbERGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE: 0.111\n"
     ]
    }
   ],
   "source": [
    "model.history[['training_loss', 'validation_loss']].plot()\n",
    "plt.ylabel('Loss')\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "y_pred = model.predict(x)\n",
    "plt.scatter(y, y_pred)\n",
    "plt.plot((y.min(), y.max()), (y.min(), y.max()), 'k-')\n",
    "plt.xlabel('True')\n",
    "plt.ylabel('Predicted')\n",
    "plt.show()\n",
    "plt.close()\n",
    "print('MAE: {:.3f}'.format(np.mean(np.abs(y_pred - y))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train the Model with P_Fun\n",
    "\n",
    "Here we train the model with loss weights (0.2, 0.8) which still weights the mean absolute error of y_predicted vs. y_noise but also gives much more weight to the p_fun calculation. We can see by the results that supplementing the pure NN training with a physical estimate of the true value helps the overall model prediction capabilities. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "PhysicsGuidedNeuralNetwork.seed(0)\n",
    "model = PhysicsGuidedNeuralNetwork(p_fun=p_fun_pythag, \n",
    "                                   hidden_layers=hidden_layers, \n",
    "                                   loss_weights=(0.2, 0.8),\n",
    "                                   n_features=2, n_labels=1)\n",
    "model.fit(x, y_noise, p, n_batch=4, n_epoch=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU1d348c93luwLWViSsCQBFAhLWF1Q1LoBorjgXiso8nOvPo8LfWytWq222upj5RE3rFqsCIq1ioooilYREmRHFtkSEiAESEIWssz5/XEnYQiTyWSZDCTf9+t1XzNz77l3vrlJ5jvnnnPPEWMMSimlVH22YAeglFLq+KQJQimllFeaIJRSSnmlCUIppZRXmiCUUkp55Qh2AK0pMTHRpKamBjsMpZQ6YWRnZ+8zxnT2tq1dJYjU1FSysrKCHYZSSp0wRGRHQ9v0EpNSSimvNEEopZTyShOEUkopr9pVG4RSqu1UVVWRm5tLRUVFsENRfggLC6N79+44nU6/99EEoZRqltzcXKKjo0lNTUVEgh2O8sEYQ2FhIbm5uaSlpfm9n15iUko1S0VFBQkJCZocTgAiQkJCQpNre5oglFLNpsnhxNGc35UmCFcNfPMX2PJFsCNRSqnjiiYImx3+87/w08fBjkQppY4rmiAA4tNh/9ZgR6GUaoKDBw/yf//3f03eb/z48Rw8eNBnmYcffphFixY1NzSvoqKiWvV4bUETBFgJ4sC2YEehlGqChhJETU2Nz/0WLFhAp06dfJZ57LHHOO+881oUX3ug3VwB4tJg3XyorgRHSLCjUeqE8+i/17E+r7hVjzkgOYbfX5zR4Pbp06fz888/k5mZidPpJCoqiqSkJFauXMn69eu59NJLycnJoaKigl//+tdMmzYNODJm26FDhxg3bhxnnHEG3333HSkpKfzrX/8iPDycyZMnM2HCBCZNmkRqaio33ngj//73v6mqqmLu3Ln069ePgoICrrvuOgoLCxk5ciSffvop2dnZJCYm+vy5jDE88MADfPLJJ4gIv/3tb7n66qvJz8/n6quvpri4mOrqal588UVOP/10br75ZrKyshARbrrpJu69995WPc++aA0CrBqEcUFRTrAjUUr56amnnqJ3796sXLmSp59+mmXLlvHEE0+wfv16AGbNmkV2djZZWVk8//zzFBYWHnOMzZs3c8cdd7Bu3To6derEe++95/W9EhMTWbFiBbfddhvPPPMMAI8++ii/+MUvWLFiBZdddhk7d+70K+7333+flStXsmrVKhYtWsT9999Pfn4+b7/9NhdeeGHdtszMTFauXMmuXbtYu3Yta9asYcqUKc08W82jNQiwEgTA/m2Q0Du4sSh1AvL1Tb+tjBo16qibwJ5//nnmz58PQE5ODps3byYhIeGofdLS0sjMzARg+PDhbN++3euxL7/88roy77//PgDffvtt3fHHjh1LXFycX3F+++23XHvttdjtdrp27cpZZ53F8uXLGTlyJDfddBNVVVVceumlZGZmkp6eztatW7nrrru46KKLuOCCC/w/Ia1AaxDgkSC0oVqpE1VkZGTd86+++opFixbx/fffs2rVKoYOHer1JrHQ0NC653a7nerqaq/Hri3nWcYY06w4G9pvzJgxLFmyhJSUFG644QbefPNN4uLiWLVqFWeffTYzZsxg6tSpzXrP5tIEARDVBZyRmiCUOoFER0dTUlLidVtRURFxcXFERETw008/sXTp0lZ//zPOOIN3330XgIULF3LgwAG/9hszZgxz5syhpqaGgoIClixZwqhRo9ixYwddunThlltu4eabb2bFihXs27cPl8vFFVdcwR/+8AdWrFjR6j+HL3qJCUAE4tM0QSh1AklISGD06NEMHDiQ8PBwunbtWrdt7NixzJw5k8GDB3PyySdz6qmntvr7//73v+faa69lzpw5nHXWWSQlJREdHd3ofpdddhnff/89Q4YMQUT485//TLdu3XjjjTd4+umn6xrc33zzTXbt2sWUKVNwuVwAPPnkk63+c/giza0mHY9GjBhhmj2j3JxfQsFGuHN56walVDu1YcMG+vfvH+wwgubw4cPY7XYcDgfff/89t912GytXrgx2WD55+52JSLYxZoS38lqDqBWXBps+s4besNmDHY1S6ji3c+dOrrrqKlwuFyEhIbzyyivBDqnVaYKoFZ8ONZVQvAs69Qx2NEqp41zfvn358ccfj1pXWFjIueeee0zZL7744pgeVCcCTRC1PLu6aoJQSjVDQkLCcX+ZqSm0F1Mt7eqqlFJH0QRRKyYZ7CGaIJRSyk0TRC2bHeJSNUEopZSbJghP8elwYHuwo1BKqeOCJghPce6b5drRvSFKKUvtfAx5eXlMmjTJa5mzzz6bxu6leu655ygrK6t77c/8Ek0xefJk5s2b12rHa4mAJggRGSsiG0Vki4hM97L9ehFZ7V6+E5EhHtu2i8gaEVkpIs28+62J4tOhqgwO7WmTt1NKtb3k5OQWfQDXTxD+zC9xogpYN1cRsQMzgPOBXGC5iHxojFnvUWwbcJYx5oCIjANeBk7x2H6OMWZfoGI8hmdX1+hubfa2Sp3wPpkOu9e07jG7DYJxTzW4+cEHH6RXr17cfvvtADzyyCOICEuWLOHAgQNUVVXx+OOPM3HixKP22759OxMmTGDt2rWUl5czZcoU1q9fT//+/SkvL68rd9ttt7F8+XLKy8uZNGkSjz76KM8//zx5eXmcc845JCYmsnjx4rr5JRITE/nrX//KrFmzAJg6dSr33HMP27dvb3DeicZ88cUX3HfffVRXVzNy5EhefPFFQkNDmT59Oh9++CEOh4MLLriAZ555hrlz5/Loo49it9uJjY1lyZIlzTnrRwlkDWIUsMUYs9UYUwm8Axz1mzLGfGeMqR3hainQPYDxNC7ePVSwNlQrddy75pprmDNnTt3rd999lylTpjB//nxWrFjB4sWL+e///m+fo66++OKLREREsHr1ah566CGys7Prtj3xxBNkZWWxevVqvv76a1avXs3dd99NcnIyixcvZvHixUcdKzs7m9dff50ffviBpUuX8sorr9TdSOfvvBOeKioqmDx5MnPmzGHNmjV1kwjt37+f+fPns27dOlavXs1vf/tbwJoF77PPPmPVqlV8+OGHTTqXDQnkjXIpgOcMPLkcXTuo72bgE4/XBlgoIgZ4yRjzsredRGQaMA2gZ88W3uDWqSeIXROEUk3l45t+oAwdOpS9e/eSl5dHQUEBcXFxJCUlce+997JkyRJsNhu7du1iz549dOvm/YrAkiVLuPvuuwEYPHgwgwcPrtv27rvv8vLLL1NdXU1+fj7r168/ant93377LZdddlndsOOXX34533zzDZdcconf80542rhxI2lpaZx00kkA3HjjjcyYMYM777yTsLAwpk6dykUXXcSECRMAGD16NJMnT+aqq66qm7+ipQKZIMTLOq+pXETOwUoQZ3isHm2MyRORLsDnIvKTMeaYOpM7cbwM1mB9LYrY7oROPXR+aqVOEJMmTWLevHns3r2ba665htmzZ1NQUEB2djZOp5PU1FSv80B4Ejn2o2rbtm0888wzLF++nLi4OCZPntzocXzVVOrPO+F5Kaupx3M4HCxbtowvvviCd955hxdeeIEvv/ySmTNn8sMPP/Dxxx/XzUbX0uE9AnmJKRfo4fG6O5BXv5CIDAZeBSYaY+rmBDTG5Lkf9wLzsS5ZBV58utYglDpBXHPNNbzzzjvMmzePSZMmUVRURJcuXXA6nSxevJgdO3b43H/MmDHMnj0bgLVr17J69WoAiouLiYyMJDY2lj179vDJJ0cubjQ0D8WYMWP44IMPKCsro7S0lPnz53PmmWc2+2fr168f27dvZ8uWLQC89dZbnHXWWRw6dIiioiLGjx/Pc889Vze0x88//8wpp5zCY489RmJiIjk5LZ9COZA1iOVAXxFJA3YB1wDXeRYQkZ7A+8ANxphNHusjAZsxpsT9/ALgsQDGekRcGuRmW11dvXyzUEodPzIyMigpKSElJYWkpCSuv/56Lr74YkaMGEFmZib9+vXzuf9tt93GlClTGDx4MJmZmYwaZX0PHTJkCEOHDiUjI4P09HRGjx5dt8+0adMYN24cSUlJR7VDDBs2jMmTJ9cdY+rUqQwdOtSvy0nehIWF8frrr3PllVfWNVLfeuut7N+/n4kTJ1JRUYExhmeffRaA+++/n82bN2OM4dxzz2XIkCGNvEPjAjofhIiMB54D7MAsY8wTInIrgDFmpoi8ClwB1Kb5amPMCBFJx6o1gJXE3jbGPNHY+7VoPoha370ACx+CB7ZBRHzLjqVUO9bR54M4ER1X80EYYxYAC+qtm+nxfCpwzCSrxpitQMvTX3N4dnXVBKGU6sD0Tur6tKurUqoN3HHHHWRmZh61vP7668EO6yg6H0R9canWoyYIpRpljPHaC0g1bsaMGW36fs1pTtAaRH3OcIhJ0a6uSjUiLCyMwsLCZn3wqLZljKGwsJCwsLAm7ac1CG+0q6tSjerevTu5ubkUFBQEOxTlh7CwMLp3b9pgFZogvIlLhU2fBjsKpY5rTqeTtLS0YIehAkgvMXkTnw6lBXD42JthlFKqo9AE4Y1nV1ellOqgNEF4o11dlVJKE4RXcZoglFJKE4Q3YTEQ2Vm7uiqlOjRNEA2JT9c2CKVUh6YJoiFxaXqJSSnVoWmCaEh8OhTvgqrGJ/ZQSqn2SBNEQ2q7uh7wPeGIUkq1V5ogGqJdXZVSHZwmiIbU1SC0oVop1TFpgmhIeByExWoNQinVYWmCaIiIjuqqlOrQNEH4ol1dlVIdmCYIX+LT4WAO1FQFOxKllGpzmiB8iU8HUwMHdwY7EqWUanOaIHyp6+qqPZmUUh2PJghftKurUqoD0wThS1RXcEZoQ7VSqkPSBOGLdnVVSnVgmiAaE5eqbRBKqQ5JE0Rj4tOtNghXTbAjUUqpNqUJojHx6VBTCcV5wY5EKaXalCaIxuiorkqpDkoTRGO0q6tSqoPSBNGYmBSwh2gNQinV4QQ0QYjIWBHZKCJbRGS6l+3Xi8hq9/KdiAzxd982Y7NDp16aIJRSHU7AEoSI2IEZwDhgAHCtiAyoV2wbcJYxZjDwB+DlJuzbduLTYf/2oL29UkoFQyBrEKOALcaYrcaYSuAdYKJnAWPMd8aYA+6XS4Hu/u7bpmpvljMmaCEopVRbC2SCSAFyPF7nutc15Gbgk6buKyLTRCRLRLIKCgpaEK4P8elQVQqH9gbm+EopdRwKZIIQL+u8fgUXkXOwEsSDTd3XGPOyMWaEMWZE586dmxVoo2q7umpPJqVUBxLIBJEL9PB43R045m4zERkMvApMNMYUNmXfNlPb1VUbqpVSHUggE8RyoK+IpIlICHAN8KFnARHpCbwP3GCM2dSUfdtUbA8QuyYIpVSH4gjUgY0x1SJyJ/AZYAdmGWPWicit7u0zgYeBBOD/RASg2n25yOu+gYq1UY4QiO2uCUIp1aEELEEAGGMWAAvqrZvp8XwqMNXffYMqPl1HdVVKdSh6J7W/dF4IpVQHownCX/HpUHEQyvYHOxKllGoTmiD8pV1dlVIdjCYIf9V1ddUEoZTqGDp8gqiucbFkUwE/7S72XTAu1XrUdgilVAfR4ROEiHDrP7J5Z1mO74LOcIhO1hqEUqrD6PAJwm4T+ifFsC6vqPHC2pNJKdWBdPgEAZCRHMP6vGJcrkZGa41P0wShlOowNEEAA5NjKa2sYcf+Mt8F49OhdC8cLmmbwJRSKog0QQADkmMAWLurkctMdV1dtwc2IKWUOg5oggBO6hqN0y6sy2ukJ5OO6qqU6kA0QQAhDhsndY1uvKE6zl2D0AShlOoANEG4ZSTHsC6vGONrWtGwGIhI1K6uSqkOQROEW0ZyLPtLK9ldXOG7oHZ1VUp1EJog3Aam1DZU+9EOoTUIpVQHoAnCrV+3GERovB0iPh2Kd0FVIzUNpZQ6wWmCcIsMdZCeGOlHT6Y0wMDBHW0Sl1JKBYsmCA8ZybGsa/ReCO3qqpTqGDRBeMhIjiGvqIIDpZUNF9IEoZTqIDRBeBiYEgvg+zJTeByExmpDtVKq3dME4SGjdsgNXw3VIjpon1KqQ9AE4aFTRAgpncL9G3JDE4RSqp3TBFFPRnKMHw3VaVCUAzVVbROUUkoFgSaIejKSY9lWWErp4eqGC8Wng6vaShJKKdVO+ZUgRCRSRGzu5yeJyCUi4gxsaMExMCUGY2BDvo/LTNqTSSnVAfhbg1gChIlICvAFMAX4e6CCCqaMZKsnk8+5IeoShPZkUkq1X/4mCDHGlAGXA38zxlwGDAhcWMHTNSaUhMgQ3w3VUV3BGaEJQinVrvmdIETkNOB64GP3OkdgQgouESEjJdZ3ghCx5obQS0xKqXbM3wRxD/AbYL4xZp2IpAOLAxdWcGUkx7BpTwmHq2saLqT3Qiil2jm/EoQx5mtjzCXGmD+5G6v3GWPuDnBsQZORHEO1y7B5z6GGC8WnWXNTu1xtFpdSSrUlf3sxvS0iMSISCawHNorI/YENLXgGJtcOudFIQ3XNYSjJa6OolFKqbfl7iWmAMaYYuBRYAPQEbghYVEHWMz6CqFCH78mDtKurUqqd8zdBON33PVwK/MsYUwX4mLzZIiJjRWSjiGwRkeletvcTke9F5LCI3Fdv23YRWSMiK0Uky884W4XNJgxIimm8BgHak0kp1W75myBeArYDkcASEekF+BywSETswAxgHFaX2GtFpH7X2P3A3cAzDRzmHGNMpjFmhJ9xtpqMlBg25JdQ42ogD8akgM2pNQilVLvlbyP188aYFGPMeGPZAZzTyG6jgC3GmK3GmErgHWBivePuNcYsB467QY0ykmMpr6ph274GGqptdohL1QShlGq3/G2kjhWRv4pIlnv5C1ZtwpcUwHOwolz3On8ZYKGIZIvINB+xTauNq6CgoAmH921gijX0t8/7IeLT9RKTUqrd8vcS0yygBLjKvRQDrzeyj3hZ12i7hYfRxphhWJeo7hCRMd4KGWNeNsaMMMaM6Ny5cxMO71vvzlGEOGyNDLmRBge2gWnKj6WUUicGf++G7m2MucLj9aMisrKRfXKBHh6vuwN+9wk1xuS5H/eKyHysS1ZL/N2/pZx2G/26RTdeg6g8BKUFENWlrUJTSqk24W8NolxEzqh9ISKjgfJG9lkO9BWRNBEJAa4BPvTnzdyjx0bXPgcuANb6GWuryUi2htwwDdUQtKurUqod87cGcSvwpojEul8fAG70tYMxplpE7gQ+A+zALPcwHbe6t88UkW5AFhADuETkHqweT4nAfBGpjfFtY8ynTfvRWi4jOYZ/LttJ7oFyesRHHFvAs6trz1PbNjillAowvxKEMWYVMEREYtyvi90f5qsb2W8B1o11nutmejzfjXXpqb5iYIg/sQVS7RzV6/KKvSeI2B4gNq1BKKXapSbNKGeMKXbfUQ3wXwGI57jSPykGu00avmHOEWIlCU0QSql2qCVTjnrrpdSuhDnt9O4c2XhD9QHt6qqUan9akiA6RN9Oq6G6ka6uWoNQSrVDPhOEiJSISLGXpQRIbqMYgyojOYY9xYcpKDnsvUB8OpQfgLL9bRuYUkoFmM8EYYyJNsbEeFmijTHtcka5+jIaG/q7tidTCy4z/XtVHmtyfdRSlFIqCFpyialDGJDcyJAbLRzVNb+onLvf+ZFrX1nq+1KWUkq1MU0QjYgNd9IzPqLhD++4VOux8OdmHf/9FbswBsJD7Ex+fTk5+8uaF6hSSrUyTRB+yEiOabgG4QyHpExY+Q+oqmjScY0xzMvOZVRaPP+85RSqalz8atYyCg810N6hlFJtSBOEHwamxLKjsIziigZGJT//UTi4E354sUnHzd5xgG37SrlyeHf6dInmtRtHkl9Uzk1/X07p4epWiFwppZpPE4Qfatsh1jdUi0g/G04aB0v+Aof2+n3cuVm5RITYGT8oCYDhveKYcd0w1uYVc+s/sqmsdrUwcqWUaj5NEH7IaKyhGuCCx6G6HBY/4dcxyyqr+Wh1HuMHJREZeqRD2Ln9u/Lk5YP4ZvM+Hpi3CldDM9oppVSAaYLwQ5foMLpEh7LO19wQiX1g5C2w4k3Ys67RY36yZjellTVcOfzYoaiuGtGD+y88mQ9W5vHkJxtaErpSSjWbJgg/+WyornXWAxAaA5/9T6OTCM3NzqFXQgSj0uK9br/97N7ceFovXvlmG68s0Tu1lVJtTxOEnzKSY9lScIiKqpqGC0XEw9m/ga1fwabPGiy2s7CMpVv3M2lYd9xDmh9DRHj44gwuGpTEEws2MP/H3Bb+BEop1TSaIPw0MCWGGpfhp90lvguOvBkS+sLC30KN915P81bkIgJXeLm85MluE/569RBOS0/g/rmr+XpT6825rZRSjdEE4adGh9yoZXdaDdaFmyFr1jGbXS7De9m5nNEnkeRO4Y2+b6jDzku/Gk7frtHc9o9sVuUcbFb8SinVVJog/NQ9LpyYMEfj7RAAJ10IaWfBV09aA/l5WLq1kF0Hy5nUSO3BU0yYkzemjCQhKoQpf1/O1oJDTQ1fKaWaTBOEn0TEGvrbV0+mI4Xhwj9CRRF8/eejNs3NziU6zMGFGd2a9P5dYsJ486ZTEOBXs5axt7hpd20rpVRTaYJogozkGDbsLqGqxo8b2LoNhKE3wLKXYd8WAIorqvhkbT4XD0kmzGlv8vunJUby+pSR7C+t5MbXlzd8Z7dSSrUCTRBNMDAllspqFz/7e4nnF78FRzh8/jsAPl6dT0WVy+u9D/4a3L0TM385nM17Spj2ZhaHq330qlJKqRbQBNEEdXdU7/KjHQIgqguc+V+wcQFs/Zq5WTn06RJFZo9OLYpjzEmdeebKISzdup9756ykRu+2VkoFgCaIJkjvHEWY0+ZfQ3WtU2+HTj05/PF0Vu7cz5XDG773oSkuHZrCQ+P7s2DNbh799zpMIzfmKaVUU2mCaAK7TeifFMPapkzs4wyD8x4ltHA9VzmWcNmwlFaL55Yx6dxyZhpvfr+Dv325pdWOq5RSoAmiyTKSY9iQV9ykQfSq+01klZzMb0Lm0iWkdRuWfzOuP5cPTeGvn2/ije+2t+qxlVIdmyaIJspIjqXkcDU5B/yf+e2bLYU8XHE9sa4D8O2zrRqPzSb8edJgzh/Qld9/uI73snVIDqVU69AE0UQD3XdUr/W3oRqYl53LzvD+1Ay8Er57wZpcqBU57Db+du1QTu+dwAPvreazdbtb9fhKqY5JE0QTndQtCodNGh9yw+1gWSWfr9/DxMwU7Oc/AmKDRY+0elxhTjsv/2oEg1JiuevtH/nPln2t/h5KqY5FE0QThTrs9O0a7XdPpn+tzKOyxsWVI7pDbHc4/S5Y+x7kLGv12KJCHfx9ykjSEiO55c0sVuw80PhOSinVAE0QzWDNDVHkV9fSudk5DEiKqRvsj9G/hqhu8OlvwNX6U4p2igjhrZtH0Tk6lMmzlrEhvwldcpVSyoMmiGbISI5h36FK9pYc9lluQ34xa3cVW7WHWqFRcO7DsCsL1r0fkPi6xITxj5tPISLEwQ2vLWP7vtKAvI9Sqn3TBNEMA1NqG6p9t0PMzcrFaRcmZta792HItZA0BD7/PVSVByTGHvER/GPqKGpcLq5/9QfyiwLzPkqp9iugCUJExorIRhHZIiLTvWzvJyLfi8hhEbmvKfsGU/+kGETw2Q5RWe3ig5W7OK9/V+IjQ47eaLNZo70W58L3LwQszj5donnzplMoKq/il6/+QOEh3zUepZTyFLAEISJ2YAYwDhgAXCsiA+oV2w/cDTzTjH2DJirUQWpCpM+eTF/+tJf9pZVHX17ylHoG9JsA3zwLJYHrljqoeyyv3TiC3APl3Pj6Mh0BVinlt0DWIEYBW4wxW40xlcA7wETPAsaYvcaY5UD9T61G9w22jOQYn/dCzMvOpXN0KGP6dm74IOc/BjWV8OUfAhDhEaekJzDzl8P5Kb+EqX/PorxSR4BVSjUukAkiBcjxeJ3rXteq+4rINBHJEpGsgoK2m7M5IzmWXQfLOVhWecy2gpLDLN64l8uHpuCw+zjFCb3hlP8HP86GTQuh+thjtZZz+nXh2aszWb5jP7fNzqayuvV7UCml2pdAJghvQ5b6O4CR3/saY142xowwxozo3NnHt/VWVjv093ov7RAf/LiLGpdp+PKSpzH3Q1RXePtKeKonvD4eFj0KGz+Fsv2tGvPFQ5L542WD+GpjAfe+q8OEK6V8cwTw2LlAD4/X3YG8Nti3TdQmiLV5RZzeJ7FuvTGGudk5ZPboRJ8u0Y0fKLwT3P49bFti3TyXsxS+ex5c1db2xJOgxynW0vNUSOhjTWnaTNeO6klJRRV/XPAT0aEOnrx8UKsMP66Uan8CmSCWA31FJA3YBVwDXNcG+7aJhKhQkmLDjunJtDq3iE17DvHEZQP9P1hEPGRcai0AlWWQ96OVLHb+ABv+DT++ZW0Lj3cnjFFWwkgeCs7wJsU+bUxvisureWHxFqLDHPzP+P6aJJRSxwhYgjDGVIvIncBngB2YZYxZJyK3urfPFJFuQBYQA7hE5B5ggDGm2Nu+gYq1uTKSY49JEHOzcwh12Lh4SHLzDxwSAamjrQWsO64LN0POD1bCyFkKmz6xttmc1j0VGZfBiCkQEunXW/z3BSdRXFHFK99sIzbcyZ2/6Nv8eJVS7VIgaxAYYxYAC+qtm+nxfDfW5SO/9j3eZCTH8MVPeyirrCYixEFFVQ0frsxj7MBuxIQ5W++NbDbofLK1DPuVta600EoYOT/A9m9g4UPWUOKn3wkjp0Ko78tbIsIjF2dQUlHNMws3ERcZwvWn9Gq9mJVSJzy9k7oFMpJjMAY25JcAsHD9HoorqrlyeI9G9mwFkQnQbzyc/yjc8iXc9JlVk1j0CDw3CJY8DRW+7/S22YSnJw3m7JM788iH61idezDwcSulThiaIFqgdsiN2hvm5mblkBwbxmm9E9o+mJ6nwg3vw9QvrTaKLx+3EsXiJ6G84VFdHXYbz16VSeeoUO54ewVF5XojnVLKogmiBZJiw4iLcLJuVzH5ReV8u2UfVwzvjt0WxAbf7sPhujkw7avo6j0AABcDSURBVCtIPRO+fgqeHQRf/KHBbrNxkSH87bph5B+s4MF5q/0apVYp1f5pgmgBEbEaqvOLeH/FLoyBScP9uPehLSQPhWtmw63/gT7nwjd/gWcHwucPw6Fjbygc3iuOB8f249N1u/m7zm2tlEITRItlpMSwcXcJ72blMCotnl4J/vUiajPdBsJVb8DtS602i+/+Zl16+uyhY8aAmnpmGuf178IfF2xgVY62RyjV0WmCaKGM5Fiqagw7Csu48nipPXjTpR9c8Srcscy632Lpi/DcYFjwABTtAqwa0TNXDqFLdJi2RyilNEG0VO0d1REhdsYPSgpyNH5I7AuXzYS7smDwVZD1GjyfCT+8BFgz0v3tuqHsLqrggXmrtD1CqQ5ME0QLpSVEEh8ZwsTMZCJDA3pbSeuKT4eJL8BdKyD9bPh0Omz9GoBhPeOYPq4fn63bw+v/2R7MKJVSQaQJooVsNuGju87g4QkZwQ6leeJ6waTXIaEvvHczFOcDcPMZaZzXvytPfrKBldoeoVSHpAmiFSR3Cic8xB7sMJovNAqufssaA2reFKipcrdHDLbaI2avoKhM2yOU6mg0QShL55Ph4v+Fnd/DF48BVnvEC9cNZU9xBfdpe4RSHY4mCHXE4CthxM3WcOMbPgJgqLs94vP1e5il7RFKdSiaINTRxj5p3WT3we2wfytgtUecP6ArTy7YwI87Gx62QynVvmiCUEdzhMKVb1iTEr37K6gqt9ojJg2ha0wYd779o9dpVpVS7Y8mCHWsuF5w2Uuwew188iAAsRFOZlw/jL0lFdw3V8drUqoj0AShvDt5LJzxX7DiDVj5NgCZPTrxm3H9WbRhD699uy3IASqlAk0ThGrYOQ9ZI8J+9F+wx5rQb8roVC4Y0JWnPvmJFdoeoVS7pglCNczugCteg7AYmHMDVBQjIjw9aQjdYsO4S9sjlGrXNEEo36K7WndaH9gOH94FxljtEdfVtkfo/RFKtVeaIFTjUkfDuQ/D+g/qBvUb0qMT/zO+P4s27OXVb5rfHuFyGcorazTJKHUcOoFGl1NBdfrdkPMDLHwIUoZBj1FMPj2VH7bu50+f/kR5VQ02gbLKGsqraqioqqG8subo11XW6wr3Omu9C4BzTu7Mi78cTpjzBB6yRKl2RtrTN7cRI0aYrKysYIfRfpUfgJfOAlc1/L9vIDKBovIqrnjxO7bsPQSA3SZEOO2EhdiJCLET7rQT5rQeI0Ks9bXPa7eVVVbzyjfb+EW/Lsz85XBCHFqxVaqtiEi2MWaE122aIFST5K2E1y6wLjtdPw9sdqprXJRW1hARYsdpb96H++wfdvDQ/LWMH9SN568ZiqOZx1FKNY2vBKH/happkjNh3J/g5y9hyTMAOOw2YsOdzU4OANef0ovfTRjAgjW7uX/ealyu9vPFRakTlbZBqKYbPhl2LoWvnoQeI6H3L1rlsDefkUZFVQ1Pf7aRMKedP142EBFplWMrpZpOaxCq6URgwl+hcz94b2rdnNat4Y5z+nDHOb3557KdPPbReu3dpFQQaYJQzRMSaU0yVH24bpKh1nLfBSdz0+g0Xv/Pdp5ZuLHVjquUahpNEKr5EvvCJc9b3V/n3FA3PHhLiQi/m9Cf607pyYzFP/PCl5tb5bhKqabRNgjVMgOvsOaxXvwEvDASRtwEYx6AqM4tOqyI8PjEgVRU1vDMwk2EOe1MPTO9lYJWSvlDaxCq5U6/E+7+EYbeAMtfg+eHwtdPQ2Vpiw5rswl/njSYiwYl8fjHG3hr6Y5WClgp5Q9NEKp1RHeDi5+D25dC+lmw+HF4fhhkvQ411c0+rMNu49mrMzmvfxd+98Fa5mXntmLQSilfNEGo1tX5JLhmNtz0mTXx0Ef3wIunwU8fQzN7JIU4bLxw3TDO7JvIA/NW8e9Vea0ctFLKm4AmCBEZKyIbRWSLiEz3sl1E5Hn39tUiMsxj23YRWSMiK0VEb48+0fQ81UoSV8+2EsM718GssZCzrFmHC3PaeemG4YzoFc+9c1by+fo9rRywUqq+gCUIEbEDM4BxwADgWhEZUK/YOKCve5kGvFhv+znGmMyGbgNXxzkR6D/Buuw04Vk4sA1eOx/euR72Nb1nUkSIg9cmjyAjJZY7Zq9gyaaCAAStlKoVyBrEKGCLMWarMaYSeAeYWK/MROBNY1kKdBKRpADGpILB7rB6N939ozVL3davYMYp8NG9UNK0mkB0mJM3p4yiT5copr2VxdKthYGJWSkV0G6uKUCOx+tc4BQ/yqQA+YABFoqIAV4yxrzs7U1EZBpW7YOePXu2TuQqMEIi4awHYPgUWPJnyJoFq+ZYvaAGXAp2J9jsYHPUW2rXOcHmIDbCyVs3j+Lql5dy89+X89bUUxjWM85qDD9cDIdLPJbiox8riuttK8YcLkE694fRv4Yu/YJ9lpQ6bgQyQXgbRKd+K6WvMqONMXki0gX4XER+MsYsOaawlTheBms015YErNpIVGcY/zSccit88Rh8/Sdr8ZuQYHPwuc1BmQ2qZ9mostXgdFU0umcNNsoIp1QiKDHhlJhwDpkQRua9R9iqf1J10kWEnPMAJA1p/s+nVDsRyASRC/TweN0dqN/9pMEyxpjax70iMh/rktUxCUKdwBJ6w1VvwO41ULgFXDXWXBNHLfXX1VjDeriqEVc1pvwwi9bkUFgh1oc94ZQSTqUjikp7FNXOSGqc0bicUdSExGALjSA8xGHNReGekyLEYePVzVsZuWcOkzd+Rsimjyjqfg4xF0xHep4a7LOkVNAEbD4IEXEAm4BzgV3AcuA6Y8w6jzIXAXcC47EuPz1vjBklIpGAzRhT4n7+OfCYMeZTX++p80F0TIerazhYVlU3MZHTLs0aBXZ9XjHvf7eWmDVv8Es+Il4OsTtuBNEX/g+RJ//CanQ/HhXtAmc4RMQHOxJ1AgrahEEiMh54DrADs4wxT4jIrQDGmJli/Re/AIwFyoApxpgsEUkH5rsP4wDeNsY80dj7aYJQraH0cDUfZ2+h6NtXuKR0Hl3lIDvDB1B9xn2knXYZYgvy7UPVlbDze9jyOWz+HAp+AmcEnPlfcNqdVrLwU1F5FQvW5FNcXsXwXnEMTInVaV87GJ1RTqlmWrN9N5sXvsyoXW/SXQrYak8jb9DtDL3wRiLDQ5t0LGMMe0sOs7WglG37Stm27xDb9pWxbd8h4iJCGD8oiYsGJ9E1JuzYnYt2HUkIW7+CykNWo32v06HPeZC7HDZ8CJ16wgWPQ/9LGqzxuFyG734u5N2sHD5bt5vD1a66bSF2G4O6xzIiNY4RveIZ3iuO+MiQJv2cqhUZY11SrS6HqoqGH0Xg5HHNegtNEEq1UHFpGWs+eZWe62fSw7WLrSaZFT2n0P/Cm8jonnhU2aKyKrbuO+ROAkcvZZU1deVCHTbSEiNJTYhkx/4yNuQXIwIjU+O5eGAiF8fl0mnXYti8CPa6r8zGdIe+51tL2hgIjT7yxtuWwKe/gT1rodcZMO4p6DaobvPOwjLmZefw3opd7DpYTmy4k4mZyVw1ogfdYsPI3nGA7B0HyNq+nzW7iqiqsT4beneOtJJFahwjU+NJTYjQiZxaQ0Wxley3fA57f4Kqcu8JwLgaPVRVWCLO6T83KwxNEEq1ElNTzbZv/kno98+RcngLOa7OfBRzFbuTz2frgWo2FVaxp8xFbQc9u03oERdOWmIkaYlRpCVGWI+dI0mKCcNmO/JBu337FrZ8O5+Q7V+SWfUjMVJONXYK44cRPWg8ERnjrEmafH0411TDijfgy8eh4iBVmb/ik843M3tNKT9s248IjOnbmStHdOe8/l0bvJxUUVXD6twisnbsJ3v7AbJ2HKCo3JrzIyEyhOG94qxaRmo8A5NjCXHoqD2NMgb2rHPXBBdBzlJwVeMKjaEkLoPDEkYFIVQQQrkJodTlpNTl5FCNk+JqO8XVDoqq7ByscnCwym6VNU4qCCE0Ipr3fje5WWFpglCqtRnDoTUfU7boKboUrzlmc40tBByh2JxhiD0UHB6Lt9cFG2GP+zjRSRxMOZslJpOXc3uyttBgtwmj+yQyYXASFw7oRmyE00dohpWbtnNo4ROcuu99ygnlzdBrsJ8yjUtHpJIU638bRS2Xy/BzwSGydhxg+fb9ZO84wI7CMsCqCaXEhYOx+qi7jMEYMLgfzZG4jHXqjmxzvw6xCzHhTmLdi+fz2GPWO+q2hzqO8/YSz1rC5kVQYnXkLOnUj5Vho/igZAD/Kkym2kuH0ugwBzFh1s8cE+ZwP9aeh6O3dYoIYVRa8zopaIJQKlCMgR3/gb0brNn1qiugptJ6rK6s9/qwtdQcPvK8dp+YZKstoe8F0DWjrpZgjGF9fjEfrc7no9V55Owvx2kXzuzbmQmDkzh/QFeiw6xksbuogvdW5PJedi5b95USGWJn8smV3FL2Kp12fQ0JfWHsk9blqVawt6Sirnaxp7gCEUGwQrcerdcICOKx3uO1u0BltYviiiqKyqsoLrcei8qrjrok502Y00ZsuJPIEAd2m+Cw23DYBIddcNgEu01w2m3WNpvgsNmw2wWnTbDbjpSNCXfSu3MUfbpE0btzZN05bTJjYO96q61o8+d1tYRqZzQ/R4/k86rBzC7sS74rjlCHjZGp8ZzeJ4Eh3TsdlQSjQq2fpy1oglCqHTDGsDq3iI9W5/Hx6nzyiioIcdg4+6TOHK528c3mAlwGRqXFc+Xw7owflERkqPub6aaF8NlvrPtN+pwPF/7RGnn3OFdV4zqSMMoOU34gn5oDOzEHc3AU5xBauouIsjwcNeXWJRoJo1zCKCecckIpw3pdZkIpNaGUEkapK4RDJowSl7WuxBXKvgqDq6YGGwbB0C06hN6J4aQlhJGWEE6qe0kIdyDG5a4GuY4shVusWsKWL6DYmqN9X+RJfG8byj8P9GNZdW+Mzcng7rGM7p3I6X0SGNYz7rjoMaYJQql2xuUy/JhzgI9W57NgTT4Om43Lh6UwaXh3eiVEet+puhKWvwJfPQVVZTBqGpz1IIR3atvgG+KqgZJ8OLgTDuZYj0Wez3Ot2pen8DiI7QGhMVBVak1SVVlm9fKqLAVX682V3phKRxRrQofyr9IMPq0YxF7iOLlrNKf3SWB070RGpccT09yaSQBpglBKHXGowJrQKfsN6+a6X/wWht1ojXkVSDVV1rfroxKA+/HgTmubq97kUpFdoFMPKwl06nlkie1hrffsxeVNdaVH4qi3eK6vqQSxg9iOLDYbBhtFh2vYW1LJ3kNV7C6pZHdxJfnFVRRXVFODDRc29pkYVpo+dO0Uzeg+CYzuk8hpvRPoEu2ly/JxRhOEUupY+autbrE7vrVutAuPh4g461t5eJz1OjzOSiLe1oV1AofHPRLVh61v+bUf+EclgByrgfaoLpsC0UnuD/0eR3/wd+oFsd2bdNNfWysqq2JLQQlb9h7CGDitdwI940+8LsCaIJRS3hljzfa383so2w/lB6C89tG91P9W7ykkykoaNZVwaPfR28QOMSleEoD7dUz3oxOMCgpfCSKQg/UppY53tZM69Z/gfbsx1tDo9RNH2X4oP3hkvc0OsT2PTgbRydZcIOqEpb89pVTDRCAsxlriegU7GtXG9PZHpZRSXmmCUEop5ZUmCKWUUl5pglBKKeWVJgillFJeaYJQSinllSYIpZRSXmmCUEop5VW7GmpDRAqAHc3cPRHY14rhtDaNr2U0vpbR+FrmeI6vlzGms7cN7SpBtISIZDU0HsnxQONrGY2vZTS+ljne42uIXmJSSinllSYIpZRSXmmCOOLlYAfQCI2vZTS+ltH4WuZ4j88rbYNQSinlldYglFJKeaUJQimllFcdKkGIyFgR2SgiW0RkupftIiLPu7evFpFhbRxfDxFZLCIbRGSdiPzaS5mzRaRIRFa6l4fbOMbtIrLG/d7HzO8azHMoIid7nJeVIlIsIvfUK9Om509EZonIXhFZ67EuXkQ+F5HN7se4Bvb1+fcawPieFpGf3L+/+SLSqYF9ff4tBDC+R0Rkl8fvcHwD+wbr/M3xiG27iKxsYN+An78WM8Z0iAWwAz8D6UAIsAoYUK/MeOATQIBTgR/aOMYkYJj7eTSwyUuMZwMfBfE8bgcSfWwP6jms9/vejXUTUNDOHzAGGAas9Vj3Z2C6+/l04E8NxO/z7zWA8V0AONzP/+QtPn/+FgIY3yPAfX78/oNy/upt/wvwcLDOX0uXjlSDGAVsMcZsNcZUAu8AE+uVmQi8aSxLgU4iktRWARpj8o0xK9zPS4ANQEpbvX8rCeo59HAu8LMxprl31rcKY8wSYH+91ROBN9zP3wAu9bKrP3+vAYnPGLPQGFPtfrkU6N7a7+uvBs6fP4J2/mqJiABXAf9s7fdtKx0pQaQAOR6vczn2w9efMm1CRFKBocAPXjafJiKrROQTEclo08DAAAtFJFtEpnnZfrycw2to+B8zmOcPoKsxJh+sLwVAFy9ljpfzeBNWjdCbxv4WAulO9yWwWQ1cojsezt+ZwB5jzOYGtgfz/PmlIyUI8bKufh9ff8oEnIhEAe8B9xhjiuttXoF12WQI8DfggzYOb7QxZhgwDrhDRMbU2x70cygiIcAlwFwvm4N9/vx1PJzHh4BqYHYDRRr7WwiUF4HeQCaQj3UZp76gnz/gWnzXHoJ1/vzWkRJELtDD43V3IK8ZZQJKRJxYyWG2Meb9+tuNMcXGmEPu5wsAp4gktlV8xpg89+NeYD5WVd5T0M8h1j/cCmPMnvobgn3+3PbUXnZzP+71Uiao51FEbgQmANcb9wXz+vz4WwgIY8weY0yNMcYFvNLA+wb7/DmAy4E5DZUJ1vlrio6UIJYDfUUkzf0N8xrgw3plPgR+5e6JcypQVHspoC24r1m+Bmwwxvy1gTLd3OUQkVFYv8PCNoovUkSia59jNWaurVcsqOfQrcFvbsE8fx4+BG50P78R+JeXMv78vQaEiIwFHgQuMcaUNVDGn7+FQMXn2aZ1WQPvG7Tz53Ye8JMxJtfbxmCevyYJdit5Wy5YPWw2YfVueMi97lbgVvdzAWa4t68BRrRxfGdgVYNXAyvdy/h6Md4JrMPqlbEUOL0N40t3v+8qdwzH4zmMwPrAj/VYF7Tzh5Wo8oEqrG+1NwMJwBfAZvdjvLtsMrDA199rG8W3Bev6fe3f4Mz68TX0t9BG8b3l/ttajfWhn3Q8nT/3+r/X/s15lG3z89fSRYfaUEop5VVHusSklFKqCTRBKKWU8koThFJKKa80QSillPJKE4RSSimvNEEodRxwjzL7UbDjUMqTJgillFJeaYJQqglE5Jcissw9hv9LImIXkUMi8hcRWSEiX4hIZ3fZTBFZ6jGvQpx7fR8RWeQeMHCFiPR2Hz5KROa552KYXXvHt1LBoglCKT+JSH/gaqxB1jKBGuB6IBJr7KdhwNfA7927vAk8aIwZjHXnb+362cAMYw0YeDrWnbhgjd57DzAA607b0QH/oZTywRHsAJQ6gZwLDAeWu7/ch2MNtOfiyKBs/wDeF5FYoJMx5mv3+jeAue7xd1KMMfMBjDEVAO7jLTPusXvcs5ClAt8G/sdSyjtNEEr5T4A3jDG/OWqlyO/qlfM1fo2vy0aHPZ7XoP+fKsj0EpNS/vsCmCQiXaBubuleWP9Hk9xlrgO+NcYUAQdE5Ez3+huAr401v0euiFzqPkaoiES06U+hlJ/0G4pSfjLGrBeR32LNAmbDGsHzDqAUyBCRbKAIq50CrKG8Z7oTwFZginv9DcBLIvKY+xhXtuGPoZTfdDRXpVpIRA4ZY6KCHYdSrU0vMSmllPJKaxBKKaW80hqEUkoprzRBKKWU8koThFJKKa80QSillPJKE4RSSimv/j+dFJhhh30LHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXyU9bX48c+Zyb4nJCCCARQEN9yiuAsq4r7iShdXRGvb21avWmut7bW1t+2t/toKIlLXqnUpxYqVRQEFWWUT3BAVEkWWyb7PzPn9McFGMpNMJvNkZjLn/Xr5IpPnPDOnaWZOvruoKsYYY5KXK9YJGGOMiS0rBMYYk+SsEBhjTJKzQmCMMUnOCoExxiS5lFgn0F3FxcU6dOjQWKdhjDEJZfXq1btUtSTYtYQrBEOHDmXVqlWxTsMYYxKKiHwe6pp1DRljTJKzQmCMMUnOCoExxiQ5KwTGGJPkrBAYY0ySS7hZQ8YYk2xmrangd69/yBdVjexbkMntE0Zy0ZGDovb8VgiMMSaOzVpTwV0vb6Cx1QdARVUjd728ASBqxcAKgTHGxKGfzdrAs8u34QtyVEBjq4/fvf6hFQJjjOmrfjZrA08v29ppzBdVjVF7PSsExhgTJyY9+g5LPvGEFbtvQWbUXtdmDRljTBzoThHITHVz+4SRUXttaxEYY0yMhNMFtLdBNmvIGGP6hkiKwLeOK+V/Ljos6rlYITDGmF7UnS6gPdwiXDVmP0eKAFghMMaYXhNJETjxgCKeufF4hzIKsEJgjDEOG/9/C/l4R3237+uNIgBWCIwxxlGRFAGnxgJCsUJgjDEOGHP/PL6qbenWPQJMClIE2q8ydmK8wAqBMcZE0aw1FfzX82u7fd/erYBQ4wk+1a9nG0WrGNiCMmOMiZJRd89xtAi09+zybd1+nVCsRWCMMT0UyZoA6FgAOhtPUPUDIBL4+z3YZnSRcqxFICIzRWSHiLzXRdwxIuITkYlO5WKMMU4Zdfccx4tA8/bNbH/qNurWvR5xnp1xskXwOPBn4MlQASLiBn4LOPO/zhhjHBLJmgD4ZgHo6jl8TXVULX6KujVzcGXl48rIjTjfzjhWCFR1sYgM7SLs+8BLwDFO5WGMMdE27M5XiaRjZk8RmLWmgttfWEurP3icqlL/3htULpyJv7GW3KPPo+CkSbgycr6OGRTF3UdjNkYgIoOAi4HT6KIQiMhkYDJAaWmp88kZY0wQkS4Me/CKI77eJG70vf+mptkXMrZlx6d45k2luXwTafuOpN/lvyRtwAHfiOlLu48+CNyhqj4R6TRQVacD0wHKysqiN0JijDFhGnrnq92+R4A/XnEEAEfcN5eqxtaQsf7mBqrefoba1a/gysih6KwfkDP6jK8Hh/foa7uPlgHPtRWBYuAcEfGq6qwY5mSMMV+btaaCn768noZQfTid+NZxpZQNKeK+VzZS2RC6AKgqDe8vovKNx/DVV5FzxAQKTvkO7sy8oM/Zp3YfVdVhe74WkceBf1kRMMbEi0gHgyHwgQ3wo+fXdjqW0LJrK55502jeup60fYZTcsnPSN+3Y5dPwu4+KiLPAmOBYhEpB+4FUgFUdZpTr2uMMT3VVT9+KA+2dQN11QrwtzRSvfQ5albOwpWaQdGZt5Bz+ATE5e4Qm9C7j6rqVd2IvcapPIwxJlyRLgzbMxbwwqqtnbYiVJWGD5dQ+cYMfLW7yD7sDArHXos7Kz9ovO0+aowxvWTWmgp+8ve1+CKYivLZA+fys1kbuuwGavVU4Jk3jabP1pBaMpTiC/6bjMEHB4213UeNMaYXRbJL6B6DCjK7nE3kb22i+p0XqFnxEuJOo/D0yeQedW7QbqBQu486zQqBMSYp9WQweET/bD7eUU9FVWOncQ0fL8ezYDq+6q/IPngsBeOuIyWnqEPcnq6laE4J7Q4rBMaYpDP8rlfxRtANJIAIXS4qa63aTuX8R2j8ZCWp/UopvurXZJSODhqbnuLit5eOjlkRACsExpgkEulZARD4wG72+uls00/1tlC9/CVqlr0A4qJg7HXklV2AuDt+1A7ITWP53eMjyiXarBAYY5LCqLvn0BTJaHCbZm/ni8oat6zGM38a3sovyRp1MoXjriclrzhobG/NBgqXFQJjTJ/Wk7EAgS43l/PW7KBywQwaPlpKStEg+l/+KzKHHRk0dkT/bOb9eGxEuTjJCoExps+KZH+g9jorAuprpWblLKqXPgcKBad8h7xjLkZSUoPG9/aU0O6wQmCM6XMiXRgWrsbP1+GZOxWvp5zMEcdRdPqNpOQPCBobb91AwVghMMb0GT0ZDA6Ht3YXlW/OpOH9xaQU7EPJxHvJOiD4LvqpLvjdZbGbEtodVgiMMX1CpPsDhUN9XmpXv0LVkr+hPi/5J15F3piJuFLTg8YnQiugPSsExpiE5nQroGnbe3jmTqV11+dk7l9G4Rk3kVo4MGhsPI8DdMYKgTEmITldAHz1lVS+OZP6jW/iziuh5OK7yRxxHKEO0krUIgBWCIwxCSjSIyPDoX4ftWvmUPXW02hrM3nHX07+cZfjSssIGv/ZA+c6kkdvskJgjEkokR4cH47mig/YPfdhWndsIWPIERSNn0Jqv8FBY1MENv8m8YsAWCEwxiSIniwM64qvoZqqRU9Qt34u7pwiii+4g6xRJ4XsBnowhhvEOcEKgTEmrjm5JkDVT92616la9AT+lkbyjr2E/BOuxJWeFTQ+wy18cP85juQSS04eVTkTOA/YoaqHBrk+Cbij7WEdcLOqrnMqH2NM4nFySmjzlx/jmfcwLV9+TPp+h1I0/mbSSoaEjO9rrYD2nGwRPA78GXgyxPVPgVNVtVJEzgamA2MczMcYk0B6uj1EKL6mOqoWP0ndmtdwZefT77yfkH3w2JDdQH21FdCek2cWLxaRoZ1cX9ru4TIg+IiMMSap9OTEsM6o+qnfsIDKhX/F31RH7tHnUXDyt3ClZ4e8py+3AtqLlzGC64HXQl0UkcnAZIDS0tLeyskY04ucKgAALTu24Jk7jeaKTaTvO4qiM28hbcD+IePj6ayA3hDzQiAi4wgUgpNCxajqdAJdR5SVlTk1c8wYEwNOLgzzN9dT9dYz1L77L1wZOfQ7+4dkH3Y6Iq6g8Wlu4X8nHp4UrYD2YloIRGQ0MAM4W1V3xzIXY0zvc6oIqCr1mxZS9eZMfPVV5BxxFgWnfAd3Zm7IexJ5ZXBPxawQiEgp8DLwbVX9KFZ5GGNio6cnhoXSsmsrnnlTad66gbR9RlBy6T2kDzwwZHy8HhbTm5ycPvosMBYoFpFy4F4gFUBVpwE/B/oBD7eN1ntVtcypfIwx8cGp7SH8LY1UL3mWmlX/xJWWSdGE75Ez+kzE5Q55T7IMBnfFyVlDV3Vx/QbgBqde3xgTX5xaGayqNHy4hMoFj+Kr2032YeMpHHsN7qz8kPckczdQMDEfLDbG9H1OLQxr3V2OZ/4jNH22htT++1Ny0Z2kDzooZHxeupv1950V9TwSnRUCY4yjnFgY5m9tovqdv1Oz/GUkJY3CM24i98hzOu0GSrTDYnqTFQJjjCOcWBegqjRuXo5n/nR8NTvIPmQchWOvw51TGPKe7DQ39198mI0FdMIKgTEmqpzaJK61ajuV86bRuGUVqcWlFF/9ABn7ddjG7GvJtiisJ6wQGGOiwrE1Ad4Wqpe9SPWyFxB3CoXjriP36AsQd/CPL7fAHy632UDdYYXAGNNjTrUCGj9ZiWf+I3irtpM16mQKT7uelNzikPE2GygyVgiMMT3ixFiAt3oHngXTafx4GSlFg+l/xf+QOfSIkPG2KKxnrBAYYyLiRCtAfa3UrPgH1UufB4GCU79L3jEXIe7UkPdYK6DnrBAYY7rFsW6gz9bimTcNr6eczAOPp+j0G0nJ6x8yPivVxa8vGW1jAVFghcAYEzYntofw1u6i8o3HaPjgLVIKBtJ/4i/IPCD0bjN96dD4eGGFwBgTlmF3vko0t4hTn5fa1bOpWvIs+H3knzSJ/DGXIilpIe+xvYGcYYXAGNOlaK8Obtq6Ac+8qbTu2krmAcdQeMZNpBbsEzI+GY6LjCUrBMaYoJxYF+Crq6Ry4UzqN76JO68/JZfcQ+bwY0OeFwzWCugNVgiMMR0Mv+tVvFHsB1K/j9o1c6ha/BTqayHv+CvIP/4yXKkZIe+xvYF6jxUCY8zXnJgR1FzxPrvnTqV1xxYyhh5J0fgppBaF/gvftobofVYIjDFA9BeG+RqqqVz4OPUb5uHO6UfxhXeSNfLETruBbE1AbFghMCbJRXssQP0+6ta9TtXiJ/G3NJJ37CXkn3gVrrTMkPfYOEBsOXlU5UzgPGCHqnbYIlACfxY8BJwDNADXqOq7TuVjjOko2lNCm7/8CM/cqbRs/5j00sMoOmMKaSVDQsbbQTHxwckWwePAn4EnQ1w/GxjR9t8YYGrbv8YYh0V7LMDXWEvV4ieoW/s67uwCis+/jayDTrVuoATh5JnFi0VkaCchFwJPqqoCy0SkQEQGquqXTuVkTLKLejeQ+qnfMJ/KhY/jb6ojt+wCCk66Gld6dsh7UlzC7y873LqC4kgsxwgGAdvaPS5v+16HQiAik4HJAKWlpb2SnDF9TbRbAS1fbcEz92Gav/iA9EEHU3TmFNL67x8yPtUFv7vMxgLiUSwLQbA2Y9DuSlWdDkwHKCsri2aXpjFJIZorg/3N9VS99TS1776KKzOXfuf8F9mHnoaIK+Q9Nhgc32JZCMqB/do9Hgx8EaNcjOmTojklVFWp37SQyjcfw19fTc6R51BwyrdxZ+SEvMfGARJDLAvBbOBWEXmOwCBxtY0PGBMd0R4LaNn5GZ5502je9h5pAw+kaOIvSN9neMh4OygmsTg5ffRZYCxQLCLlwL1AKoCqTgPmEJg6upnA9NFrncrFmGQSza2i/c0NVC95lppV/8SVnk3RhFvJOfzMTruBbGuIxOPkrKGruriuwPecen1jklG0xgJUlYYP3qbyjUfx1XnIGX0mBad+F3dWfsh77KCYxGUri43pA6I5FtC6uxzPvGk0fb6WtAEHUHLRT0kfNCpkvB0Uk/isEBiTwKI5FuBvaaL6neepWfEPJDWdovFTyDnibMTlDnmPzQbqG6wQGJOAZq2p4Cd/X4svCpOpVZXGj5fhWTAdX81Osg89jcKx1+LOLgx5jx0U07dYITAmwUSzFdBa+SWe+dNo2rKa1OIhFF/9ABn7ddga7GvZaW7uv/gwawX0MZ0WAhEp6uy6qnqim44xJpRotgL8rc3ULH+R6mUvIu4UCk+7gdyjzkPcoT8SrBuo7+qqRbCawGpfAUqByravC4CtwDBHszPGADDq7jk0RaMCAA2frKRy/iN4q7aTddCpFI67jpTcfiHjbU1A39dpIVDVYQAiMg2Yrapz2h6fDZzhfHrGJLdo7g/krf4Kz4JHafx4GSlFg+l/5f1kDjm803tsZXByCHeM4BhVnbLngaq+JiK/cignYwzRmxKq3lZqVrxM9Tt/B4GCU68h75gLEXdqyHusGyi5hFsIdonIz4CnCXQVfQvY7VhWxiSxaA4GN366Bs/8aXg9FWQdeAKFp99ASl7/kPHWDZScwi0EVxHYIuIfBArB4rbvGWOiKFonhnlrdlH5xgwaPnyblMKB9L/sPjL3PzpkvJ0RkNzCKgRts4N+KCI5qlrncE7GJJ1Jj77Dkk96PglPfV5qVv2T6iXPgvrJP2kS+WMuRVLSQt5j4wAmrEIgIicAM4AcoFREDgduUtVbnEzOmL4umt1ATVvX45k7jdbdW8kcfiyFp08mtWCfkPF2UIzZI9yuoT8CEwhsHY2qrhORUxzLypgkEK1dQr11HqrenEn9poW48wdQcuk9ZA0Pffx3mlv434nWDWT+I+yVxaq6ba+DqH3RT8eY5BCNXULV76P23X9R9dYzqK+F/OOvIO/4y3ClZoS8x2YDmWDCLQTb2rqHVETSgB8A7zuXljF9U7RaAU3lm/DMfZjWnZ+RMewois64idSi0B/wNg5gOhNuIZgCPETgcPlyYC5g4wPGhClaYwG++ioqFz5O/XvzcecWU3zRXWQdeAJ7tda/wVoBpivhFoKRqjqp/TdE5ERgSfRTMqZviUYrQP0+6ta9TtWiJ/C3NpE3ZiL5J1yBKy0z5D1WAEy4wi0EfwKOCuN73yAiZxFoSbiBGar6wF7X8wksUitty+X3qvrXMHMyJu5FYyyg+cuP8MydSsv2j0kvHU2/8TeTWrxfyPi8dDfr7zurx69rkkdXu48eD5wAlIjIj9tdyiPw4d7ZvW7gL8B4At1JK0Vktqpuahf2PWCTqp4vIiXAhyLyjKpG56glY2IkGusCfI01VC16krp1r+POKaT4/NvJOuiUkN1AdkaAiVRXLYI0AmsHUoDcdt+vASZ2ce+xwGZV3QIgIs8BFwLtC4ECuRL4zc4BPIA37OyNiTPR2CRO1U/d+vlULXocf1MduWUXUHDSJFzpWSHvsVaA6Ymudh9dBCwSkcdV9fNuPvcgYFu7x+XA3pOb/0xgbcIXBArNFarq7+brGBNzs9ZU8KPn1/Z4e4iWrz5h99yHafniQ9IHH0zR+JtJ6x96t3c7KMZEQ7hjBDNE5DJVrQIQkULgOVWd0Mk9wdqve79PJgBrgdOAA4B5IvKWqtZ844lEJgOTAUpLS8NM2ZjeEY0ZQf6mOqreepraNXNwZebR79wfkX3IaTYbyPSKcAtB8Z4iAKCqlSISegvDgHKg/YjWYAJ/+bd3LfCAqiqwWUQ+BUYBK9oHqep0YDpAWVlZdE7nMCYKetoVpKrUb3yDyjf/ir+xhtwjz6bg5G/jysgJeY8VABNt4RYCv4iUqupWABEZQse/7ve2EhghIsOACuBK4Oq9YrYCpwNvicgAYCSwJdzkjYmVaIwFtOz8DM/cqTSXbyRt4EiKLvsF6fsMDxl/4gFFPHPj8T16TWOCCbcQ3A28LSKL2h6fQltXTSiq6hWRW4HXCcwwmqmqG0VkStv1acCvgMdFZAOBrqQ7VHVXBP87jOkVUekGam6gasnfqF01G1dGDkVnfZ+c0eMRcQWNF+CP1gowDpJAr0wYgSLFwHEEfi/fidUHdllZma5atSoWL22SXE+nhKoqDe8vpvLNx/DVVZJz+JkUnPpd3Jl5Ie+xrSFMtIjIalUtC3atq3UEo1T1AxHZs3BsTx9/aVtX0bvRTNSYeBSNbqDW3dvwzJtK0+frSRtwACUX3036viNDxls3kOlNXXUN/QS4EfhDkGtKYLaPMX3W6Hv/TU1z5Bvt+luaqH7nOWpWzMKVmk7R+JvJOeIsxBV8Paa1AEwsdLWO4Ma2f8f1TjrGxIeejgWoKo0fvYNnwaP4aneSfegZFI69Bnd2Qch7bDaQiZWuuoYu6ey6qr4c3XSMib2enhvc6qnAM386TZ+uJrVkKMUX3EbG4ENCxlsBMLHWVdfQ+W3/9iew59AbbY/HAQsBKwSmz+jpYLC/tZmaZS9QvfxFxJ1K4ek3knvUeSG7gVIENv/m3Ihfz5ho6apr6FoAEfkXcLCqftn2eCCBDeWM6RNG3T2HJl/k7YCGzSuonP8I3uqvyDr4VArHXU9KTlHIeGsFmHgS7jqCoXuKQJuvgAMdyMeYXtXTsYDWqu1ULphO4+YVpPbbjwFX/pqMIaNDxg/ITWP53eMjfj1jnBBuIVgoIq8DzxKYLXQl8KZjWRnjsB4PBntbqV7xEjXv/B3ERcHYa8kruwBxpwaNT3EJv7/MDow38SmsQqCqt4rIxQRWFANMV9V/OJeWMc6YtaaCn/x9LT3oBaLx03fxzJuGt/ILskaeSOFpN5CSVxIy3qaEmngXbosA4F2gVlXni0iWiOSqaq1TiRkTbT1tBXhrdlK54FEaPlpKSuG+9L/8l2QOC31InxUAkyjCKgQiciOBvYWKCGwXPQiYRmDDOGPiXk9WB6uvlZpV/6R6yXOgSsHJ3ybv2EuQlODdQG6BP1xug8EmcYTbIvgegRPHlgOo6sdhbENtTMz1tBXQ9Pl6PPOm0rp7G5nDx1B0xmRS8geEjLfZQCYRhVsImlW1Zc8hGSKSQtfbUBsTMz09Mcxb56HyjcdoeH8RKfkDKLn052QNPzZkvHUDmUQWbiFYJCI/BTJFZDxwC/CKc2kZE7metALU76N29b+oevtp1Ocl/4SryDtuIq7U9KDx1g1k+oJwC8EdwA3ABuAmYA4ww6mkjInU+P9byMc76iO6t6l8I565U2nd+RkZw46maPxNpBbuGzQ2wy18cP85PUnVmLjRZSGQwGkZ61X1UOBR51Mypvt60grw1VdRufCv1L+3AHduCSUX/5TMEceHPC/YuoFMX9NlIVBVv4isa39UpTHxoidjAer3Ubf231QtfhJ/azN5x00k//grcaVlBI23AmD6qnC7hgYCG0VkBfB1u1tVL3AkK2PC0JNuoOYvPsQz92FavvqEjCGjKRp/M6n99gsZb7OBTF8WbiG4L5InF5GzgIcInFk8Q1UfCBIzFngQSAV2qeqpkbyWSR492SXU11hD1aInqFs3F3dOIcUX/DdZo04O2Q1kJ4WZZNDVeQQZwBRgOIGB4sdU1RvOE4uIm8AOpeOBcmCliMxW1U3tYgqAh4GzVHWrrU0wXYn0xDBVP3Xr5lK16An8zfXkHXMR+SdehSs9K2i8HRhvkklXLYIngFbgLeBs4GDgh2E+97HAZlXdAiAizwEXApvaxVwNvLxn7EFVd4SfukkmPRkMbt6+Gc/cqbR8+SHpgw+h6MybSSsZGjQ2O83N/RcfZgXAJJWuCsHBqnoYgIg8BqzoxnMPAra1e1wOjNkr5kAgVUQWArnAQ6r65N5PJCKTCWxxQWlpaTdSMH3BmPvn8VVtS7fv8zXVUbX4KerWzMGVlU+/c39M9iHjgnYDWQvAJLOuCkHrni9U1RuqHzWEYMF7T+5IAY4msGdRJvCOiCxT1Y++cZPqdGA6QFlZma1oThKRDgarKvXvvUHlwpn4G2vJPfo8Ck6ahCsjJ2i8DQSbZNdVIThcRGravhYCK4tr2r5WVc3r5N5yoP00jMHAF0FidqlqPVAvIouBw4GPMElt+F2v4o2g5Lfs+BTPvKk0l28ibd+R9Lv8l6QNOCBorA0EGxPQ1VGVwQ9bDc9KYISIDAMqCBxmc/VeMf8E/ty2d1Eaga6jP/bgNU2Ci3SXUH9zA1VvP0Pt6ldwZeRQdNYPyBl9BoH1kB1ZK8CY/+jOeQTd0taVdCvwOoHpozNVdaOITGm7Pk1V3xeRfwPrAT+BKabvOZWTiW+RzAhSVRreX0zlm4/hq6sk54gJFJzyXdyZuUHjbVGYMR2JamJ1uZeVlemqVatinYaJkp6cGNa6axu7502leet60vYZTtGZt5A+MPRR2tYKMMlMRFaralmwa461CIzpSsTdQC2NVC99jpqVs3ClZVJ05i3kHD4BcQXvybQCYEznrBCYmIhkSqiq0vDRUioXPIqvdhfZh51B4dhrcWflB41PT3Hx20tHWxEwpgtWCEyvinRhWKunAs/8R2j69F1S+w+j+IL/JmPwwUFjR/TPZt6Px/YwU2OShxUC0ysi7gZqbaLmnReoXvES4k6j8PTJ5B51btBuoMKsVO49/xBrARjTTVYIjOMiLQINm5fjmT8dX/VXZB8yjoKx15KSU9QhLi/dzfr7zopGqsYkJSsExjGRFoDWqu1Uzn+Exk9WktqvlOKrfkNGafApnzYd1Jies0JgHBHJ9hDqbaF6+UvULHsBxEXB2OvIK7sAcXf8NbVVwcZEjxUCE1WRtgIat6zGM38a3sovyRp1MoXjriclr7hD3IDcNJbfPT4aqRpj2lghMFETSRHw1uygcsEMGj5aSkrRIPpf/isyhx0ZNNa6gYxxhhUC02ORFAD1tVKzchbVS58DhYJTvkPeMRcjKakdYq0AGOMsKwSmRyLZH6jx83V45k7F6yknc8RxFJ1+Iyn5AzrE2TiAMb3DCoGJSCTnBntrd1P55mM0vL+YlIJ9KJl4L1kHHNMhzgqAMb3LCoHplsi6gbzUrn6FqiV/Q31e8k+8irwxE3GlpneItW4gY3qfFQITlkhnAzVtew/P3Km07vqczP3LKDzjJlILB3aIs20hjIkdKwSmS5EUAV99JZUL/0r9e2/gziuh5OK7yRxxXIfzglNd8LvLbHdQY2LJCoEJKZKzAtTvo3bNHKreehptbSbv+MvJP+5yXGkZHWJte2hj4oMVAhNUJIPBzRUf4Jk3lZavPiFjyBEUjZ9Car/BHeI+e+DcaKVpjIkCRwuBiJwFPETgqMoZqvpAiLhjgGXAFar6opM5mc5FUgB8DdVULXqCuvVzcecUUXzBHWSNOqlDN1CKwObfWBEwJt44VghExA38BRgPlAMrRWS2qm4KEvdbAmcbmxiZtaaCO15aT7PXH/Y9qn7q1s2latET+FsayDv2EvJPuBJXelaHWJsNZEz8crJFcCywWVW3AIjIc8CFwKa94r4PvAR0nFBuekUkp4U1b9+MZ+7DtHz5Een7HUrR+JtJKxnSIS7DLXxw/znRStUY4wAnC8EgYFu7x+XAmPYBIjIIuBg4jU4KgYhMBiYDlJaWRj3RZBXJaWG+pjqqFj9J3ZrXcGXn0++8n5B98NgO3UBgYwHGJAonC0HHTwbYe/7Jg8AdquoL9kHy9U2q04HpAGVlZd2Yw2JC6W4rQNVP/XtvULnwr/gba8k9+jwKTv4WrvTsDrF2UIwxicXJQlAO7Nfu8WDgi71iyoDn2opAMXCOiHhVdZaDeSW1SKaEtuzYgmfuNJorNpG+7yiKLv8VaQP27xBni8KMSUxOFoKVwAgRGQZUAFcCV7cPUNVhe74WkceBf1kRcE53ZwT5m+upeusZat/9F66MHPqd/UOyDzsdEVeHWFsTYEzicqwQqKpXRG4lMBvIDcxU1Y0iMqXt+jSnXtt01J0Tw1SV+k0LqXpzJr76KnKOOIuCU76DOzO3Q6xtEGdM4nN0HYGqzgHm7PW9oAVAVa9xMpdk1p0i0LJrK555U2neuoG0gSMoufTnpA8cETTWBoON6RtsZXEfN+ruOTSFMSDgb2mkesmz1Kz6J660TIomfI+c0WciLrbon1EAAA4FSURBVHeHWFsTYEzfYoWgjwp3PEBVafhwCZULHsVXt5uc0WdScOp3cWfld4i12UDG9E1WCPqgcFsBrZ4KPPOm0fTZGlL770/JRXeSPuigoLHWCjCm77JC0IeEu0DM39pE9TsvULPiJSQlncIzbiL3yHOCdgMJ8KmNBRjTp1kh6CPCWSCmqjRuXo5n/nR8NTvIPmQcheOuw51dGDTepoQakxysECS4WWsq+NHzazss2d5ba9V2Kuc/QuMnK0ktLqX46gfI2O/QoLE2JdSY5GKFIIGFMy1UvS1UL3+JmmUvgMtN4bjryT36fMQd/P96mxJqTPKxQpCght/1Kt4umgGNn6zCM/8RvFVfkjXqZApPu56U3OKgsdYKMCZ5WSFIQMPufLXTriBvzQ48Cx6l8aN3SCkaTP8r/ofMoUcEjbXDYowxVggSRDibxamvlZoV/6B66fMgUHDqd8k75iLEnRo03jaJM8aAFYKE8LNZG3h62dZOYxo/W4tn3jS8nnIyDzyeotNvJCWvf8h4GwswxuxhhSDOdTUt1Fu7i8o3HqPhg7dIKRhI/4m/IPOAspDxLoH/uzx4N5ExJjlZIYhTXU0LVZ+X2tWzqVryLPh95J80ifwxlyIpaSGfc1BBJrdPGGlrA4wx32CFIM6Esy6gadt7eOY+TOuurWQecAyFZ9xEasE+IeNtLMAY0xkrBHGkq7EAX10llQtnUr/xTdx5/Sm55B4yhx8b9LxgsE3ijDHhsUIQB2atqeC+VzZS2dAa9Lr6fdSumUPV4qdQXwt5x19B/vGX4UrNCPmcti7AGBMuKwQx1lUroLnifXbPnUrrji1kDD2SovFTSC0K3cdfkJnKLy44xMYBjDFhc7QQiMhZwEMEjqqcoaoP7HV9EnBH28M64GZVXedkTvFk9L3/pqbZF/Sar6GayoWPU79hHu6cfhRfeCdZI08M2Q0kwB9tkzhjTAQcKwQi4gb+AowHyoGVIjJbVTe1C/sUOFVVK0XkbGA6MMapnOLBrDUV3PbCOrz+4MPB6vdRt34uVYuewN/SSN6YS8k/4UpcaZlB49NTXPz20tFWAIwxEXOyRXAssFlVtwCIyHPAhcDXhUBVl7aLXwYMdjCfmOvq1LDmLz/GM/dhWrZ/THrpYRSNv5m04tKQ8TYOYIyJBicLwSBgW7vH5XT+1/71wGsO5hMzXRUAX2MtVYufpG7tv3FnF1B8/m1kHXSqzQYyxvQKJwtBsE+xoP0hIjKOQCE4KcT1ycBkgNLS0H8hx5uuTgxT9VO/YQGVC/+Kv6mO3LILKDhpEq70rJD3WCvAGBNtThaCcmC/do8HA1/sHSQio4EZwNmqujvYE6nqdALjB5SVlXV9GG+MzVpTwY+fX4u/k5iWr7bgmTeV5or3SR90MEVnTiGt//4h49Pcwv9OPNzGAowxUedkIVgJjBCRYUAFcCVwdfsAESkFXga+raofOZhLrwhnczh/cz1Vbz1N7buv4srMpd85/0X2oach4goaL8AkOzjeGOMgxwqBqnpF5FbgdQLTR2eq6kYRmdJ2fRrwc6Af8HBbf7hXVUPvmBanwmkBqCr1mxZS+eZj+OuryTnyHApO+TbujJyQ93zLCoAxpheIatz3tHxDWVmZrlq1KtZpAOEdFQnQsvPzQDfQtvdIG3ggRWfeQvo+w0PGZ6a6+M0lNiXUGBM9IrI61B/atrK4m2atqeCnL6+nobWzv/8D/M0NVC95lprVs3GlZVE04VZyDj8zZDeQDQQbY2LBCkGYupoB1J6q0vDB21S+MQNf3W5yRp9JwanfxZ2VHzTeCoAxJpasEIShq8Nh2mvdXY5n3jSaPl9L2oADKLnoLtIHjQoaawPBxph4YIUghHD7//fwtzRR/c7z1Kz4B5KaTtH4KeQccTbicneItRaAMSaeWCFop7sf/hDoBmr8eBmeBdPx1ewk+9DTKBx7Le7swg6xdkCMMSYeJX0hmLWmgt+9/iEVVY3dvre18ksq5z9C45ZVpBYPofjqB8jY79CgsTYV1BgTr5KyEIy6ew5Nvsinzfpbm6lZ/iLVy15E3CkUnnYDuUedh7j/8+MUgUlj7MPfGBP/kqIQRNLlE0rjJyvxzH8Eb9V2sg46lcJx15GS2w+AFIHNvzk3Kq9jjDG9pc8XgmgVAW/1DjwLptP48TJSigbT/8r7yRxyOGB9/8aYxNbnC0FPi4B6W6lZ+Q+qlz4PAgWnXkPeMRci7lQ7FtIY0yf0+ULQE42frsEzfxpeTwVZB55A4ek3MKR0CLdPGGkf/saYPsMKQRDeml1UvjGDhg/fJqVwIAdd8xt+/YNv24e/MaZP6vOFYET/7LC7h9TnpWbVbGqW/o0UUX75y19y++23k5GR4XCWxhgTO8F3P+tD5v14LCP6Z3/jewNy0xhUkIkAgwoyefCKI3ji7Bxy5vyUqoUzOXfCeD54/33uueceKwLGmD6vz7cIgE5n9Gzfvp3bbruNZ555hqFDhzJ79mzOP//83kvOGGNirM+3CELxer089NBDjBw5khdeeIF77rmHjRs3WhEwxiSdpGgR7G3p0qXccsstrFu3jgkTJvCnP/2JESNGxDotY4yJiaRqEezcuZPrrruOE088kd27d/Piiy/y2muvWREwxiQ1RwuBiJwlIh+KyGYRuTPIdRGR/9d2fb2IHOVULnPmzGHkyJE89dRT3HHHHXzwwQdceumltJ2VbIwxScuxriERcQN/AcYD5cBKEZmtqpvahZ0NjGj7bwwwte3fqDvwwAM57rjj+MMf/sBBBx3kxEsYY0xCcnKM4Fhgs6puARCR54ALgfaF4ELgSVVVYJmIFIjIQFX9MtrJDB8+nDlz5kT7aY0xJuE52TU0CNjW7nF52/e6G4OITBaRVSKyaufOnVFP1BhjkpmThSBY5/vehwCEE4OqTlfVMlUtKykpiUpyxhhjApwsBOXAfu0eDwa+iCDGGGOMg5wsBCuBESIyTETSgCuB2XvFzAa+0zZ76Dig2onxAWOMMaE5Nlisql4RuRV4HXADM1V1o4hMabs+DZgDnANsBhqAa53KxxhjTHCOrixW1TkEPuzbf29au68V+J6TORhjjOlcUq0sNsYY05EVAmOMSXIS6J1JHCKyE/i8G7cUA7scSscJiZRvIuUKiZVvIuUKiZVvIuUK0ct3iKoGnX+fcIWgu0RklaqWxTqPcCVSvomUKyRWvomUKyRWvomUK/ROvtY1ZIwxSc4KgTHGJLlkKATTY51ANyVSvomUKyRWvomUKyRWvomUK/RCvn1+jMAYY0znkqFFYIwxphNWCIwxJsn1mUIQT8didiWMXCe15bheRJaKyOGxyLNdPp3m2y7uGBHxicjE3sxvrxy6zFVExorIWhHZKCKLejvHvXLp6nchX0ReEZF1bfnGbD8uEZkpIjtE5L0Q1+PpPdZVrvH2Hus033ZxzrzHVDXh/yOwqd0nwP5AGrAOOHivmHOA1wicgXAcsDyOcz0BKGz7+uxY5Rpuvu3i3iCwt9TEeM0VKCBwSl5p2+P+8fyzBX4K/Lbt6xLAA6TFKN9TgKOA90Jcj4v3WJi5xs17LJx82/2+OPIe6ystgq+PxVTVFmDPsZjtfX0spqouAwpEZGBvJ0oYuarqUlWtbHu4jMA5DbESzs8W4PvAS8CO3kxuL+HkejXwsqpuBVDVeM9XgVwRESCHQCHw9m6abYmoLm57/VDi5T3WZa5x9h4L52cLDr7H+kohiNqxmL2gu3lcT+CvrFjpMl8RGQRcDEwjtsL52R4IFIrIQhFZLSLf6bXsOgon3z8DBxE4sGkD8ENV9fdOet0WL++x7or1e6xLTr/HHN2GuhdF7VjMXhB2HiIyjsAv6UmOZtS5cPJ9ELhDVX2BP1xjJpxcU4CjgdOBTOAdEVmmqh85nVwQ4eQ7AVgLnAYcAMwTkbdUtcbp5CIQL++xsMXJeywcjr7H+kohSKRjMcPKQ0RGAzOAs1V1dy/lFkw4+ZYBz7X9ghYD54iIV1Vn9U6KXwv392CXqtYD9SKyGDgciEUhCCffa4EHNNBJvFlEPgVGASt6J8VuiZf3WFji6D0WDmffY7EcIIniQEsKsAUYxn8G3Q7ZK+ZcvjmQtSKOcy0lcGrbCYnws90r/nFiN1gczs/2IGBBW2wW8B5waBznOxX4RdvXA4AKoDiGvw9DCT0AGxfvsTBzjZv3WDj57hUX9fdYn2gRaAIdixlmrj8H+gEPt/0F4NUY7ZYYZr5xIZxcVfV9Efk3sB7wAzNUtdMpe7HMF/gV8LiIbCDwAXuHqsZkC2UReRYYCxSLSDlwL5DaLte4eI9BWLnGzXsMwsrX2ddvqzDGGGOSVF+ZNWSMMSZCVgiMMSbJWSEwxpgkZ4XAGGOSnBUCY4xJcn1i+qgxThKRfgTWHgDsA/iAnW2Pj9XAPkHGJCybPmpMN4jIL4A6Vf19u++lqGpMNoIzJhqsRWBMBETkcQK7RR4JvCsitbQrEG37yp+nqp+JyLeAHxBYPbwcuEVVfbHJ3JiObIzAmMgdCJyhqj8JFSAiBwFXACeq6hEEupUm9VJ+xoTFWgTGRO6FMP6yP53Abqcr27YyyCS2ZzYY04EVAmMiV9/uay/fbGFntP0rwBOqelevZWVMN1nXkDHR8RmBowZpO6t3WNv3FwATRaR/27UiERkSkwyNCcEKgTHR8RJQJCJrgZtpO99AVTcBPwPmish6YB4Qk+MbjQnFpo8aY0ySsxaBMcYkOSsExhiT5KwQGGNMkrNCYIwxSc4KgTHGJDkrBMYYk+SsEBhjTJL7//ACXIqJ9fipAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE: 0.023\n"
     ]
    }
   ],
   "source": [
    "model.history[['training_loss', 'validation_loss']].plot()\n",
    "plt.ylabel('Loss')\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "y_pred = model.predict(x)\n",
    "plt.scatter(y, y_pred)\n",
    "plt.plot((y.min(), y.max()), (y.min(), y.max()), 'k-')\n",
    "plt.xlabel('True')\n",
    "plt.ylabel('Predicted')\n",
    "plt.show()\n",
    "plt.close()\n",
    "print('MAE: {:.3f}'.format(np.mean(np.abs(y_pred - y))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bad Physics Function Test\n",
    "\n",
    "Here we see that a p_fun input with numpy operations instead of tensorflow operations cannot be used in phygnn. This is because of how the stochastic gradient descent algorithm works in tensorflow. SGD finds the gradient of the loss with respect to the change in node weights. The loss function must be automatically differentiable for this to work. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_fun_bad(y_predicted, y_true, p):\n",
    "    \"\"\"This is an example of a poorly formulated p_fun() that uses numpy operations.\"\"\"\n",
    "    \n",
    "    y_physical = p[:, 0]**2 + p[:, 1]**2\n",
    "    p_loss = np.mean(np.abs(y_predicted.numpy() - y_physical))\n",
    "    p_loss = tf.convert_to_tensor(p_loss, dtype=tf.float32)\n",
    "    \n",
    "    return p_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The input p_fun was not differentiable! Please use only tensor math in the p_fun.\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "The input p_fun was not differentiable! Please use only tensor math in the p_fun.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-6886f135bbff>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      3\u001b[0m                                    \u001b[0mloss_weights\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0.5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0.5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m                                    n_features=2, n_labels=1)\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my_noise\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mc:\\sandbox\\phygnn\\phygnn\\phygnn.py\u001b[0m in \u001b[0;36mfit\u001b[1;34m(self, x, y, p, n_batch, n_epoch, shuffle, validation_split, p_kwargs, run_preflight, return_diagnostics)\u001b[0m\n\u001b[0;32m    691\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    692\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_loss_weights\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mrun_preflight\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 693\u001b[1;33m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpreflight_p_fun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx_val\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my_val\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_val\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp_kwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    694\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    695\u001b[0m         \u001b[0mt0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\sandbox\\phygnn\\phygnn\\phygnn.py\u001b[0m in \u001b[0;36mpreflight_p_fun\u001b[1;34m(self, x, y_true, p, p_kwargs)\u001b[0m\n\u001b[0;32m    431\u001b[0m                         'Please use only tensor math in the p_fun.')\n\u001b[0;32m    432\u001b[0m                 \u001b[0mlogger\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merror\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0memsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 433\u001b[1;33m                 \u001b[1;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0memsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    434\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    435\u001b[0m         \u001b[0mlogger\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdebug\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'p_fun passed preflight check.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mRuntimeError\u001b[0m: The input p_fun was not differentiable! Please use only tensor math in the p_fun."
     ]
    }
   ],
   "source": [
    "model = PhysicsGuidedNeuralNetwork(p_fun=p_fun_bad, \n",
    "                                   hidden_layers=hidden_layers, \n",
    "                                   loss_weights=(0.5, 0.5),\n",
    "                                   n_features=2, n_labels=1)\n",
    "model.fit(x, y_noise, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:mlclouds]",
   "language": "python",
   "name": "conda-env-mlclouds-py"
  },
  "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
}