preprocessing/profiling/profile_cathub_feature_space.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both CGCNN and GP modeling were too good at making predictions on CatHub's validation set. We double-checked the splitting and don't see any obvious leakage into the training set. This notebook documents our journey into figuring this out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "\n",
    "# Load data\n",
    "splits = {}\n",
    "for dataset in ['cathub', 'gasdb']:\n",
    "    with open('../splits_%s.pkl' % dataset, 'rb') as file_handle:\n",
    "        _splits = pickle.load(file_handle)\n",
    "\n",
    "    # Parse data\n",
    "    docs_train = _splits['docs_train']\n",
    "    docs_val = _splits['docs_val']\n",
    "    fingerprints_train = _splits['fingerprints_train']\n",
    "    fingerprints_val = _splits['fingerprints_val']\n",
    "    splits[dataset] = _splits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lipschitz distributions\n",
    "First, make a function to calculate the lipschitz distribution of a dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "import itertools\n",
    "import multiprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import h5py\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "\n",
    "def calculate_lipschitz(_docs, _features, suffix, force=False, processes=8):\n",
    "    '''\n",
    "    Calculates the Lipschitz values of a series of data.\n",
    "    \n",
    "    Args:\n",
    "        _docs       A list of dictionaries of parsed information\n",
    "        _features   A numpy array of features. Should map to `docs`\n",
    "        suffix      A string indicating the suffix you want to use\n",
    "                    for the cache of results that this function will use\n",
    "        force       A Boolean indicating whether or not you want to force\n",
    "                    this function to recalculate the Lipschitz values\n",
    "                    (as opposed to just loading a cache).\n",
    "        processes   An integeter indicating how many processes you want\n",
    "                    to use whem multithreading\n",
    "    Returns:\n",
    "        lips    A pandas.Series of all the Lipschitz values for each\n",
    "                combination of data points\n",
    "    '''\n",
    "    # If we're not forcing a re-do, then try to just load the data\n",
    "    if force is False:\n",
    "        try:\n",
    "            with h5py.File('lips_%s.hdf5' % suffix, 'r') as h5_handle:\n",
    "                lips = h5_handle['lips']\n",
    "                lips = pd.Series(lips)\n",
    "            return lips\n",
    "\n",
    "        # If the file is not there, then move on to the rest of this code,\n",
    "        # which will do the calculation\n",
    "        except OSError:\n",
    "            pass\n",
    "\n",
    "\n",
    "    # Push these big objects to global space for faster process sharing\n",
    "    global docs\n",
    "    global features\n",
    "    docs = _docs\n",
    "    features = _features\n",
    "\n",
    "    def __calculate_lipschitz(indices):\n",
    "        '''\n",
    "        Calculates the Lipschitz value between document with index i\n",
    "        and document with index j, where indices=(i, j). We make a\n",
    "        function so that we can multimap it.\n",
    "\n",
    "        Args:\n",
    "            indices 2-tuple contaning the indices of the documents you\n",
    "                    want to get the Lipschitz value for\n",
    "        Returns:\n",
    "            lipschitz   The Lipschitz value between document i and j\n",
    "        '''\n",
    "        i, j = indices\n",
    "        energy_diff = docs[i]['energy'] - docs[j]['energy']\n",
    "        feature_diff = features[i, :] - features[j, :]\n",
    "        lipschitz = abs(energy_diff) / np.linalg.norm(feature_diff, 2)\n",
    "        return lipschitz\n",
    "\n",
    "    # Figure out all the 2-pair combinations\n",
    "    index_combinations = itertools.combinations(range(len(docs)), 2)\n",
    "    n_docs = len(docs)\n",
    "    n_combinations = int((n_docs * (n_docs-1))/2)\n",
    "\n",
    "    # Ignore divide-by-zero warnings\n",
    "    with warnings.catch_warnings(): \n",
    "        warnings.simplefilter('ignore', category=RuntimeWarning)\n",
    "        # Multithreaded calculations\n",
    "        with multiprocess.Pool(processes=processes, maxtasksperchild=1000) as pool:\n",
    "            iterator = pool.imap(__calculate_lipschitz, index_combinations, chunksize=10000)\n",
    "            lips = list(tqdm_notebook(iterator, total=n_combinations, unit_scale=True))\n",
    "\n",
    "    # Filter out things with infinity\n",
    "    lips = [lip for lip in lips if lip != np.inf]\n",
    "\n",
    "    # Save the data\n",
    "    with h5py.File('lips_%s.hdf5' % suffix, 'w') as h5_handle:\n",
    "        h5_handle.create_dataset('lips', data=lips, dtype='f', compression='gzip')\n",
    "        lips = pd.Series(lips)\n",
    "\n",
    "    # Clean up and output\n",
    "    del globals()['docs']\n",
    "    del globals()['features']\n",
    "    return lips"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "\n",
    "# Load data\n",
    "lips_val = {}\n",
    "lips_train = {}\n",
    "for dataset in ['cathub', 'gasdb']:\n",
    "\n",
    "    # Calculate Lipschitz distribution\n",
    "    _lips_val = calculate_lipschitz(docs_val, fingerprints_val, '%s_val' % dataset)\n",
    "    _lips_train = calculate_lipschitz(docs_train, fingerprints_train, '%s_train' % dataset)\n",
    "\n",
    "    # We accidentally forgot to drop NaN values. Do so here.\n",
    "    _lips_val = _lips_val.dropna()\n",
    "    _lips_train = _lips_train.dropna()\n",
    "\n",
    "    lips_val[dataset] = _lips_val\n",
    "    lips_train[dataset] = _lips_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeyUlEQVR4nO3debgV1Znv8e9PRFE4IoMmCgTQeA0amTwBg0ZJNN0ZHBOUmGjU2CHmxtyY1nQGu1ucEtsx3TGtwRG9OA+PaGKM18hD0Gg8GGTSjhqJYggYUAQEZHjvH7vOcYPnFLU5p/auzfl9nmc/e++qtareKjfnddVatUoRgZmZWVu2q3UAZmZWbE4UZmaWyonCzMxSOVGYmVkqJwozM0vlRGFmZqmcKKzTkHStpH/r6LLtJWmBpCOqsS+zrSHfR2H1QNIC4J8i4v/VOpaOVsmxSQpgn4h4KffAzBJuUdg2QdL2tY7BbFvlRGGFJ+lW4EPAg5JWSvoXSYMkhaTTJb0K/DYpe7ekv0laLmm6pP3LtnOzpIuSz2MlLZR0tqQlkhZJOm0ry/aR9KCktyU9I+kiSTNSjudkSX+RtFTSuZutGyXp95LeSvZztaQdknXTk2LPJedhvKRekh6S9IakN5PP/dt7zs3KOVFY4UXEycCrwFER0SMiLi1bfRgwBPjH5PvDwD7A7sCzwJSUTX8Q6An0A04Hfi6p11aU/TmwKilzSvJqlaT9gGuAk4E9gT5A+R/2DcB3gb7Ax4HDgf+dnIdDkzLDkvNwJ6V/wzcBAykl09XA1SnHbFYxJwqrdxMjYlVErAaIiBsjYkVErAUmAsMk9Wyj7jrggohYFxG/AlYC+1ZSVlIX4IvAeRHxTkTMByanxDsOeCgipicx/huwsXllRMyMiKciYn1ELAB+QSkZtioilkbEvcm+VwAXp5U32xq+rmv17rXmD8kf7YuB44HdeO8PcF9geSt1l0bE+rLv7wA92thPW2V3o/Tv6LWydeWfN7dn+fqIWCVpadkx/C/gSqAR2DnZ9sy2NiZpZ+Aq4DNAcwunQVKXiNiQEodZZm5RWL1oa3he+fIvA8cAR1C6TDQoWa78wuINYD2bXj4akFJ+Ufn65A99n7L11wAvUBrZtAvwI9LjP5tSK2h0Ur758lSex2ydjBOF1YvFwF5bKNMArAWWUvq/8R/nHVTyf+33ARMl7SzpI8BXU6rcAxwp6ZCkk/oCNv132AC8DaxMtvXNzepvfh4aKPVLvCWpN3Beuw7IrBVOFFYvfgL8azIa6Jw2ytwC/AV4HZgPPFWl2M6k1IL5G3ArcDulhPU+ETEP+BZwG6XWxZvAwrIi51BqGa0ArgPu3GwTE4HJyXk4AfgpsBPwd0rH++sOOSKzMr7hzqyDSfoP4IMR0eboJ7N64haFWTtJ+oikoSoZRWn47P21jsuso3jUk1n7NVC63LQnsAS4AnigphGZdSBfejIzs1S+9GRmZqnq4tJT3759Y9CgQbUOw8ysrsycOfPvEbFbe7dTF4li0KBBNDU11ToMM7O6IukvHbEdX3oyM7NUThRmZpbKicLMzFLVRR+FmRnAunXrWLhwIWvWrKl1KIXSrVs3+vfvT9euXXPZvhOFmdWNhQsX0tDQwKBBg5A8QS5ARLB06VIWLlzI4MGDc9mHLz2ZWd1Ys2YNffr0cZIoI4k+ffrk2spyojCzuuIk8X55nxMnCjMzS+U+CjOrW7c9/WqHbu/Loz/UodvLauzYsVx++eU0NjbSo0cPVq5cWZM42lIXiWLZqncr+kHU6j+2mdm2qC4ShZlZkVx44YVMmTKFAQMG0LdvXw488EB69uzJpEmTePfdd/nwhz/Mrbfeys4778zdd9/N+eefT5cuXejZsyfTp09n9erVnHbaacyfP58hQ4awevXqTbZ/9tln8/jjj9OrVy/uuOMOdtut3dM1tYv7KMzMKtDU1MS9997LH//4R+67776Weei+8IUv8Mwzz/Dcc88xZMgQbrjhBgAuuOACHnnkEZ577jmmTp0KwDXXXMPOO+/M7NmzOffcc5k5c2bL9letWsXIkSN59tlnOeywwzj//POrf5CbcaIwM6vAjBkzOOaYY9hpp51oaGjgqKOOAmDu3Ll84hOf4IADDmDKlCnMmzcPgIMPPphTTz2V6667jg0bNgAwffp0TjrpJACGDh3K0KFDW7a/3XbbMX78eABOOukkZsyYUc3Da5UThZlZBdp62Nupp57K1VdfzZw5czjvvPNa7mu49tprueiii3jttdcYPnw4S5cuBbIPaS3CcGAnCjOzChxyyCE8+OCDrFmzhpUrV/LLX/4SgBUrVrDHHnuwbt06pkyZ0lL+5ZdfZvTo0VxwwQX07duX1157jUMPPbSlzNy5c5k9e3ZL+Y0bN3LPPfcAcNttt3HIIYdU8eha585sM6tbtRjh+LGPfYyjjz6aYcOGMXDgQBobG+nZsycXXngho0ePZuDAgRxwwAGsWLECgO9973u8+OKLRASHH344w4YNY9999+W0005j6NChDB8+nFGjRrVsv3v37sybN6+lg/zOO++s+jFuri6emb3XkKFx0c0PZS7v4bFm26bnn3+eIUOG1DoMVq5cSY8ePXjnnXc49NBDmTRpEiNHjqxpTK2dG0kzI6Kxvdt2i8LMrEITJkxg/vz5rFmzhlNOOaXmSSJvThRmZhW67bbbah1CVbkz28zMUjlRmJlZKicKMzNL5URhZmap3JltZvWr6aaO3V7jaR26uWnTprHDDjswZswYoHT39pFHHsm4ceMyb2PQoEE0NTXRt2/fDo2tEm5RmJnlZNq0aTz55JO1DqPdnCjMzCp0yy23MHToUIYNG8bJJ5/Mgw8+yOjRoxkxYgRHHHEEixcvZsGCBVx77bVcddVVDB8+nN/97ndAaULAMWPGsNdee7VM1TFt2jSOPPLIlu2feeaZ3HzzzS3fL7vsMkaNGsWoUaN46aWXqnqskOOlJ0ndgOnAjsl+7omI8yT1Bu4EBgELgBMi4s284jAz60jz5s3j4osv5oknnqBv374sW7YMSTz11FNI4vrrr+fSSy/liiuu4IwzzqBHjx6cc845ANxwww0sWrSIGTNm8MILL3D00Udnugy1yy678Ic//IFbbrmFs846i4ceyj5TRUfIs49iLfCpiFgpqSswQ9LDwBeAxyLiEkk/AH4AfD/HOMzMOsxvf/tbxo0b19Jn0Lt3b+bMmcP48eNZtGgR7777LoMHD26z/rHHHst2223Hfvvtx+LFizPt88QTT2x5/+53v9v+g6hQbpeeoqT5wa9dk1cAxwCTk+WTgWPzisHMrKNFxPum/v72t7/NmWeeyZw5c/jFL37RMsV4a3bcccdNtgWw/fbbs3Hjxpblm9cv318tph3PtY9CUhdJs4AlwKMR8TTwgYhYBJC8795G3QmSmiQ1rXhrWZ5hmplldvjhh3PXXXe1PFdi2bJlLF++nH79+gEwefLklrINDQ0ts8imGThwIPPnz2ft2rUsX76cxx57bJP1zTPI3nnnnXz84x/vqEPJLNfhsRGxARguaVfgfkkfraDuJGASlGaPzSlEM6tnHTycNYv999+fc889l8MOO4wuXbowYsQIJk6cyPHHH0+/fv046KCDeOWVVwA46qijGDduHA888AA/+9nP2tzmgAEDOOGEExg6dCj77LMPI0aM2GT92rVrGT16NBs3buT222/P9fhaU7VpxiWdB6wCvg6MjYhFkvYApkXEvml1Pc24mUFxphkvojynGc/t0pOk3ZKWBJJ2Ao4AXgCmAqckxU4BHsgrBjMza788Lz3tAUyW1IVSQrorIh6S9HvgLkmnA68Cx+cYg5mZtVNuiSIiZgMjWlm+FDg8r/2a2battVFHnV3eXQi+M9vM6ka3bt1YunRp7n8Y60lEsHTpUrp165bbPjwpoJnVjf79+7Nw4ULeeOONWodSKN26daN///65bd+JwszqRteuXVPverZ8+NKTmZmlcqIwM7NUThRmZpbKicLMzFI5UZiZWSonCjMzS+VEYWZmqZwozMwslROFmZmlcqIwM7NUThRmZpbKicLMzFI5UZiZWSonCjMzS+VEYWZmqZwozMwslROFmZmlcqIwM7NUThRmZpYqt0QhaYCkxyU9L2mepO8kyydKel3SrOT1ubxiMDOz9ts+x22vB86OiGclNQAzJT2arLsqIi7Pcd9mZtZBcksUEbEIWJR8XiHpeaBfXvszM7N8VKWPQtIgYATwdLLoTEmzJd0oqVcbdSZIapLUtOKtZdUI08zMWpF7opDUA7gXOCsi3gauAfYGhlNqcVzRWr2ImBQRjRHR2LBr77zDNDOzNuSaKCR1pZQkpkTEfQARsTgiNkTERuA6YFSeMZiZWfvkOepJwA3A8xFxZdnyPcqKHQfMzSsGMzNrvzxHPR0MnAzMkTQrWfYj4ERJw4EAFgDfyDEGMzNrpy0mCkm3RsTJW1q2uYiYAaiVVb+qLEQzM6ulLJee9i//IqkLcGA+4ZiZWdG0mSgk/VDSCmCopLeT1wpgCfBA1SI0M7OaajNRRMRPIqIBuCwidkleDRHRJyJ+WMUYzcyshrbYRxERP5TUDxhYXj4ipucZmJmZFUOWzuxLgC8B84ENyeIAnCjMzDqBLMNjjwP2jYi1eQdjZmbFk2XU05+BrnkHYmZmxZSlRfEOMEvSY0BLqyIi/k9uUZmZWWFkSRRTk5eZmXVCWUY9Ta5GIGZmVkxZRj29QmmU0yYiYq9cIjIzs0LJcumpsexzN+B4wA+IMDPrJLY46ikilpa9Xo+InwKfqkJsZmZWAFkuPY0s+7odpRZGQ24RmZlZoWS59FT+qNL1lJ4hcUIu0ZiZWeFkGfX0yWoEYmZmxbTFPgpJPSVdKakpeV0hqWc1gjMzs9rLMoXHjcAKSpebTgDeBm7KMygzMyuOLH0Ue0fEF8u+n1/2DGwzM9vGZWlRrJZ0SPMXSQcDq/MLyczMiiRLi+KbwOSyfok3gVNzi8jMzAoly6inWcAwSbsk39/OPSozMyuMLKOefixp14h4OyLeltRL0kUZ6g2Q9Lik5yXNk/SdZHlvSY9KejF579URB2JmZvnI0kfx2Yh4q/lLRLwJfC5DvfXA2RExBDgI+Jak/YAfAI9FxD7AY8l3MzMrqCyJooukHZu/SNoJ2DGlPAARsSgink0+rwCeB/oBxwDNU5dPBo6tNGgzM6ueLJ3Z/xd4TNJNlKYb/xrv/aHPRNIgYATwNPCBiFgEpWQiafc26kwAJgD0/WC/SnZnZmYdKEtn9qWSZgNHAAIujIhHsu5AUg/gXuCspI8jU72ImARMAthryND3PQ/DzMyqI0uLgoj4NfDrSjcuqSulJDElIu5LFi+WtEfSmtgDWFLpds3MrHqy9FFsFZWaDjcAz0fElWWrpgKnJJ9PAR7IKwYzM2u/TC2KrXQwcDIwp2zKjx8BlwB3STodeJXSE/PMzKygcksUETGDUp9Gaw7Pa79mZtax2kwUkuZQGuXUqogYmktEZmZWKGktiiOT928l77cm718B3sktIjMzK5Q2E0VE/AVKs8VGxMFlq34g6QnggryDMzOz2ssy6qn7ZtOMjwG65xeSmZkVSZbO7NOBG5NpxgNYTunubDMz6wSy3Jk9k/emGVdELM8/LDMzK4os04x/QNINwJ0RsVzSfsk9EGZm1glk6aO4GXgE2DP5/ifgrLwCMjOzYsmSKPpGxF3ARoCIWA9syDUqMzMrjCyJYpWkPiQ330k6iFKHtpmZdQJZRj2dTWkiv72T+yd2w/MzmZl1GplGPUk6DNiX0txN/xMR63KPzMzMCiHLqKeXgX+KiHkRMTci1kl6qAqxmZlZAWS59LQO+KSk0cA3IuJdSs++Lqzbnn41c9kvj/5QjpGYmdW/LJ3Z70TEeOB54HeSBpIyq6yZmW1bsrQoBC3Pzp5J6Z6K3rlGZWZmhZElUfx784eIeEzSP/Leo0zNzGwbl/bgoo9ExAvA65JGbrbandlmZp1EWovibODrwBWtrAvgU7lEZGZmhZL24KKvJ++frF44ZmZWNGmXnr6QVjEi7uv4cMzMrGjSLj0dlbIuACcKM7NOIO3S02nt2bCkG4EjgSUR8dFk2URK/R5vJMV+FBG/as9+zMwsX1mGxyLp88D+QLfmZRFxwRaq3QxcDdyy2fKrIuLyCmI0M7MayjLX07XAeODblG6+Ox4YuKV6ETEdWNbeAM3MrLayTOExJiK+CrwZEecDHwcGtGOfZ0qaLelGSb3aKiRpgqQmSU0r3nK+MTOrlSyJYnXy/o6kPSlNEjh4K/d3DbA3MBxYROv3aAAQEZMiojEiGht29YwhZma1kqWP4iFJuwKXAc9SGvF0/dbsLCIWN3+WdB2+w9vMrPCyPLjowuTjvclzKLpFxFY9ClXSHhGxKPl6HDB3a7ZjZmbVs8VEIakL8HlgUHN5SUTElVuodzswFugraSFwHjBW0nBKrZIFwDfaEbuZmVVBlktPDwJrgDnAxqwbjogTW1l8Q9b6ZmZWDFkSRf+IGJp7JCl2fPdN9n717szlX/7Q8TlGY2bWuWQZ9fSwpH/IPRIzMyukLC2Kp4D7JW1HaWisgIiIXXKNzMzMCiFLoriC0k12cyLCz8o2M+tkslx6ehGY6yRhZtY5ZWlRLAKmSXoYWNu8cEvDY83MbNuQJVG8krx2SF5mZtaJpCaK5Ga7HhHxvSrFY2ZmBZOaKCJig6SR1Qqmo1RyzwVdekNju57RZGa2Tcty6WmWpKnA3cCq5oV+ZraZWeeQJVH0BpYCnypb5mdmm5l1Ellmj/V1GTOzTizLo1D7S7pf0hJJiyXdK6l/NYIzM7Pay3LD3U3AVGBPoB+l2WRvyjMoMzMrjiyJYreIuCki1ievm4Hdco7LzMwKIkui+LukkyR1SV4nUercNjOzTiBLovgacALwN0rTeYxLlpmZWSeQZdTTq8DRVYjFzMwKqM1EIenfU+pFRFyYQzxmZlYwaS2KVa0s6w6cDvQBnCjMzDqBNhNFRFzR/FlSA/Ad4DTgDkoPMzIzs05gS7PH9gb+GfgKMBkYGRFvViMwMzMrhjZHPUm6DHgGWAEcEBETK0kSkm5M7uaeW7ast6RHJb2YvPdqV/RmZpY7tfWEU0kbKT3Rbj2lSQBbVlHqzN4ldcPSocBK4JaI+Giy7FJgWURcIukHQK+I+P6Wghyy94C4+ZKzshxPxUYP7l15JU9LbmZ1QNLMiGhs73bS+iiy3GPRpoiYLmnQZouPAcYmnycD04AtJgozM6uddiWDrfCBiFgEkLzvXuX9m5lZhaqdKDKTNEFSk6Smt95ubaSumZlVQ7UTxWJJewAk70vaKhgRkyKiMSIad92le9UCNDOzTVU7UUwFTkk+nwI8UOX9m5lZhXJLFJJuB34P7CtpoaTTgUuAT0t6Efh08t3MzAosyzOzt0pEnNjGqsPz2qeZmXW8wnZmm5lZMThRmJlZKicKMzNL5URhZmapnCjMzCxVbqOe6sXTryyrqPzowb2h6abKduJJBM2sjrlFYWZmqZwozMwslROFmZmlcqIwM7NUThRmZpbKicLMzFI5UZiZWSonCjMzS+VEYWZmqZwozMwslROFmZmlcqIwM7NUThRmZpbKicLMzFJ1+mnGq8LTkptZHXOLwszMUjlRmJlZqppcepK0AFgBbADWR0RjLeIwM7Mtq2UfxScj4u813L+ZmWXgS09mZpaqVokigN9ImilpQmsFJE2Q1CSp6a23V1U5PDMza1arS08HR8RfJe0OPCrphYiYXl4gIiYBkwCG7D0gahGkmZnVqEUREX9N3pcA9wOjahGHmZltWdVbFJK6A9tFxIrk8z8AF1Q7jq319CvLKio/enDvnCIxM6uOWlx6+gBwv6Tm/d8WEb+uQRxmZpZB1RNFRPwZGFbt/dYVT/lhZgXi4bFmZpbKicLMzFI5UZiZWSonCjMzS+VEYWZmqZwozMwslZ9wty3wcFozy5FbFGZmlsotipx5yg8zq3duUZiZWSonCjMzS+VEYWZmqZwozMwslTuzOyMPpzWzCrhFYWZmqdyiKJhCDqd1C8SsU3OLwszMUjlRmJlZKl96qnOVXKqq2l3flV6qqpQvbZlVlVsUZmaWyi0Kqz/uXDerKieKTqSQI6qqwYnFrF1qkigkfQb4T6ALcH1EXFKLOCydE0tBOHFZjVU9UUjqAvwc+DSwEHhG0tSImF/tWKxjVZpYKrHNJKGtUe+Jyy26uleLFsUo4KWI+DOApDuAYwAnCmtTnkmo3lU9ieaduIqWGK0miaIf8FrZ94XA6M0LSZoATEi+rj3ohHPmViG29uoL/L3WQWTgODtOPcQIjrOj1Uuc+3bERmqRKNTKsnjfgohJwCQASU0R0Zh3YO3lODtWPcRZDzGC4+xo9RRnR2ynFvdRLAQGlH3vD/y1BnGYmVkGtUgUzwD7SBosaQfgS8DUGsRhZmYZVP3SU0Ssl3Qm8Ail4bE3RsS8LVSblH9kHcJxdqx6iLMeYgTH2dE6VZyKeF/3gJmZWQvP9WRmZqmcKMzMLFVNE4Wkz0j6H0kvSfpBK+sl6b+S9bMljcxat0BxLpA0R9Ksjhqq1o44PyLp95LWSjqnkroFirNI5/MryX/v2ZKelDQsa90CxVmV85khxmOS+GZJapJ0SNa6BYqzML/NsnIfk7RB0rhK624iImryotSR/TKwF7AD8Byw32ZlPgc8TOnei4OAp7PWLUKcyboFQN+CnM/dgY8BFwPnVFK3CHEW8HyOAXolnz9b4N9nq3FW63xmjLEH7/WZDgVeKOi5bDXOov02y8r9FvgVMK4957OWLYqWqTwi4l2geSqPcscAt0TJU8CukvbIWLcIcVbTFuOMiCUR8QywrtK6BYmzmrLE+WREvJl8fYrSPUGZ6hYkzmrJEuPKSP6SAd157ybcop3LtuKspqzn5NvAvcCSrai7iVomitam8uiXsUyWuh2lPXFC6Yf0G0kzVZqWJC/tOSdFO59pino+T6fUqtyauu3RnjihOuczU4ySjpP0AvBL4GuV1C1AnFCg36akfsBxwLWV1m1NLZ9HkWUqj7bKZJoGpIO0J06AgyPir5J2Bx6V9EJETO/QCLccQ551K9XefRXufEr6JKU/wM3Xqwt5PluJE6pzPrNO23M/cL+kQ4ELgSOy1u0g7YkTivXb/Cnw/YjYIG1SfKvOZy1bFFmm8mirTDWnAWlPnERE8/sS4H5KTb9axZlH3Uq1a19FO5+ShgLXA8dExNJK6hYgzmqdz4rOR/LHdW9JfSut207tibNov81G4A5JC4BxwH9LOjZj3ffLu+MlpUNme+DPwGDe61TZf7Myn2fTTuI/ZK1bkDi7Aw1ln58EPlOrOMvKTmTTzuxCnc+UOAt1PoEPAS8BY7b2GGscZ1XOZ8YYP8x7ncQjgdeTf09FO5dtxVmo3+Zm5W/mvc7srTqfHX4QFR7w54A/UeqFPzdZdgZwRvJZlB5y9DIwB2hMq1u0OCmNLHguec0rQJwfpPR/FG8DbyWfdyng+Ww1zgKez+uBN4FZyaupoL/PVuOs5vnMEOP3kxhmAb8HDinouWw1zqL9NjcrezNJotja8+kpPMzMLJXvzDYzs1ROFGZmlsqJwszMUjlRmJlZKicKMzNL5URhdUXSylaWnSHpqx24j2mSGltZ3ijpv5LPYyWN6ah9trKvsZIeymv7ZpWo5RQeZh0iIjafzyav/TQBzdNHjwVWUrqxymyb5haF1T1JE5ufW5G0Bn6aPHdhrqRRyfLDkucEzJL0R0kNyfJ/SZ4h8JykS8o2e7ykP0j6k6RPJGXHSnpI0iBKNzd9N9neJ8q2PUvSakmHbRbj05L2L/s+TdKBkkYlsf4xed837fiS73OTGJB0UhLnLEm/kNSlQ06qWRm3KGxb1D0ixiSTtt0IfBQ4B/hWRDwhqQewRtJngWOB0RHxjqTeZdvYPiJGSfoccB7vTfxGRCyQdC2wMiIuTxYPB5B0FPAvvL+lcQdwAnBeMgX9nhExU9IuwKERsV7SEcCPgS9mOUhJQ4DxlCajWyfpv4GvALdkPE9mmThR2LbodihN2iZpF0m7Ak8AV0qaAtwXEQuTP8w3RcQ7SfllZdu4L3mfCQzKslNJ+wCXAZ+KiM2fpXEX8CilpHMCcHeyvCcwOakbQNcKjvNw4EDgmWSG0J3Y9NkDZh3CicK2RZvPSxMRcYmkX1Ka5+apJEmolbLN1ibvG8jw70RSd0rJ4OuRzCK6WQCvS1qazOI6HvhGsupC4PGIOC65nDStlc2vZ9PLxN2adwtMjogfbik+s/ZwH4Vti8YDqPQ84+URsVzS3hExJyL+g1KH9EeA3wBfk7RzUr53m1t8vxVAQ9n3myi1Tn6XUucOSpelekbEnGRZT0ozkAKc2ka9BZRmKkWl57EPTpY/BoxLnn+ApN6SBlZwDGaZOFFYvdlZ0sKy1z+3UuZNSU9SerrX6cmys5JO4OeA1cDDEfFrYCrQJGkWpX6MrB4Ejks6kQ+jNOf/18o6tN83vBa4B/gSpZZHs0uBn0h6gtLzjFtzL9A7ifGblGb+JCLmA/9K6alqsyld2qr2I3itE/DssbZNkTSN0jMsmrZU1syycYvCzMxSuUVhZmap3KIwM7NUThRmZpbKicLMzFI5UZiZWSonCjMzS/X/Af5y3dFs/pkNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "# Plot the distributions\n",
    "ax = sns.distplot(lips_train['gasdb'], kde=False,\n",
    "                  bins=1000,\n",
    "                  norm_hist=True,\n",
    "                  label='gasdb',\n",
    "                 )\n",
    "ax = sns.distplot(lips_train['cathub'], kde=False,\n",
    "                  bins=100,\n",
    "                  norm_hist=True,\n",
    "                  label='cathub',\n",
    "                 )\n",
    "\n",
    "# Format\n",
    "plot_lims = (0, 0.4)\n",
    "_ = ax.set_xlim(plot_lims)\n",
    "_ = ax.set_xlabel('Lipschitz value')\n",
    "_ = ax.set_ylabel('Normalized count')\n",
    "_ = ax.set_title('training data')\n",
    "_ = ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgdVZ3/8feHsASyEgJMTGICyIMByWabIGsEVGTHCSAKQuRnxJ8w4oAOyIysKiMgjuKAYQ1I2AQeAqjIE8kPw96BkA1GQCIEMgkmBBIgYcn390dVNzdNd3X1UvdWpz+v57nPvXWqzqnvLZr7zalTdUoRgZmZWUs2qnUAZmZWbk4UZmaWyYnCzMwyOVGYmVkmJwozM8vkRGFmZpmcKGyDIGmCpMUVywskTcizbTv2dYWk/2hv/Tbsp0NxmnUWJwrbIEXELhExs6PtSDpB0qwmbZ8UEed3tO3O1FycZp3FicLMzDI5UVhpSDpD0u+alP2XpF+mnydJekbSKkl/k/StjLYWSdo//by5pOskvS5pIfCZZvb7QtruQklHpOUjgCuAz0paLWllWn6dpAsq6n9T0vOSVkiaLuljFetC0kmSnkv3/2tJaiHmzo7zIElPSXpT0suSzsn8D2DWAicKK5ObgAMl9QWQ1AM4CpiWrl8GHAz0BSYBl0oam6Pds4Ed0tcXgeObrH8B2AvoB5wL/FbSoIh4BjgJeCQiekdE/6YNS9oX+Gka5yDg78DNTTY7mORHf1S63RerFOdbwNeB/sBBwLclHd7Cvs1a5ERhpRERfweeBBp+zPYF3o6IR9P190bEC5H4f8CfSH44W3MU8OOIWBERLwO/bLLf2yLi1YhYFxG3AM8B43KG/TXgmoh4MiLWAmeS/Mt+eMU2F0bEyoh4CXgAGF2NOCNiZkTMS7efS5KI98n5vcwaOVFY2UwDjkk/f5UPexNI+pKkR9NTPCuBA4GBOdr8GPByxfLfK1dK+rqkOZJWpu1+Kme7DW03thcRq4HlwOCKbf634vPbQO9qxClpvKQHJL0m6Q2SXkfe72XWyInCyuY2YIKkIcARpIlC0mbA7cDFwLbp6ZXfA82e729iCTC0YvnjDR8kDQOuBE4GtkrbnV/RbmvTK78KDKtorxewFfBKjriKjnMaMB0YGhH9SMYx8hwvs/U4UVipRMRrwEzgWuDF9Pw7wKbAZsBrwPuSvgR8IWeztwJnStoyTUCnVKzrRfIj+xokA+Yk/1JvsBQYImnTFtqeBkySNDpNZj8BHouIRTljKzLOPsCKiFgjaRxJD82szZworIymAftTcdopIlYB/0LyY/o6yY/e9JztnUtyGudFknGNGyraXQhcAjxC8mO7K/BQRd0/AwuA/5X0j6YNR8QM4D9IejtLSAaiv5IzrqLj/L/AeZJWAT8iOXZmbSY/uMjMzLK4R2FmZpmcKMzMLJMThZmZZXKiMDOzTBvXOoA8Bg4cGMOHD691GGZmXcrs2bP/ERFbd7SdLpEohg8fTn19fa3DMDPrUiT9vfWtWudTT2ZmlsmJwszMMjlRmJlZpi4xRmFmBvDee++xePFi1qxZU+tQSqVnz54MGTKETTbZpJD2nSjMrMtYvHgxffr0Yfjw4bTwoMBuJyJYvnw5ixcvZrvttitkHz71ZGZdxpo1a9hqq62cJCpIYquttiq0l+VEYWZdipPERxV9TJwozMwsk8cozKzLmvbYS53a3lfHf7z1jQowYcIELr74Yurq6ujduzerV6+uSRwt6RKJYsVb72b+QdTqP66ZWXdQeKKQ1AOoB16JiIMlDQBuAYYDi4CjIuL1ouMwM+ss559/PjfeeCNDhw5l4MCBfPrTn6Zfv35MmTKFd999l0984hPccMMNbLHFFtx2222ce+659OjRg379+vHggw/yzjvvMGnSJBYuXMiIESN455131mv/tNNO44EHHmDLLbfk5ptvZuutOzxdU4dUY4ziu8AzFctnADMiYkdgRrpsZtYl1NfXc/vtt/PUU09xxx13NM5D9+Uvf5knnniCp59+mhEjRnD11VcDcN5553Hffffx9NNPM3168vTeyy+/nC222IK5c+dy1llnMXv27Mb233rrLcaOHcuTTz7JPvvsw7nnnlv9L9lEoYkifUD8QcBVFcWHAVPTz1OBw4uMwcysM82aNYvDDjuMzTffnD59+nDIIYcAMH/+fPbaay923XVXbrzxRhYsWADAHnvswQknnMCVV17JBx98AMCDDz7IscceC8DIkSMZOXJkY/sbbbQRRx99NADHHnsss2bNqubXa1bRPYpfAD8A1lWUbRsRSwDS922aqyhpsqR6SfWrVq4oOEwzs3wiotnyE044gcsuu4x58+Zx9tlnN97XcMUVV3DBBRfw8ssvM3r0aJYvXw7kv6S1DJcDF5YoJB0MLIuI2a1u3IyImBIRdRFR16f/gE6Ozsysffbcc0/uvvtu1qxZw+rVq7n33nsBWLVqFYMGDeK9997jxhtvbNz+hRdeYPz48Zx33nkMHDiQl19+mb333rtxm/nz5zN37tzG7detW8fvfvc7AKZNm8aee+5ZxW/XvCIHs/cADpV0INAT6Cvpt8BSSYMiYomkQcCyAmMwsw1YLa54/MxnPsOhhx7KqFGjGDZsGHV1dfTr14/zzz+f8ePHM2zYMHbddVdWrVoFwPe//32ee+45IoL99tuPUaNGsdNOOzFp0iRGjhzJ6NGjGTduXGP7vXr1YsGCBY0D5LfcckvVv2NTaqkb1ak7kSYAp6dXPV0ELI+ICyWdAQyIiB9k1d9+xMi44Lp7Wlzvy2PNuodnnnmGESNG1DoMVq9eTe/evXn77bfZe++9mTJlCmPHjq1pTM0dG0mzI6Kuo23X4j6KC4FbJZ0IvAQcWYMYzMzabfLkySxcuJA1a9Zw/PHH1zxJFK0qiSIiZgIz08/Lgf2qsV8zsyJMmzat1iFUled6MjOzTE4UZmaWyYnCzMwyOVGYmVmmLjF7rJlZs+qv7dz26iZ1anMzZ85k0003ZffddweSu7cPPvhgJk6cmLuN4cOHU19fz8CBAzs1trZwj8LMrCAzZ87k4YcfrnUYHeZEYWbWRtdffz0jR45k1KhRHHfccdx9992MHz+eMWPGsP/++7N06VIWLVrEFVdcwaWXXsro0aP5y1/+AiQTAu6+++5sv/32jVN1zJw5k4MPPrix/ZNPPpnrrruucfmiiy5i3LhxjBs3jueff76q3xV86snMrE0WLFjAj3/8Yx566CEGDhzIihUrkMSjjz6KJK666ip+9rOfcckll3DSSSfRu3dvTj/9dACuvvpqlixZwqxZs3j22Wc59NBDc52G6tu3L48//jjXX389p556Kvfc0/JMFUVwojAza4M///nPTJw4sXHMYMCAAcybN4+jjz6aJUuW8O6777Lddtu1WP/www9no402Yuedd2bp0qW59nnMMcc0vn/ve9/r+JdoI596MjNrg4j4yNTfp5xyCieffDLz5s3jN7/5TeMU483ZbLPN1msLYOONN2bdug+fxtC0fuX+ajHtuBOFmVkb7Lffftx6662Nz5VYsWIFb7zxBoMHDwZg6tSpjdv26dOncRbZLMOGDWPhwoWsXbuWN954gxkzZqy3vmEG2VtuuYXPfvaznfVVcvOpJzPrujr5ctY8dtllF8466yz22WcfevTowZgxYzjnnHM48sgjGTx4MLvtthsvvvgiAIcccggTJ07krrvu4le/+lWLbQ4dOpSjjjqKkSNHsuOOOzJmzJj11q9du5bx48ezbt06brrppkK/X3OqMs14R3macTOD8kwzXkZFTjPuU09mZpbJicLMzDI5UZhZl9IVTpdXW9HHpLBEIamnpMclPS1pgaRz0/JzJL0iaU76OrCoGMxsw9KzZ0+WL1/uZFEhIli+fDk9e/YsbB9FXvW0Ftg3IlZL2gSYJekP6bpLI+LiAvdtZhugIUOGsHjxYl577bVah1IqPXv2ZMiQIYW1X1iiiCTlr04XN0lf/meAmbXbJptsknnXsxWj0DEKST0kzQGWAfdHxGPpqpMlzZV0jaQtW6g7WVK9pPpVK1cUGaaZmWUoNFFExAcRMRoYAoyT9CngcmAHYDSwBLikhbpTIqIuIur69B9QZJhmZpahKlc9RcRKYCZwQEQsTRPIOuBKYFw1YjAzs/Yp8qqnrSX1Tz9vDuwPPCtpUMVmRwDzi4rBzMw6rsirngYBUyX1IElIt0bEPZJukDSaZGB7EfCtAmMwM7MOKvKqp7nAmGbKjytqn2Zm1vl8Z7aZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCxTkY9C7SnpcUlPS1og6dy0fICk+yU9l75vWVQMZmbWcUX2KNYC+0bEKGA0cICk3YAzgBkRsSMwI102M7OSKixRRGJ1urhJ+grgMGBqWj4VOLyoGMzMrOMKHaOQ1EPSHGAZcH9EPAZsGxFLANL3bVqoO1lSvaT6VStXFBmmmZllKDRRRMQHETEaGAKMk/SpNtSdEhF1EVHXp/+A4oI0M7NMrSYKSTfkKcsSESuBmcABwFJJg9J2BpH0NszMrKTy9Ch2qVyQ1AP4dGuVJG0tqX/6eXNgf+BZYDpwfLrZ8cBdbQnYzMyqa+OWVkg6E/ghsLmkNxuKgXeBKTnaHgRMTRPLRsCtEXGPpEeAWyWdCLwEHNmRL2BmZsVqMVFExE+Bn0r6aUSc2daGI2IuMKaZ8uXAfm1tz8zMaqPFRNEgIs6UNBgYVrl9RDxYZGBmZlYOrSYKSRcCXwEWAh+kxQE4UZiZdQOtJgrgCGCniFhbdDBmZlY+ea56+hvJXdVmZtYN5elRvA3MkTSDZP4mACLiXwqLyszMSiNPopievszMrBvKc9XT1Na2MTOzDVeeq55eJLnKaT0RsX0hEZmZWankOfVUV/G5J8md1J6lz8ysm2j1qqeIWF7xeiUifgHsW4XYzMysBPKcehpbsbgRSQ+jT2ERmZlZqeQ59XRJxef3gUXAUYVEY2ZmpZPnqqfPVSMQMzMrpzwPLuon6ecNjyWVdImkftUIzszMai/PFB7XAKtITjcdBbwJXFtkUGZmVh55xih2iIh/rlg+V9KcogIyM7NyydOjeEfSng0LkvYA3ikuJDMzK5M8ieLbwK8lLZK0CLgMOKm1SpKGSnpA0jOSFkj6blp+jqRXJM1JXwd26BuYmVmh8lz1NAcYJalvuvxmK1UavA+cFhFPSuoDzJZ0f7ru0oi4uF0Rm5lZVeW56uknkvpHxJsR8aakLSVd0Fq9iFgSEU+mn1cBzwCDOx6ymZlVU55TT1+KiJUNCxHxOtCm00WShgNjgMfSopMlzZV0jaQtW6gzueGS3FUrV7Rld2Zm1onyJIoekjZrWJC0ObBZxvbrkdQbuB04NT1tdTmwAzAaWML6d343iogpEVEXEXV9+nsOQjOzWslzeexvgRmSriWZbvwbQK5nVEjahCRJ3BgRdwBExNKK9VcC97Q16KamPfZS5vqvjv94R3dhZtZt5RnM/pmkucD+gIDzI+K+1upJEnA18ExE/LyifFBELEkXjwDmtytyMzOrijw9CiLij8Af29j2HsBxwLyKG/R+CBwjaTRJ72QR8K02tmtmZlWUK1G0R0TMIumBNPX7ovZpZmadL89gtpmZdWNOFGZmlqnFU0+S5pGMIzQrIkYWEpGZmZVK1hjFwen7d9L3G9L3rwFvFxaRmZmVSouJIiL+DslssRGxR8WqMyQ9BJxXdHBmZlZ7ecYoejWZZnx3oFdxIZmZWZnkuTz2ROCa9PGnAbxBcne2mZl1A3nuzJ7Nh9OMKyLeKD4sMzMrizzTjG8r6Wrgloh4Q9LOkk6sQmxmZlYCecYorgPuAz6WLv8VOLWogMzMrFzyJIqBEXErsA4gIt4HPig0KjMzK408ieItSVuR3nwnaTeSAW0zM+sG8lz1dBowHdghvX9ia+DIQqMyM7PSyHXVk6R9gJ1IZoP9n4h4r/DIzMysFPJc9fQC8H8iYkFEzI+I9yR1+Kl0ZmbWNeQZo3gP+JykayVtmpYNLjAmMzMrkTyJ4u2IOBp4BviLpGFkzCrbQNJQSQ9IekbSAknfTcsHSLpf0nPp+5Yd+wpmZlakPIlCkDw7m+RRpvcBQ3LUex84LSJGALsB35G0M3AGMCMidgRmpMtmZlZSeRLFjxo+RMQM4IvAZa1VioglEfFk+nkVSY9kMHAYMDXdbCpweBtjNjOzKsp6cNEnI+JZ4BVJY5usbtNgtqThwBjgMWDbiFgCSTKRtE0LdSYDkwEG/pOHRMzMaiXr8tjTgG8ClzSzLoB98+xAUm/gduDUiHhTUq7AImIKMAVg+xEjWx0TMTOzYmQ9uOib6fvn2tu4pE1IksSNEXFHWrxU0qC0NzEIWNbe9s3MrHhZp56+nFWx4oe/pfoCrgaeiYifV6yaDhwPXJi+35U7WjMzq7qsU0+HZKwLIDNRAHsAxwHzJM1Jy35IkiBuTacqfwlPB2JmVmpZp54mdaThiJhFemltM/brSNtmZlY9eSYFRNJBwC5Az4ayiDivqKDMzKw88sz1dAVwNHAKSQ/hSGBYwXGZmVlJ5LnhbveI+DrwekScC3wWGFpsWGZmVhZ5EsU76fvbkj5GMkngdsWFZGZmZZJnjOIeSf2Bi4AnSa54uqrQqMzMrDTyPLjo/PTj7elzKHpGhB+FambWTbSaKCT1AA4ChjdsL4kmN9GZmdkGKs+pp7uBNcA8YF2x4ZiZWdnkSRRDImJk4ZGYmVkp5bnq6Q+SvlB4JGZmVkp5ehSPAndK2ojk0lgBERF9C43MzMxKIU+iuITkJrt5EeHnQpiZdTN5Tj09B8x3kjAz657y9CiWADMl/QFY21Doy2PNzLqHPInixfS1afoyM7NuJDNRpDfb9Y6I71cpHjMzK5nMMYqI+AAY256GJV0jaZmk+RVl50h6RdKc9HVge9o2M7PqyXPqaY6k6cBtwFsNha09Mxu4DrgMuL5J+aURcXFbgjQzs9rJkygGAMuBfSvKWn1mdkQ8KGl4uyMzM7NSyDN7bIeend2MkyV9HagHTouI1zu5fTMz60R5HoU6RNKd6XjDUkm3SxrSzv1dDuwAjCa57PaSjP1OllQvqX7VyhXt3J2ZmXVUnlNP1wLTSJ6VDXBsWvb5tu4sIpY2fJZ0JXBPxrZTgCkAI3YYGju8dFvu/bzw8SNb38jMzHLJc2f21hFxbUS8n76uA7Zuz84kDapYPAKY39K2ZmZWDnl6FP+QdCxwU7p8DMngdiZJNwETgIGSFgNnAxMkjSYZDF8EfKsdMZuZWRXlSRTfILnM9VKSH/iH07JMEXFMM8VXtyk6MzOruTxXPb0EHFqFWAoz7bGXMtd/dfzHqxSJmVnX02KikPSjjHoREecXEI+ZmZVMVo/irWbKegEnAlsBThRmZt1Ai4kiIhrvcZDUB/guMAm4mYz7H8zMbMPS2uyxA4B/Bb4GTAXG+k5qM7PuJWuM4iLgyyQ3ve0aEaurFpWZmZVGVo/iNJIn2v07cJakhnKRDGb3LTi2dmvLXdwAjD+tmEDMzDYAWWMUee7aNjOzDZyTgZmZZXKiMDOzTE4UZmaWyYnCzMwyOVGYmVkmJwozM8vkRGFmZpnyPI9iw1d/bdu2r5tUTBxmZiXkHoWZmWUqLFFIukbSMknzK8oGSLpf0nPp+5ZF7d/MzDpHkT2K64ADmpSdAcyIiB2BGemymZmVWGGJIiIeBFY0KT6MZLpy0vfDi9q/mZl1jmqPUWwbEUsA0vdtWtpQ0mRJ9ZLqV77Z3MP2zMysGko7mB0RUyKiLiLq+vftVetwzMy6rWoniqWSBgGk78uqvH8zM2ujaieK6cDx6efjgbuqvH8zM2ujIi+PvQl4BNhJ0mJJJwIXAp+X9Bzw+XTZzMxKrLA7syPimBZW7VfUPs3MrPOVdjDbzMzKwYnCzMwyeVLA9vAkgmbWjbhHYWZmmZwozMwsk089AY+92HRKqo8av92AKkRiZlY+7lGYmVkmJwozM8vkRGFmZpmcKMzMLJMThZmZZXKiMDOzTE4UZmaWyYnCzMwyOVGYmVkm35ldDZ5E0My6sJokCkmLgFXAB8D7EVFXizjMzKx1texRfC4i/lHD/ZuZWQ4eozAzs0y1ShQB/EnSbEmTm9tA0mRJ9ZLqV775VpXDMzOzBrU69bRHRLwqaRvgfknPRsSDlRtExBRgCsCIHYZGLYKs1NpU5J6G3Mw2VDXpUUTEq+n7MuBOYFwt4jAzs9ZVPVFI6iWpT8Nn4AvA/GrHYWZm+dTi1NO2wJ2SGvY/LSL+WIM4zMwsh6onioj4GzCq2vs1M7P28Z3ZZeQ7uc2sRHwfhZmZZXKiMDOzTE4UZmaWyWMUncQ35JnZhso9CjMzy+QexYbAV0mZWYHcozAzs0xOFGZmlsmnnqrEg91m1lW5R2FmZpnco+iOPPhtZm3gHoWZmWVyj6IkSj2G4R6IWbfmHoWZmWVyj6KLKHWPoyn3QMw2KE4UG4jWEgmULJlUcmIxK7WaJApJBwD/BfQAroqIC2sRR3fTpXolWdqaWMrGic66mKonCkk9gF8DnwcWA09Imh4RC6sdi60vT68kS5dJNLXW3RJd0T1G90gLV4sexTjg+fTZ2Ui6GTgMcKLo4jqaaMrAyS6HohNdV29/A1SLRDEYeLlieTEwvulGkiYDk9PFtbsddfr8KsTWUQOBf9Q6iBwcZ+fpCjGC4+xsXSXOnTqjkVokCjVTFh8piJgCTAGQVB8RdUUH1lGOs3N1hTi7QozgODtbV4qzM9qpxX0Ui4GhFctDgFdrEIeZmeVQi0TxBLCjpO0kbQp8BZhegzjMzCyHqp96ioj3JZ0M3Edyeew1EbGglWpTio+sUzjOztUV4uwKMYLj7GzdKk5FfGR4wMzMrJHnejIzs0xOFGZmlqmmiULSAZL+R9Lzks5oZr0k/TJdP1fS2Lx1SxTnIknzJM3prEvVOhDnJyU9ImmtpNPbUrdEcZbpeH4t/e89V9LDkkblrVuiOKtyPHPEeFga3xxJ9ZL2zFu3RHGW5m+zYrvPSPpA0sS21l1PRNTkRTKQ/QKwPbAp8DSwc5NtDgT+QHLvxW7AY3nrliHOdN0iYGBJjuc2wGeAHwOnt6VuGeIs4fHcHdgy/fylEv99NhtntY5nzhh78+GY6Ujg2ZIey2bjLNvfZsV2fwZ+D0zsyPGsZY+icSqPiHgXaJjKo9JhwPWReBToL2lQzrpliLOaWo0zIpZFxBPAe22tW5I4qylPnA9HxOvp4qMk9wTlqluSOKslT4yrI/0lA3rx4U24ZTuWLcVZTXmPySnA7cCydtRdTy0TRXNTeQzOuU2eup2lI3FC8of0J0mzlUxLUpSOHJOyHc8sZT2eJ5L0KttTtyM6EidU53jmilHSEZKeBe4FvtGWuiWIE0r0tylpMHAEcEVb6zanls+jyDOVR0vb5JoGpJN0JE6APSLiVUnbAPdLejYiHuzUCFuPoci6bdXRfZXueEr6HMkPcMP56lIez2bihOocz7zT9twJ3Clpb+B8YP+8dTtJR+KEcv1t/gL4t4j4QFpv83Ydz1r2KPJM5dHSNtWcBqQjcRIRDe/LgDtJun61irOIum3VoX2V7XhKGglcBRwWEcvbUrcEcVbreLbpeKQ/rjtIGtjWuh3UkTjL9rdZB9wsaREwEfhvSYfnrPtRRQ+8ZAzIbAz8DdiODwdVdmmyzUGsP0j8eN66JYmzF9Cn4vPDwAG1irNi23NYfzC7VMczI85SHU/g48DzwO7t/Y41jrMqxzNnjJ/gw0HiscAr6f9PZTuWLcVZqr/NJttfx4eD2e06np3+Jdr4hQ8E/koyCn9WWnYScFL6WSQPOXoBmAfUZdUtW5wkVxY8nb4WlCDOfyL5F8WbwMr0c98SHs9m4yzh8bwKeB2Yk77qS/r32Wyc1TyeOWL8tzSGOcAjwJ4lPZbNxlm2v80m215Hmijaezw9hYeZmWXyndlmZpbJicLMzDI5UZiZWSYnCjMzy+REYWZmmZworEuRtLqZspMkfb0T9zFTUl0z5XWSfpl+niBp987aZzP7miDpnqLaN2uLWk7hYdYpIqLpfDZF7aceaJg+egKwmuTGKrMNmnsU1uVJOqfhuRVpb+AX6XMX5ksal5bvkz4nYI6kpyT1Sct/kD5D4GlJF1Y0e6SkxyX9VdJe6bYTJN0jaTjJzU3fS9vbq6LtOZLekbRPkxgfk7RLxfJMSZ+WNC6N9an0faes75cuz09jQNKxaZxzJP1GUo9OOahmFdyjsA1Rr4jYPZ207RrgU8DpwHci4iFJvYE1kr4EHA6Mj4i3JQ2oaGPjiBgn6UDgbD6c+I2IWCTpCmB1RFycFo8GkHQI8AM+2tO4GTgKODudgv5jETFbUl9g74h4X9L+wE+Af87zJSWNAI4mmYzuPUn/DXwNuD7ncTLLxYnCNkQ3QTJpm6S+kvoDDwE/l3QjcEdELE5/mK+NiLfT7VdUtHFH+j4bGJ5np5J2BC4C9o2Ips/SuBW4nyTpHAXclpb3A6amdQPYpA3fcz/g08AT6Qyhm7P+swfMOoUThW2Ims5LExFxoaR7Sea5eTRNEmpm2wZr0/cPyPH/iaReJMngm5HOItokgFckLU9ncT0a+Fa66nzggYg4Ij2dNLOZ5t9n/dPEPRt2C0yNiDNbi8+sIzxGYRuiowGUPM/4jYh4Q9IOETEvIv6TZED6k8CfgG9I2iLdfkCLLX7UKqBPxfK1JL2Tv2TUuZnktFS/iJiXlvUjmYEU4IQW6i0imakUJc9j3y4tnwFMTJ9/gKQBkoa14TuY5eJEYV3NFpIWV7z+tZltXpf0MMnTvU5My05NB4GfBt4B/hARfwSmA/WS5pCMY+R1N3BEOoi8D8mc/9+oGND+yOW1wO+Ar5D0PBr8DPippIdInmfcnNuBAWmM3yaZ+ZOIWAj8O8lT1eaSnNqq9iN4rRvw7LG2QZE0k+QZFvWtbWtm+bhHYWZmmdyjMDOzTO5RmJlZJu2dlnsAAAAcSURBVCcKMzPL5ERhZmaZnCjMzCyTE4WZmWX6/72cjk7iBIzhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "# Plot the distributions\n",
    "ax = sns.distplot(lips_val['gasdb'], kde=False,\n",
    "                  bins=1000,\n",
    "                  norm_hist=True,\n",
    "                  label='gasdb',\n",
    "                 )\n",
    "ax = sns.distplot(lips_val['cathub'], kde=False,\n",
    "                  bins=100,\n",
    "                  norm_hist=True,\n",
    "                  label='cathub',\n",
    "                 )\n",
    "\n",
    "# Format\n",
    "plot_lims = (0, 0.4)\n",
    "_ = ax.set_xlim(plot_lims)\n",
    "_ = ax.set_xlabel('Lipschitz value')\n",
    "_ = ax.set_ylabel('Normalized count')\n",
    "_ = ax.set_title('validation data')\n",
    "_ = ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like Cathub data is actually better looking than GASdb. This means that it should be easier for models to make predictions from Cathub data than GASdb data, which is consistent with what we saw before. Let's keep doing with some other distance metrics.\n",
    "# RBF distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculation and plot timing for val set of gasdb:\n",
      "CPU times: user 1.51 s, sys: 160 ms, total: 1.67 s\n",
      "Wall time: 1.67 s\n",
      "CPU times: user 9.08 s, sys: 220 ms, total: 9.3 s\n",
      "Wall time: 9.3 s\n",
      "Calculation and plot timing for val set of cathub:\n",
      "CPU times: user 560 ms, sys: 51.8 ms, total: 612 ms\n",
      "Wall time: 611 ms\n",
      "CPU times: user 2.41 s, sys: 64 ms, total: 2.48 s\n",
      "Wall time: 2.47 s\n",
      "Calculation and plot timing for train set of gasdb:\n",
      "CPU times: user 38.2 s, sys: 4.46 s, total: 42.6 s\n",
      "Wall time: 42.6 s\n",
      "CPU times: user 4min 6s, sys: 12.5 s, total: 4min 18s\n",
      "Wall time: 4min 18s\n",
      "Calculation and plot timing for train set of cathub:\n",
      "CPU times: user 16 s, sys: 1.8 s, total: 17.8 s\n",
      "Wall time: 17.8 s\n",
      "CPU times: user 1min 9s, sys: 2.1 s, total: 1min 11s\n",
      "Wall time: 1min 11s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7xVdZ3/8dcbvKCCmEJOgVwyxtAEpBN4v6SZOF7KwbxP0oWxGWsqtZrs5wWr8Zc5TZMmoSJqouYt0fDSWAyiowkmKmiJSHLSEQVTUEHBz/yx1oHNce1z1uHstfc5e7+fj8d5nL2+67I/34Puz/5+v2t9v4oIzMzMWutR6wDMzKxrcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYdYJkg6U1NyB42dJ+mKRMZlVihOEWRclaYmkQ2odhzUuJwgzM8vkBGENT9K3Jd3cquwnkv4zfT1B0lOSVkpaLOkfO3DtT0p6WtJrki4BVLJvZ0m/lbRc0iuSrpO0XbrvWmAQcIekVZK+mZbfJOl/0+vNlrRbBf4EZpmcIMzgeuBwSdsCSOoJfBaYnu5fBhwBbAtMAH4saXR7F5XUD7gF+C7QD3gW2Kf0EODfgA8Cw4GdgPMAIuIU4HngyIjoHRE/TM+5CxgGvB94FLhuk2psloMThDW8iPgzyYftp9OiTwBvRsRD6f5fR8Szkfhv4F5gvxyXPhxYGBE3R8Q7wH8A/1vyvosi4jcRsSYiXgb+HTignVinRsTKiFhDkkxGSurboQqb5eQEYZaYDpyQvj6RDa0HJI2T9JCkFZL+SvLB3y/HNT8ILG3ZiGRmzPXbkt4v6QZJf5H0OvCLtq4rqaekCyU9mx6/JN2VJxazDnOCMEvcBBwoaSDwGdIEIWlLkm6iHwE7RsR2wExKxhLa8CJJtxHptVS6TdK9FMCIiNgWOLnVdVtPtXwicDRwCNAXGNJy6RyxmHWYE4QZkHbxzAKuAp6LiKfSXVsAWwIvA2sljQMOzXnZXwO7STpG0mbAV4G/KdnfB1gF/FXSAOCsVue/BHyo1fFrgOXA1sAPcsZhtkmcIMw2mE7y7Xx991JErCT5YP8l8CrJt/gZeS4WEa8AxwIXknyoDwMeKDnkfGA08BpJMrm11SX+DfiupL9KOhO4Bvgz8BdgIfBQx6pn1jHygkFmZpbFLQgzM8vkBGFmZpmcIMzMLJMThJmZZdqs1gFUUr9+/WLIkCG1DsPMrNuYN2/eKxHRP2tfXSWIIUOGMHfu3FqHYWbWbUj6c7l97mIyM7NMThBmZpbJCcLMzDLV1RhElnfeeYfm5mZWr15d61C6jF69ejFw4EA233zzWodiZl1Y3SeI5uZm+vTpw5AhQ0gm02xsEcHy5ctpbm5m6NChtQ7HzLqwuu9iWr16NTvssIOTQ0oSO+ywg1tUZtauuk8QgJNDK/57mFkeDZEgzMys4+p+DKK16Q8/X9HrnTh2UEWvl8eBBx7Ij370I5qamujduzerVq2qegxmVv8KSxCSpgJHAMsi4qMZ+88CTiqJYzjQPyJWSFoCrATWAWsjoqmoOM2qau5V2eVNE6obh1kORbYgpgGXkKyC9R4RcRFwEYCkI4GvR8SKkkMOSlfk6vYuuOACrrvuOnbaaSf69evHxz72Mfr27cuUKVN4++23+fCHP8y1117L1ltvzU033cT5559Pz5496du3L7Nnz+att95iwoQJLFy4kOHDh/PWW29tdP0zzjiD3/3ud7zvfe/jhhtuoH//zGlVzMw6pLAxiIiYDaxo98DECcD1RcVSS3PnzuWWW27hD3/4A7feeuv6uaKOOeYYHnnkEebPn8/w4cO58sorAZg0aRL33HMP8+fPZ8aMZGXLyy67jK233prHH3+cs88+m3nz5q2//htvvMHo0aN59NFHOeCAAzj//POrX0kzq0s1H6SWtDVwGHBLSXEA90qaJ2liO+dPlDRX0tyXX365yFA3yZw5czj66KPZaqut6NOnD0ceeSQATz75JPvttx+777471113HQsWLABgn3324dRTT+Xyyy9n3bp1AMyePZuTTz4ZgBEjRjBixIj11+/RowfHHXccACeffDJz5sypZvXMrI7VPEEARwIPtOpe2iciRgPjgH+WtH+5kyNiSkQ0RURTV+xaKbfm96mnnsoll1zCE088wbnnnrv+uYTJkyfzve99j6VLlzJq1CiWL18O5L811bewmlmldIUEcTytupci4oX09zLgNmBMDeKqiH333Zc77riD1atXs2rVKn79618DsHLlSj7wgQ/wzjvvcN11160//tlnn2Xs2LFMmjSJfv36sXTpUvbff//1xzz55JM8/vjj649/9913ufnmmwGYPn06++67bxVrZ2b1rKa3uUrqCxwAnFxStg3QIyJWpq8PBSZV6j2rfVvqxz/+cY466ihGjhzJ4MGDaWpqom/fvlxwwQWMHTuWwYMHs/vuu7Ny5UoAzjrrLJ555hkigoMPPpiRI0eyyy67MGHCBEaMGMGoUaMYM2ZDvtxmm21YsGDB+oHvG2+8sar1M7P6pXJdIJ2+sHQ9cCDQD3gJOBfYHCAiJqfHnAocFhHHl5z3IZJWAyQJbHpEfD/PezY1NUXrBYOeeuophg8f3pmqdNqqVavo3bs3b775Jvvvvz9Tpkxh9OjRNY2pK/xdGpJvc7UuRtK8co8SFNaCiIgTchwzjeR22NKyxcDIYqKqjYkTJ7Jw4UJWr17N5z73uZonBzOzPBruSepamD59eucv8kaZR0K26df5a1vllWspmOXRRVqaXWGQ2szMuiC3IMw6wy0Fq2NuQZiZWSYnCDMzy9R4XUyV7hKo4KDRrFmz2GKLLdh7772B5GnrI444gvHjx+e+xpAhQ5g7dy79+nnw2sw6xy2ILmTWrFk8+OCDtQ7DzAxwgqiKa665hhEjRjBy5EhOOeUU7rjjDsaOHcsee+zBIYccwksvvcSSJUuYPHkyP/7xjxk1ahT3338/kEzUt/fee/OhjzZx823J7K6zZj/AEeNPXH/9008/nWnTpq3fvuiiixgzZgxjxoxh0aJFVa2rmdWPxutiqrIFCxbw/e9/nwceeIB+/fqxYsUKJPHQQw8hiSuuuIIf/vCHXHzxxZx22mn07t2bM888E4Arr7ySF198kTlz5vD0o//DUcedzPjPHNXue2677bb8/ve/55prruFrX/sad955Z9HVNLM65ARRsN/+9reMHz9+/ZjA9ttvzxNPPMFxxx3Hiy++yNtvv83QoUPLnv/pT3+aHj16sOvwXXhpWb7pzE844YT1v7/+9a93vhJm1pCcIAoWEe+ZgvsrX/kK3/jGNzjqqKOYNWsW5513Xtnzt9xyy42uBbDZZj15991315e3TBXeovT9PP13hfh5B2tAHoMo2MEHH8wvf/nL9es6rFixgtdee40BAwYAcPXVV68/tk+fPutndW3L4EE7sfDpP7FmzRpee+017rvvvo32t8zoeuONN7LXXntVqipm1mAarwVR5blMdtttN84++2wOOOAAevbsyR577MF5553Hsccey4ABA9hzzz157rnnADjyyCMZP348t99+Oz/96U/LXnOngQP47DFHM2LPAxj2tx9hjz322Gj/mjVrGDt2LO+++y7XX1+XK7maWRUUNt13LXTV6b4rosKT9dXN36Vaiu5i8nTfVqqKk/XVZLpvs27JYw1m63kMwszMMjVECyLrTqJGVk/dimYNpcrrRNR9C6JXr14sX77cH4qpiGD58uX06tWr1qGYWRdX9y2IgQMH0tzczMsv53vIrMtasyq7fMuO16tXr14MHDiwkwGZWb2r+wSx+eabt/mkcrdRrmk5yne/mFkxCutikjRV0jJJT5bZf6Ck1yQ9lv6cU7LvMEl/lLRI0reLitHMzMorcgxiGnBYO8fcHxGj0p9JAJJ6ApcC44BdgRMk7VpgnGZmlqGwBBERs4EVm3DqGGBRRCyOiLeBG4CjKxqcmZm1q9Z3Me0lab6kuyTtlpYNAJaWHNOclmWSNFHSXElzu/1AtJlZF1LLBPEoMDgiRgI/BX6Vlmc9sFD2HtWImBIRTRHR1L9//wLCNDNrTDVLEBHxekSsSl/PBDaX1I+kxbBTyaEDgRdqEKKZWUOrWYKQ9DdKH2+WNCaNZTnwCDBM0lBJWwDHAzNqFaeZWaMq7DkISdcDBwL9JDUD5wKbA0TEZGA88GVJa4G3gOMjedx5raTTgXuAnsDUiFhQVJxmZpatsAQRESe0s/8S4JIy+2YCM4uIy8zM8qn7J6mtk6o8OZiZdR21vs3VzMy6KLcgujt/wzezgrgFYWZmmZwgzMwskxOEmZll8hiE1QePxZhVnBOEWVfmxGc15ARh9a1eP2DrtV7WpXgMwszMMjlBmJlZJncx1YmHn9t48b5n1z2//vWJYwdVOxwzqwNOEA1g+sPPZ5Y7cXQh5cYUKnUdj03YJnCCMDOrlUp9MSiIxyDMzCyTE4SZmWVygjAzs0xOEGZmlqndBCHp2jxlZmZWX/LcxbRb6YaknsDHigmn8fgWVDPrqsomCEn/CnwH2ErS6y3FwNvAlPYuLGkqcASwLCI+mrH/JOBb6eYq4MsRMT/dtwRYCawD1kZEU94K1bvWD8SZmRWlbBdTRPxbRPQBLoqIbdOfPhGxQ0T8a45rTwMOa2P/c8ABETECuID3Jp2DImKUk4OZWW2028UUEf8qaQAwuPT4iJjdznmzJQ1pY/+DJZsPAQPbi8XMzKqn3QQh6ULgeGAhSZcPQABtJogO+gJwV8l2APdKCuDnEVG2S0vSRGAiwKBB7rc3M6uUPIPUnwF2iYg1RQQg6SCSBLFvSfE+EfGCpPcDv5H0dLkWS5o8pgA0NTVFETGamTWiPAliMbA5UPEEIWkEcAUwLiKWt5RHxAvp72WSbgPGUNkWi7Wh9M6qnZ/fMCg+duj2tQjHzGokT4J4E3hM0n2UJImI+Gpn3ljSIOBW4JSI+FNJ+TZAj4hYmb4+FJjUmfcyM7OOy5MgZqQ/HSLpeuBAoJ+kZuBckpYIETEZOAfYAfiZJNhwO+uOwG1p2WbA9Ii4u6Pvb2ZmnZPnLqarN+XCEXFCO/u/CHwxo3wxMHJT3tPMzConz11Mz5HcVbSRiPhQIREZ8N4nrEvHAop6DzOzUnm6mEofVOsFHAt4tNLMrM7l6WJa3qroPyTNIRlDsJz8bd3Mups8XUyjSzZ7kLQo+hQWkZmZdQl5upguLnm9FlgCfLaQaMzMrMvI08V0UDUCMTOzriVPF1NfkmcY9k+L/huYFBGvFRmYdc7Oz9+UWf7soGOrHIlZ1+C1VzouTxfTVOBJNnQrnQJcBRxTVFBmZpvKN4RUTp4EsXNE/H3J9vmSHisqIKuNci2ORldugSbPS2WNIE+CeEvSvhExB0DSPsBbxYbVfdXzt5fSD8tn122op5voZvUpT4L4MnB1OhYB8CpwamERWWOZe1V2edOE6sZh3U49fxnrKvLcxfQYMFLStun26+2cYmYNpLsP/raVaLpLHYqS5y6mHwA/jIi/ptvvA86IiO8WHZxZtbR8SBQx55VZd5Wni2lcRHynZSMiXpV0OOAEYQ3Lg9fV466k2smTIHpK2rJlyVFJWwFbFhtW4yr6bqJGv1up5YO9dJDdzLLlSRC/AO6TdBXJtN+fBzZpjYhG5AfWqq/ct3vrHH+Tbzx5Bql/KOlx4BBAwAURcU/hkXVx/p/FzOpdnhYE6ZKfXvbTzHLr7nc3QX3UoTNyJQgzqy+1vLWz0T90uxMnCLM6timD8rXqPnW3bddTWIKQNBU4AlgWER/N2C/gJ8DhwJvAqRHxaLrvsHRfT+CKiLiwqDjNKsm3v1o9KZsgJD1BctdSpogY0c61pwGXANeU2T8OGJb+jAUuA8ZK6glcCnwSaAYekTQjIha2835mZlZBbbUgjkh//3P6+9r090kk3/jbFBGzJQ1p45CjgWsiIoCHJG0n6QPAEGBRRCwGkHRDeqwThJlZFZVNEBHxZ0hmb42IfUp2fVvSA8CkTr73AGBpyXZzWpZVPrbcRSRNBCYCDBrkQS4zs0rJMwaxTavpvvcGtqnAeyujLNoozxQRU4ApAE1NTWWPsxorN2trjTT6E+XWOY1yJ1aeBPEFYGo63XcAr5E8Td1ZzcBOJdsDgReALcqUm1kZlX563DMAGOR7knoeG6b7VgXXop4BnJ6OMYwFXouIFyW9DAyTNBT4C3A8cGKF3tPqSOm3OM/CalZ5eab73hH4AfDBiBgnaVdgr4i4sp3zrgcOBPpJagbOBTYHiIjJwEySW1wXkQx6T0j3rZV0OnAPyW2uUyNiwaZVz6qhUs1tr1hn1rXk6WKaBlwFnJ1u/wm4EWgzQUTECe3sDzbcIdV630ySBGJdlLsgzOpfjxzH9IuIXwLvQvINH1hXaFRmZlZzeVoQb0jagfROIkl7kgxUm22Sjg6oegoGs9rIkyDOIBlQ3jl9/qE/4H4EM7O8utht3nnluotJ0gHALiTPKPwxIt4pPDIzew8vhmTVlOcupmeBi9I7j1rK7oyII9o4zdrhB7XMrKvL08X0DnCQpLHAP0bE2yTTYZhZQRq1peC747qWPHcxvRkRxwFPAfdLGkwbU1+YmVl9yNOCEKxfm3oeyQNsntze2uUnnTfwOhHWHeVJEOe0vIiI+yR9CvhccSGZNY5G7Uqy7qGtBYM+EhFPA3+RNLrV7juLDcvMzGqtrRbEGcCXgIsz9gXwiUIiMjOzLqGtBYO+lP4+qHrhmLXNtwebVU9bXUzHtHViRNxa+XDMzKyraKuL6cg29gXgBGHWTfj5guooO/V9zyoHUiFtdTFNqGYgZqXclWRWe3luc0XS3wG7Ab1ayiJiUlFBmZlZ7bX7JLWkycBxwFdIHpo7FhhccFxmZlZjeaba2Dsi/gF4NSLOB/YCdio2LDMzq7U8CeKt9Pebkj5IMnnf0OJCMjOzriDPGMSdkrYDLgIeJbmD6YpCozIzs5rLs2DQBenLWyTdCfSKiFxLjko6DPgJ0BO4IiIubLX/LOCkkliGA/0jYoWkJcBKkvWv10ZEU573NDOzysizYFBP4O+AIS3HSyIi/j3HeZcCnwSagUckzYiIhS3HRMRFJC0TJB0JfD0iSmcvOygiXulQjczMrCLydDHdAawGngDe7cC1xwCLImIxgKQbgKOBhWWOPwG4vgPXNzOzAuVJEAMjYsQmXHsAsLRkuxkYm3WgpK2Bw4DTS4oDuFdSAD+PiCllzp0ITAQYNGjQJoRpZmZZ8tzFdJekQzfh2sooK7cS3ZHAA626l/aJiNHAOOCfJe2fdWJETImIpoho6t+//yaEaWZmWfIkiIeA2yS9Jel1SSslvZ7jvGY2fl5iIPBCmWOPp1X3UkS8kP5eBtxG0mVlZmZVkidBXEzycNzWEbFtRPSJiG1znPcIMEzSUElbkCSBGa0PktQXOAC4vaRsG0l9Wl4DhwJP5nhPMzOrkDxjEM8AT0ZEue6hTBGxVtLpJGtY9wSmRsQCSael+yenh34GuDci3ig5fUeSVktLjNMj4u6OvL+Z1T/PUlusPAniRWCWpLuANS2F7d3mmh4zE5jZqmxyq+1pwLRWZYuBkTliM7NuyLP1dg95EsRz6c8W6Y+ZmTWANhNE+rBb74g4q0rxmFkV+Zu8taXNBBER6ySNrlYwZmb16OHnVmSWjx26fZUj6Zg8XUyPSZoB3ASsH0j2mtRmZvUtT4LYHlgOfKKkzGtSWyZ3WdQ33zXUWPLM5uq1qc3MGlCeJUcHSrpN0jJJL0m6RdLAagRnZma1k+dJ6qtInoD+IMkEfHekZWZmVsfyjEH0j4jShDBN0teKCsjMrFF09bub8rQgXpF0sqSe6c/JJIPWZmZWx/IkiM8DnwX+l2TajfFpmZmZ1bE8dzE9DxxVhVjMzDL59unaKJsgJJ3TxnkRERcUEI+ZWaf5eY3KaKsF8UZG2TbAF4AdACcIM7M6VjZBRMTFLa/TxXv+BZgA3ECyiJCZWbfilkXHtDeb6/bAN4CTgKuB0RHxajUCMzOz2mprDOIi4BhgCrB7RKyqWlRmZlZzbbUgziBZQe67wNnp8p8AIhmkzrMutZk1AHfd1Ke2xiDyPCNhZmZ1Ks9UG5aD79M2s3pTaIKQdBjwE6AncEVEXNhq/4HA7SRrXgPcGhGT8pxbLdMffr4Wb2tm3VC9fVEsLEGk61lfCnwSaAYekTQjIha2OvT+iDhiE881M7OCFDnOMAZYFBGLI+Jtkucnjq7CuWZmVgFFJogBwNKS7ea0rLW9JM2XdJek3Tp4rpmZFaTIMQhllEWr7UeBwRGxStLhwK+AYTnPTd5EmghMBBg0aNCmR2tmZhspMkE0AzuVbA8EXig9ICJeL3k9U9LPJPXLc27JeVNIHuajqakpM4mYmXUnXWUhoSK7mB4BhkkaKmkL4HiSpUvXk/Q3Sp/AkzQmjWd5nnPNzKxYhbUgImKtpNOBe0huVZ0aEQsknZbun0yy+NCXJa0F3gKOj4gAMs8tKlYzs46ot9tZyyn0OYiImAnMbFU2ueT1JcAlec81M7Pq8ZPUZmbdRNmxiaZi3s/zLZmZWSa3IMysMI3SV1+v3IIwM7NMThBmZpbJXUwd5CazWf3xgkfZ3IIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwsU6HrQUg6DPgJ0BO4IiIubLX/JOBb6eYq4MsRMT/dtwRYCawD1kZEQctym5lla/T1XwpLEJJ6ApcCnwSagUckzYiIhSWHPQccEBGvShoHTAHGluw/KCJeKSpGMzMrr8gupjHAoohYHBFvAzcAR5ceEBEPRsSr6eZDwMAC4zEzsw4oMkEMAJaWbDenZeV8AbirZDuAeyXNkzSx3EmSJkqaK2nuyy+/3KmAzcxsgyLHIJRRFpkHSgeRJIh9S4r3iYgXJL0f+I2kpyNi9nsuGDGFpGuKpqamzOubmVnHFdmCaAZ2KtkeCLzQ+iBJI4ArgKMjYnlLeUS8kP5eBtxG0mVlZmZVUmSCeAQYJmmopC2A44EZpQdIGgTcCpwSEX8qKd9GUp+W18ChwJMFxmpmZq0U1sUUEWslnQ7cQ3Kb69SIWCDptHT/ZOAcYAfgZ5Jgw+2sOwK3pWWbAdMj4u6iYjUzs/cq9DmIiJgJzGxVNrnk9ReBL2actxgYWWRsZmbWtkITRHcy/eHnN9pueUBm51oEY2bWBXiqDTMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZllcoIwM7NMThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYWZmmZwgzMwsU6EJQtJhkv4oaZGkb2fsl6T/TPc/Lml03nPNzKxYhSUIST2BS4FxwK7ACZJ2bXXYOGBY+jMRuKwD55qZWYGKbEGMARZFxOKIeBu4ATi61TFHA9dE4iFgO0kfyHmumZkVaLMCrz0AWFqy3QyMzXHMgJznAiBpIknrA2CVpD9uYrz9gFc28dzuqNHqC65zo2jAOp/ZmToPLrejyAShjLLIeUyec5PCiCnAlI6F9l6S5kZEU2ev0100Wn3BdW4UrnPlFJkgmoGdSrYHAi/kPGaLHOeamVmBihyDeAQYJmmopC2A44EZrY6ZAfxDejfTnsBrEfFiznPNzKxAhbUgImKtpNOBe4CewNSIWCDptHT/ZGAmcDiwCHgTmNDWuUXFmup0N1U302j1Bde5UbjOFaKIzK59MzNrcH6S2szMMjlBmJlZpoZKEJ2Z+qO7ylHnk9K6Pi7pQUkjaxFnJeWdpkXSxyWtkzS+mvEVIU+dJR0o6TFJCyT9d7VjrLQc/233lXSHpPlpnSfUIs5KkTRV0jJJT5bZX/nPr4hoiB+Swe5ngQ+R3EY7H9i11TGHA3eRPIexJ/BwreOuQp33Bt6Xvh7XCHUuOe63JDdKjK913FX4d94OWAgMSrffX+u4q1Dn7wD/P33dH1gBbFHr2DtR5/2B0cCTZfZX/POrkVoQnZn6o7tqt84R8WBEvJpuPkTyzEl3lnealq8AtwDLqhlcQfLU+UTg1oh4HiAiunu989Q5gD6SBPQmSRBrqxtm5UTEbJI6lFPxz69GShDlpvXo6DHdSUfr8wWSbyDdWbt1ljQA+AwwuYpxFSnPv/PfAu+TNEvSPEn/ULXoipGnzpcAw0kesn0C+JeIeLc64dVExT+/inySuqvpzNQf3VXu+kg6iCRB7FtoRMXLU+f/AL4VEeuSL5fdXp46bwZ8DDgY2Ar4H0kPRcSfig6uIHnq/CngMeATwM7AbyTdHxGvFx1cjVT886uREkRnpv7ornLVR9II4ApgXEQsr1JsRclT5ybghjQ59AMOl7Q2In5VnRArLu9/269ExBvAG5JmAyOB7pog8tR5AnBhJB30iyQ9B3wE+H11Qqy6in9+NVIXU2em/uiu2q2zpEHArcAp3fjbZKl26xwRQyNiSEQMAW4G/qkbJwfI99/27cB+kjaTtDXJ7MhPVTnOSspT5+dJWkxI2hHYBVhc1Sirq+KfXw3TgohOTP3RXeWs8znADsDP0m/Ua6Mbz4SZs851JU+dI+IpSXcDjwPvAldERObtkt1Bzn/nC4Bpkp4g6X75VkR022nAJV0PHAj0k9QMnAtsDsV9fnmqDTMzy9RIXUxmZtYBThBmZpbJCcLMzDI5QZiZWSYnCDMzy+QEYd2OpEmSDtnEc0+rxDQTkqZJei6dHXW+pINL9s1KZxl9TNJTkiaW7Fsi6Yl032OS9m7nfc6TdGb6us16S/q0pF07WzezFg3zHITVj4g4J6tcUs+IWNfOuZV8DuKsiLg5naZkCjCsZN9JETFX0vbAs5KmpZPKARy0Kffjl6t3iU8Dd5LM2mrWaW5BWE1JGiLpaUlXp3PY35w+6YukcyQ9IulJSVPSWTlbvr2PT18vSY+bA5wkaV5aPlJSpE+KI+lZSVu3+kb+VUkL0/e9IS3bRsm8+49I+oOkrJlgW/sfyn3i6/8AAAQQSURBVE+K1ht4A2gzcbX6m5ydtkD+i+Tp35by0npfWBL7j9KWyFHARWnLZGdJX0rrMV/SLSV/12lK1g14UNJilayHIembaQtnvqQL07KdJd2tZJK/+yV9JC0/Nv23ma9k6g6rM25BWFewC/CFiHhA0lTgn4AfAZdExCQASdcCRwB3ZJy/OiL2TY/7lqRtgf2AuSTTS8wBlkXEm9p4cr5vA0MjYo2k7dKys4HfRsTn07LfS/qvdA6jcg4DWk/VcZ2kNSStiq+1atn8TtI6YE1EjC09SdLHSKaN2IPk/89HgXmtjtmeZDbaj0RESNouIv4qaQZwZ0TcnB7314i4PH39PZLJGH+aXuYDJBMzfoRkioabJY0jaYWMTf9W26fHTgFOi4hnJI0FfkYyAd45wKci4i8lfz+rI25BWFewNCIeSF//gg0zyh4k6eF0qoRPALuVOf/GktcPAvuQLK7yg/T3fsD9Gec9TvJBfjIb1gk4FPi2pMeAWUAvYFCZ971I0uI05h+02ndSRIxIzz1T0uCSfQdFxKjWySG1H3BbRLyZzjraen4hgNeB1cAVko4hmVYhy0fTb/xPACex8d/vVxHxbkQsBHZMyw4BroqINwEiYoWk3iSLSt2U/k1+TpJcAB4gmcriSyTTXVidcYKwrqD1fC8hqRfJN9XxEbE7cDnJh3WW0m/395N8yA4mmaBuJEnCyeoC+TvgUpJpsOdJ2oxkzp6/Tz/AR0XEoHQeo6vSrpuZJeefBXwY+C5wdWbFIl4maQVkJYNy2pz/JiLWkiyYcwvJN/67yxw6DTg9/fudz8Z/vzUlr1Xyu/V79wD+WvL3GBURw9M4TiOp+07AY5J2aK9i1r04QVhXMEjSXunrE4A5bPgweyX9Fpt33ejZwMnAM+niMCtIJjB7oPQgST2AnSLid8A3SZbk7E0y+dtXSsY79gCIiAnph+PhpddJ3+MnQA9Jn2odTNrvvwfJ8ph54/+MpK0k9QGOzLhmb6BvRMwEvgaMSnetBPqUHNoHeFHS5iQtiPbcC3y+ZKxi+7QV85ykY9MyKV23XNLOEfFwOnj+ChtPNW11wGMQ1hU8BXxO0s+BZ4DL0j7wy0lWAltCMr1zuyJiSfrZ3tJimAMMLFlWtUVP4BeS+pJ8c/5x2o9/AcmCQo+nSWIJydhHW+8ZaR//N0kSDCRdV28BWwLTImJe2QtsfK1HJd1IstDNn8nuGusD3J62sgR8PS2/Abhc0ldJEur/Ax5Or/MEGyePrPe+W9IoYK6kt0lmB/0OSXK5TNJ3SWYPvYFkDeiLJA1LY7gvLbM64tlcraYkDSEZWP1ojUMxs1bcxWRmZpncgjAzs0xuQZiZWSYnCDMzy+QEYWZmmZwgzMwskxOEmZll+j+JHXeNxwJ8BAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwV1Z338c8XxKBCcGnjGPYYHuMGaDp03DWaRB3RxMG4TyQmPGYeMpM8amJixgUziY/GJ5OJRkRF1Igat4gOLhkTB9HRCEZWTUQg0sFxARdQQJHf/FHVeGnrdlfTt3q59/t+vfrV9546VfU7Ld7fPedUnVJEYGZm1lyPzg7AzMy6JicIMzPL5ARhZmaZnCDMzCyTE4SZmWVygjAzs0xOEGabQdJESf9coWNdKOlXlTiWWSU5QVjNkbRU0uHtOUZEnBkRF1cqprwkTZH0o44+r9UmJwizZiRt0dkxmHUFThBWUyTdBAwC7pW0WtJ3JQ2RFJLOkPQi8Lu07u2S/lvSm5JmSNqj5Dgbv8lLOkRSo6SzJL0i6SVJY1uIYaik/5S0StJvgbpm2zPPK2kccArw3TT2e9PycyW9kB5voaQvV/avZrXKCcJqSkScBrwIjI6IPhFxacnmg4HdgC+m7+8HhgEfA54Gbm7h0H8D9AP6A2cAV0rarkzdqcBsksRwMfDVZtszzxsRk9LXl6axj07rvwAcmJ7/IuBXknZuIVazXJwgzD5wYUS8HRFrACJickSsioh1wIXACEn9yuz7HjAhIt6LiOnAamDX5pUkDQI+A/xzRKyLiBnAvaV12nheIuL2iFgeERsi4jbgeWBUG9tu9iFOEGYfWNb0QlJPSZekQzdvAUvTTXWZe8KKiFhf8v4doE9GvY8Dr0fE2yVlf2nHeZH095KekfSGpDeAPVuqb5aXE4TVonJLGJeWnwwcCxxOMnQzJC1XO8/9ErCdpG1Kyga14bybxC5pMHANMB7YISK2BeZXIE4zJwirSS8Dn2ilTl9gHbAC2Br4cSVOHBF/AWYBF0naUtIBwOiSKq2dt3ns25AkjVcB0snxPSsRq5kThNWinwA/TIdkzi5T50aSoZ+/AguBJyp4/pOBBmAlcEF6rrznvQ7YPY39NxGxELgc+C+S5LEX8FgFY7UaJj8wyMzMsrgHYWZmmZwgzMwskxOEmZllcoIwM7NMVbUoWV1dXQwZMqSzwzAz6zZmz579WkTsmLWtqhLEkCFDmDVrVmeHYWbWbUj6S7ltHmIyM7NMThBmZpbJCcLMzDJV1RyEmRnAe++9R2NjI2vXru3sULqM3r17M2DAAHr16pV7HycIM6s6jY2N9O3blyFDhiB5YduIYMWKFTQ2NjJ06NDc+3mIycyqztq1a9lhhx2cHFKS2GGHHdrco3KCMLOq5OSwqc35ezhBmJlZJs9BmFnVm/rkixU93skNg1qvVGGHHHIIP/3pT6mvr6dPnz6sXr268HM6QVh1m3V9dnn92I6Nw6wbcoIwMyvAxRdfzM0338zAgQOpq6vj05/+NP369WPSpEm8++67fPKTn+Smm25i66235vbbb+eiiy6iZ8+e9OvXjxkzZrBmzRrGjh3LwoUL2W233VizZs0mxz/rrLP4/e9/z3bbbcett97KjjtmLqfULp6DMDOrsFmzZnHnnXfyxz/+kbvuumvjGnHHHXccTz31FHPmzGG33XbjuuuuA2DChAk8+OCDzJkzh2nTpgFw1VVXsfXWWzN37lzOO+88Zs+evfH4b7/9Nvvssw9PP/00Bx98MBdddFEh7XCCMDOrsJkzZ3Lsscey1VZb0bdvX0aPHg3A/PnzOfDAA9lrr724+eabWbBgAQD7778/p59+Otdccw3vv/8+ADNmzODUU08FYPjw4QwfPnzj8Xv06MEJJ5wAwKmnnsrMmTMLaYcThJlZhUVEZvnpp5/OFVdcwbx587jgggs23pcwceJEfvSjH7Fs2TJGjhzJihUrgPyXphZ1SW9hCULSZEmvSJpfZvs5kp5Jf+ZLel/S9um2pZLmpdu8freZdSsHHHAA9957L2vXrmX16tX8+7//OwCrVq1i55135r333uPmm2/eWP+FF16goaGBCRMmUFdXx7JlyzjooIM21pk/fz5z587dWH/Dhg3ccccdAEydOpUDDjigkHYUOUk9BbgCuDFrY0RcBlwGIGk08J2IWFlS5dCIeK3A+MysRnT0Zamf+cxnOOaYYxgxYgSDBw+mvr6efv36cfHFF9PQ0MDgwYPZa6+9WLVqFQDnnHMOzz//PBHBYYcdxogRI9h1110ZO3Ysw4cPZ+TIkYwaNWrj8bfZZhsWLFiwceL7tttuK6QdKtcVqsjBpSHAfRGxZyv1pgK/j4hr0vdLgfq2Joj6+vrwA4NsE77MtSY9++yz7Lbbbp0aw+rVq+nTpw/vvPMOBx10EJMmTWKfffbp1Jiy/i6SZkdEfVb9Tr/MVdLWwBHA+JLiAB6SFMDVETGphf3HAeMABg3q+JtXzMyyjBs3joULF7J27Vq++tWvdnpy2BydniCA0cBjzYaX9o+I5ZI+BvxW0nMRMSNr5zR5TIKkB1F8uGZmrZs6dWpnh9BuXeEqphOBW0oLImJ5+vsV4G5gVMZ+ZmZWoE5NEJL6AQcD95SUbSOpb9Nr4AtA5pVQZmZWnMKGmCTdAhwC1ElqBC4AegFExMS02peBhyLi7ZJddwLuTq/r3QKYGhEPFBWnmZllKyxBRMRJOepMIbkctrRsMTCimKjMzCyvrjBJbWZWrHKXO2+uCl4m/cgjj7Dllluy3377Acnd1kcffTRjxozJfYwhQ4Ywa9Ys6urqKhYXdI1JajOzmvXII4/w+OOPd3YYmZwgzMwKcOONNzJ8+HBGjBjBaaedxr333ktDQwN77703hx9+OC+//DJLly5l4sSJ/OxnP2PkyJE8+uijQLJQ33777ccnPvGJjUtqPPLIIxx99NEbjz9+/HimTJmy8f1ll13GqFGjGDVqFIsWLapIGzzEZGZWYQsWLOBf/uVfeOyxx6irq2PlypVI4oknnkAS1157LZdeeimXX345Z555Jn369OHss88G4LrrruOll15i5syZPPfccxxzzDG5hps++tGP8oc//IEbb7yRb3/729x3333tbocThJlZhf3ud79jzJgxG+cEtt9+e+bNm8cJJ5zASy+9xLvvvsvQoUPL7v+lL32JHj16sPvuu/Pyyy/nOudJJ5208fd3vvOd9jcCDzGZmVVcRHxoCe5vfetbjB8/nnnz5nH11VdvXOo7y0c+8pFNjgWwxRZbsGHDho3lzfcvPV+llv92gjAzq7DDDjuMX//61xuf67By5UrefPNN+vfvD8ANN9ywsW7fvn03ruraksGDB7Nw4ULWrVvHm2++ycMPP7zJ9qYVXW+77Tb23XffirTDQ0xmVv06ePXePfbYg/POO4+DDz6Ynj17svfee3PhhRdy/PHH079/fz772c+yZMkSAEaPHs2YMWO45557+MUvflH2mAMHDuQrX/kKw4cPZ9iwYey9996bbF+3bh0NDQ1s2LCBW265pcxR2qbQ5b47mpf7tg/xct81qSss990VtXW5bw8xmZlZJicIMzPL5DkIs1IekqoaWVcS1bLNmU5wD8LMqk7v3r1ZsWLFZn0oVqOIYMWKFfTu3btN+7kHYWZVZ8CAATQ2NvLqq692dihdRu/evRkwYECb9nGCMLOq06tXrxbvVLZ8nCDMOpLnOKwb8RyEmZllcoIwM7NMThBmZpapsDkISZOBo4FXImLPjO2HAPcAS9KiuyJiQrrtCODnQE/g2oi4pKg4rYN47N2s2ymyBzEFOKKVOo9GxMj0pyk59ASuBI4EdgdOkrR7gXGamVmGwhJERMwAVm7GrqOARRGxOCLeBW4Fjq1ocGZm1qrOnoPYV9IcSfdL2iMt6w8sK6nTmJZlkjRO0ixJs3xTjJlZ5XRmgngaGBwRI4BfAL9Jy7MWTyl7v3xETIqI+oio33HHHQsI08ysNnVagoiItyJidfp6OtBLUh1Jj2FgSdUBwPJOCNHMrKZ12p3Ukv4GeDkiQtIokmS1AngDGCZpKPBX4ETg5M6K0wzwVVhWk4q8zPUW4BCgTlIjcAHQCyAiJgJjgG9KWg+sAU6MZOnF9ZLGAw+SXOY6OSIWFBWnmZllKyxBRMRJrWy/AriizLbpwPQi4rIqVe4bvpltNi/WZ7XJCcWsVU4Q3YXHwM2sg3X2fRBmZtZFOUGYmVkmJwgzM8vkBGFmZpmcIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYmVkmL7Vhm8dLf5hVPfcgzMwsU6sJQtJNecrMzKy65OlB7FH6RlJP4NPFhGNmZl1F2TkISd8HfgBsJemtpmLgXWBSB8RmtaC7z2V09/jNWlA2QUTET4CfSPpJRHy/A2My8wevWRfQ6lVMEfF9Sf2BwaX1I2JGkYGZZfKT4Mw6TKsJQtIlwInAQuD9tDiAFhOEpMnA0cArEbFnxvZTgO+lb1cD34yIOem2pcCq9HzrI6I+T2OsAP5ANqtZee6D+DKwa0Ssa+OxpwBXADeW2b4EODgiXpd0JMm8RkPJ9kMj4rU2ntM6mxOKWdXIcxXTYqBXWw+cDkGtbGH74xHxevr2CWBAW89hZmbFydODeAd4RtLDwMZeRET8YwXjOAO4v+R9AA9JCuDqiPBVU9a9uCdlVSBPgpiW/hRC0qEkCeKAkuL9I2K5pI8Bv5X0XLlJcUnjgHEAgwYNKipMM7Oak+cqphuKOrmk4cC1wJERsaLknMvT369IuhsYRZlJ8bR3MQmgvr4+iorVzKzW5LmKaQnJkM8mIuIT7TmxpEHAXcBpEfHnkvJtgB4RsSp9/QVgQnvOVQumPvli2W0nN7hnZWZtl2eIqfQS097A8cD2re0k6RbgEKBOUiNwAelkd0RMBM4HdgB+KQk+uJx1J+DutGwLYGpEPJCzPWZmViF5hphWNCv6V0kzST7gW9rvpFa2fx34ekb5YmBEa3GZmVmx8gwx7VPytgdJj6JvYRGZmVmXkGeI6fKS1+uBpcBXConGrFZ57SnrgvIMMR3aEYGYmVnXkmeIqR/JBPNBadF/AhMi4s0iA7PKKXeFk69uMrOW5BlimgzM54NhpdOA64HjigrKzKyJv+B0njwJYpeI+LuS9xdJeqaogKxlTy7ZdHmrF94vf/+DmXVTXWROKk+CWCPpgIiYCSBpf2BNsWHVDn87Mts8Ld0c2hb+f628PAnim8AN6VwEwOvA6YVFZGZmXUKeq5ieAUZI+mj6/q1WdrEKaP7taJcXy66cbtatVOqbvxUvz1VMPwYujYg30vfbAWdFxA+LDs6sszWf82nSMLTV1WZqhodJq1eeIaYjI+IHTW/SJ8AdBThBdAG7vHh7ZvkLg45vdV//j21Fck+h+8uTIHpK+kjTI0clbQV8pNiwrKNtkmh6+tuxmeVLEL8CHpZ0Pcmy318DCntGhJl1TdXaI/BS+eXlmaS+VNJc4HBAwMUR8WDhkVmn8bi7WRfVwfdH5OlBkD6Pwc9ksKpVLim2tX41JNFq7SlY2/Xo7ADMzKxrytWDMLPq4l5CPrV+pZ8TRAfx/5BWpKZ/X80ve9445OXnSthmKJsgJM0juWopU0QMLySibs6JoPO1NJ9Q9BxBh81NNJusLHenfVM8XtTRNkdLPYij09//J/19U/r7FOCdwiKyLqsaJmbbOhltVsvKTlJHxF8i4i/A/hHx3YiYl/6cC3yxtQNLmizpFUnzy2yXpH+TtEjS3NJnX0s6QtKf0m3nbk7DzMysffLMQWzTbLnv/YBtcuw3BbgCuLHM9iOBYelPA3AV0CCpJ3Al8HmgEXhK0rSIWJjjnGZmhauVyes8CeIMYHK63HcAb5LcTd2iiJghaUgLVY4FboyIAJ6QtK2knYEhwKKIWAwg6da0rhOE1bwPJqM9VGbFy3Mn9Ww+WO5bFXwWdX9gWcn7xrQsq7yh3EEkjQPGAQwaVF3Z22qXJ5etK8iz3PdOwI+Bj0fEkZJ2B/aNiOvaeW5llEUL5ZkiYhIwCaC+vr5svVrTnlVereN40ty6sjxDTFOA64Hz0vd/Bm4D2psgGoGBJe8HAMuBLcuUm9WcconerCPkSRB1EfFrSd8HiIj1kt6vwLmnAePTOYYG4M2IeEnSq8AwSUOBvwInAidX4HxGdX3g+Nu3dXvlFt/rIvIkiLcl7UA6zCPpsyQT1S2SdAtwCFAnqRG4AOgFEBETgenAUcAikvsqxqbb1ksaDzwI9AQmR8SCtjXLqokTgVnnyJMgziL5tr+LpMeAHYFWB7Ij4qRWtgcf3ITXfNt0kgTSZfmO6cqr5kTQ3drmOSyDnFcxSToY2JVkAvlPEfFe4ZGZWZdTdOJwYupa8lzF9AJwWTos1FR2X0Qc3cJutpmqaY7AzLq3PENM7wGHSmoA/ndEvEtyr4IZUB1rNJlVQrXdYZ3ngUHvRMQJwLPAo5IG08J9CWZmVh3y9CAEG59NPZvk6iJ/NTQzq3J5EsT5TS8i4mFJXwS+WlxIVu262xU9ZrWqpQcGfSoingP+WroUd+q+YsMyM7PO1lIP4izgG8DlGdsC+FwhEZmZWZdQNkFExDfS34d2XDjVx9d1m1l31dIQ03Et7RgRd1U+HDMz6ypaGmIa3cK2AJwgrEWejLauxj36tmlpiGlsRwZiZsXxB6NtjjyXuSLpb4E9gN5NZRExoaigzMyqSdk7rHt2cCBt1Oqd1JImAicA3yK5ae54YHDBcZmZWSfL04PYLyKGS5obERdJuhzPP5hZF9bWRS89BJctz1pMa9Lf70j6OMnifUOLC8nMzLqCPD2I+yRtC1wGPE1yBdO1hUZlZpaDl8cvVp4HBl2cvrxT0n1A74ho9ZGjZmbl+IO9e8jzwKCewN8CQ5rqSyIi/n+xoZmZWWfKM8R0L7AWmAdsaMvBJR0B/BzoCVwbEZc0234OcEpJLLsBO0bESklLgVXA+8D6iKhvy7nNzKx98iSIARExvK0HTnseVwKfBxqBpyRNi4iFTXUi4jKSuQ0kjQa+ExGlt98eGhGvtfXcZmbdQVd/GmOeq5jul/SFzTj2KGBRRCxOH1N6K3BsC/VPAm7ZjPOYmVkB8iSIJ4C7Ja2R9JakVZLeyrFff2BZyftGyjzLWtLWwBHAnSXFATwkabakcTnOZ2ZmFZRniOlyYF9gXkS05VnUyigrt/9o4LFmw0v7R8RySR8DfivpuYiY8aGTJMljHMCgQd3zweBmZl1RngTxPDC/jckBkh7DwJL3A4DlZeqeSLPhpYhYnv5+RdLdJENWH0oQETEJmARQX1/f1hg7jS/zM7OuLk+CeAl4RNL9wLqmwhyXuT4FDJM0FPgrSRI4uXklSf2Ag4FTS8q2AXpExKr09RcALw5oZtaB8iSIJenPlulPLhGxXtJ44EGSy1wnR8QCSWem2yemVb8MPBQRb5fsvhPJvEdTjFMj4oG85zaz6tLde9zdNf4WE0R6qWqfiDhncw4eEdOB6c3KJjZ7PwWY0qxsMTBic85pZmaV0eJVTBHxPrBPB8ViZmZdSJ4hpmckTQNuBzYOA/mZ1GZm1S1PgtgeWAF8rqTMz6Q2M6tyeVZz9bOpzapUd508tY6RZzXXAcAvgP1Jeg4zgX+KiMaCYzOzbsJPZKtOeZbauB6YBnycZKmMe9MyMzOrYnkSxI4RcX1ErE9/pgA7FhyXmZl1sjwJ4jVJp0rqmf6cSjJpbWZmVSxPgvga8BXgv0mW3RiTlpmZWRXLcxXTi8AxHRCLmZl1IWUThKTzW9gvIuLiAuIxM7MuoqUexNsZZdsAZwA7AE4QZmZVrGyCiIjLm15L6gv8EzCW5NGhl5fbz8zMqkNrq7luD/xf4BTgBmCfiHi9IwIzs+7Pd2p3by3NQVwGHEfytLa9ImJ1h0VlZmadrqUexFkkT5D7IXBe+vAeSJ41HRHx0YJjMzPrkmqlZ9TSHESeeyTMzKxKOQmYmVmmPM+DMDOzDvTkkpWZ5Q1Dt+/QONyDMDOzTIX2ICQdAfwc6AlcGxGXNNt+CHAPsCQtuisiJuTZ18ys1pTtWdQXc77CEoSknsCVwOeBRuApSdMiYmGzqo9GxNGbua+ZmRWkyCGmUcCiiFgcEe+S3IF9bAfsa2ZmFVBkgugPLCt535iWNbevpDmS7pe0Rxv3RdI4SbMkzXr11VcrEbeZmVHsHIQyyqLZ+6eBwRGxWtJRwG+AYTn3TQojJpHc7U19fX1mnfaY+uSLlT6kmVm3UGQPohEYWPJ+ALC8tEJEvNW0hEdETAd6SarLs6+ZmRWryATxFDBM0lBJWwInAtNKK0j6G6VreEgalcazIs++ZmZWrMKGmCJivaTxwIMkl6pOjogFks5Mt08keXzpNyWtB9YAJ0ZEAJn7FhVrJdTK2ixmVjsKvQ8iHTaa3qxsYsnrK4Ar8u5rZmYdx3dSm5lZJicIMzPL5ARhZmaZvJqrmVkZtX7xiXsQZmaWyQnCzMwyOUGYmVkmJwgzM8vkSeo2qvVJKzOrHe5BmJlZJicIMzPL5ARhZmaZnCDMzCyTE4SZmWVygjAzs0xOEGZmlskJwszMMjlBmJlZJicIMzPLVGiCkHSEpD9JWiTp3Iztp0iam/48LmlEybalkuZJekbSrCLjNDOzDytsLSZJPYErgc8DjcBTkqZFxMKSakuAgyPidUlHApOAhpLth0bEa0XFaGZm5RXZgxgFLIqIxRHxLnArcGxphYh4PCJeT98+AQwoMB4zM2uDIhNEf2BZyfvGtKycM4D7S94H8JCk2ZLGldtJ0jhJsyTNevXVV9sVsJmZfaDI5b6VURaZFaVDSRLEASXF+0fEckkfA34r6bmImPGhA0ZMIhmaor6+PvP4ZmbWdkX2IBqBgSXvBwDLm1eSNBy4Fjg2IlY0lUfE8vT3K8DdJENWZmbWQYpMEE8BwyQNlbQlcCIwrbSCpEHAXcBpEfHnkvJtJPVteg18AZhfYKxmZtZMYUNMEbFe0njgQaAnMDkiFkg6M90+ETgf2AH4pSSA9RFRD+wE3J2WbQFMjYgHiorVzMw+rNBHjkbEdGB6s7KJJa+/Dnw9Y7/FwIjm5WZm1nF8J7WZmWVygjAzs0xOEGZmlskJwszMMjlBmJlZJicIMzPLVOhlrt3J1Cdf7OwQzMy6FPcgzMwskxOEmZll8hBTGbu8eHtnh2Bm1qncgzAzs0xOEGZmlskJwszMMjlBmJlZJicIMzPL5ARhZmaZnCDMzCyTE4SZmWVygjAzs0xOEGZmlqnQBCHpCEl/krRI0rkZ2yXp39LtcyXtk3dfMzMrVmEJQlJP4ErgSGB34CRJuzerdiQwLP0ZB1zVhn3NzKxARfYgRgGLImJxRLwL3Aoc26zOscCNkXgC2FbSzjn3NTOzAhW5mmt/YFnJ+0agIUed/jn3BUDSOJLeB8BqSX/azHjrgNc2c9/uqtbaXGvtBbe5RpzdnjYPLrehyAShjLLIWSfPvklhxCRgUttC+zBJsyKivr3H6U5qrc211l5wm2tFUW0uMkE0AgNL3g8Aluess2WOfc3MrEBFzkE8BQyTNFTSlsCJwLRmdaYBf59ezfRZ4M2IeCnnvmZmVqDCehARsV7SeOBBoCcwOSIWSDoz3T4RmA4cBSwC3gHGtrRvUbGm2j1M1Q3VWptrrb3gNteKQtqsiMyhfTMzq3G+k9rMzDI5QZiZWaaaShDtWfqju8rR5lPSts6V9LikEZ0RZyXlXaZF0mckvS9pTEfGV4Q8bZZ0iKRnJC2Q9J8dHWOl5fi33U/SvZLmpG0e2xlxVoqkyZJekTS/zPbKf35FRE38kEx2vwB8guQy2jnA7s3qHAXcT3IfxmeBJzs77g5o837AdunrI2uhzSX1fkdyocSYzo67A/47bwssBAal7z/W2XF3QJt/APy/9PWOwEpgy86OvR1tPgjYB5hfZnvFP79qqQfRnqU/uqtW2xwRj0fE6+nbJ0juOenO8i7T8i3gTuCVjgyuIHnafDJwV0S8CBAR3b3dedocQF9JAvqQJIj1HRtm5UTEDJI2lFPxz69aShDllvVoa53upK3tOYPkG0h31mqbJfUHvgxM7MC4ipTnv/P/AraT9Iik2ZL+vsOiK0aeNl8B7EZyk+084J8iYkPHhNcpKv75VeSd1F1Ne5b+6K5yt0fSoSQJ4oBCIypenjb/K/C9iHg/+XLZ7eVp8xbAp4HDgK2A/5L0RET8uejgCpKnzV8EngE+B+wC/FbSoxHxVtHBdZKKf37VUoJoz9If3VWu9kgaDlwLHBkRKzootqLkaXM9cGuaHOqAoyStj4jfdEyIFZf33/ZrEfE28LakGcAIoLsmiDxtHgtcEskA/SJJS4BPAX/omBA7XMU/v2ppiKk9S390V622WdIg4C7gtG78bbJUq22OiKERMSQihgB3AP/QjZMD5Pu3fQ9woKQtJG1Nsjrysx0cZyXlafOLJD0mJO0E7Aos7tAoO1bFP79qpgcR7Vj6o7vK2ebzgR2AX6bfqNdHN14JM2ebq0qeNkfEs5IeAOYCG4BrIyLzcsnuIOd/54uBKZLmkQy/fC8iuu0y4JJuAQ4B6iQ1AhcAvaC4zy8vtWFmZplqaYjJzMzawAnCzMwyOUGYmVkmJwgzM8vkBGFmZpmcIKzbkTRB0uGbue+ZlVhmQtIUSUvS1VHnSDqsZNsj6Sqjz0h6VtK4km1LJc1Ltz0jab9WznOhpLPT1y22W9KXJO3e3raZNamZ+yCsekTE+VnlknpGxPut7FvJ+yDOiYg70mVKJgHDSradEhGzJG0PvCBpSrqoHMChm3M9frl2l/gScB/Jqq1m7eYehHUqSUMkPSfphnQN+zvSO32RdL6kpyTNlzQpXZWz6dv7mPT10rTeTOAUSbPT8hGSIr1THEkvSNq62Tfyf5S0MD3vrWnZNkrW3X9K0h8lZa0E29x/UX5RtD7A20CLiavZ3+S8tAfyHyR3/zaVl7b7kpLYf5r2RI4BLkt7JrtI+kbajjmS7iz5u05R8tyAxyUtVsnzMCR9N+3hzJF0SVq2i6QHlCzy96ikT6Xlx6f/beYoWbrDqox7EOc3w8EAAAOtSURBVNYV7AqcERGPSZoM/APwU+CKiJgAIOkm4Gjg3oz910bEAWm970n6KHAgMItkeYmZwCsR8Y42XZzvXGBoRKyTtG1adh7wu4j4Wlr2B0n/ka5hVM4RQPOlOm6WtI6kV/HtZj2b30t6H1gXEQ2lO0n6NMmyEXuT/P/5NDC7WZ3tSVaj/VREhKRtI+INSdOA+yLijrTeGxFxTfr6RySLMf4iPczOJAszfopkiYY7JB1J0gtpSP9W26d1JwFnRsTzkhqAX5IsgHc+8MWI+GvJ38+qiHsQ1hUsi4jH0te/4oMVZQ+V9GS6VMLngD3K7H9byevHgf1JHq7y4/T3gcCjGfvNJfkgP5UPnhPwBeBcSc8AjwC9gUFlznuZpMVpzD9utu2UiBie7nu2pMEl2w6NiJHNk0PqQODuiHgnXXW0+fpCAG8Ba4FrJR1HsqxClj3Tb/zzgFPY9O/3m4jYEBELgZ3SssOB6yPiHYCIWCmpD8lDpW5P/yZXkyQXgMdIlrL4BslyF1ZlnCCsK2i+3ktI6k3yTXVMROwFXEPyYZ2l9Nv9oyQfsoNJFqgbQZJwsoZA/ha4kmQZ7NmStiBZs+fv0g/wkRExKF3H6Pp06GZ6yf7nAJ8EfgjckNmwiFdJegFZyaCcFte/iYj1JA/MuZPkG/8DZapOAcanf7+L2PTvt67ktUp+Nz93D+CNkr/HyIjYLY3jTJK2DwSekbRDaw2z7sUJwrqCQZL2TV+fBMzkgw+z19JvsXmfGz0DOBV4Pn04zEqSBcweK60kqQcwMCJ+D3yX5JGcfUgWf/tWyXzH3gARMTb9cDyq9DjpOX4O9JD0xebBpOP+e5M8HjNv/F+WtJWkvsDojGP2AfpFxHTg28DIdNMqoG9J1b7AS5J6kfQgWvMQ8LWSuYrt017MEknHp2VS+txySbtExJPp5PlrbLrUtFUBz0FYV/As8FVJVwPPA1elY+DXkDwJbCnJ8s6tioil6Wd7U49hJjCg5LGqTXoCv5LUj+Sb88/ScfyLSR4oNDdNEktJ5j5aOmekY/zfJUkwkAxdrQE+AkyJiNllD7DpsZ6WdBvJg27+QvbQWF/gnrSXJeA7afmtwDWS/pEkof4z8GR6nHlsmjyyzv2ApJHALEnvkqwO+gOS5HKVpB+SrB56K8kzoC+TNCyN4eG0zKqIV3O1TiVpCMnE6p6dHIqZNeMhJjMzy+QehJmZZXIPwszMMjlBmJlZJicIMzPL5ARhZmaZnCDMzCzT/wAao4liRty4BwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy\n",
    "from scipy.spatial.distance import cdist\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "# Calculate RBF distances\n",
    "for partition in ['val', 'train']:\n",
    "    fig = plt.figure()\n",
    "    for dataset in ['gasdb', 'cathub']:\n",
    "        _splits = splits[dataset]\n",
    "        fps = _splits['fingerprints_%s' % partition]\n",
    "        scale = 55.  # This was tuned until there was a relatively flat distribution within [0, 1]\n",
    "        print('Calculation and plot timing for %s set of %s:' % (partition, dataset))\n",
    "        %time distances = scipy.exp(-cdist(fps, fps, 'sqeuclidean') / scale**2)\n",
    "\n",
    "        # Plot\n",
    "        %time ax = sns.distplot(distances.flatten(), kde=False, norm_hist=True, label=dataset)\n",
    "\n",
    "    # Format\n",
    "    _ = ax.set_xlabel('pairwise-RBF distances')\n",
    "    _ = ax.set_ylabel('Normalized count')\n",
    "    _ = ax.set_title('%s data' % partition)\n",
    "    _ = ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like GASdb might be more diverse than Cathub by a bit. Let's verify using another distance metric.\n",
    "# Cosine similarity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculation and plot timing for val set of gasdb:\n",
      "CPU times: user 411 ms, sys: 188 ms, total: 599 ms\n",
      "Wall time: 443 ms\n",
      "CPU times: user 9.3 s, sys: 216 ms, total: 9.52 s\n",
      "Wall time: 9.52 s\n",
      "Calculation and plot timing for val set of cathub:\n",
      "CPU times: user 182 ms, sys: 99.9 ms, total: 282 ms\n",
      "Wall time: 155 ms\n",
      "CPU times: user 2.28 s, sys: 99.9 ms, total: 2.38 s\n",
      "Wall time: 2.37 s\n",
      "Calculation and plot timing for train set of gasdb:\n",
      "CPU times: user 8.35 s, sys: 3.86 s, total: 12.2 s\n",
      "Wall time: 11.2 s\n",
      "CPU times: user 4min 12s, sys: 12.3 s, total: 4min 24s\n",
      "Wall time: 4min 24s\n",
      "Calculation and plot timing for train set of cathub:\n",
      "CPU times: user 3.31 s, sys: 1.74 s, total: 5.05 s\n",
      "Wall time: 4.53 s\n",
      "CPU times: user 1min 9s, sys: 2.21 s, total: 1min 11s\n",
      "Wall time: 1min 11s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5gV1Znv8e9PRFFBFFsNctfwGHUEJB1Q8a7JUeIlOnghapRkhnFOTMYcTWLGHMfLk0wejTEXMyLeUAPe70qiHhMG0YgCclcjCMEWBhQUQQW5vOePqm62Te3uoum9d9P9+zzPfnbtqlVV7y5xv73WqlpLEYGZmVl921U6ADMza5mcIMzMLJMThJmZZXKCMDOzTE4QZmaWyQnCzMwyOUGYbQVJx0iq2YLyEyT9UyljMmsuThBmLZSkhZJOqHQc1nY5QZiZWSYnCGvzJF0u6aF6634j6bfp8ghJr0taJeltSf+yBcf+qqQ3JK2UdBOggm37SfqzpOWS3pc0VtJu6bZ7gJ7Ak5JWS/pRuv5BSf+THm+ipIOa4RKYZXKCMIN7gaGSdgWQ1A44CxiXbl8GnAzsCowAbpQ0sLGDSqoCHgZ+ClQB84EhhUWA/wT2AQ4AegBXAUTE+cAi4JSI6BgR16X7/BHoC+wFTAPGNukbm+XgBGFtXkT8neTH9hvpquOATyLi5XT70xExPxL/DTwLHJnj0EOBuRHxUESsA34N/E/BeedFxHMRsTYi3gN+BRzdSKx3RMSqiFhLkkz6S+q8RV/YLCcnCLPEOGB4uvxNNtUekHSSpJclrZD0IckPf1WOY+4DvFP7IZKRMes+S9pL0n2S3pX0EfCHho4rqZ2kX0ian5ZfmG7KE4vZFnOCMEs8CBwjqTtwOmmCkLQjSTPRL4G9I2I3YDwFfQkNWELSbER6LBV+JmleCqBfROwKnFfvuPWHWv4mcBpwAtAZ6F176ByxmG0xJwgzIG3imQDcCSyIiNfTTTsAOwLvAeslnQR8LedhnwYOknSGpO2B7wNfKNjeCVgNfCipG/DDevsvBfatV34tsBzYGfh5zjjMmsQJwmyTcSR/ndc1L0XEKpIf9geAD0j+in8iz8Ei4n3gTOAXJD/qfYEXC4pcDQwEVpIkk0fqHeI/gZ9K+lDSZcDdwN+Bd4G5wMtb9vXMtow8YZCZmWVxDcLMzDI5QZiZWSYnCDMzy+QEYWZmmbavdADNqaqqKnr37l3pMMzMthlTp059PyL2zNrWqhJE7969mTJlSqXDMDPbZkj6e7FtbmIyM7NMThBmZpbJCcLMzDK1qj4IMzOAdevWUVNTw5o1ayodSovRoUMHunfvTvv27XPv4wRhZq1OTU0NnTp1onfv3iSD6LZtEcHy5cupqamhT58+ufcrWROTpB6S/pJO1ThH0r+l67tIek7SW+n77kX2P1HSm5LmSbq8VHGaWeuzZs0a9thjDyeHlCT22GOPLa5RlbIPYj1waUQcABwKfFfSgcDlwPMR0Rd4Pv38OemUj78HTgIOBIan+5qZ5eLk8HlNuR4lSxARsSQipqXLq4DXgW4kE57clRa7i03TPBYaBMyLiLcj4jPgvnQ/MzMrk7L0QUjqDRwCTCaZlWsJJElE0l4Zu3SjYGpGoAYYXOIwzayVGjd5UbMe75uDezbr8fI45phj+OUvf0l1dTUdO3Zk9erVJT9nyROEpI4kUzZeEhEf5azmZBXKnLhC0khgJEDPnlvxH23KnZuWq0c0/ThmZq1ESROEpPYkyWFsRNTOlrVUUte09tAVWJaxaw2fn7u3O7A46xwRMRoYDVBdXe3Zj8ysRbj22msZO3YsPXr0oKqqii9/+ct07tyZ0aNH89lnn/HFL36Re+65h5133pkHH3yQq6++mnbt2tG5c2cmTpzIp59+yogRI5g7dy4HHHAAn3766eeOf+mll/KXv/yF3Xffnfvuu48998wcTmmrlPIuJgG3A69HxK8KNj0BXJAuXwA8nrH7q0BfSX0k7QCcQ85pHs3MKm3KlCk8/PDDvPbaazzyyCN1Y8SdccYZvPrqq8yYMYMDDjiA22+/HYBrrrmGZ555hhkzZvDEE8lP3c0338zOO+/MzJkzueKKK5g6dWrd8T/++GMGDhzItGnTOProo7n66qtL8j1KeRfTEOB84DhJ09PXUJL5eb8q6S3gq+lnJO0jaTxARKwHLgaeIencfiAi5pQwVjOzZjNp0iROO+00dtppJzp16sQpp5wCwOzZsznyyCM5+OCDGTt2LHPmJD9rQ4YM4cILL+TWW29lw4YNAEycOJHzzjsPgH79+tGvX7+642+33XacffbZAJx33nlMmjSpJN+jZE1METGJ7L4EgOMzyi8GhhZ8Hg+ML010ZmalE5Hd2n3hhRfy2GOP0b9/f8aMGcOECRMAGDVqFJMnT+bpp59mwIABTJ8+Hch/a2qpbun1WExmZs3siCOO4Mknn2TNmjWsXr2ap59+GoBVq1bRtWtX1q1bx9ixY+vKz58/n8GDB3PNNddQVVXFO++8w1FHHVVXZvbs2cycObOu/MaNG3nooYcAGDduHEcccURJvoeH2jCzVq/ct6V+5Stf4dRTT6V///706tWL6upqOnfuzLXXXsvgwYPp1asXBx98MKtWrQLghz/8IW+99RYRwfHHH0///v3Zf//9GTFiBP369WPAgAEMGjSo7vi77LILc+bMqev4vv/++0vyPVSsKrQtqq6ujiZPGOTbXM1ajddff50DDjigojGsXr2ajh078sknn3DUUUcxevRoBg4cWNGYsq6LpKkRUZ1V3jUIM7MSGDlyJHPnzmXNmjVccMEFFU8OTeEEYWZWAuPGjat0CFvNndRmZpbJNYjmVtiXAe7PMLNtlmsQZmaWyQnCzMwyuYnJzFq/+k2/W6sZm44nTJjADjvswOGHHw4kT1uffPLJDBs2LPcxevfuzZQpU6iqqmq2uMA1CDOzipowYQIvvfRSpcPI5ARhZlYCd999N/369aN///6cf/75PPnkkwwePJhDDjmEE044gaVLl7Jw4UJGjRrFjTfeyIABA3jhhReAZKC+ww8/nH333bduSI0JEyZw8skn1x3/4osvZsyYMXWfr7/+egYNGsSgQYOYN29es3wHNzGZmTWzOXPm8LOf/YwXX3yRqqoqVqxYgSRefvllJHHbbbdx3XXXccMNN3DRRRfRsWNHLrvsMgBuv/12lixZwqRJk3jjjTc49dRTczU37brrrrzyyivcfffdXHLJJTz11FNb/T2cIMzMmtmf//xnhg0bVtcn0KVLF2bNmsXZZ5/NkiVL+Oyzz+jTp0/R/b/xjW+w3XbbceCBB7J06dJc5xw+fHjd+w9+8IOt/xK4icnMrNlFxGZDcH/ve9/j4osvZtasWdxyyy2sWbOm6P477rjj544FsP3227Nx48a69fX3Lzxfcw3/7QRhZtbMjj/+eB544AGWL18OwIoVK1i5ciXdunUD4K677qor26lTp7pRXRvSq1cv5s6dy9q1a1m5ciXPP//857bXjuh6//33c9hhhzXL93ATk5m1fmUe0eCggw7iiiuu4Oijj6Zdu3YccsghXHXVVZx55pl069aNQw89lAULFgBwyimnMGzYMB5//HF+97vfFT1mjx49OOuss+jXrx99+/blkEMO+dz2tWvXMnjwYDZu3Mi9997bLN+jZMN9S7oDOBlYFhH/kK67H9g/LbIb8GFEDMjYdyGwCtgArC82FG19LWK4bw+1YVZxLWG475aoJQ33PQa4Cbi7dkVEnF0Q1A3Aygb2PzYi3i9ZdGZm1qBSzkk9UVLvrG1KelDOAo4r1fnNzGzrVKqT+khgaUS8VWR7AM9KmippZBnjMrNWojXNltkcmnI9KpUghgMN9aIMiYiBwEnAdyUdVaygpJGSpkia8t577zV3nGa2DerQoQPLly93kkhFBMuXL6dDhw5btF/Z72KStD1wBvDlYmUiYnH6vkzSo8AgYGKRsqOB0ZB0Ujd7wGa2zenevTs1NTX4j8ZNOnToQPfu3bdon0rc5noC8EZE1GRtlLQLsF1ErEqXvwZcU84AzWzb1r59+wafVLZ8StbEJOle4K/A/pJqJH0n3XQO9ZqXJO0jaXz6cW9gkqQZwCvA0xHxp1LFaWZm2Up5F9PwIusvzFi3GBiaLr8N9C9VXGZmlo+H2jAzs0xOEGZmlskJwszMMjlBmJlZJicIMzPL5ARhZmaZnCDMzCyTE4SZmWVygjAzs0xOEGZmlskJwszMMjlBmJlZJicIMzPL5ARhZmaZnCDMzCyTE4SZmWVygjAzs0ylnHL0DknLJM0uWHeVpHclTU9fQ4vse6KkNyXNk3R5qWI0M7PiSlmDGAOcmLH+xogYkL7G198oqR3we+Ak4EBguKQDSxinmZllKFmCiIiJwIom7DoImBcRb0fEZ8B9wGnNGpyZmTWqEn0QF0uamTZB7Z6xvRvwTsHnmnRdJkkjJU2RNOW9995r7ljNzNqscieIm4H9gAHAEuCGjDLKWBfFDhgRoyOiOiKq99xzz+aJ0szMypsgImJpRGyIiI3ArSTNSfXVAD0KPncHFpcjPjMz26SsCUJS14KPpwOzM4q9CvSV1EfSDsA5wBPliM/MzDbZvlQHlnQvcAxQJakG+A/gGEkDSJqMFgL/kpbdB7gtIoZGxHpJFwPPAO2AOyJiTqniNDOzbI0mCEn3RMT5ja2rLyKGZ6y+vUjZxcDQgs/jgc1ugTUzs/LJ08R0UOGH9DmFL5cmHDMzaymK1iAk/QT4d2AnSR/VrgY+A0aXIbbWYcqdm5arR1QuDjOzLVS0BhER/xkRnYDrI2LX9NUpIvaIiJ+UMUYzM6uARvsgIuInkroBvQrLp09Km5lZK5Wnk/oXJLeazgU2pKsDcIIwM2vF8tzmejqwf0SsLXUwZmbWcuS5i+ltoH2pAzEzs5YlTw3iE2C6pOeBulpERHy/ZFGZmVnF5UkQT+ChLszM2pw8dzHdVY5AzMysZclzF9MCMobbjoh9SxKRmZm1CHmamKoLljsAZwJdShOOmZm1FI3exRQRywte70bEr4HjyhCbmZlVUJ4mpoEFH7cjqVF0KllEZmbWIuRpYiqcFnQ9yTwOZ5UkGjMzazHy3MV0bDkCMTOzlqXRPghJnSX9StKU9HWDpM7lCM7MzConz1AbdwCrSJqVzgI+Au5scA9A0h2SlkmaXbDueklvSJop6VFJuxXZd6GkWZKmS5qS76uYmVlzypMg9ouI/4iIt9PX1UCeZyDGACfWW/cc8A8R0Q/4G9DQvBLHRsSAiKhuoIyZmZVIngTxqaQjaj9IGgJ82thO6XwRK+qtezYi1qcfXwa6b0GsZmZWRnnuYvpX4K6CfocPgAub4dzfBu4vsi2AZyUFcEtEFJ3iVNJIYCRAz549myEsMzODfHcxTQf6S9o1/fxRI7s0StIVJLfMji1SZEhELJa0F/CcpDeKzWCXJo/RANXV1ZsNCWJmZk2T50G5nwPXRcSH6efdgUsj4qdNOaGkC4CTgeMjIvMHPSIWp+/LJD0KDKIlz2A3pdE+ezOzbU6ePoiTapMDQER8AAxtyskknQj8GDg1Ij4pUmYXSZ1ql4GvAbOzypqZWenkSRDtJO1Y+0HSTsCODZSvLXcv8Fdgf0k1kr4D3EQyTMdz6S2so9Ky+0gan+66NzBJ0gzgFeDpiPjTFn0rMzPbank6qf8APC/pTpLO428Djc4RERHDM1bfXqTsYtJaSUS8DfTPEZeZmZVQnk7q6yTNBE4ABFwbEc+UPDIzM6uoPDUI0iYeN/OYmbUhefogzMysDXKCMDOzTLmamKyZFD4vUT2icnGYmeVQNEFImkVy11KmdMA9MzNrpRqqQZycvn83fb8nfT8XyHzIzczMWo+iCSIi/g7J6K0RMaRg0+WSXgSuKXVwrZqbm8yshcvTSb1LveG+Dwd2KV1IZmbWEuTppP4OcEc63HcAK0mepjYzs1Ysz5PUU9k03LciYmXpwzIzs0prtIlJ0t6Sbgfuj4iVkg5MB94zM7NWLE8fxBjgGWCf9PPfgEtKFZCZmbUMeRJEVUQ8AGwESOeU3lDSqMzMrOLyJIiPJe1B+tCcpENJOqrNzKwVy3MX06XAE8B+6fMPewJnljQqMzOruFx3MUk6GtifZD6INyNiXckjMzOzispzF9N84J8iYk5EzI6IdZKeyrHfHZKWSZpdsK6LpOckvZW+715k3xMlvSlpnqTLt+gbmZlZs8jTB7EOOFbSnZJ2SNd1y7HfGODEeusuB56PiL7A8+nnz5HUDvg9cBJwIDBc0oE5zmdmZs0oTx/EJxFxtqQfAS9IOosGRnmtFRETJfWut/o04Jh0+S5gAvDjemUGAfPSuamRdF+639wcsZbNuMmL6pb3W7Sibnlwny6VCMfMrNnlSRCCurmpp5I8E9HUX8G9I2JJerwlkvbKKNMNeKfgcw0wuGhw0khgJEDPnj2bGJaZmdWXJ0FcWbsQEc9L+l/ABaULKUlI9TQ0L8VoYDRAdXV1ozWbLVVYUzAza0samjDoSxHxBvCupIH1NjfaSV3EUkld09pDV2BZRpkaoEfB5+7A4iaez8zMmqihGsSlwD8DN2RsC+C4JpzvCZLaxy/S98czyrwK9JXUB3gXOAf4ZhPOZWZmW6GhCYP+OX0/tikHlnQvSYd0laQa4D9IEsMD6WB/i0gfuJO0D3BbRAyNiPWSLibp62gH3BERc5oSg5mZNV1DTUxnNLRjRDzSyPbhRTYdn1F2MTC04PN4YHxDx2+pJi9YkbnedzeZ2bamoSamUxrYFkCDCcLMzLZtDTUxeaJkM7M2LM9trkj6OnAQ0KF2XURcU6qgzMys8vKMxTQKOBv4HskzCmcCvUocl5mZVViesZgOj4hvAR9ExNXAYXz+OQUzM2uF8jQxfZq+f5Lejroc6FO6kNq4KXduWq52N5CZVU6eBPGUpN2A64FpJHcw3VbSqFqhBm9/LUwKZmYtRJ4Jg65NFx9O54HoEBGectTMrJVrNEGk8zN8HehdW14SEfGr0oZmZmaVlKeJ6UlgDTAL2FjacMzMrKXIkyC6R0S/kkdiZmYtSp7bXP8o6Wslj8TMzFqUPDWIl4FHJW1HMj+1gIiIXUsamZmZVVSeBHEDycNxsyKi2WdsMzOzlilPE9NbwGwnBzOztiVPDWIJMEHSH4G1tSt9m6uZWeuWJ0EsSF87pC8zM2sDGkwQ6UNyHSPih811Qkn7A/cXrNoXuDIifl1Q5hiS+aoXpKse8fDiZmbl1WCCiIgNkgY25wkj4k1gANQloHeBRzOKvhARJzfnuVsiT1FqZi1Vniam6ZKeAB4EPq5d2dic1DkdD8yPiL83w7HMzKwZ5UkQXUiG+D6uYF1zzUl9DnBvkW2HSZoBLAYui4g5WYUkjQRGAvTs2bMZQjIzM8g3mmtJJiWQtANwKvCTjM3TgF4RsVrSUOAxoG+R+EYDowGqq6tLeivufoserFue3/PMUp7KzKzi8kw52l3So5KWSVoq6WFJ3Zvh3CcB0yJiaf0NEfFRRKxOl8cD7SVVNcM5zcwspzxNTHcC40jmogY4L1331a0893CKNC9J+gKwNCJC0iCSRLZ8K8+3TZm8YAXzNyzabP03B7sZzczKI8+T1HtGxJ0RsT59jQH23JqTStqZJME8UrDuIkkXpR+HAbPTPojfAuf4SW4zs/LKU4N4X9J5bPprfzhb+dd8RHwC7FFv3aiC5ZuAm7bmHGZmtnXyJIhvk/xY30hy99JL6TqrgHGTN292Ajc9mVnzy3MX0yKSu43MzKwNKZogJF3ZwH4REdeWIJ6WYcqdBR+Or1gYZmaV1FAN4uOMdbsA3yHpP2i9CcLMzIoniIi4oXZZUifg34ARwH0kkwi1aYUPzZmZtUaNjebaBfg/wLnAXcDAiPigHIGZmVllNdQHcT1wBskwFgfXPtlsZmZtQ0MPyl0K7AP8FFgs6aP0tUrSR+UJz8zMKqWhPog8T1mbmVkr5SRgZmaZ8jxJ3SYUm9ltPyp3t1LJhhcvfM6juiSjuZtZK+AE0UoUG4IDPAyHmTWNm5jMzCyTE4SZmWVygjAzs0zug9hGeD5sMys31yDMzCxTRRKEpIWSZkmaLmlKxnZJ+q2keZJmShpYiTjNzNqySjYxHRsR7xfZdhLQN30NBm5O383MrExaahPTacDdkXgZ2E1S10oHZWbWllSqBhHAs5ICuCUiRtfb3g14p+BzTbpuSf0DSRoJjATo2dMPhGWp/xDdfouSp8YH9+lSiXDMbBtRqQQxJCIWS9oLeE7SGxExsWC7MvaJrAOlyWU0QHV1dWaZ1qz+xEW+w8nMmktFEkRELE7fl0l6FBgEFCaIGqBHwefuwOLyRdiyeTY7MyuHsvdBSNolncIUSbsAXwNm1yv2BPCt9G6mQ4GVEbFZ85KZmZVOJWoQewOPSqo9/7iI+JOkiwAiYhQwHhgKzAM+IZkL28zMyqjsCSIi3gb6Z6wfVbAcwHfLGZfV4yHBzdo8D7XRhk1esIL5GzYfJnyz4cELkwU4YZi1ES31OQgzM6sw1yBamWKD+vnOJzPbUk4QtuXcP2HWJjhBtGKuNZjZ1nAfhJmZZXINwjYzbvKiuvGaCnnsJrO2xTUIMzPL5ARhZmaZ3MRkuU1esHmz0/wNizZ/sM7MWgXXIMzMLJMThJmZZXITk22V/RY9CO3Su5v80JxZq+IahJmZZXINoo0rNnbTlqjrvF5ww+eO5c5rs22baxBmZpbJCcLMzDKVvYlJUg/gbuALwEZgdET8pl6ZY4DHgQXpqkci4ppyxmlbb9zkzScjgowJicysRapEH8R64NKImCapEzBV0nMRMbdeuRci4uQKxGdmZlSgiSkilkTEtHR5FfA60K3ccZiZWcMq2gchqTdwCDA5Y/NhkmZI+qOkgxo4xkhJUyRNee+990oUqZlZ21Ox21wldQQeBi6JiI/qbZ4G9IqI1ZKGAo8BfbOOExGjgdEA1dXVUcKQrZm4b8Js21CRGoSk9iTJYWxEPFJ/e0R8FBGr0+XxQHtJVWUO08ysTavEXUwCbgdej4hfFSnzBWBpRISkQSSJbHkZw2yTPEWpmRWqRBPTEOB8YJak6em6fwd6AkTEKGAY8K+S1gOfAudEhJuPzMzKqOwJIiImAWqkzE3ATeWJyMzMsngsJiuJpozx5M5rs5bFQ22YmVkmJwgzM8vkBGFmZpmcIMzMLJM7qa3Fc+e1WWW4BmFmZplcg7AWIe9T3IW3zLpmYVZaThBWcs0x77WZlZ8ThLU6xWoW4NqF2ZZwgrCKqcTggG6WMsvPndRmZpbJNQgrKw8pbrbtcIIww01PZlmcIGyb4juizMrHCcKsAa5ZWFvmBGGtQp6aRXPWPpw4rC2oSIKQdCLwG6AdcFtE/KLedqXbhwKfABdGxLSyB2otWrEO70o2QxUmjjxxOKFYS1b2BCGpHfB74KtADfCqpCciYm5BsZOAvulrMHBz+m5Wclt6p9XWJKGGHurL4oRi5VSJGsQgYF5EvA0g6T7gNKAwQZwG3B0RAbwsaTdJXSNiSfnDtW1ZnlpGuc9XzoTSnJyc2p5KJIhuwDsFn2vYvHaQVaYbsFmCkDQSGJl+XC3pzSbGVQW838R9WwtfgyZdg8u28BRbWr7sMq/BuRUIpILa0v8LvYptqESCUMa6aEKZZGXEaGD0VgclTYmI6q09zrbM18DXAHwNwNegViWG2qgBehR87g4sbkIZMzMroUokiFeBvpL6SNoBOAd4ol6ZJ4BvKXEosNL9D2Zm5VX2JqaIWC/pYuAZkttc74iIOZIuSrePAsaT3OI6j+Q21xFlCG2rm6laAV8DXwPwNQBfAwCU3ChkZmb2eR7u28zMMjlBmJlZpjaVICSdKOlNSfMkXZ6xXZJ+m26fKWlgJeIspRzX4EuS/ippraQWf8N+U+S4Buem//1nSnpJUv9KxFlKOa7Baen3ny5piqQjKhFnKTV2DQrKfUXSBknDyhlfixARbeJF0iE+H9gX2AGYARxYr8xQ4I8kz2EcCkyudNwVuAZ7AV8BfgZcVumYK3QNDgd2T5dPaqP/DjqyqY+yH/BGpeMu9zUoKPdnkhtnhlU67nK/2lINom6Ij4j4DKgd4qNQ3RAfEfEysJukruUOtIQavQYRsSwiXgXWVSLAMshzDV6KiA/Sjy+TPIfTmuS5Bqsj/YUEdqHIg6rbsDy/BwDfAx4GlpUzuJaiLSWIYsN3bGmZbVlr/355bOk1+A5JrbI1yXUNJJ0u6Q3gaeDbZYqtXBq9BpK6AacDo8oYV4vSlhJEsw7xsY1q7d8vj9zXQNKxJAnixyWNqPxyXYOIeDQivgR8A7i25FGVV55r8GvgxxGxoQzxtEhtacIgD/HR+r9fHrmugaR+wG3ASRGxvEyxlcsW/TuIiImS9pNUFRGtZQC7PNegGrgvmZ6GKmCopPUR8Vh5Qqy8tlSD8BAf+a5Ba9foNZDUE3gEOD8i/laBGEstzzX4YjpxF+ndfDsArSlRNnoNIqJPRPSOiN7AQ8D/bkvJAdpQDSJa7hAfZZPnGkj6AjAF2BXYKOkSkrs7PqpY4M0o57+DK4E9gP9KfyPXRysa2TPnNfhHkj+W1gGfAmcXdFpv83JegzbPQ22YmVmmttTEZGZmW8AJwszMMjlBmJlZJicIMzPL5ARhZmaZnCCsxZJ0jaQTmrjvRZK+1dwxNYWk8ZJ2K8FxF0qqSpdfaqTsvzf3+a31822uts2R1K4tD39QS9JCoDrP082SVkdEx9JHZa2JaxBWFpJ6S3pD0l3pPAMPSdo53XalpFclzZY0uuAJ3jG1Y/Cnfy1fKWkScK6kqen6/pIiffoZSfMl7Szpqtr5LCR9X9Lc9Lz3pet2kXRHet7XJGWN5Fn7RPH/kzRD0rR0yAlJuj6Nd5aks9OyXSVNTOdQmC3pyILYq9Jr8LqkWyXNkfSspJ3SMvtJ+pOkqZJekPSljFj2SPd5TdItFIwnJGl1sRgk/QLYKV03Ni33WHquOZJGFh5H0s/S7/uypL3T9XtLejRdP0PS4en68yS9kh77Fknt0teYguvzgyb+s7FKq/R44361jRfQmza5w04AAAPdSURBVGQwtCHp5ztI55sAuhSUuwc4JV0eQzoGP7AQ+FFBuTkkT3tfTDJswrlAL+Cv6farCo6/GNgxXd4tff85cF7tOuBvwC4ZcU8GTk+XOwA7kzxl/BzJE7h7A4uArsClwBVp2XZAp4LYq9JrsB4YkK5/oCCG54G+6fJg4M8ZsfwWuDJd/np6PavSz6vT92IxrK53rC7p+07AbGCP9HMUXP/rgJ+my/cDlxQctzNwAPAk0D5d/1/At4AvA88VnGu3Sv/786tprzYz1Ia1CO9ExIvp8h+A7wO/BI6V9COSH98uJD/+T2bsf3/B8kvAEOAokh/7E0n+on4hY7+ZwFhJjwG1Y+l8DThVm2bN6wD0BF6v3UlSJ6BbRDwKEBFr0vVHAPdG0sy1VNJ/k0yy9Cpwh6T2wGMRMT0jlgUF66cCvSV1JJmk6MG08gSwY8a+RwFnpLE8LemDjDJ5YgD4vqTT0+UeQF+SsZY+A54qiO+r6fJxJD/+pN97paTzSZLBq2ncO5HMm/AksK+k35EMFf5skRishXMTk5VT/Q6vkNSB5C/PYRFxMHAryY91lo8Lll8AjiSpNTwO9AeOACZm7Pd14PckP2ZTJW1Pkkz+MSIGpK+eEfG6pDvT5pLxZA8JTbH1ETGR5Ef8XeAeZXeSry1Y3kAyHtp2wIcFsQyIiAOKnLvBTsM8MUg6BjgBOCwi+gOvsemar4uI2nPUxleMgLsKYt4/Iq6KZLKl/sAE4Lsko+LaNsgJwsqpp6TD0uXhwCQ2/TC9n/4lnXfe34nAecBbEbERWEEy0OKLhYUkbQf0iIi/AD8iaU7qSDJI2/cK+jsOAYiIEemP3dBIBiiskfSNtMyOSvpNJgJnp23te5L8IL8iqRewLCJuBW4Hcs1pnp5ngaQz0/NI2fNgTyRpSkPSScDu9Qs0EMO6tFYBSfPQBxHxSdrXcWiOMJ8H/jU9RztJu6brhknaK13fRVIvJXdWbRcRDwP/l5zXwVoeJwgrp9eBCyTNJGlKujkiPiSpNcwiaf55Nc+BImJhulhbY5hE8ld4/WaXdsAfJM0i+Uv5xvSc1wLtgZmSZlN8QpzzSZpjZpI0a30BeJSk2WoGyXzFP4qI/wGOAaZLeo2kn+I3eb5L6lzgO5JmkDSxZXWaXw0cJWkaSRPZoowyxWIYnX7XscCfgO3T73QtybSqjfk3kqbAWSRNTwdFxFzgp8Cz6bGeI+mL6QZMkDSdpB/pJzmOby2Qb3O1spDUG3gqIv6hwqGYWU6uQZiZWSbXIMzMLJNrEGZmlskJwszMMjlBmJlZJicIMzPL5ARhZmaZ/j8vB1nveMSB5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZRcVbn38e+PMISQEMigYqYOyMIwZMCGMIOA96IyCAYRCQKybha+FxRfEOXiZYp6XSCiV7zGQCAMYR6UoAi8YG4ICNKBzKAQiSEYISQYEkhCSJ73j3MSiqaH0911qrr6/D5r1epTp06d/eyu5Knd++yztyICMzMrji2qHYCZmVWWE7+ZWcE48ZuZFYwTv5lZwTjxm5kVjBO/mVnBOPGbNSJpgqT/LNO5LpV0SznOZVYuTvzWpUhaJOnIjpwjIs6KiPHliikrSZMlfb/S5VrxOPFboUjastoxmFWbE791GZJuBgYDUyWtlnSBpDpJIelMSYuBx9Jj75L0D0krJU2XtEfJeTa3vCUdJmmJpPMkvS5pqaQzWohhqKT/lbRK0iNAv0avN1mupHHAKcAFaexT0/3flbQwPd8CSceX97dmReTEb11GRJwKLAaOiYieEXFFycuHAsOAf02fPwjsCnwEeBaY0sKpPwb0BgYAZwK/kLRjM8feCswkSfjjgdMavd5kuRExMd2+Io39mPT4hcDBafmXAbdI2qmFWM1a5cRvRXFpRLwdEWsAIuL6iFgVEeuAS4ERkno38971wOURsT4ifgesBnZrfJCkwcA+wH9GxLqImA5MLT2mjeUSEXdFxN8jYmNE3AG8COzbxrqbfYATvxXFK5s2JHWT9KO0C+UtYFH6Ur8m3wnLI+K9kufvAD2bOO7jwJsR8XbJvr91oFwkfVXSLEn/lPRPYM+WjjfLwonfuprmppst3f8V4DjgSJIulLp0vzpY9lJgR0nblewb3IZyPxC7pCHAtcDZQN+I2AGYV4Y4reCc+K2reQ3YuZVjegHrgOVAD+CH5Sg4Iv4GNACXSdpa0kHAMSWHtFZu49i3I/kyWAaQXlTesxyxWrE58VtX81/A99KukfObOeYmki6YV4EFwFNlLP8rwGhgBXBJWlbWcicBu6ex/zoiFgBXAX8k+VLYC3iijLFaQckLsZiZFYtb/GZmBePEb2ZWME78ZmYF48RvZlYwNTFhVb9+/aKurq7aYZiZ1ZSZM2e+ERH9G++vicRfV1dHQ0NDtcMwM6spkv7W1H539ZiZFYwTv5lZwTjxm5kVTE308ZuZAaxfv54lS5awdu3aaofSqXTv3p2BAwey1VZbZTreid/MasaSJUvo1asXdXV1SJ6kFCAiWL58OUuWLGHo0KGZ3uOuHjOrGWvXrqVv375O+iUk0bdv3zb9FeTEb2Y1xUn/w9r6O3HiNzMrGPfxm1nNuvXpxWU931dGD279oBwcdthh/PjHP6a+vp6ePXuyevXqXMvLrcUv6XpJr0ua18Rr50sKSbWxdmjDDe8/zMxqXJ4t/snANXxwBSIkDQI+A5T3q9rMrELGjx/PlClTGDRoEP369eNTn/oUvXv3ZuLEibz77rt84hOf4Oabb6ZHjx7cddddXHbZZXTr1o3evXszffp01qxZwxlnnMGCBQsYNmwYa9as+cD5zzvvPP7whz+w4447cvvtt9O//4em2+mQ3Fr8ETGdZPm5xq4GLqD5RbHNzDqthoYG7rnnHp577jnuvffezfOInXDCCTzzzDPMnj2bYcOGMWnSJAAuv/xyHnroIWbPns39998PwC9/+Ut69OjBnDlzuOiii5g5c+bm87/99tvsvffePPvssxx66KFcdtllZa9DRS/uSjoWeDUiZmc4dpykBkkNy5Ytq0B0ZmatmzFjBscddxzbbrstvXr14phjjgFg3rx5HHzwwey1115MmTKF+fPnA3DggQdy+umnc+2117JhwwYApk+fztixYwEYPnw4w4cP33z+LbbYgpNOOgmAsWPHMmPGjLLXoWKJX1IP4CLg4izHR8TEiKiPiPpy/5ljZtZeza1Tfvrpp3PNNdcwd+5cLrnkks3j6idMmMD3v/99XnnlFUaOHMny5cuB7EMw8xi+WskW/y7AUGC2pEXAQOBZSR+rYAxmZh1y0EEHMXXqVNauXcvq1av57W9/C8CqVavYaaedWL9+PVOmTNl8/MKFCxk9ejSXX345/fr145VXXuGQQw7ZfMy8efOYM2fO5uM3btzI3XffDcCtt97KQQcdVPY6VGw4Z0TMBT6y6Xma/Osj4o1cCy4diVN/Rq5FmVllVWP45T777MOxxx7LiBEjGDJkCPX19fTu3Zvx48czevRohgwZwl577cWqVasA+Pa3v82LL75IRHDEEUcwYsQIdtttN8444wyGDx/OyJEj2XfffTeff7vttmP+/PmbLxjfcccdZa+DmvuzpcMnlm4DDgP6Aa8Bl0TEpJLXF5Ex8dfX10e7F2IpR+L3l4dZp/D8888zbNiwaofB6tWr6dmzJ++88w6HHHIIEydOZO+9965qTE39biTNjIj6xsfm1uKPiJNbeb0ur7LNzPI0btw4FixYwNq1aznttNOqnvTbynfutpVb/2aFd+utt1Y7hA7xXD1mZgXjxG9mVjBO/GZmBePEb2ZWML64a2a1q9wz5pZ5wMa0adPYeuutOeCAA4Dk7t6jjz6aMWPGZD5HXV0dDQ0N9OtXvsmM3eI3M8vJtGnTePLJJ6sdxoc48ZuZtdFNN93E8OHDGTFiBKeeeipTp05l9OjRjBo1iiOPPJLXXnuNRYsWMWHCBK6++mpGjhzJ448/DiQTtB1wwAHsvPPOm6dmmDZtGkcfffTm85999tlMnjx58/Mrr7ySfffdl3333ZeXXnqpw/G7qycPHutv1mXNnz+fH/zgBzzxxBP069ePFStWIImnnnoKSVx33XVcccUVXHXVVZx11ln07NmT888/H4BJkyaxdOlSZsyYwQsvvMCxxx6bqdtn++23509/+hM33XQT5557Lg888ECH6uDEb2bWBo899hhjxozZ3Ofep08f5s6dy0knncTSpUt59913GTp0aLPv/8IXvsAWW2zB7rvvzmuvvZapzJNPPnnzz29961sdroO7eszM2iAiPjRV8jnnnMPZZ5/N3Llz+dWvfrV5SuambLPNNh84F8CWW27Jxo0bN+9v/P7S8soxTbMTv5lZGxxxxBHceeedm+fVX7FiBStXrmTAgAEA3HjjjZuP7dWr1+ZZOlsyZMgQFixYwLp161i5ciWPPvroB17fNEPnHXfcwf7779/hOrirx8xqVxWuoe2xxx5cdNFFHHrooXTr1o1Ro0Zx6aWXcuKJJzJgwAD2228/Xn75ZQCOOeYYxowZw29+8xt+/vOfN3vOQYMG8aUvfYnhw4ez6667MmrUqA+8vm7dOkaPHs3GjRu57bbbOlyH3KZlLqdONS1zqebO54u7ZrnoLNMyd0ZtmZbZXT1mZgXjrp5yKfcdhGZmOXGL38xqSi10T1daW38nTvxmVjO6d+/O8uXLnfxLRATLly+ne/fumd/jrh4zqxkDBw5kyZIlLFu2rNqhdCrdu3dn4MCBmY934jezmrHVVlu1eFesZZNbV4+k6yW9Lmleyb4rJb0gaY6k+yTtkFf5ZmbWtDz7+CcDRzXa9wiwZ0QMB/4CXJhj+WZm1oTcEn9ETAdWNNr3cES8lz59CsjeKWVmZmVRzVE9XwMebO5FSeMkNUhq8IUcM7PyqcrFXUkXAe8BU5o7JiImAhMhmbKhQqG1jW/aMrMaVPHEL+k04GjgiPBgXDOziqto4pd0FPAd4NCIeKeSZZuZWSLP4Zy3AX8EdpO0RNKZwDVAL+ARSbMkTcirfDMza1puLf6IOLmJ3ZPyKs/MzLLxXD1mZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgVTlRW4OrNbn178oX27LF7B6KF9qhCNmVn5ucVvZlYwTvxmZgVT2K6eprp0zMyKwC1+M7OCabXFL+nmiDi1tX1dzS6L76p2CGZmucjS1bNH6RNJ3YBP5RNOdZQm+YWDT6xiJGZm+Wu2q0fShZJWAcMlvZU+VgGvA79p7cSSrpf0uqR5Jfv6SHpE0ovpzx3LUgszM8us2cQfEf8VEb2AKyNi+/TRKyL6RsSFGc49GTiq0b7vAo9GxK7Ao+lzMzOroFa7eiLiQkkDgCGlx0fE9FbeN11SXaPdxwGHpds3AtOA72SO1szMOizLxd0fAV8GFgAb0t0BtJj4m/HRiFgKEBFLJX2khXLHAeMABg8e3I6izMysKVku7h4P7BYR6/IOplRETAQmAtTX10clyzYz68qyjOP/K7BVmcp7TdJOAOnP18t03s6r4YbkYWbWSWRp8b8DzJL0KLC51R8R32hHefcDpwE/Sn+2Ojqos3j65RXNvuYJ3MyslmRJ/PenjzaRdBvJhdx+kpYAl5Ak/DslnQksBjxo3syswrKM6rmxPSeOiJObeemI9pzPzMzKI8uonpdJRvF8QETsnEtEZmaWqyxdPfUl291JumfcqW1mVqNaHdUTEctLHq9GxE+BwysQm5mZ5SBLV8/eJU+3IPkLoFduEZmZWa6ydPVcVbL9HrAI+FIu0ZiZWe6yjOr5dCUCMTOzymi1j19Sb0k/kdSQPq6S1LsSwZmZWfllmbLhemAVSffOl4C3AM9BYGZWo7L08e8SEV8seX6ZpFl5BWRmZvnK0uJfI+mgTU8kHQisyS8kMzPLU5YW/9eBG0v69d8ETs8tIjMzy1WWUT2zgBGStk+fv5V7VGZmlpsso3p+KGmHiHgrIt6StKOk71ciODMzK78sXT2fjYj/2PQkIt6U9Dnge/mFVVuam6vf8/SbWWeU5eJuN0nbbHoiaVtgmxaONzOzTixLi/8W4FFJN5BMz/w1oF1z9JuZWfVlubh7haQ5wJGAgPER8VDukZmZWS6ytPiJiN8Dv885FjMzq4AsffxmZtaFOPGbmRVMVRK/pG9Jmi9pnqTbJHWvRhxmZkXUbB+/pLk0scj6JhExvD0FShoAfAPYPSLWSLoT+DIwuT3nMzOztmnp4u7R6c9/T3/enP48BXinDOVuK2k90AP4ewfPZ2ZmGTWb+CPib5DMxhkRB5a89F1JTwCXt6fAiHhV0o+BxSSzfD4cEQ83Pk7SOGAcwODBg9tTlJmZNSFLH/92jaZlPgDYrr0FStoROA4YCnw8Pf/YxsdFxMSIqI+I+v79+7e3ODMzayTLOP4zgevTaZkDWEly9257HQm8HBHLACTdCxxAcoewmZnlLMuduzN5f1pmRcTKDpa5GNhPUg+Srp4jgIYOnrNsdll8V7VDMDPLVZZpmT8qaRJwR0SslLS7pDPbW2BEPA3cDTwLzE1jmNje85mZWdtk6eOfDDxE0h8P8Bfg3I4UGhGXRMQnI2LPiDg1ItZ15HxmZpZdlsTfLyLuBDYCRMR7wIZcozIzs9xkubj7tqS+pDdzSdqP5AKvtUXDDe9v159RvTjMrPCyJP7zgPuBXdLx+/2BE3ONyszMcpNpVI+kQ4HdSObj/3NErM89sjIpXRZx4YbFVYzEzKxzaDXxS1oIXBkRE0r2PRARR7fwNqOFtXhxt4+ZVU+Wi7vrgU9LukHS1um+ATnGZGZmOcqS+N+JiJOA54HHJQ2hhVk7zcysc8tycVewee3dmSRj+vvkGpWZmeUmS+K/eNNGRDwq6V+B0/ILyczM8tTSQiyfjIgXgFcl7d3o5QfyDcvMzPLSUov/PODfgKuaeC2Aw3OJqACaG2L6ldFed8DM8tfSQiz/lv78dOXCqRzPwmlmRdVSV88JLb0xIu4tfzhmZpa3lrp6jmnhtQCc+M3MalBLXT2+pdTMrAvKMpwTSZ8H9gC6b9oXEe1abN3MzKorywpcE4CTgHNIbuY6ERiSc1xmZpaTLFM2HBARXwXejIjLgP2BQfmGZWZmecmS+NekP9+R9HGSSduG5heSmZnlKUsf/wOSdgCuJFkgPYDrco3KzMxyk2UhlvHp5j2SHgC6R0SHll5Mv0iuA/Yk+SL5WkT8sSPn7ApufbrphWJ8R6+ZlVOWhVi6AZ8H6jYdL4mI+EkHyv0Z8PuIGJPO8d+jA+cyM7M2yNLVMxVYC8wFNna0QEnbA4cApwNExLvAux09bxaepsHMLFviHxgRw8tY5s7AMuAGSSOAmcA3I+Lt0oMkjQPGAQwe7K4OM7NyyTKq50FJ/1LGMrcE9gZ+GRGjgLeB7zY+KCImRkR9RNT379+/jMWbmRVblhb/U8B9krYgGcopICJi+3aWuQRYEhFPp8/vponEX0SlXVELB59YxUjMrCvLkvivIrlpa25EdHit3Yj4h6RXJO0WEX8GjgAWdPS8taq56w4f2D/6vApFY2ZFkCXxvwjMK0fSL3EOMCUd0fNXwBPCmZlVSJbEvxSYJulBYN2mnR0ZzhkRs4D69r7fzMzaL0vifzl9bJ0+zMyshrWY+NObt3pGxLcrFI+ZmeWsxeGcEbGBZOilmZl1EVm6emZJuh+4i2TMPeA1d83MalWWxN8HWA4cXrLPa+5WW8MN72/Xe1CUmWWXZXZOZxUzsy4ky9KLAyXdJ+l1Sa9JukfSwEoEZ2Zm5Zdlrp4bgPuBjwMDSGbrvKHFd5iZWaeVpY+/f0SUJvrJks7NKyD7sKYWaNll8QpGD+1ThWjMrNZlafG/IWmspG7pYyzJxV4zM6tBWVr8XwOuAa4mGc3zZLrPqsCLyZhZR2UZ1bMYOLYCsZiZWQU0m/glXdzC+6JkEXYzM6shLbX4325i33bAmUBfwInfzKwGNZv4I+KqTduSegHfJJk3/3aSxVnMzKwGtTY7Zx/g/wKnADcCe0fEm5UIzMzM8tFSH/+VwAnARGCviFhdsajMzCw3LY3jP4/kbt3vAX+X9Fb6WCXprcqEZ2Zm5dZSH3+Wm7usip5+eQUACzd88M7er4weXI1wzKxGOLmbmRVM1RJ/Ov3Dc5IeqFYMZmZFVM0W/zeB56tYvplZIVUl8afz+X8euK4a5ZuZFVmWSdry8FPgAqBXlcqvKZ6YzczKqeItfklHA69HxMxWjhsnqUFSw7JlyyoUnZlZ11eNrp4DgWMlLSKZ/uFwSbc0PigiJkZEfUTU9+/fv9Ixmpl1WRVP/BFxYUQMjIg64MvAYxExttJxmJkVVbX6+K3SGkpWz6w/o3pxmFnVVTXxR8Q0YFo1Y+gKSi/+Lhx8YhUjMbNa4BZ/F7RpcfbGo4G8OLuZgRN/Mbnbx6zQPFePmVnBOPGbmRWME7+ZWcE48ZuZFYwTv5lZwXhUT4FsWrGr1MINi71il1nBuMVvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF4+Gctnk2z8Y8zNOsa3KL38ysYJz4zcwKxom/4HZZfNeHFmwxs67Nid/MrGB8cbeLKWfr3Rd9zbomt/jNzAqm4olf0iBJf5D0vKT5kr5Z6RjMzIqsGl097wHnRcSzknoBMyU9EhELqhCLpUq7iBYOPrGKkZhZ3ire4o+IpRHxbLq9CngeGFDpOMzMiqqqF3cl1QGjgKebeG0cMA5g8GBfTOxMfNHXrLZV7eKupJ7APcC5EfFW49cjYmJE1EdEff/+/SsfoJlZF1WVFr+krUiS/pSIuLcaMVjz3N9v1rVVY1SPgEnA8xHxk0qXb2ZWdNXo6jkQOBU4XNKs9PG5KsRhZlZIFe/qiYgZgCpdrpmZJXznrplZwXiuHiub5oZ5god6mnUmbvGbmRWME7+ZWcE48ZuZFYz7+C2zjtzY5WkezDoPJ35rkZdlNOt6nPitqvyXgFnluY/fzKxgnPjNzArGXT3WKbkLyCw/TvzWLq1d9PV0zmadl7t6zMwKxonfzKxg3NVjuchrFS/3/Zt1nFv8ZmYF4xa/5a4Sa/j6LwGz7Jz4rUvzF4LZhznxW9Vk+Utg0zEeHmpWPk78VlGdZdI3/yVgRVaVxC/pKOBnQDfguoj4UTXisM7DXwhmlVPxxC+pG/AL4DPAEuAZSfdHxIJKx2K1o61fDG3pGsrS5dTSesJt5S8Rq7ZqtPj3BV6KiL8CSLodOA5w4reyae6LYlNib+71ao5AKhd/sVhrqpH4BwCvlDxfAoxufJCkccC49OlqSX9uZ3n9gDfa+d5aVKT6tqOu5+d0bEVkqu8pFQikQvxvueOGNLWzGolfTeyLD+2ImAhM7HBhUkNE1Hf0PLWiSPUtUl3B9e3KKl3Xaty5uwQYVPJ8IPD3KsRhZlZI1Uj8zwC7ShoqaWvgy8D9VYjDzKyQKt7VExHvSTobeIhkOOf1ETE/xyI73F1UY4pU3yLVFVzfrqyidVXEh7rXzcysC/PsnGZmBePEb2ZWMF0m8Us6StKfJb0k6btNvC5J/52+PkfS3tWIsxwy1PWTkv4oaZ2kTjcYva0y1PeU9DOdI+lJSSOqEWe5ZKjvcWldZ0lqkHRQNeIsh9bqWnLcPpI2SBpTyfjKLcNne5iklelnO0vSxbkEEhE1/yC5SLwQ2BnYGpgN7N7omM8BD5LcR7Af8HS1486xrh8B9gF+AJxf7ZgrUN8DgB3T7c/W6mfbhvr25P3rc8OBF6odd151LTnuMeB3wJhqx53zZ3sY8EDesXSVFv/maSAi4l1g0zQQpY4DborEU8AOknaqdKBl0GpdI+L1iHgGWF+NAMssS32fjIg306dPkdwbUquy1Hd1pFkC2I4mboCsEVn+3wKcA9wDvF7J4HKQtb656yqJv6lpIAa045ha0FXqkVVb63smyV92tSpTfSUdL+kF4LfA1yoUW7m1WldJA4DjgQkVjCsvWf8t7y9ptqQHJe2RRyBdJfFnmQYi01QRNaCr1COrzPWV9GmSxP+dXCPKV9YpTe6LiE8CXwDG5x5VPrLU9afAdyJiQwXiyVuW+j4LDImIEcDPgV/nEUhXSfxZpoHoKlNFdJV6ZJWpvpKGA9cBx0XE8grFloc2fb4RMR3YRVK/vAPLQZa61gO3S1oEjAH+R9IXKhNe2bVa34h4KyJWp9u/A7bK47PtKok/yzQQ9wNfTUf37AesjIillQ60DIo25UWr9ZU0GLgXODUi/lKFGMspS30/IUnp9t4kFwpr8cuu1bpGxNCIqIuIOuBu4P9ERC6t4ArI8tl+rOSz3ZckR5f9s+0SSy9GM9NASDorfX0CyYiAzwEvAe8AZ1Qr3o7IUldJHwMagO2BjZLOJRk98FbVAm+njJ/txUBfktYgwHtRo7M6ZqzvF0kaMeuBNcBJJRd7a0bGunYZGes7Bvi6pPdIPtsv5/HZesoGM7OC6SpdPWZmlpETv5lZwTjxm5kVjBO/mVnBOPGbmRWME79VnKTLJR3ZzveeJemr5Y6pPST9TtIOOZx30aabdiQ92cqx/1Hu8q3r83BO6zQkdesit+Z3SHqXan1EvJHh2NUR0TP/qKwrcYvfOkRSnaQXJN2YzhF/t6Qe6WsXS3pG0jxJE0vuSJy8aV71tHV7saQZwCmSZqb7R0iK9K5cJC2U1EPSpUrXGJD0DUkL0nJvT/dtJ+n6tNznJDU5+2F69+v/SyfDelbSLuld3Vem8c6VdFJ67E6Spqfzo8+TdHBJ7P3S38Hzkq6VNF/Sw5K2TY/ZRdLvJc2U9LikTzYRS9/0Pc9J+hUlc7pIWt1cDJJ+BGyb7puSHvfrtKz5ksaVnkfSD9L6PiXpo+n+j0q6L90/W9IB6f6xkv6UnvtXkrqlj8klv59vtfOfjVVbteeo9qO2H0AdyURTB6bPryddAwDoU3LczcAx6fZk0nnVgUXABSXHzSe54/hsklvcTwGGAH9MX7+05Px/B7ZJt3dIf/4QGLtpH/AXYLsm4n4aOD7d7g70ILkj9hGSuyo/CiwGdgLOAy5Kj+0G9CqJvV/6O3gPGJnuv7MkhkeBXdPt0cBjTcTy38DF6fbn099nv/T56vRnczGsbnSuPunPbYF5QN/0eZT8/q8Avpdu3wGcW3Le3sAwYCqwVbr/f4CvAp8CHikpa4dq//vzo32PLjFlg1XdKxHxRLp9C/AN4MfApyVdQJJU+5Ak9alNvP+Oku0ngQOBQ0iS+FEkLeDHm3jfHGCKpF/z/iyG/wIcq/dXHusODAae3/QmSb2AARFxH0BErE33HwTcFkl302uS/pdkQZtngOslbQX8OiJmNRHLyyX7ZwJ1knqSLBJzV/rHDsA2Tbz3EOCENJbfSnqziWOyxADwDUnHp9uDgF1J5np5F3igJL7PpNuHkyR10nqvlHQqSZJ/Jo17W5K58KcCO0v6Ocl00A83E4N1cu7qsXJofKEoJHUnaSmOiYi9gGtJknBT3i7Zfhw4mKSV/xtgBHAQML2J930e+AVJkpopaUuSL4kvRsTI9DE4Ip6XdEPabfE7mp4el+b2RzID5iHAq8DNavri8rqS7Q0k82BtAfyzJJaRETGsmbJbvNiWJQZJhwFHAvtHMq3vc7z/O18fEZvK2BRfcwTcWBLzbhFxaSSL3YwApgH/TjIbqtUgJ34rh8GS9k+3TwZm8H7CeSNt+WZdK3U6MBZ4MSI2AitIJtd7ovQgSVsAgyLiD8AFJN06PUkmwDqn5HrCKICIOCNNYp+LZLK6JUqn95W0jZLrEtOBk9K+7P4kifZPkoYAr0fEtcAkINN6zWk5L0s6MS1Hano94OkkXVpI+iywY+MDWohhffpXACTdNG9GxDvptYT9MoT5KPD1tIxukrZP942R9JF0fx9JQ5SMNNoiIu4B/pOMvwfrfJz4rRyeB06TNIekS+eXEfFPklb+XJJumGeynCgiFqWbm1r4M0hazY27P7oBt0iaS9KyvTotczywFTBH0jyaX6TkVJJukTkk3UsfA+4j6T6aTbLG6wUR8Q+SdVBnSXqO5DrAz7LUJXUKcKak2SRdXU1dbL4MOETSsyRdVYubOKa5GCamdZ0C/B7YMq3TeJJlKFvzTZIuubkkXUB7RMQC4HvAw+m5HiG51jEAmCZpFsl1mgsznN86IQ/ntA6RVEeyOPSeVQ7FzDJyi9/MrGDc4jczKxi3+M3MCiikT78AAAAdSURBVMaJ38ysYJz4zcwKxonfzKxgnPjNzArm/wO25ohACvZnowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import scipy\n",
    "from sklearn.metrics.pairwise import cosine_similarity\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "# Calculate cosine distances\n",
    "for partition in ['val', 'train']:\n",
    "    fig = plt.figure()\n",
    "    for dataset in ['gasdb', 'cathub']:\n",
    "        _splits = splits[dataset]\n",
    "        fps = _splits['fingerprints_%s' % partition]\n",
    "        print('Calculation and plot timing for %s set of %s:' % (partition, dataset))\n",
    "        %time distances = 1 - cosine_similarity(fps, fps)\n",
    "\n",
    "        # Plot\n",
    "        %time ax = sns.distplot(distances.flatten(), kde=False, norm_hist=True, label=dataset)\n",
    "\n",
    "    # Format\n",
    "    _ = ax.set_xlabel('pairwise-cosine distances')\n",
    "    _ = ax.set_ylabel('Normalized count')\n",
    "    _ = ax.set_title('%s data' % partition)\n",
    "    _ = ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like GASdb is a bit more diverse wrt cosine distances, but not by that much."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unique & Duplicate sites\n",
    "Let's calculate the fraction of \"unique calculations\" there are in both Cathub and GASdb. How we define \"unique\" is subjective. Let's start with identifying whatever is within the 0th and 1st coordination shell---i.e., the adsorbate and its binding atoms.\n",
    "\n",
    "And while we're at it, let's also calculate the fraction of sites in the validation set that show up in the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For gasdb dataset:\n",
      "    92% of training set is duplicated\n",
      "    84% of validation set is duplicated\n",
      "    96% of validation is also in training\n",
      "For cathub dataset:\n",
      "    99% of training set is duplicated\n",
      "    97% of validation set is duplicated\n",
      "    100% of validation is also in training\n"
     ]
    }
   ],
   "source": [
    "import json\n",
    "\n",
    "\n",
    "def get_site(doc):\n",
    "    ads = doc['adsorbate']\n",
    "\n",
    "    # Turn the coordination string into a dictionary whose\n",
    "    # keys are elements and whose values are counts\n",
    "    coord = {}\n",
    "    for el in doc['coordination'].split('-'):\n",
    "\n",
    "        # EAFP to instantiate or increment the value\n",
    "        try:\n",
    "            coord[el] += 1\n",
    "        except KeyError:\n",
    "            coord[el] = 1\n",
    "\n",
    "    # Convert the dictionary to a string so that we can throw it\n",
    "    # in a tuple to hash\n",
    "    coord_string = json.dumps(coord, sort_keys=True)\n",
    "    site = (ads, coord_string)\n",
    "    return site\n",
    "\n",
    "\n",
    "def calculate_percent_duplicates(docs):\n",
    "    sites = {get_site(doc) for doc in docs}\n",
    "    n_repeats = len(docs) - len(sites)\n",
    "    percent_duplicates = (n_repeats / len(docs)) * 100\n",
    "    return percent_duplicates\n",
    "\n",
    "\n",
    "# Parse the data\n",
    "for dataset in ['gasdb', 'cathub']:\n",
    "    print('For %s dataset:' % dataset)\n",
    "    docs_train = splits[dataset]['docs_train']\n",
    "    docs_val = splits[dataset]['docs_val']\n",
    "\n",
    "    # Report duplicates within sets\n",
    "    print('    %i%% of training set is duplicated'\n",
    "          % calculate_percent_duplicates(docs_train))\n",
    "    print('    %i%% of validation set is duplicated'\n",
    "          % calculate_percent_duplicates(docs_val))\n",
    "\n",
    "    # Calculate the overlap between validation and training\n",
    "    sites_train = {get_site(doc) for doc in docs_train}\n",
    "    n_duplicates = 0\n",
    "    for doc in docs_val:\n",
    "        site = get_site(doc)\n",
    "        if site in sites_train:\n",
    "            n_duplicates += 1\n",
    "    percent_duplicates = n_duplicates / len(docs_val) * 100\n",
    "\n",
    "    # Report\n",
    "    print('    %i%% of validation is also in training' % percent_duplicates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like there are a lot of repeat \"sites\" between the training and validation sets within both GASdb and CatHub. And Cathub seems to have a bit more repeats.\n",
    "\n",
    "Let's expand the definition of \"site\" to include the coordination of adsorbate neighbors', as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For gasdb dataset:\n",
      "    46% of training set is duplicated\n",
      "    20% of validation set is duplicated\n",
      "    60% of validation is also in training\n",
      "For cathub dataset:\n",
      "    99% of training set is duplicated\n",
      "    96% of validation set is duplicated\n",
      "    100% of validation is also in training\n"
     ]
    }
   ],
   "source": [
    "import json\n",
    "\n",
    "\n",
    "def get_site(doc):\n",
    "    ads = doc['adsorbate']\n",
    "\n",
    "    # Turn the coordination string into a dictionary whose\n",
    "    # keys are elements and whose values are counts\n",
    "    coord = {}\n",
    "    for el in doc['coordination'].split('-'):\n",
    "\n",
    "        # EAFP to instantiate or increment the value\n",
    "        try:\n",
    "            coord[el] += 1\n",
    "        except KeyError:\n",
    "            coord[el] = 1\n",
    "\n",
    "    # Convert the dictionary to a string so that we can throw it\n",
    "    # in a tuple to hash\n",
    "    coord_string = json.dumps(coord, sort_keys=True)\n",
    "\n",
    "    # Sort the neighbors coordinations so that they hash correctly\n",
    "    neighbor_coords = doc['neighborcoord']\n",
    "    neighbor_coords.sort()\n",
    "    neighbor_coords = tuple(neighbor_coords)\n",
    "    site = (ads, coord_string, neighbor_coords)\n",
    "    return site\n",
    "\n",
    "\n",
    "def calculate_percent_duplicates(docs):\n",
    "    sites = {get_site(doc) for doc in docs}\n",
    "    n_repeats = len(docs) - len(sites)\n",
    "    percent_duplicates = (n_repeats / len(docs)) * 100\n",
    "    return percent_duplicates\n",
    "\n",
    "\n",
    "# Parse the data\n",
    "for dataset in ['gasdb', 'cathub']:\n",
    "    print('For %s dataset:' % dataset)\n",
    "    docs_train = splits[dataset]['docs_train']\n",
    "    docs_val = splits[dataset]['docs_val']\n",
    "\n",
    "    # Report duplicates within sets\n",
    "    print('    %i%% of training set is duplicated'\n",
    "          % calculate_percent_duplicates(docs_train))\n",
    "    print('    %i%% of validation set is duplicated'\n",
    "          % calculate_percent_duplicates(docs_val))\n",
    "\n",
    "    # Calculate the overlap between validation and training\n",
    "    sites_train = {get_site(doc) for doc in docs_train}\n",
    "    n_duplicates = 0\n",
    "    for doc in docs_val:\n",
    "        site = get_site(doc)\n",
    "        if site in sites_train:\n",
    "            n_duplicates += 1\n",
    "    percent_duplicates = n_duplicates / len(docs_val) * 100\n",
    "\n",
    "    # Report\n",
    "    print('    %i%% of validation is also in training' % percent_duplicates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within the first two shells, every single item in the CatHub validation split already showed up in the training split. But only 60% of them happened in the GASdb dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ktran",
   "language": "python",
   "name": "ktran"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}