molecular-property-prediction/kernel-methods/fchl-in-one-notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FCHL in One Notebook\n",
    "This notebook is designed to explain the [\"FCHL\" machine learning technique](https://aip.scitation.org/doi/full/10.1063/1.5020710) from 2018.\n",
    "FCHL is a\n",
    "- [ ] kernel method\n",
    "- [ ] created with atomistic, scalable kernels\n",
    "- [ ] that use \"alchemical\" many-body similarity functions\n",
    "- [ ] to regress atomic properties with ridge regression.\n",
    "\n",
    "Our goal in this notebook is to explain what each of those bullets means.\n",
    "We will presume that you are familiar, in general, with regression methods and how to use the PyData stack (e.g., numpy, pandas, sklearn)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "from qml.fchl import get_local_kernels, get_local_symmetric_kernels\n",
    "from qml import Compound\n",
    "from sklearn.model_selection import GridSearchCV, ShuffleSplit, train_test_split\n",
    "from sklearn.metrics.pairwise import rbf_kernel\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from sklearn.base import BaseEstimator\n",
    "from sklearn.pipeline import Pipeline\n",
    "from scipy.stats import linregress\n",
    "from time import perf_counter\n",
    "from io import StringIO\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernel Methods\n",
    "The [\"kernel\" trick](https://en.wikipedia.org/wiki/Kernel_method) is make a supervised learning problem easier by creating new input features based on pairwise similarities between a new input point and each point from the training set. The kernel creates a new set of $N$ features where each feature is the similarity between a certain entry and each of $N$ points from the training set. The similiarity function used to compute the new features is called _the kernel_.\n",
    "\n",
    "Here, we will demonstrate how the kernel can drastically simplify the learning problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x = np.arange(-2, 2.1, 1)\n",
    "y = (-1 <= x) & (x <= 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAACsCAYAAABM8oFkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPdUlEQVR4nO3de3BW5Z0H8O8vEEpCbuYCCATCNXJJFkpAwVbigihMQooRBI2g1gnFCxeZHVvE1RWqHRjXluFSaHVAxEqGai2z7ACxSrKwosCEnboRtgMbyBAGAlFiArn+9o8gW8z1vZycPM/5fmaYIXnfnOd34Pd9z3nO874noqogIjuEuF0AEQUPA01kEQaayCIMNJFFGGgii3T35cnx8fGalJTkUCm+KSkpAQAMGjTI5UrIdl2x144dO1auqgnf/75PgU5KSsLRo0eDVxUR+UVESlr6Pk+5iSxibKBzc3ORm5vrdhnkASb1mk+n3F3JqVOn3C6BPMKkXjP2CE1EzTHQRBZhoIksYuwceuzYsW6XQB5hUq+JLx+fTEtLU65DE7lPRI6patr3v89TbiKLGBvonJwc5OTkuF0GeYBJvWbsHLq0tNTtEsgjTOo1Y4/QRNQcA01kEQaayCLGzqEnTZrkdgnkESb1GtehiQzk9zq0iOSKyFEROXrp0iVnqiOioGg30Kq6VVXTVDUtIaHZHU9ck52djezsbLfLIA8wqdeMnUNfvnzZ7RLII0zqNV7lJrIIA01kEQaayCLGzqGnTp3qdgnkESb1GtehiQzEz0MTeYCxgZ4xYwZmzJjhdhnkASb1mrFz6GvXrrldAnmESb1m7BGaiJpjoIkswkATWcTYOXRGRobbJZBHmNRrXIcmMhDXoYk8wNhAp6enIz093e0yyANM6jVjA01EzTHQRBZhoIkswkATWcTYdei5c+e6XQJ5hEm91u46tIjkAsgFgIEDB44vKSnpjLqIqA1+r0N31dv4VldXo7q62u0yyANM6jVjT7lnzpwJAPj000/dLYSsZ1Kv8aIYkUUYaCKLMNBEFmGgiSxi7EWxxx9/3O0SyCNM6jV+HprIQNZ9Hrq8vBzl5eVul0EeYFKvGXvK/dBDDwEwY22QzGZSrxl7hCai5hhoIosw0EQWYaCJLGLsRbHFixe7XQJ5hEm9xnVoIgO1tg4d1CN0cXExPv/8cwBAWloaRo8eHczN3+LcuXMAgMTEREe2X1tbi/z8fJSVlSEmJgbTp09HZGSkI2N1hrNnz6KgoAB1dXUYNWoUJk6cCBFxuyy/NDY24uDBgzhz5gzCw8MxdepUOPlZfad7rby8HPn5+aiurkZSUhLS09MREuLnbFhVO/xn/Pjx2pITJ07onePv1MiwKB0cMUIHR4zQyLAoTRubpsePH2/xZwI1ZcoUnTJlStC329jYqOvW/Ur79onRH90Vrwvn9dHp9ybobbeF6/PPP6vXr18P+phOOnv2rGbNuk9jbwvTOVl99LG5fXTYkGhNTRmq+/btc7s8n+3Y8Y4OGXy7jk2J1YXz+mjWzD4aHd1TFy6YqxUVFY6M6VSvVVRU6Py58zXsB+E6KHKoDumVrL0j+2q/Pv10+/btbf4sgKPaQkYDPkIXFRVhyo/T0b9qGCboNIRI0yvLYB2DCyfOIv3H9yL/kwOYMGFCoEN1iuXLn8Whg+/hwK4ojEr+wc3vl57vhSWr3seDs7/ER3/ej+7du/7lh9LSUvzo7jT8dJ5ix5v90Cu86f+msTESez+uwoLHHsSWrTuRlZXlcqUds379r/GbN1/G9vVRmJQWe/MM4/KVXvjntZ/g3vS7cLDgc0RFRblcafsqKysx+c67cb2kHhNqpqJHbVOvqSq+qbyMpYuX4/Lly1i+fLlP2w3oKreq4tF5OUisGoH+GHwzzAAQIiHohyQMrLoDjzz8KNSHubpbDh8+jI8+fBf73o+9JcwAMKBfKHZtiUFlxQls377dpQp9s+L5xVg4R7Hq+ZibYQaAkBBBxn0R+NO2WDz11GNG/ELz0tJSvPLyShzIi8XkCWG3TBfiYrthw+vRGD28HK+/vtrFKjvu1VdXo6qkBkNqUtBD/r/XRAQxEo/R1Xdi1cqXbp7ud1RAgf7iiy9QVlqGvjqw1ef0wQBcuViBQ4cOBTJUp9i08Q0892RPxER3a/Hx0FDBz5/riY0b1nZyZb4rKyvD/gMfY/mi1uf9E8f1xPjUnsjLy+vEyvyzZcsmzJ8diaTE0BYfFxG8uDQCb7+1FTU1NZ1cnW9qamrw+y2/Q/+aYa1exwiTXuiridi8abNP2w4o0IWFhYip693mxRURQfT1eBQWFgYyVKcoKCjArAfC23zO9PRwFH91GlVVVZ1UlX+OHDmCyROiEB3V8ovTd7LuFxQU7O+kqvz3H4X7Mev+tqc5ycN6ID62O06ePNlJVfnn1KlTCEUP9JK2L7JG1yQgf9/HPm273Yng927je8tjDQ0NQGP7g4gK6uvrfSqsPStWrAjq9oCm/Qlt518kJETQvXtI0Pcn2BoaGhDa8sHsFqGhQEN9nfMFBaihvgGhoe1flQ8N7fq91tDQgBBp+4UWAELge58FdBvflJQUXAu72u4g13pVIiUlxafC2pOZmYnMzMygbjM1dQwKPmt7PnnsxHXExcZ0+QsvY8aMwZFj36K2tu1rFwWfAWNSuv4FyzGp41H4WdsvPBfL61FSWo2hQ4cGdexg99qQIUNQWXsVNXq9zed90/0Kfpg2zqdtB3TKPX36dKCnokIvtV6UXkZ9aA0yMjICGaqZkydPBv3U6meLV2D97+tQX996CN7ceg2Lfraky6/hJicnY+TI0Xjvw8pWn1N6vg579lfiiSee7MTK/LN48VJsfbcald+2fkq48e1KZGdnIzo6OqhjB7vXoqKikJ2djfPdTrf6nHqtx8Ue57Bk2RKfth1QoLt164bf/m4z/iesCFf1SrPHK/VrnAw/jg2bNwR9mWfRokVYtGhRULeZkZGBuN5j8PjSr1FdfWvj1NcrXln3DYr+OwJPP/1sUMd1ytp1G/HC6irs/bj5fL/kXB1mPlqBF154EXFxcS5U55uUlBRkZD6E2U9U4EpFwy2PqSreeq8Sb/2hAS+9tCboYzvRa6t/+Sq+jrqI8/K/zVaA6rQWX4UfxezsLKSmpvq24ZYWp1v709obS/Ly8jQ6Mkb7RwzU4UjV4UjVARGDNCoiWnfu3BnI2nurnFrsr6qq0kfmz9a42DB97ql4Xf/LBP3Fkngd0C9C06dM1LKysqCP6aTCwkIdNLCPTvxhnP5qVZz+enWCzsnqrTEx4bp27eva2NjodokdVl9fr8uWPa0xMWH65CMJ+ps1Cfov/xSndwyP0ZQxw7S4uNiRcZ3qteLiYh0xdITGRSToUBmtIzBWk3qO0PCe4frM4me0rq6u1Z9FK28sCdp7uWtqarB7924cKjwEVWDS3Xdhzpw5CAsL8+0VpoPS09MBOHcXidOnT+Pdd9/BhQvnEB0dh7lz52PcON/mM11FQ0MD9u7di7/8ZR/q6moxcuQ/ICcnJ+inpp3lwoUL2LHjHZw5cxLh4VHIzPwJ7rnnHsemQU72mqqisLAQH37wJ3x7tRLD7xiOBQsWoG/fvm3+XGvv5Tb2wxlOB5roO12x16y7SSARNdf135DcilWrVrldAnmESb1mbKCnTZvmdgnkESb1mrGn3EVFRSgqKnK7DPIAk3rN2CP0smXLAHStCxVkJ5N6zdgjNBE1x0ATWYSBJrIIA01kEWMvir322mtul0AeYVKvGRvoyZMnu10CeYRJvWbsKffhw4dx+PBht8sgDzCp14w9Qq9cuRKAGWuDZDaTes3YIzQRNcdAE1mEgSaySLuBFpFcETkqIkcvXWr9ZoBE5D5j71jy3adfxo4d62odZL+u2Gud8utkO1NX+sclu5nUa8bOofPz85Gfn+92GeQBJvWasUfoNWua7r9s0t0kyEwm9ZqxR2giao6BJrIIA01kEQaayCLGXhTbsmWL2yWQR5jUa8YGOjk52e0SyCNM6jVjT7n37NmDPXv2uF0GeYBJvWbsEfqNN94AAGRmZrpcCdnOpF4z9ghNRM0x0EQWYaCJLMJAE1nE2ItiO3bscLsE8giTes3YQCcmJrpdAnmESb1m7Cn3rl27sGvXLrfLIA8wqdeMPUJv3rwZAPDwww+7XAnZzqReM/YITUTNMdBEFuFtfIks0m6gVXWrqqapalpCQkJn1EREfjL2otju3bvdLoE8wqReMzbQ8fHxbpdAHmFSrxl7UWzbtm3Ytm2b22WQB5jUaww0UTtM6jVjA01EzTHQRBZhoIkswkATWcTYZau9e/e6XQJ5hEm9Zmygw8PD3S6BPMKkXjP2lHvTpk3YtGmT22WQB5jUa8YGOi8vD3l5eW6XQR5gUq8ZG2giao6BJrIIA01kEQaayCKiqh1/ssglACUdeGo8gHJ/i+qCbNofm/YFsGt/fNmXQara7I4jPgW6o0TkqKqmBX3DLrFpf2zaF8Cu/QnGvvCUm8giDDSRRZwK9FaHtusWm/bHpn0B7NqfgPfFkTk0EbmDp9xEFmGgiSziWKBFZJ2IfCUi/yUiH4pIjFNjOU1E5ojIlyLSKCLGLpGIyAMiclJE/iYiP3e7nkCIyNsiclFE/up2LYESkUQR+UREim/02VJ/t+XkEfoAgDGqmgrgFIBfODiW0/4K4EEABW4X4i8R6QZgI4AZAEYBmC8io9ytKiDbADzgdhFBUg9ghaqOBHAXgGf8/b9xLNCqul9V6298+RmAAU6N5TRVLVbVk27XEaCJAP6mqqdVtRbA+wCyXK7Jb6paAOCK23UEg6qWqerxG3+vBFAMoL8/2+qsOfSTAP69k8ailvUHcO7vvi6Fn01DzhGRJADjABzx5+cDugWRiOQD6NvCQy+q6kc3nvMimk4pdgYyltM6si+Gkxa+xzXLLkREIgD8EcAyVb3qzzYCCrSqTmvrcRFZCCADwFTt4gve7e2LBUoBJP7d1wMAnHepFvoeEQlFU5h3quoH/m7HyavcDwB4AcAsVa12ahzqsC8ADBeRwSLSA8A8AH92uSYCICIC4C0Axar6r4Fsy8k59AYAkQAOiEiRiPzWwbEcJSKzRaQUwCQA/yYi+9yuyVc3LlA+C2Afmi665Knql+5W5T8R+QOA/wSQLCKlIvJTt2sKwN0AHgPwjzeyUiQiM/3ZEN/6SWQRvlOMyCIMNJFFGGgiizDQRBZhoIkswkATWYSBJrIIA+0RNz5ve9+Nv68RkfVu10TBZ+zvhyafvQzgVRHpjaZP88xyuR5yAN8p5iEichBABID0G5+7JcvwlNsjRCQFwO0AahhmezHQHiAit6Pp8+hZAKpE5H6XSyKHMNCWE5FwAB+g6Z5VxQBWA3jF1aLIMZxDE1mER2giizDQRBZhoIkswkATWYSBJrIIA01kEQaayCL/B8vAPCKFsD0CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 252x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5, 2.5))\n",
    "\n",
    "ax.scatter(x, [0]*5, c=y, s=100, ec='k')\n",
    "ax.set_ylim(ax.get_ylim())\n",
    "for i in [-1.5, 1.5]:\n",
    "    ax.plot([i]*2, ax.get_ylim(), 'k--')\n",
    "\n",
    "ax.set_yticklabels([])\n",
    "ax.set_xlabel('$x$')\n",
    "fig.tight_layout()\n",
    "fig.savefig('figures/without-kernel.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how this problem is not linearly separable. You need 2 rules to separate the two classes.\n",
    "\n",
    "Now, we apply a kernel trick to create a new set of features. We use the [\"RBF Kernel\"](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.rbf_kernel.html)\n",
    "to generate new features. The RBF kernel uses a Gaussian function to create new features\n",
    "\n",
    "$f_i = k(x, x_i) = \\exp(-(x - x_i)^2/\\sigma))$\n",
    "\n",
    "Applying the kernel function to all of the data yields the so-called, \"kernel matrix\" where\n",
    "\n",
    "$K_{ij} = k(x_i, y_j)$\n",
    "\n",
    "We will use this term in the remainder of the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f = rbf_kernel(x[:, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAACsCAYAAABM8oFkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASWElEQVR4nO3deXRUZZrH8e9DWLKxCAFElk5QEB2FKEHAVkH2TZBNQAFBp7FxAw6MqGDbp0E9PS3TNDoIqDQiiCCLgoJCtJElIAQaWcQgQgNpGNkNSTAa88wfVUJCKkUFk1tVN8/nnJyk7vtW8QD1y626dd/niqpijHGHcsEuwBhTcizQxriIBdoYF7FAG+MiFmhjXKR8sAv4teLi4jQ+Pj7YZZQpOTk5/PDDD4gIsbGxlCtn+wWnbdu27aSq1rx0e9gHOj4+ntTU1GCXUSZs3LiRiRPGkH7kIIk3x5KdlcfutH/zwAODeeGFv1ClSpVgl1hmiMghX9vtV6sJyKpVq+h9b2eG9zvMv7bWYeW8KqxdVo0dybXIOrWUu9u2IiMjI9hllnmOBVpEZovIcRHZXcS4iMg0EdkvIjtF5FanajP+ZWdnM3ToAJbOvorB/apQsaJcGKt3TQVen1KNZk1OMHHiU0Gs0oCze+g5QBc/412BRt6vEcBrDtRkArBgwQJaNY/k9hZRPsdFhOfHxjJ//jyysrIcrs7k51igVXUdcNrPlF7AXPXYDFQTkTrOVGf8WbP6ffp28/9UqV+3Ajc0imbLli0OVWV8CaX30HWBI/lup3u3FSIiI0QkVURST5w44UhxZdmPOTnERMtl50VFlSMnJ8eBikxRQinQvp4xPleOqOosVU1S1aSaNQsduTclrHGTpnyx/We/c3Jy8ti5J5PGjRs7VJXxJZQCnQ7Uz3e7HnA0SLWYfH73u5G8tSiLc5l5Rc5Z8H4mzZo1o2HDhg5WZi4VSoFeDgz1Hu1uBXyvqseCXZSBa6+9ln79B9D3oTNknCu8p1636TxPT87iT5OmBKE6k59jJ5aIyAKgLRAnIunA80AFAFWdAawEugH7gWxguFO1mct75ZVZPPmk0Kj1uwztH0OLWyLIys5j8YfK9p0/8s6CZbRq1SrYZZZ5Eu4NDpKSktTOFHPOgQMHeOONGexL20mlipF07tqH/v37ExXl+yMtUzpEZJuqJhXaboE2JvwUFehQeg9tjPmVLNDGuIgF2hgXsUAb4yIWaGNcxAJtjItYoI1xEQu0MS5igTbGRSzQxriIBdoYF7FAG+MiFmhjXMQCbYyLOBpoEekiImne3ttP+xivKiIrRORLEdkjItbkwJhicLLRfgTwv3j6b98IDBKRGy+Z9hjwlao2w9PdZIqIVHSqRmPCnZN76NuA/ap6QFV/BN7F04s7PwUqi4gAsXj6eOc6WKMxYc3JQAfSd/tV4AY83T53AaNUtVCrSevLbYxvTgY6kL7bnYEdwDVAIvCqiBS6pKH15TbGNycDHUjf7eHAUu/lcPYDB4EmDtVnTNhzMtBbgUYikuA90DUQTy/u/A4D7QFEpDZwPXDAwRqNCWuO9eVW1VwReRz4BIgAZqvqHhH5vXd8BjAJmCMiu/C8RB+vqiedqtGYcOdYoAFUdSWehvr5t83I9/NRoJOTNRnjJnammDEuYoE2xkUs0Ma4iAXaGBexQBvjIhZoY1zEAm2Mi1igjXERC7QxLmKBNsZFLNDGuIgF2hgXsUAb4yIWaGNcJKTa+HrntBWRHd42vp87WZ8x4c6x9dD52vh2xNOOaKuILFfVr/LNqQZMB7qo6mERqeVUfca4Qai18b0fT0+xwwCqetzB+owJe6HWxrcxcJWIrBWRbSIy1NcDWRtfY3wLtTa+5YHmQHc8LX2fE5HGhe5kbXyN8cnJnmKBtPFNB06qahaQJSLrgGbAPmdKNCa8hVob3w+AO0WkvIhEAy2BvQ7WaExYC6k2vqq6V0Q+BnYCecAbqrrbqRqNCXeieunb2PCSlJSkqampwS7DGEeJyDZVTbp0u50pZoyLWKCNcRELtDEuYoE2xkWKHWgR6Sgir4tIovf2iBKvyhhzRa7kY6tH8VzHeaKIVMdzYXZjTAi47B5aRP7m/R7l3XRCVc+q6jg8V4psUYr1GWOKIZCX3O293zd4v3/0y4CqPg3MLemijDFXJpBAfywim4CrReQhIF1EIn8ZVNVXSq06Y0yxXPY9tKqOE5GGwFogAegJ/IeI/AjsVtUBpVuiMSZQAR0UU9UDItJBVS+sehKRWOCmUqvMGFNsAR/lzh9m7+1MYHOJV2SMuWJ2YokxLmKBNsZFLNDGuEjI9eX2zmshIj+LSD8n6zMm3IVUX+588/6Mp7NJmbB161befHM6hw5+Q3RMDD3uGcjAgQOJioq6/J1Nkb799ltmznyVPbu3EVEugjvbdOWhhx6mRo0awS6t1IRaX26AJ4AlgOt7cmdmZtKrZyf6921Pg+oreWzoQXp3+JL33hlPQvw1rF+/PtglhqW8vDzGjHmcVi2bkpc5n0cG7Wdo7zR2bn2Z665rwJw5fw92iaXGya6fvvpyt8w/QUTqAr2Bdvg5R9y7wmsEQIMGDUq8UCeoKgPu60n12F2kbaxNhQoXuxzf3wfWfJ5F3z7d+PSzFG6++eYgVhp+Jkx4ii0b3+HrDVdzVbWIC9v7dIe9+6LpMuhJqlatRu/evYNYZekItb7cU4HxqvqzvwdyQ1/uDRs2sP+b7bwxpVqBMP+iY5sYxj8WxYsvPBeE6sLX8ePHmTFjOktnVysQ5l/c0Lgis6dWZeKEMYR7Pz1fnAx0IH25k4B3ReRfQD9guojc60h1Dnt91jRGDq3kM8y/GD6oMqs+/oTTp087WFl4mzv3LXp3jaVmXNEvPtvdEQV5GaSkpDhYmTNCqi+3qiaoaryqxgOLgUdV9X0Ha3TM/v1fc2vTin7nVKsawW/qxXD48GGHqgp/33yzh+ZN/c8REW5tWpH9+/c7U5SDQqovt1O1hILIyEgys/L8zlFVzmX+RGRkpN955qLIyBgysy7/UvpcFq78d3X0c2hVXamqjVX1WlV9wbtthq8wq+owVV3sZH1O6tylLwuX5/qds3VHDkg0jRo1cqiq8NetW08WLs/z+/74zNmf+Twlg7vvvtvBypxhZ4oFycMP/ycrVmfyz10/+BzPzVX+8N9Z/H7kKCIiCh/cMb517NiRzOxo3luRVeScyVPP0aN7N2rVct/lxy3QQRIXF8esWW/RffBp5i/JICfn4svvnV/l0PPB01SMTmTMmLFBrDL8lCtXjoWLljNqYhZTXjvL9xkXPzA5nP4TI8efYfW6WKb+bWYQqyw9Fugg6tevHwsXfcRbSxOIb/F/tO3zPYntT9F9cCa/bfsEy97/mAoVKgS7zLBzyy238Pm6L0jd25KGLY/Rpvf3tOp2hls7naBilT5s2LjNtWeL2bWtQsTBgwc5dOgQMTExJCYmWpBLyHfffUdaWhrly5enadOmxMbGBrukElHUta2cPFPM+JGQkEBCQkKwy3Cd2rVrU7t27WCX4Rh7yW2Mi1igjXERC7QxLmKBNsZFLNDGuIgF2hgXsUAb4yIWaGNcxAJtjIuEVBtfEXlARHZ6v1JEpJmT9RkT7hwLdL42vl2BG4FBInLjJdMOAm1UtSkwCZjlVH3GuEFItfFV1RRVPeO9uRlP3zFjTIBCqo3vJR4GVl3uQdPS0mjbtm2Bbffddx+PPvoo2dnZdOvWrdB9hg0bxrBhwzh58iT9+hW+OMfIkSMZMGAAR44cYciQIYXGx44dyz333ENaWhqPPPJIofGJEyfSoUMHduzYwejRowuNv/jii9x+++2kpKTw7LPPFhqfOnUqiYmJJCcnM3ny5ELjM2fO5Prrr2fFihVMmTKl0Pjbb79N/fr1WbhwIa+99lqh8cWLFxMXF8ecOXOYM2dOofGVK1cSHR3N9OnTWbRoUaHxtWvXAvDyyy/z4YcfFhiLiopi1SrPf9ukSZP49NNPC4zXqFGDJUuWAPDMM8+wadOmAuP16tVj3rx5AIwePZodO3YUGG/cuDGzZnleuI0YMYJ9+wpcFJXExESmTp0KwODBg0lPTy8w3rp1a1566SUA+vbty6lTpwqMt2/fnuee83Ra7dq1K+fPny8w3qNHD8aNGwdQ6HkHzj33iuJkoANp4+uZKHI3nkDfUcT4hb7clSpVKqn6jAl7jq2HFpHWwB9VtbP39jMAqvrSJfOaAsuArpdek9oXt6yHNqY4iloPHVJtfEWkAbAUGBJImI0xBYVaG98/ADXwNNgHyPX1W8gY45u1IDImDIXCS25jTCmzQBvjIhZoY1zEAm2Mi1igjXERC7QxLmKBNsZFLNDGuIgF2hgXsUAb4yKuvVhdRkYGS5cu5ejRo1SpUoVevXpRv379YJdljE+qymeffca2bdsAaNmyJXfddRfeNQ3Fe6Bw/mrevLnml5ubq0/913iNiYrVBrENtWHEDRof1UijK8Voz+499cyZM2pMKElOTtYGdRtozdjamlC+iTYs30RrxNbUhAYNdd26dT7vA6Sqjzy4anGGqjL8weF8siSZ67KbESnRF+blai6HKu0lOqECX6RuJiYmJlglG3NBcnIyfXr24drzTanB1Rf2yKrKCY5xIHoXKz/+iDvvvLPA/crE4oyUlBSWL11Bk+wWBcIMUF7K0zDnJjIOZTFt2rQgVWjMRXl5eQwbMpzrzjcjTuoUeHktItSSa2iYfRPDhgwn0B2vqwL91ylTqXW+AeXF96EBEeHq8wlM++sr5OXlOVydMQWtWbOGn7JyqU7RF6SvyTV8fyqD9evXB/SYodaXW0Rkmnd8p4jcWpzH37xpM9XzavmdU1Wqc+7cOU6ePFnM6o0pWVu2bCEm6yq/B75EhKo/1mDr1q0BPWao9eXuCjTyfo0ACresNMYUKaT6cntvz/UeyNsMVBOROgH/Abfdxmk57ndOhp4mNjaWuLi4YpZvTMlKSkoiK/qs3/fHqkpGxVM0b948oMd0MtC++nLXvYI5iMgIEUkVkdQTJ05c2D5m3GhORB/hZ831WYCqcizqIE+Mepxy5Vx1+MCEoU6dOhERLZym6J3QSY5R+apY2rRpE9BjOvmsDqQvd0C9u1V1lqomqWpSzZo1L2y/44476HJPF76OTiVHfyhwn581l4OVviKmXiRPjnryCso3pmRFRETw97mz2R/1Jaf0uwJ7alXlpB7jQPRuZs+dHfAJJk4GOh3If6pWPeDoFcwpkogwd95bDBpxH9si/8G+6H9ygD18W2knWyKTadK2IRs3b6By5cpX/JcwpiR17tyZ95Yt4rvaB9hVOYUD5fbwbbk97IzdwOm6R1j+0Qc+r9BRFCcb7ZcH9gHtgX/j6dN9v6ruyTenO/A40A3PZXKmqept/h63qK6fZ8+eZfHixRdO/bz33nuJj48vsb+PMSUpLy+P1atXs337dsBz6me7du2K3DMXdWKJo2eKiUg3YCoX+3K/kL8vt3iqfxXoAmQDw1XVb49eETkBHLrMHx0HhMvnVFZr6XBbrb9R1ZqXbgz7Uz8DISKpvn6bhSKrtXSUlVrtUK8xLmKBNsZFykqgZwW7gGKwWktHmai1TLyHNqasKCt7aGPKBAu0MS5SJgItIn8Rka+9SzKXiUi1YNdUFBHpLyJ7RCRPRELyY5bLLYMNJSIyW0SOi8juYNdyOSJSX0T+ISJ7vc+BUcV9jDIRaGANcJOqNsVzttozQa7Hn91AH2BdsAvxJcBlsKFkDp4TlcJBLjBWVW8AWgGPFffftkwEWlVXq15YgrUZzzniIUlV96pqWrDr8COQZbAhQ1XXAaeDXUcgVPWYqm73/nwO2IuP1Yb+lIlAX+IhYFWwiwhjAS1xNb+OiMQDtwBfFOd+runLLSLJwNU+hiao6gfeORPwvKyZ72Rtlwqk1hAW0BJXc+VEJBZYAoxW1Yzi3Nc1gVbVDv7GReRBoAfQXoP84fvlag1xv2qJq/FPRCrgCfN8VV1a3PuXiZfcItIFGA/0VNXsYNcT5rYCjUQkQUQqAgOB5UGuyRW8qw3fBPaq6v9cyWOUiUDjWZJZGVgjIjtEZEawCyqKiPQWkXSgNfCRiHwS7Jry8x5cfBz4BM9Bm0X517SHGhFZAGwCrheRdBF5ONg1+fFbYAjQzvs83eFdchwwO/XTGBcpK3toY8oEC7QxLmKBNsZFLNDGuIgF2hgXsUAb4yIWaGNcxAJtis27Zrej9+fJIjIt2DUZD9ecy20c9TzwJxGphWdFUM8g12O87Ewxc0VE5HMgFmjrXbtrQoC95DbFJiI3A3WAHAtzaLFAm2IRkTp41pP3ArJEpHOQSzL5WKBNwEQkGliKp+/VXmAS8MegFmUKsPfQxriI7aGNcRELtDEuYoE2xkUs0Ma4iAXaGBexQBvjIhZoY1zk/wF9bsAS5ydxxAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 252x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5, 2.5))\n",
    "\n",
    "ax.scatter(x, f[:, 2], c=y, s=100, ec='k')\n",
    "ax.set_xlim(ax.get_xlim())\n",
    "ax.plot(ax.get_xlim(), [0.2]*2, 'k--')\n",
    "\n",
    "ax.set_ylabel('$f_2$')\n",
    "ax.set_xlabel('$x$')\n",
    "fig.tight_layout()\n",
    "fig.savefig('figures/with-kernel.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now my features are linear separable based on $f_2$. We've made the learning problem easier by transforming the data to a new space.\n",
    "\n",
    "FCHL operates based on this idea of \"apply a similiarty function to make learning easier.\" The secret is finding a function that maps \"materials\" into a similiarly easy space for ML to learn a simple function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Atomistic Kernel\n",
    "The first concept used by FCHL is the notion of a atomistics kernel, which leverages the fact that many molecular properties can be represented as sums over atomic properties.\n",
    "Energy is a good example due to most quantum effects being local: two atoms with similar _local_ environemnts should have similiar energies.\n",
    "An atomistic kernel definies the kernel (read: similarity) function between two molecules as a sum over kernels describe each atom:\n",
    "\n",
    "$K(x, y) = \\sum_i \\sum_j k(x_i, y_j)$\n",
    "\n",
    "Here, we show how that simplifies molecular machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider an (artificial) case where we have a single molecule where the atom's local environment can be described by a feature vector: `[1]`\n",
    "\n",
    "We assign it an energy of -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "molA = np.ones((1, 1))\n",
    "molA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume an energy model based on a kernel regression: $E(x) = \\sum_i K(x, x_i)$ where $x_i$ is a molecule from the training set and $\\alpha_i$ is a learned parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def scalable_kernel(x: np.ndarray, y: np.ndarray) -> float:\n",
    "    \"\"\"Compute a scalable kernel\n",
    "    \n",
    "    Args:\n",
    "        x: 2D representation for a molecule\n",
    "        y: 2D representation for a second molecule\n",
    "    \"\"\"\n",
    "    \n",
    "    k = rbf_kernel(x, y)\n",
    "    return k.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving for the only learned parameter, $\\alpha_0$: $E(molA) = \\alpha_0 \\times K(molA, molA)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha_0 = -1.00\n"
     ]
    }
   ],
   "source": [
    "alpha_0 = -1 / scalable_kernel(molA, molA)\n",
    "print(f'alpha_0 = {alpha_0:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now predict the properties of several molecules which have atoms with similar environemnts but with different numbers of atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAACsCAYAAABM8oFkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAL90lEQVR4nO3dbYxcZRnG8f9ly8sqgVVbIi3UxSigVGxhRbTBF1qFIBZjAI1RiyY2oOFNbbXUD/pBQIv4DqapEI0ICtYKKrRFxC8GyJYWSoUi4U22EtoPLSK1tPX2w5xtl+3uzJl2ps+ZZ65fMunMmZ1z7m32ynnOmXOfRxGBmeXhVakLMLPWcaDNMuJAm2XEgTbLiANtlpHxqQtoxoQJE6Kvry91GWbJrVq1alNETBy5vKMC3dfXx8DAQOoyzJKT9PRoyz3kNstIR+2h61m2epBFy9ezYfNWJvX2MO/0Y/no9MmpyzLbr7II9LLVgyxYupat23cCMLh5KwuWrgVwqK2rZDHkXrR8/a4wD9m6fSeLlq9PVJFZGlkEesPmrU0tN8tVFoGe1NvT1HKzXGUR6HmnH0vPAeNesazngHHMO/3YRBWZpZHFSbGhE18+y23dLotAQy3UDrB1uyyG3GZW40CbZcSBNsuIA22WEQfaLCMOtFlGkgRa0iJJj0p6SNLvJPWmqMMsN6n20CuBqRFxAvAYsCBRHQ0tWz3IjKvu5uiv/ZEZV93NstWDqUsyG1OSQEfEiojYUby8FzgyRR2NDLVlDm7eSrC7LdOhtqqqwjH054A7UhcxGrdlWqdp26Wfku4C3jDKWwsj4vfFzywEdgA31lnPXGAuwJQpU9pQ6djclmmdpm2BjohZ9d6XNAc4C5gZdSbYiojFwGKA/v7+/ToR16TeHgZHCa/bMq2qUp3lPgP4KjA7Il5KUUMZbsu0TpOq2+rHwEHASkkA90bEBYlqGZPbMq3TJAl0RLw5xXb3htsyrZNU4Sy3mbWIA22WEQfaLCMOtFlGHGizjDjQZhnJ5q6fncaT61k7ONAJeHI9axcPuRNwF5e1iwOdgLu4rF0c6AQ8uZ61iwOdgLu4rF18UiwBd3FZuzjQibiLy9rBQ26zjDjQZhlxoM0y4kCbZcSBNsuIA22WkaSBlvQVSSFpQso6zHKR7HtoSUcBHwSeSVVDjtyW2d1S7qG/B8wH9utsGDnz5HpWKtCSrpZ0fKs2Kmk2MBgRD7Zqnea2TCs/5H4UWCxpPHADcFNEbKn3gXqT1QGXAx8qs+GUk9V1GrdlWqk9dEQsiYgZwGeAPuAhSb+S9IE6n5kVEVNHPoAngKOBByU9RW1u6AckjRZ+ImJxRPRHRP/EiROb++26jNsyrfQxtKRxwHHFYxPwIPAlSTc3s8GIWBsRh0dEX0T0Ac8CJ0bEc82sx/bktkwrNeSWdA0wG/gzcEVE3F+89W1JPkCrCLdlWtlj6IeBr48x9evJ+1JAsZe2FnFbZncrG+g1wHHF1K9DtgBPNzo5Zmb7T9lAXwucCDwECJhaPH+9pAsiYkWb6jOzJpQ9KfYUML0423wSMJ3aMHwW8J021WZmTSob6OMiYt3Qi4j4O7WAP9Gessxsb5Qdcj8m6Tpg6CuqjxfLDgK2t6UyM2ta2T30HOBx4FLgMmoXh5xPLcxjXlxiZvtXwz10cUHJ7RExC/juKD/yYsursspxF1dnaBjoiNgp6SVJh/krqu7kyfU6R9lj6P8CayWtBP4ztDAiLm5LVVYp9bq4HOhqKRvoPxYP60Lu4uocpQIdET+X1ANMiQhfu91lJvX2MDhKeN3FVT1lb3DwEWqXf95ZvJ4m6bY21mUV4i6uzlH2a6tvUGvC2AwQEWuo9TRbF/jo9Mlc+bG3M7m3BwGTe3u48mNv9/FzBZU9ht4REVtGNGf4XmBdxF1cnaF0+6SkTwLjJL0FuBj4W/vKMrO9UXbIfRFwPLANuAl4gdpVY2ZWIWXPcr9E7eZ+C9tbjpnti7K3IDoG+Aq1GwTu+kxEnNaessxsb5Q9hr4F+CmwBNjZ4GfNLJFmznJf19ZKzGyflT0pdrukL0g6QtLrhh77smFJF0laL2mdJN/1xKwFyu6h5xT/zhu2LIA37c1Gixv0nw2cEBHbJB2+N+uxzue2zNYqe5a71VeFXQhcFRHbivU/3+L1WwdwW2br1R1yS5o/7Pm5I967Yh+2ewxwqqT7JP1V0jv3YV3WoTy5Xus1Oob+xLDnC0a8d0a9D0q6S9LDozzOpjYyeC1wCrVh/G804rrSYeuZK2lA0sDGjRsblGudxG2ZrddoyK0xno/2+hWKWxaNvlLpQmBpRARwv6T/AROAPRIbEYuBxQD9/f2+fjwjbstsvUZ76Bjj+Wivm7EMOA12XbRyILUJ8KyLuC2z9Rrtod8h6QVqe+Oe4jnF64P3YbvXA9dLehh4GZhT7K2ti3hyvdarG+iIGFfv/b0VES8Dn2rHuq2zuC2ztUrPD21m1edAm2XEgTbLiANtlhEH2iwjDrRZRsp2W5l1nG7s5HKgLUvd2snlIbdlqVs7uRxoy1K3dnI50JalsTq2cu/kcqAtS93ayeWTYpalbu3kcqAtW93YyeUht1lGHGizjDjQZhlxoM0y4kCbZcSBNstIkkBLmibpXklripvon5yiDrPcpPoe+jvANyPiDklnFq/fn6gWs6ZUuS0zVaADOLR4fhiwIVEdZk2peltmqmPoS4FFkv4JXM2e82aZVVLV2zLbtoeWdBfwhlHeWgjMBC6LiN9KOg/4GTDqXFiS5gJzAaZMmdKmas3KqXpbZtsC3WCyul8AlxQvbwGW1FmPJ6uzyqj6BHuphtwbgPcVz08D/pGoDrOmVL0tM9VJsc8DP5A0HvgvxZDarOqq3papTpr0sb+/PwYGBlKXYZacpFUR0T9yua8UM8uIA22WEQfaLCMOtFlGHGizjDjQZhnxXT/NKqIVXVwOtFkFtKqLy0NuswpoVReXA21WAa3q4nKgzSqgVZPrOdBmFdCqLi6fFDOrgFZ1cTnQZhXRisn1POQ2y0hH9UNL2gg8nWjzE4BNiba9N1xv+6Ws+Y0RMXHkwo4KdEqSBkZrKK8q19t+VazZQ26zjDjQZhlxoMtbnLqAJrne9qtczT6GNsuI99BmGXGgzTLiQNch6ShJf5H0iKR1ki5p/Kn0JI2TtFrSH1LXUoakXkm3Snq0+L9+d+qa6pF0WfH38LCkmyQdnLqmIQ50fTuAL0fEW4FTgC9Kelvimsq4BHgkdRFN+AFwZ0QcB7yDCtcuaTJwMdAfEVOBccAn0la1mwNdR0T8KyIeKJ7/m9ofWjXmPBmDpCOBD1NnAsAqkXQo8F5qM5ASES9HxOakRTU2HugppnJ6NRWa39yBLklSHzAduC9xKY18H5gP/C9xHWW9CdgI3FAcJiyR9JrURY0lIgapzWn+DPAvYEtErEhb1W4OdAmSDgF+C1waES+krmcsks4Cno+IValracJ44ETguoiYDvwH+FraksYm6bXA2cDRwCTgNZI+lbaq3RzoBiQdQC3MN0bE0tT1NDADmC3pKeBm4DRJv0xbUkPPAs9GxNDI51ZqAa+qWcCTEbExIrYDS4H3JK5pFwe6Dkmidmz3SERck7qeRiJiQUQcGRF91E7U3B0Rldl7jCYingP+KWno1hwzgb8nLKmRZ4BTJL26+PuYSYVO4vkGB/XNAD4NrJW0plh2eUT8KV1JWboIuFHSgcATwGcT1zOmiLhP0q3AA9S+BVlNhS4B9aWfZhnxkNssIw60WUYcaLOMONBmGXGgzTLiQNsrSFpYdBI9JGmNpHcVl2N2QlNK1/PXVrZL0bZ4DfD+iNgmaQJwYERUpvnA6vMe2oY7AtgUEdsAImJTRGyQdI+kfkmzi732GknrJT0JIOkkSX+VtErScklHJP0tupgDbcOtAI6S9JikayW9b/ibEXFbREyLiGnAg8DVxbXuPwLOiYiTgOuBb+3vwq3Gl37aLhHxoqSTgFOBDwC/lrRH55Ok+cDWiPiJpKnAVGBl7dJmxlFrK7QEHGh7hYjYCdwD3CNpLTBn+PuSZgLnUrspAYCAdRFR6dsGdQsPuW0XScdKesuwRdMYNpeYpDcC1wLnRcTWYvF6YOLQfcAkHSDp+P1Uso3gPbQNdwjwI0m91DqJHgfmUutRBjgfeD3wu2J4vSEizpR0DvBDSYdR+5v6PrBuv1ZugL+2MsuKh9xmGXGgzTLiQJtlxIE2y4gDbZYRB9osIw60WUb+D0ESegnuRddOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 252x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5, 2.5))\n",
    "\n",
    "n = np.arange(1, 10, 1)\n",
    "e = [alpha_0 * scalable_kernel(np.ones((i, 1)), molA) for i in n]\n",
    "\n",
    "ax.scatter(n, e)\n",
    "\n",
    "ax.set_xlabel('Size')\n",
    "ax.set_ylabel('Energy')\n",
    "fig.tight_layout()\n",
    "fig.savefig('figures/atomic-kernel.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scalable kernel captures the scaling between molecule size and energy by construction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alchemical Kernels\n",
    "The other ingredient of FCHL is how local environments are compared using a \"many-body atomistic kernel.\"\n",
    "The atomistic kernels are based on a simple concept: similar atoms have similar elements\n",
    "and neighboring atoms at similar positions with similar elements.\n",
    "This principle of \"similar environments\" is captured by representing the local environment as a series of \n",
    "functions that capture 1-body, 2-body and 3-body interactions in the environment.\n",
    "The functions are based on Gaussians in a coordinate space where some directions spatial (e.g., distances\n",
    "and angles between atoms) and others are chemical (e.g., position on periodic table). \n",
    "The similarity between atomic environments are then expressed as integrals over the difference between \n",
    "the two functions.\n",
    "\n",
    "\n",
    "For example, the 1-body representation can be written as\n",
    "\n",
    "$A_1(I, x_0, x_1) = \\exp(-\\frac{(P_I - x_0)^2}{2\\sigma_P^2}-\\frac{(G_I - x_1)^2}{2\\sigma_G^2})$\n",
    "\n",
    "where $P_I$ and $G_I$ are the period and group of atom $I$, and $x_0$ and $x_1$ are the inputs to the function.\n",
    "We then compute the difference between two atomic environments by integrating out $x_0$ and $x_1$\n",
    "\n",
    "$k_1(I, J) = \\iint (A_1(I, x_0, x_1) - A_2(J, x_0, x_1))^2 \\,dx_0 \\,dx_1\n",
    "= \\frac{1}{2} \\exp(-\\frac{(P_I - P_J)^2}{4\\sigma_P^2}-\\frac{(G_I - G_J)^2}{4\\sigma_G^2})$\n",
    "\n",
    "Choosing Gaussians as the basis of the representation function leads to these integrals having analytical solutions.\n",
    "The technique for the 2- and 3-body terms is similar, with integrals over not just the elemental similarity variables\n",
    "but also over spatial terms.\n",
    "\n",
    "<center>\n",
    "<img width=300px src=\"https://aip.scitation.org/na101/home/literatum/publisher/aip/journals/content/jcp/2018/jcp.2018.148.issue-24/1.5020710/20180320/images/medium/1.5020710.figures.online.f1.jpg\"/>\n",
    "</center>\n",
    "\n",
    "Here, we demonstrate how to compute such kernels using [QML](https://github.com/qmlcode/qml/tree/develop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Downsample to only a few thousand points. As we will demonstrate\n",
    "#  FCHL can be costly to run on large datasets\n",
    "data = pd.read_json('../datasets/qm9.json.gz', lines=True).sample(2000, random_state=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QML is primary designed for providing the low-level operations needed to build ML tools, \n",
    "but does offer [a scikit-learn-compatible interface](https://github.com/qmlcode/qml/blob/develop/qml/qmlearn/__init__.py). \n",
    "We will eschew using it for teaching purposes and, instead, build these tools ourselves.\n",
    "\n",
    "The first step of using FCHL from QML is to pre-compute the representation used to compute the FCHL kernel.\n",
    "For this, we must first parse the molecular structure from an XYZ into QML's internal [Compound object](http://www.qmlcode.org/qml.html#qml-compound-class).\n",
    "We then use the [`generate_fchl_representation`](http://www.qmlcode.org/qml.html#qml.Compound.generate_fchl_representation) method to compute a\n",
    "description of the molecule that is ready for using to compute the kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 835 ms, sys: 0 ns, total: 835 ms\n",
      "Wall time: 834 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# I use StringIO to turn the XYZ strings into the file object needed by Compound,\n",
    "#  `map` for parsimony, and set the `max_size` to 30 to capture all molecuels in QM9\n",
    "def make_representation(x):\n",
    "    # Step 1: Compute the representation\n",
    "    c = Compound(StringIO(x))\n",
    "    c.generate_fchl_representation(max_size=30)\n",
    "    \n",
    "    # Step 2: Return it\n",
    "    return c.representation\n",
    "data['fchl_rep'] = data['xyz'].apply(make_representation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FCHL method uses Kernel Ridge Regression (KRR) to learn atomic properties.\n",
    "\n",
    "### Theory: Kernel Ridge Regression\n",
    "Kernel ridge regression creates a model of the functional form:\n",
    "\n",
    "$f(x) = \\sum_i \\alpha_i K(x, x_i)$\n",
    "\n",
    "where the sum is over all points in the training set and $x_i$ is training point $i$.\n",
    "\n",
    "The coefficients of this equation, $\\alpha_i$, are learned by minimizing\n",
    "\n",
    "$\\sum_i (y_i - f(x_i))^2 + \\lambda \\left\\Vert{\\bf \\alpha}\\right\\Vert^2$\n",
    "\n",
    "where the squared error to the training data is $\\sum_i (y_i - f(x_i))^2$,\n",
    "the sum of squares of the weights is $\\left\\Vert{\\bf \\alpha}\\right\\Vert^2$ (i.e., $L_2$-norm), \n",
    "and $\\lambda$ controls the tradeoff between minimizing model error and\n",
    "model complexity (larger weights lead to a more complex model).\n",
    "\n",
    "This equation has a convenient [analytic solution](https://en.wikipedia.org/wiki/Tikhonov_regularization)\n",
    "based on the fact it involves solving a system of linear equations.\n",
    "The costly step in that equation is a matrix inversion, which can be solved efficiently \n",
    "using [Cholesky decomposition](https://en.wikipedia.org/wiki/Cholesky_decomposition) - a \n",
    "convenience due to symmetry properties of the kernel function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building a KRR Model using `scikit-learn`\n",
    "\n",
    "We implement the KRR model using a two-step process. First, we use the [FCHL kernel functions from QML](http://www.qmlcode.org/examples.html#generating-the-fchl-kernel) to compute the kernel matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class FCHLKernel(BaseEstimator):\n",
    "    \"\"\"Class for computing the kernel matrix using the qml utility functions\n",
    "    \n",
    "    The input `X` to all of the function is a list of FCHL representation vectors\n",
    "    \n",
    "    Follows the \"BaseEstimator\" API so that we can \n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self):\n",
    "        super(FCHLKernel, self).__init__()\n",
    "        self.train_points = None\n",
    "    \n",
    "    def fit(self, X, y=None):\n",
    "        # Store the training set\n",
    "        self.train_points = np.array(X)\n",
    "        return self\n",
    "        \n",
    "    def fit_transform(self, X, y=None):\n",
    "        self.fit(X)\n",
    "        # Uses the get_localget_local_symmetric_kernels to halve the\n",
    "        #  computational cost (as the matrix is symmetric)\n",
    "        return np.squeeze(get_local_symmetric_kernels(self.train_points))\n",
    "    \n",
    "    def transform(self, X, y=None):\n",
    "        return get_local_kernels(np.array(X), self.train_points)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fchl_kernel = FCHLKernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This class provides a method for fitting the FCHL kernel, which captures the points used as training entries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[408.41636274, 445.73433157, 309.88690898],\n",
       "       [445.73433157, 487.71323198, 337.78980536],\n",
       "       [309.88690898, 337.78980536, 236.99023028]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fchl_kernel.fit_transform(data.head(3)['fchl_rep'].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training set contains 3 entries\n"
     ]
    }
   ],
   "source": [
    "print(f'Training set contains {len(fchl_kernel.train_points)} entries')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can then compute the kernel matrix for new molecules on demand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed feature matrix of shape: (2000, 3)\n",
      "CPU times: user 24.7 s, sys: 240 ms, total: 24.9 s\n",
      "Wall time: 5.11 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "k = fchl_kernel.transform(data['fchl_rep'].tolist())\n",
    "print(f'Computed feature matrix of shape: {k.shape}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we build a [`Pipeline`](https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html) object to streamline fitting the KRR model.\n",
    "The Pipeline will automatically pipe the kernel matrix computed from the FCHL kernel as inputs into the KRR model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "model = Pipeline([\n",
    "    ('kernel', FCHLKernel()),\n",
    "    ('krr', KernelRidge(kernel='precomputed'))\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we demonstrate fitting it to a small subset of the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "train_data = data.sample(100, random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 19.8 s, sys: 62.1 ms, total: 19.9 s\n",
      "Wall time: 847 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Pipeline(steps=[('kernel', FCHLKernel()),\n",
       "                ('krr', KernelRidge(kernel='precomputed'))])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "model.fit(train_data['fchl_rep'].tolist(), train_data['u0_atom'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we can compute the kernel matrix and fit the KRR model in one step.\n",
    "\n",
    "Evaluating the model across the entire training set is also simple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12min 11s, sys: 199 ms, total: 12min 11s\n",
      "Wall time: 34.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "data['u0_atom_pred'] = model.predict(data['fchl_rep'].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE = 0.0308 Ha = 19.33 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "mae = mean_absolute_error(data['u0_atom'], data['u0_atom_pred'])\n",
    "print(f'MAE = {mae:.4f} Ha = {mae*627.50947406:.2f} kcal/mol')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We acheive an error of ~20 kcal/mol with only 100 training points. Not bad!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAAD0CAYAAACsLwv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoH0lEQVR4nO3dd3iUZdb48e+ZSSUESEihKh1FUEHsIrKAlNXFrqi7Ku+KBWGLuOriqmBZ38UVX9a26upvdcHeVxQiNuyAFaIBRJGADqEnkJBk5vz+mJk4SWbCJJlkSs7nunKZp84Z5HDfz/3cRVQVY0xicEQ7AGNM5FhCG5NALKGNSSCW0MYkEEtoYxJIUrQDaGk5OTnaq1evaIdhTERs2rSJiooKdu7cuVVVc+seT/iE7tWrFytWrIh2GMY021133cXVV1/NpZdeykMPPbQh2DlW5TYmDjz66KNcffXVnH322dx///0hz7OENibGqSoLFy7k5JNP5vHHH8fpdIY8N+Gr3MbEOxHhlVdewe12k5qa2uC5VkIbE6OWL1/O+PHj2bFjB2lpaWRkZOz3GiuhjYkhBYUulq0t4UDHDq695HQyMzMpLy8nKysrrOstoY2JEQWFLmY88RmlW3/EteBPZKY4KCgooFu3bmHfwxLamBixbG0JZTu34XrqBtyV5Zx9y7/p169fo+5hz9DGxIgR/XNJ8exDHE56njebs8ae0Oh7WAltTJT4n5dH9M9lZL8sxhycxwPTfsk7o49g5EFdGDsov9H3jImEFpGzgZuBg4GjVDVo1y4R+R4oBdxAtaoOb60YjYkk//NyeZWbpz/+nqwP/8FhA3vzwAMPNCmR/WKlyr0KOAN4N4xzR6nq4ZbMJp4tW1tCeZUbVQ/FL9/FR28v4bDDDmv2fWMioVX1a1UtinYcxrSWEf1zSUtysGPpQ+xZ/RYXTb+WK664otn3jYmEbgQFlojIShGZGu1gjGmqsYPyOWb3O5SufIXTfz2VR//vrxG5b6s9Q4vIG0CXIIdmqepLYd7meFXdLCJ5QIGIfKOq9arpvmSfCnDAAQc0OWZjWtJ5E08ivWoX9957LyISkXtKLM36KSJvAzNDNYrVOfdmoExV72zovOHDh6sNnzSxZNOmTXTv3r1Z9xCRlcHakeKmyi0iGSKS6f8dOBlvY5oxcWPRokX07duXl19+uUXuHxMJLSKni0gxcCzwqogs9u3vJiKLfKflA++JyBfAJ8Crqvp6dCI2pvGWLVvGmWeeyeDBgznppJNq9hcUurjxpVUUFLqa/RkxVeVuCVblNtHg7zSSmZZMaUUV3apdXDvlDLp27cqyZcvIzc2tOc//Pjo92cn8yUPDeg8dqsodEx1LjEkEgUn8yHvfUV7lBsC9dxc/PjKN7MxMlixZUpPM8PP7aIDyKjfL1pY0q2OJJbQxERBY0joF3AEVX2e7jnQ45mzOOv2Uem9dRvTP5ZkVxTUl9Ij+9eb9axRLaGMiILCkdSs4HULlnl24y7aTktuL/GPP4PRRQ+tdN3ZQPvMnD63p092c0hlipFHMmHg3on8u6cneub7Sk51cclQXPK/9lZ3P3US/rCSmnNA7ZLKOHZTPnEmDm53MYCW0MRERWNIefWAH5s2cwk/rC+l6xizW7ahm03vfcXjPThFJ2oZYCW1MhIwdlM9xfbKZeeVvWbp0KadMm01Sn6OAnxu8WpoltDERUlDo4qKZt/DNhwXkjZ3KCRPPqlUNb26DVzisym1MhCxbW0LakHHkpHUk/eARlFZURbTBKxyW0MY0UeCMIwUvPc3GjIGkpqYiB4+oKZHHDspvlUT2s4Q2JkyBCQzUvHf+5z8fYstr8+lw9Jnkj/4fRg3M5fyjD2zVRPazhDYmDIEdR55ZUUy/vAzKq9zs+eY9tr5+D2m9j6DTiAupdHvomd0uKskM1ihmTFjqdtFcvWk35d99xtZX7iSt+0F0P2sW4kxutcavUKyENiYMgV00nQ6h2u1mx9uPkNy5B6dfN59LRg1u1cavUGy0lTFh8j9Dby2rZPGqH6ks3UFasoP7Lh0NwMKPN1BSVklu+5QWf4a20VbGNFFBoYuFH3vXV+/iKOXJ++8lc8RFpHTIZurIvny+cSf3v7UOT8A176/bxr0XDGv10toS2pgGFBS6mLbgUyrdHtxlO/hpwZ/wVJSSfthEJKsrhZt38e7arbWSGaDS7Wn2UMimsIQ2Joi5i4t4o/AnKqrcVLo9eCrKcD39F9x7dtDl3FtJzuqK0yGUlFXi9tR/bE1xOqLSOGYJbUyAgkIX85eu4atNu2v2eaoq2PLsHKq2FZN31k2k9zgIVXB7lDU/lZLidFDp9uAQ6JHVjr65GfYe2phoC3zXHKhqWzGVWzeQc+pM0nsPrTV5QaXbw6iBufTMbhf1Fm6whDamRuC75kCpXfrR/bKHcaZn1uxLcgjVHiU92Rm10jgYS2hjfALfNasqO954gKSsbnQYPqlWMnvPzYmZUjmQ9RQzxmfsoHymnNAbp8CuZf+h9NNXcZduq3eev1SO1CwjkWQltDEBSiuq2PHJS+z68CnaH3oy2SddguJttb70xD6UVlTFXKkcyBLamABr3nuVHW8+RLsBx5E9bhrDemVTVlHFmEFduGbcwGiHt18xkdAiMhc4FagEvgUuUdWdQc4bD/wf4AQeVtU7WjNOk1jmLi7ixc+Kyc5IYcboAQB88s1G0noNJefUaxCHk8837sTtUX7Y3vCcYIFDK9t8X24RORl4U1WrReR/AVT12jrnOIE1wFigGFgOTFbVwobubX25TTBzFxdx71vrara1uhJHUgoKqMeNOJw4oFYPsN8c631urqupq180R0wvVqeqS1S12rf5EdAjyGlHAetUdb2qVgJPApNaK0aTWF78rLjm930/rWPTPy+lvHg1AOLwzgM2qHvHsOYEC7b6RbTERELXMQV4Lcj+7sDGgO1i3756RGSqiKwQkRUlJdH7wzWxqaDQxY+7KgCo2raRLU/fCA4nSR1rl6ojB+Qy5YTeDMxv3+C82nXn5G4T46HDWfBdRGYB1cCCYLcIsi/o84KqPgg8CN4qd5MCNglr/tI1eBSqd2/B9dSNIA7yz72FpMycWucVbt7FR+u3U17lbvAZOtKrXzRHqyW0qo5p6LiIXAScAozW4A/2xUDPgO0ewObIRWjagrmLi/hq027c5aW4nroRz749dDn/DpKza1f2/CVuuAvJtfZkgKHESiv3eOBaYKSq7g1x2nKgv4j0BjYB5wHnt1KIJk4Ftj5/vnEn9/kawhwp6aT1PISMQ0aRkt+n1jX+Sf6AmhI62lXpcMVEQgP3AKlAgYgAfKSql4tIN7yvpyb6WsCvAhbjfW31iKqujl7IJtYFtj4v/PgHqj2KVlfi2bcXZ0YnOo+fHvS6wEn+YqUqHa6YSGhV7Rdi/2ZgYsD2ImBRa8Vl4tvCjzfUVJmrPYp63JS8/Deqtm2k68XzcSSn1rumbkkcK1XpcMVEQhvTXHXnzF748QbeLvr5DYeqh22vzad87UdkjbksaDIDjDooL64SuC5LaBP3AqvWT36yEY8q1QGziKgqO978F3tWLaXj8efT4YhTQ95r8eqfKCh0xW1Sx+J7aGMaJbBjR6XbUyuZAco+f43SFS+RecSpdDx+coP3cns0qh1DmstKaBPXCgpdbNy+t2bCgWDaDTwe995ddDzuXHyNriE5hbhozQ7FEtrErVBTBvlVFK8mtesAnO060ilIyZzidFDt8eD/d8AhcPlJ/eK2ug1W5TZxLNSUQQDl61fieuLP7HxvYdDjGSlOBnZpT2ChPnJAblwMkWyIJbSJWyP65wbtD1xRXEjJC7eTktuLjsecFfTaPZVuin4qI8XpTQH/LCTxzqrcJm698Nmmep35K7d8R8mzs3F2yCHv7Nk4UjNCXh9rM3ZGgiW0iUtzFxex6Ksfa+1Tj5uSl/4XSUkn/5xbcGZ0avAeTofE1IydkWAJbeKGfzWL9mnJrNywo95xcTjJ/dWfwJlEUse8oPcQvEP0nAKXj+ybUMkMltAmTtSdYSSQu7yUvWs+JPOwk+sNtKhLgJN8gy8SLZnBEtrEiRc/LQ6631NZwZZnb6bS9S1pBwwmOatbg/fxUHvwRaJpdCu3iGT45vcyplUUFLpwle6rt1+rqyh54TYqf1xL7q/+FDKZhZ9nx4iXYZBNtd8SWkQceMceXwAcCewDUkWkBO/IpwdVdW2LRmnarIJCF9c+90W9XmDqcbP1v3+n4vvP6Dzhd7QbcFzIe/ir2PE0DLKpwqlyvwW8AVwPrFJVD4CIZAOjgDtE5AVV/U/LhWnaooaem/cVF7K36H2yRk2h/aFjQ97DIdQ8LydyIvvtdxpfEUlW1armnhMtNo1vfPEPg8xMSw6ZzH6VJd+Tktur3v4kB6QlO+mdk8GM0QMSMpFDTeO73xI6nESN1WQ28SWwb3aoIRS7l79EcucepPc5Imgyg7cH2aOXHNViccayRrVyi0gW0B9I8+9T1XcjHZRpmwL7ZgerN5Z+sYQdbz5ExiGjSO9zRNB7pDgdCdGFs6nCTmgR+S3wO7yzbX4OHAN8CPyiRSIzbYK/s8iYQV1qLeda196iD9i++B7Seg+j84QZQe81ML89M8cdlJBV7HA1poT+Hd5W7o9UdZSIHATMbpmwTFsQ2OhV5FrHxCFdyctMZcP22hO/ln//OSWv/I3UrgPIPe3PiDO53r1SnI42n8zQuISuUNUKEUFEUlX1GxGJ77FmJqreKPyp1nbdvtl+5es+JjmrO7ln3YQjJa3e8SQH3HvBsDafzNC4hC4WkU7Ai3in292BTXRvmmHMoC4UuUK3ZKsqIkLW6Knovj040toHPa9DerIls0/YPcVU9XRV3amqNwN/Af6FLRZnmuGacQPpmxt8eGP1LheuBddStX0TIhIymQEmH9V2G8HqCqen2D8IsYYUMAYI3kJhTBg2bNtTb597zw5cT92AZ+9u1N3wG9GJQ7rG/SwjkRROlTuwV8Zs4KZIB9GIBd+/B0oBN1Ad7MW6iR9XLviUak/tfZ59e3A9fRPusu3kn3tr0HfNiT5iqjnC6Vjyb//vIvL7wO0IKgCuD1jw/Xq8a10FM0pVt7ZADKaV+HuDLVldu1HMU1XBlmfnULX1B/LO/Aup3Q8Oev2Vo/pZqRxCY4dPtsjSrKq6JGDzIyD4RFAm7hUUupi24FMq3Z76B93VgJJzyh9DdhwxDYvFSQJDLfgO3n9QlojIShGZGuoGtuB77Fr48YZ6yazqQaurcKS1J//8O8g4+MQG71H3dZf5WTiNYqX8XDK3E5Hd/kOAqmqHcD4oAgu+AxyvqptFJA/vq7NvgnU9tQXf44eqsmPpQ1SWbCD/7NlIUv1OI3WNGRTsr5GB8J6hMyPxQRFY8N2/GiWqukVEXgCOAqwveYy7csGnfLCuhL55mfXmAtv1/hOUrnyFDkeeDs6G/zp2Sk/igmN62fNzA8IpoSVUgjXmnP1cv98F30UkA3Coaqnv95OBOU39TNM6rlzwaU0PsLrJvHvFy+x6fyEZQ8bQadSUBpepSU92Mvfsw61Vez/CeYZ+S0Smi8gBgTtFJEVEfiEi/wYuamYc9wCZeKvRn4vIA77P6CYi/vWg84H3ROQL4BPgVVV9vZmfa1rYB+uCt2HsKXyHHUsfJL3/MXQePz1kMqc4HYwamMv8yUMtmcMQTiv3eLwNVU+ISG9gJ97hk05gCTBPVT9vThDhLPiuquuBw5rzOab1dW6fys7y6nr7k/N6kzHoJDpPmIE46k9RN3FIV3LapyT8lEGRFs4zdAVwH3CfiCQDOUB5sI4fxgSau7iIb0tq9wSr3uXC2SGPlJwDyDl1ZtDrRg3M5b4LhrVGiAmnUe+hfTOTBB8SYwy1pxC6r84UQpWub/lp4fV0PPackGtOJcoaU9Fi83KbiJm7uIgH3l6HO0jzaNX2TbievhFHWgYZg0aGvMeUE3pbFbsZYrFjiYlDBYUuHnjn26DJXL17K66nbgAg/9xbSeoQel7s0gqbnq45mpXQIvL7CMVh4tz8pWtwe+pns3rcbHn2ZjwVZeSdPZvk7O4h75Hok+C3huZWuf8I3B2BOEwc8j8vby2r5KtNu4OeIw4nnU64AEdae1K71H6ZIcDg7h04cUAepRVV1qIdAc1N6NA9AUxCa3CQBaDVlez7cQ1pPQfTbsCxQc+ZMKSrtWZHWHOfoa2fdBsVbJCFn3rclLwyF9eTs6jaGXogxeJVP1JQ6GqpENuk/Sa0iJSKyG7ff2v9AKEfiExCKymrDLpfVdn2+j8oX/MhWaOmkNwp9EAKt3rn4jaR02qDM0xiyW2fUm+fqrLzrUfY89UbdDxuMh2GNzzlXIrTYY1gERbO4IyXGzquqr+KXDgmHsxdXMQXG3fW21/x/WfsXv4CmcNOoeMJ59c77hAIbAg/vl9nawSLsHAaxY4FioGFwMdYQ1ib1tCKkGm9hpJz2vW0G3BsvcEWDuCKk/rxyHvfUV7lth5hLSSchO4CjAUmA+cDrwJPqOrqlgzMxKYnPtlQb9/eNR+SlN2dlJwDyBh4fNDrOmUkc824gRzes1ObWKc5WvbbKKaqblV9XVUvwrue1TrgbRGZ3uLRmZhSUOhi+57aPbnK16+k5KX/Zee7jzV4rX/u7LGD8pkzabAlcwsJ6z20iKQCv8RbSvcC5gPPt1xYJhZNW7Cy1va+TV9T8uLtJOf0JGfi70NeN6R7B5tlpJWE0yj2b2Aw3on7ZqvqqhaPysSMgkIX85euYfWm3QS+da4s+Z4tz87G2T6b/HPmhFzZIj3ZyYzRA1onWBNWCf1rYA8wAJgR0NjRqEkCTfwpKHRxxX9WUh2kj/auD59GklLIP/dWnBlZ9Y7bBAXREc57aBuR1UYtW1sSNJkBcib+nurSbSR1rJ+s02wi/KixZDUhjeifiyPg7ZO7ooxtr83HU1GGJKWQnNU16HU2BDJ6LKFNUHMXF3HZYytqOoJ4KisoeXY2ZavfpLLk+5DXWe+v6LIZS0wto//+dr15wNRdRcmLf2Xf5iJyJl1LWs/BQa/tkJbE38+xqXajyUpoUyNoMnvcbP3vXVR8t5LscdNCdhwBb9dOS+boanJCi0hX3/tpkyDWlwRZq3nvLvZtLqLTSReTedi4Bq8/rp9VtaOtOSX048A3InJnc4MQkVtE5EvfJPtLRKRbiPPGi0iRiKwTkeua+7nmZwWFrnqD21WVpPbZdJvyDzoe3fCCoEccmGWTFcSAJj9Dq+oY8b6UHhSBOOaq6l8ARGQGcCNweeAJIuIE7sXbr7wYWC4iL6tqYQQ+v00rKHQx9bEVtfbt/uQFqrYXk33ylThSM4Je5wBG2sLrMaVJCS0ig4ADga8iMUhDVQMnpMog+EwoRwHrfCtoICJPApMAS+gmKih0cet/V7Nhe3mt/WVfFrDjrX/RbuAJDV7vAXpmt7NkjiFNrXLPxrsW1VRf19BmE5HbRGQjcAHeErqu7sDGgO1iQsyYYutD719BoYvLHl9RL5n3rvmQba//wzsU8pSrgy5T4381bbN0xp6mVrkLVPVp4OlwL9jf+tCqOguYJSLXA1cBN9W9RZBrQy07a+tD78eytSXU7QRWseFLSl7+Gyld+5N7+p9DrtU8wbp1xqymJvRxviVgtwFfq+pd+7tgf+tDB1iId8x13YQuBnoGbPcANod5T1PH++u21tun7ipS8nqRd9bNOFLSg1430WbqjGlNTehVqnqniCQBhzQ3CBHpr6prfZu/Ar4JctpyoL9vBcxNwHl4J1wwjXTm/R/Uet/sqarAkZxGep8jSOs9FJHgT2IP/Wa4lcgxrqnP0Kf4Jjjoo6pfRCCOO0RklYh8iXch999B7fWhVbUab1V8MfA18LTNmtJ4BYWuWguvV+/ewuaHr6Tsq6UAIZP5iAOzLJnjQFNL6HOBocDpItJPVS9tThCqemaI/TXrQ/u2FwGLgp1rwjN/6Zqa3917duJ66i949u0hJb93yGty26fw3BXHtUZ4ppmamtDTgIOBvcDtkQvHtKSCQhdrXKUAePbtZcszN+HevZW8c+eQktcn6DXpyU5uP+PQ1gzTNENTEzpLVc8WkRRgHt4ENzGq7vtmdVez5bk5VJZ8T94ZN5DWo34zyBEHZnFItw7Wkh1nmprQlSIyDPgKb0cQE6MKCl3eYZCBOx1O0vscQebhE0jve2TQ6wo37+bykX0tmeNMUxvFZgGj8b7rfSpy4ZhIKih0ce1zX9Qks6qH6l0uRISOx5zd4MLr5VVuW6YmDjU1oa9U1bmqegmwPpIBmcgoKHQx44nPaqbdVVV2LH2YzY/OoHr3/hPVeoHFp0ZVuUWkE95n5oNEpAL4AvgtcEnkQzPNsWxtCeVV7prtXR88SenKl8kcPglnZk6D146yARdxq1EJrao7gUtEZBxQAhwKPNcCcZlmKCh08eQnP9Rsl376X3a9t4CMwaPJ+sX/1Fumpi4bcBG/mtoodixwNN6ul59GLhzTXHXXnir/7jO2FzxAer+j6TxhRr2OI6MG5rK7orqms4lVteNbUxO6E/AR8Bjwh4hFY5ol2EJyaT0H03HEhXQ48vR6I6ccQk3VuqDQZWtOJYCmJvR237VbfL+bKCsodNVK5n0/riGpUxec6R3odNx5Qa/x+BZcHzsov+bHxLf9tnKLSL1eB6o6B3gA7xpXu1ogLtNICz/+eVXIStd6XE/ewLbX5jd4jQOsep1gwimhHweGAYjIb1X1YfD2sxaRGapaf2Y502ICq8bg7Zv9zY+lVPkGN1dt34Tr6RtxpLYje8zUkPcR4IpR/axUTjDhJHRgk+iVwMMB2+8CR0Q0IhOS/91yeZWbJz/ZSJXbU2uGh+rSrbie+guoh/xzbyGpQ17Iez1oQyETUjgJHfh3pu77DpvXuxUFvluudHvqHd++5H48FaXkT/4ryZ171jvuZ0MhE1c4Cd1FRC7G24mkbkLb9D6taET/XJ5ZUUx5lRsHUDelO4+7iqqdP5HapV+9azulJ7G30s2QHp1sKGQCCyehbwaG4+0N1kNEVuOdUeQboOEuRyaixg7KZ/7koSxbW0JmWjIPvbueffv2UfrpK94eYO2zcLavv7SrzTTSdoSznOyDgdsi0gNvD7EheJ+hTSvyv166csGnVFZVsfW/d7K36H2Sc3uT3ntovfMnDulqydyGNPo9tKoW452wz2YOiZLRf3+bdVvK2L74XvYWvU/WL34bNJkPzE63Cf3aGGvUijNzFxfxbckedr7z/yj7cgkdjz2XDkeeVu+8FKeDG05p9vyNJs7YcrJx5t8ffEf17i2UfraI9kN/SccRF9Y6npokHNc3x0ZLtVGW0HGg7hRCSR3y6HrR3SRlda03cmpAfqYlcxtmVe4YN3dxEZc+5l2yZk/R++xe+QoAydnd642ccgh8tWk3M574jIJCVzTCNVFmCR3DAkdPlX/3GVtfmcver99F3dW1zhuY355RA3Nrlrax6YPaLqtyx6jA0VP7NhdR8sJtJGf3IPesmxDnz//bUpwOZo47CICP1m+nvMptY5rbsJhIaBG5Be/SsB68QzIv9k2yX/e874FSwA1Uq+rw1oyzNflL2MqSDWx55iacGZ3IO2cOzrT2NeckOeDeC4bVPC/7O53YmOa2KyYSmjAWfA8wSlXrr7SWIEb//e1a607t2/Q1kpRC3rm3ktQ+u2Z/RoqD1XMm1LrWxjSbmEjoMBd8T3iByayqiAiZh48n4+ATcaS2q3Xu3edZhxFTX8w0ioWx4Dt4E32JiKwUkZCDfeN1wXd/MnsqynAtvI6KDV8C1EtmqL1GlTF+rZbQIvKGb4XJuj+TAFR1lqr2BBbgXWUymONVdRgwAZgmIicGO0lVH1TV4ao6PDc3fhqHkp2Cp6qCLc/NYd/mItRTHfLc7XsqWzEyEy9aLaFVdYyqDg7y81KdUxcCDa1GiapuAV4AjmrZqFvX/xx3AFtfvIN9xV+Tc8rVpPf2VquTHELf3NorDp02tEc0QjQxLiaq3CLSP2Az6ILvIpIhIpn+3/GuI72qdSJseYtX/chDt86kfP0KcidMI+PgEYD3f9BlI/uy9OqTmDaqHwPz2zNtVD+uGTcwugGbmBQTjWJ4F3wfiPe11QZ8Ldwi0g14WFUnAvnAC76ujknAQlV9PUrxRlRBoYurFqzEVVpJp5EX0e7Q8TXHPEDhZu88jNeMG2iJbBoUEwkdzoLvqroeOKw142otS7/aQJUKnX/5x2iHYuJcTFS527K7776bx/90Hu6yHYhIzWAL/5CLFKeD848+MHoBmrgSEyV0W/XYY4/xhz/8gTPOOIOKft1Z/dPemmODu3dg6AFZ1uvLNIoldCvzz6v91LPPs2bBbDJ6H86Uv8wjJSWVaQs+pdLtIcXpYMboAZbIptFENbE7ZQ0fPlxXrFgR7TAAbzJf8Z+VlP2wGteTs0jJ60X+ubfhSG3HQ7/xdku3vtgmHCKyMthYBiuhW9H8pWuo9ijJ2d3JGHg8WWOm1vQCW/jxBh695ChLZNMsltCt6MdNG1F3O5ztOpJz6sxoh2MSkLVyt4KCQhdH//lJvrj/90EXkEtyiLVkm4iwErqF+Bu/Vm/ezSdff89PC67FXbGn1gydB2an0ye3vc0BZiLGEroFBC4q59m3ly3P3Ix79xbyzplDSn7fmvNGDsxjzqTBUYzUJBqrcreAwEXlti6aR6XrW3ImXUdaz5+T16YJMi3BEroFjOifS3qyE4BOx00m59RraNfv54FhQ7p3YP7koVbNNhFnVe4WMPqgXC7svo19XQ9ja1lXvt1SSvu0ZMoqqhgzqIsNsDAtxhI6wlSVmTNnMm/ePJYsWcLYSWOjHZJpQ6zKHWG333478+bNY/r06YwZMyba4Zg2xhI6gu6//35uuOEGLrzwQu6+++56y9QY09IsoSNk/fr1TJ8+nVNOOYVHHnkEh8P+aE3rs2foCOnTpw+LFi1ixIgRJCcnRzsc00ZZQjfTBx98wDurfqA8fzAj+h9Genp6tEMybZgldDN8+eWXjJswkcrUTnS5eD7PrCi298smquxBr4m+/fZbTj75ZCQ5ldwzb0IcTlv10USdJXQTbN68mbFjx1JdXc28R58lM6crYN05TfRZlTtM/tFTI/rn8vHzj1BSUsKbb77JkUceyQF9XTbTiIkJltBhCBw99cyKYi45/ix+87eh7Mw4ALBVH03siKkqt4jMFBEVkZwQx8eLSJGIrBOR61orrmVrS9hbXsG21+az2/UD/3x3Pa9ugBlPfEZBoau1wjBmv2ImoUWkJzAW+CHEcSdwL96F6gYBk0VkUGvEdlyfbHYsuouyL5dQ9eMa3B7vxIrWCGZiTcwkNDAP+BOh14Y+ClinqutVtRJ4EpjU0kGpKi/cM5vSr5cx5uKZXHPVpTVDI60RzMSamHiGFpFfAZtU9YsG+j93BzYGbBcDR4e431RgKsABBxzQrNhmzZrFgw8+yPXXX8/tt98OwOE9O1kjmIlJrZbQIvIG0CXIoVnAn/GuJtngLYLsC1qaq+qDwIPgnZe7EWHWsm/fPt5++20uu+wybrvttpr91ghmYlWrJbSqBh1LKCJDgN6Av3TuAXwqIkep6k8BpxYDPQO2ewCbWyhcVJXU1FTeeOMNUlNTbeSUiQtRf4ZW1a9UNU9Ve6lqL7yJO6xOMgMsB/qLSG8RSQHOA15uiZief/55xo0bR2lpKe3atcPpdLbExxgTcVFP6IaISDcRWQSgqtXAVcBi4GvgaVVdHenPXLp0KZMnT6a0tNSGQJq4ExONYoF8pbT/95r1oX3bi4BFLfXZn3zyCZMmTWLAgAG8+uqrZGRktNRHGdMirAjyKSwsZMKECeTl5bFkyRKys7OjHZIxjWYJ7aOq9O7dm4KCArp27RrtcIxpkpircre2srIyMjIyOOSQQ1i+fLm1Zpu41qZL6F27djFixAiuu87bLdyS2cS7NpvQ5eXlnDB6PF9+tYqOfQ6PdjjGRESbTOiqqipGTZjEqpUfkz3xDzxenGWjpkxCaJMJfemll/LxOwVkn3wFGYNG2qgpkzDaZKPYaaedRlJ2D97POIHyKreNmjIJQ1SbPHYhLgwfPlxXrFgBwNq1a+nfv3/NscBphWywhYknIrJSVYfX3d9mqtz33HMPBx98MO+9917NvrGD8pkzabAls0kYbSKhFyxYwPTp0zn11FM55phjoh2OMS0m4avc/fv31++++44TTzyRRYsWkZaWFu2QjGm2UFXuhE9oh8Ohw4YN480336RDhw7RDseYiGizCS0iJcCGCNwqB9gagfvEC/u+se1AVa33aibhEzpSRGRFsH8RE5V93/jUJhrFjGkrLKGNSSCW0OF7MNoBtDL7vnHInqGNSSBWQhuTQCyhjUkgltCNtL8VMhOJiNwiIl+KyOciskREukU7ppYkInNF5Bvfd35BRDpFO6bGsoRuhP2tkJmA5qrqoap6OPBf4MYox9PSCoDBqnoosAa4PsrxNJoldOPsb4XMhKKquwM2M0jw762qS3wLOgB8hHe5pbjSJic4aIowV8hMOCJyG/AbYBcwKsrhtKYpwFPRDqKx7LVVgHBWyFTVXSLyPTBcVeOp729QDX1nVX0p4LzrgTRVvanVgmsB4XxfEZkFDAfO0DhLEEvoMPhWyFwK7PXt8q98WXeFzIQlIgcCr6rq4GjH0pJE5CLgcmC0qu7d3/mxxqrcYVDVr4A8/3YildANEZH+qrrWt/kr4JtoxtPSRGQ8cC0wMh6TGayEbpI2lNDPAQMBD94hqJer6qboRtVyRGQdkAps8+36SFUvj2JIjWYJbUwCsddWxiQQS2hjEogltDEJxBLamARiCW1MArGENiaBWEIbk0AsoeOYiFwmIj/6xiv7f4a00me7fZ+3WkS+EJE/ioijzjH/zxEBv/8kIpsCtlPq3DddRN4REadvu6zO8YtF5J79xJYiIu+KSJvrCdnmvnCCORS4QVX/FYXPLveNk0ZE8oCFQEfgpsBjAfzn3gyUqeqdIe47BXheVd1NDUxVK0VkKXAusKCp94lHVkLHtyHA59EOQlW3AFOBq6T5Y0svAF7a71k+IvKiiKz01RSmBhx60XevNsVK6Ph2CPCoiHh82/epar3paEUkS1V3tGQgqrreV+XOA9JF5HPfoe9U9fRw7uGrfvdR1e8DdgfeCyAbeDlge4qqbheRdGC5iDynqtuAVcCRTfs28csSOk75pkPa4psuZ3/mARe3bEQA+EvnYFXucOQAO+vsq3UvEbkY71hlvxki4v8HoyfQH9imqm4RqRSRTFUtbUIscckSOn4dSojhjCJyE96SbCfwIXCQiMwE7gVuwzuiaAdwD95ZOV7FW9p/gHfOtJtVdVVjghGRPoAb2NKE7+JXDoS93q+InASMAY5V1b0i8nad61OBimbEE3fsGTp+DSFIQotIdyAZbzIfg3dFxf/4GqGmAwtVdRpwEDAUbwPU3/A2aD0EPAMc6LvXUt/9GiQiucADwD3NmeHD91jgFJFwk7ojsMOXzAfh/b7+mDoDJapa1dR44pGV0PFrCDBSRCb4thUYAdwC/A7IxVsFPRT4wnfOIcD/+Z5V9+JteX5BRJLxVlM9IjIYeMj3PNwP2B7i8/3PtslANfA4cFcEvtcS4ATgjTDOfR24XES+BIrwTuznNwpYFIF44oqNh04wInI10AHojPcv+QbgTOAOoA9wDt5kng/MBC7Fm/RjVPVOEXlcVX/tS+wpqvrHVo5/KPBHVf11M+/zPHC9qhZFJrL4YAltYo6ITAH+3dR30b4ayHmq+lhkI4t9ltDGJBBrFDMmgVhCG5NALKGNSSCW0MYkEEtoYxKIJbQxCcQS2pgE8v8BCPmRYG3NxCUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x252 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5, 3.5))\n",
    "\n",
    "# Plot the actual vs predicted\n",
    "ax.scatter(data['u0_atom'], data['u0_atom_pred'], s=10)\n",
    "\n",
    "# Place a y=x line on the graph. We fix\n",
    "#  the x and y limits first so that they do not\n",
    "#  re-adjust to create more whitespace after \n",
    "#  adding the new line\n",
    "ax.set_xlim(ax.get_xlim())\n",
    "ax.set_ylim(ax.get_xlim())\n",
    "ax.plot(ax.get_xlim(), ax.get_xlim(), 'k--')\n",
    "\n",
    "ax.set_xlabel('$E_{atom}$, DFT (Ha)')\n",
    "ax.set_ylabel('$E_{atom}$, ML (Ha)')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting Hyperparameters\n",
    "The key hyperparamter for a KRR model is $\\lambda$, which controls the tradeoff between fitness to the training set and model complexity.\n",
    "\n",
    "Scikit-learn provides a convenient tool, [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html), for optimizing hyperparameters that we will use for this exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# GridSearchCV wraps over a model performs hyperparameter optimization when you call fit\n",
    "#  by re-running the model with a series of different parameters. For efficiency,\n",
    "#  I wrap the KRR part of the model rather than the whole pipeline so that we do \n",
    "#  not need to re-run the kernel code\n",
    "model = Pipeline([\n",
    "    ('kernel', FCHLKernel()),\n",
    "    ('krr', GridSearchCV(\n",
    "        estimator=KernelRidge(kernel='precomputed'), # Set the estimator to be KRR\n",
    "        param_grid={'alpha': np.logspace(-8, 0, 32)},  # Lambda is known as alpha in sklearn\n",
    "        cv=ShuffleSplit(test_size=0.1, random_state=1),  # Use 10, 90/10% splits\n",
    "        scoring='neg_mean_absolute_error',  # Score based on MAE\n",
    "    ))\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 23.7 s, sys: 60 ms, total: 23.8 s\n",
      "Wall time: 1.11 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Pipeline(steps=[('kernel', FCHLKernel()),\n",
       "                ('krr',\n",
       "                 GridSearchCV(cv=ShuffleSplit(n_splits=10, random_state=1, test_size=0.1, train_size=None),\n",
       "                              estimator=KernelRidge(kernel='precomputed'),\n",
       "                              param_grid={'alpha': array([1.00000000e-08, 1.81160919e-08, 3.28192787e-08, 5.94557071e-08,\n",
       "       1.07710506e-07, 1.95129342e-07, 3.53498111e-07, 6.40400427e-07,\n",
       "       1.16015530e-06, 2.10174801...78538e-06,\n",
       "       1.24960914e-05, 2.26380341e-05, 4.10112707e-05, 7.42963951e-05,\n",
       "       1.34596032e-04, 2.43835410e-04, 4.41734470e-04, 8.00250228e-04,\n",
       "       1.44974067e-03, 2.62636353e-03, 4.75794431e-03, 8.61953566e-03,\n",
       "       1.56152301e-02, 2.82886943e-02, 5.12480588e-02, 9.28414545e-02,\n",
       "       1.68192432e-01, 3.04698957e-01, 5.51995432e-01, 1.00000000e+00])},\n",
       "                              scoring='neg_mean_absolute_error'))])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "model.fit(train_data['fchl_rep'].tolist(), train_data['u0_atom']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAACsCAYAAABM8oFkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAU/0lEQVR4nO3deZQU9bXA8e+dBQaVRZZhF0RAHVBBx4AxEiPBIAkiJuaJifKUFzVGnzE5SXxZntt5LxyjWYg+cSPBaFQEoriQQDAoiBhmAAVBhIxs4zAM+8gMzNL3/VHV2Iw90zPTVV3dNfdzTp/urq1vcbhTVb/61f2JqmKMCYesoAMwxnjHEtqYELGENiZELKGNCRFLaGNCJCfoAFKhe/fuOnDgwKDDMMYzxcXFe1S1R8PpbSKhBw4cSFFRUdBhGOMZEdkWb7qdchsTIpbQxoSIJbQxaSTZnpuW0Makkbe27OXS37zBvyo+adX6ltDGpJFFG3axfV8VfTp3aNX6ltDGpAlVZfGGcsYM6UGHdtmt2oavCS0i40Vkk4hsEZE748wXEZnhzn9PRM51p+eJyD9F5F0ReV9E7olZp6uILBaRze77yX7ugzGpsr70EGUHjzCuoGert+FbQotINvAwcBlQAEwRkYIGi10GDHFfNwKPuNOPApeo6jnACGC8iIx2590JLFHVIcAS97sxGW/Rhl1kCYw9Mw0TGvgcsEVVS1S1BngOmNRgmUnAU+pYCXQRkd7u92irQK770ph1ZrufZwNX+LgPxqTM4g3lFA7sStcT27V6G34mdF9gR8z3ne60Zi0jItkishbYDSxW1XfcZXqqahmA+54f78dF5EYRKRKRooqKimT3xRhfbd9bxQe7Krk0idNt8DehJc60hjfZGl1GVetVdQTQD/iciAxvyY+r6mOqWqiqhT16fKbLqzFpZdGGXQBcWtArqe34mdA7gf4x3/sBH7d0GVU9ACwFxruTykWkN4D7vtuziI0JyKIN5ZzRqyOndDshqe34mdCrgCEicqqItAOuBhY0WGYBcJ3b2j0aOKiqZSLSQ0S6AIhIB+DLwAcx60x1P08FXvJxH4zx3b7DNRRt3Zf06Tb4+LSVqtaJyK3A34BsYJaqvi8iN7vzZwKvAROALUAVcL27em9gtttSngXMUdVX3HnTgTkiMg3YDlzl1z4YkwpLNpYTURiX5Ok2+Pz4pKq+hpO0sdNmxnxW4Htx1nsPGNnINvcCY72N1JjgLNpQTu/OeQzv2ynpbVlPMWMCVF1Tz7LNFYwr6IlIvDbilrGENiZAyzZXcKQ2knTrdpQltDEBWryhnI55OYwa1NWT7VlCGxOQuvoIf99YziVn5JOb7U0qWkIbE5DibfvZX1Xr2ek2WEIbE5hFG8ppl53FF0/3riejJbQxAairj/Dyux8zZmh3Tmrv3d1jS2hjArB8yx52Vx7l6+f283S7ltDGBGDe6lI6d8jlkjPjPizYapbQxqTYwepaFr2/i8vP6UP7nNaVGmqMJbQxKfbqe2UcrYvwjfO8Pd0GS2hjUm7e6p0Mzj+Js/t19nzbltDGpNBHew5TvG0/Xz+3nyd9txuyhDYmheav3kmWwOSRDatxecMS2pgUiUSU+atLuXBwd3p1zvPlNyyhjUmRlSV7KT1Q7UtjWFS6FtrvLyL/EJGNbqH922PWuVtESkVkrfua4Oc+GOOVuat30rF9Dl8Z5l3f7YZ8q1gSU2h/HE4xwFUiskBVN8QsFltofxROof1RQB3wQ1VdLSIdgWIRWRyz7m9U9QG/YjfGa4eP1vHX9c6957xcb+89x0rXQvtlqroaQFUrgY18tqa3MRlj4fpdVNXU83UfT7chjQvtR4nIQJz6Yu/ETL7VPUWf1djYVlZo36STecU7GdDtBAoH+DsUW9oW2gcQkZOAecD3VfWQO/kR4DScMa/KgAfj/bgV2jfpoqTiE94u2evbvedYaVtoX0RycZL5GVWdH11AVcvdUTUiwOM4p/bGpK2nV24nJ0u4+vz+iRdOUroW2hfgSWCjqv46doXoqBmuycB6/3bBmORU1dTxQvEOxg/vRX4nf+49x2pWK7eI5AMXAn2AapwkKnKPknElWWj/QuBaYJ07YB3AT9063/eLyAicU/OtwE3N3VljUu2ltR9TeaSO6y4YmJLfE6fWfSMzRb6EM/5yV2ANzjhSecBQnOvYucCDMde3aamwsFCLioqCDsO0MarKhBnLUVUW3n6Rp9fPIlKsqoUNpyc6Qk8AvqOq2+NsMAf4Gs595nmeRGlMiBRv28/GskP8z+ThvjeGRTWZ0Kr6oybm1QEveh2QMWHx1Nvb6Ng+hytGpK4LRbN7ionIV4FhOKfcAKjqvX4EZUymq6g8ysL1ZXxr1ABO9LAIYCLNauUWkZnAvwG34dw7vgoY4GNcxmS051dtp7ZeufaC1KZJc29bfV5VrwP2q+o9wAUcf//YGOOqq4/wzDvb+cLg7pzW46SU/nZzE7rafa8SkT5ALXCqPyEZk9n+vnE3ZQePpPzoDM2/hn5FRLoAvwJW49wDfsKvoIzJZH9auZU+nfMYe4a3JXqbo1kJrar3uR/nicgrQJ6qHvQvLGMy05bdlby1ZS8/+srp5Hg0AF1LNJnQInJlE/OI7WNtjIGZb5SQl5uVkn7b8SQ6Qk9s8PnlmO8KWEIb4yo9UM2La0r59ugBdDupfSAxJOpYEu1bjYisif1ujDne42+WAPCdMYMCi6ElJ/mNd/o2po3b88lRnv3ndiaP7EvfLh0Ci8OqfhrjgT+89RE19RFuvvi0QONI1Cj2Mp8emQeJyHHPM6vq5X4FZkymOHSklqdWbOOy4b1S3pGkoUSNYrGVNeOW+jGmrXt65TYqj9Zxy8WDgw4lYaPYG6kKxJhMVF1Tz5PLPmLM0B4M7+v94HMt1eQ1tIi8LCIT3fpeDecNEpF7ReQG/8IzJr3NKdrB3sM13BLwtXNUokax7wAXAR+IyCoReU1EXheREuBRoFhVZzW2sk8jZ3QVkcUistl997cuqjGNqK2P8NibJZw34GRGndo16HCAxKfcu4AfAz9262P3xnlQ40NVrWpqXR9HzrgTWKKq090/EncCP2nhfhuTtJfWfkzpgWruu2JYyiqSJNLs21aqulVV31bVtYmS2eXXyBmTgNnu59nAFc3dB2O8Ulcf4fevb+bM3p340umpfwijMZk4ckZPVS0DcN/j/mvayBnGT/PXlLJtbxV3fHlI2hydITNHzmgWGznD+KXWPToP79uJcQU9gw7nOIlauTs1Me+UBNv2ZeQMoDxabN99350gDmM8Na94Jzv2VfODcUPT6ugMiY/QS6MfRGRJg3kvJljXl5Ez3HWmup+nAi8liMMYz9TURfj961s4p3+XtLp2jkqU0LF/fhq2yzf5p8kt8xsdOWMjMCc6ckZ09AyckTNKcEbOeBy4xZ0eHTnjkjgDu08HxonIZpwW9OkJ9sEYz8wp2kHpgfQ8OkPirp/ayOd43z+7sjN0zWsNps2M+azA9+Kst5xG/mCo6l5gbKLfNsZrR2rrefgfWygccDJjhnQPOpy4EiV0voj8ACe5op9xv1tLk2lTnl+1g7KDR3jwqnPS8ugMiRP6caBjnM9gRQJNGxI9Oo86tSsXnNYt6HAalain2D2NzROR870Px5j09Mw729ldeZQZU0am7dEZWjAUDoCIFOC0Vk8BDgKfGf3OmLDZf7iGh/+xhc+f1o3Rg9L36AzNSGgRGYCTwFNw+lgPAApVdau/oRmTHqYv/IBD1bX898SCoENJKFHHkhU4rdS5wDdU9Tyg0pLZtBX//GgfzxftYNpFp3JGr0b7WaWNRPehK3Aawnryaau2FQs0bUJNXYSf/WUdfbt04PaxQ4IOp1maTGhVnQSchTP8zT0i8hFwsoh8LhXBGROkJ5aXsHn3J9w7aRgntEvdkLDJSBilO+TNLGCWiPTEGVb2tyLSX1VtBEoTSjv2VTFjyWbGD+vF2DPT6wGMprToaStVLVfVGar6eeALPsVkTKBUlV+8tJ5sEe66PP0bwmIlKuPb8GGKhqyMrwmd19btYummCn7xtQJ6dw6uaH5rJDrlvgCnAMGzOAUG0veOujEeOFhdyz0vv8+wPp2YGsD4zslKlNC9cJ5omgJcA7wKPKuq7/sdmDGppqr8dP469h2u4YmphYEMB5usRK3c9ar6V1WdCozGecxxqYjclpLojEmhOUU7eHVdGT+89HTO7tcl6HBapTk9xdoDX8U5Sg8EZmDDyJqQ2bK7krsXbOALg7tzU4CjRyYrUaPYbGA4sBC4R1XXpyQqY1LoSG09tz27lg7tsvn1N88hKytzm4oSXSRcCwwFbgdWiMgh91UpIgmL9rW20L47b5aI7BaR9Q3WuVtESuNUMjGmVaYv/ICNZYd44Kqzye+UF3Q4SUn0+GSrWwWSLLQP8EfgIeCpOJv/jao+EGe6MS2yZGM5f1yxlesvHMglZ2ROB5LG+NmM1+pC+wCq+iawz8f4TBtXfugIP5r7Hmf27sSdl50RdDieSPtC+4241T1Fn9XY2FZWaN80pT6ifP+5tVTX1PP7KSNon5MddEieSOtC+414BDgNGAGU0ci41VZo3zTlode38HbJXu6dNIzB+R0Tr5Ah/EzopArtN8btT16vqhGcOmf25JdpkZUle/ndkg+5cmRfvnFev6DD8ZSfCd3qQvtNbTR6je2aDNitNNNsez85yu3PrWFgtxO574rhaV0frDV8e8hTVetEJFpoPxuYFS20786fiVMNZQJOD7Qq4Pro+iLyLHAx0F1EdgJ3qeqTwP0iMgLn1HwrcJNf+2DCJRJRfvjCu+yvqmXWv5/Pie0z4xnnlvB1j1pbaN+dN6WR6dd6GaNpOx5fVsLSTRXcd8VwhvXpHHQ4vsi83ufGtMLq7fv51d82MeGsXnx7VKJxFjOXJbQJvf2Ha7jtz2vo1TmPX155duium2OF7yLCmBiRiHLHnLVUVB5l7ncvoHOH3KBD8pUdoU2oPfLGv5zqIxMLMvaRyJawhDahtWLLHh5ctIlJI/qE+ro5liW0CaXyQ0f4z+fWMKjHSfzv5LNCfd0cy66hTejU1ke49c+rqaqp57kbzw3l/ebGtJ09NW3GA3/bxKqt+/nd1SNC1U+7OeyU24TKM+9s49E3S/j26FOYNKI5D+6FiyW0CY1X3vuYn7+4ni+d3oO7Jg4LOpxAWEKbUHjzwwrueH4thQNO5v++dR65GViC1wttc69NqKzevp+b/lTM4PyOPDH1fDq0C0exgtawhDYZbdOuSq7/wyryO7Vn9g3nh74nWCKW0CZjbS6v5Non3yEvN4unp40iv2NmV+z0giW0yUgvFO3g8ofeIqLKUzeMon/XE4IOKS3YfWiTUapq6vj5i+uZv7qU0YO6MuPqkRlfS9tLvh6hfSq031VEFovIZvc9btVPEz6bdlVy+UNv8Zc1pdw+dgjP/MdoS+YGfEvomEL7lwEFwBQRaTh6dmyh/RtxKnpG/REYH2fTdwJLVHUIsMT9bkLsaF09Ty7/iEkPL+dAVS1PTxvFHeOGkp3BQ9b4xc9T7mOF9gFEJFpoP3bkjGOF9oGVItJFRHqrapmqvikiA+NsdxJOrTGA2cBS4Cf+7IIJUiSivPzexzywaBM79lUzZmgPZ7gaa/xqlJ8JHa+I/qhmLNMXp952Y3pGK4OqapmI5HsQq0kjqsobH1Zw/183saHsEAW9OzH7hrMYM6R7m3lqqrX8TGi/Cu0378dFbsQ5jeeUU9rGs7CZrvJILQvX7+KFoh2s2rqf/l078LurRzDx7D4ZPSJkKvmZ0L4U2gfKo6flbo3u3fEWUtXHgMcACgsLPfkjYbxXWx9h2eYK5q8uZfGGco7WRRjY7QTunljANaMG0C7H7qy2hJ8JfazQPlCKU2j/mgbLLMAZp+o5nNPxhIX23XWmAtPd95eSCXLOqh0s37InmU20SLJnjLGrR08/JWamIMd+Iyvmu/MSBMgSceaJkCVCdhZkZQnZImRnOdNysoSc7Cxys51pOdlZ5GbJ8ctlOctlSXSbQlaWG5c640fVqxJx36tq6qmoPMqug0coP3SE8sqjbN1zmIPVtZx8Qi7fLOzP5HP7MrJ/Fzu1bqVMLLQ/HZgjItOA7cBVycS580A160oPJrOJZnPa/pJY/7htRafpse+xm1d15kRUnXnRae7niJtoEXU+10U+Tbwkw0yoc4dcenZqT89OeYwf1osvF/Tki0N72NHYA5Lsf7JMUFhYqEVFRUGHkTHUTfC6eqU2EqHefa+rV+ojSkQ/fXf+EHz6hyOizh8FgWNH++ws55WXk01+p/bk5bbdhye8IiLFqlrYcLr1FDOfISLkZgu52dABS75MYuc4xoSIJbQxIWIJbUyItIlGMRGpALal+Ge7A6m7H9Y6FqM3gohxgKr2aDixTSR0EESkKF4rZDqxGL2RTjHaKbcxIWIJbUyIWEL757GgA2gGi9EbaROjXUMbEyJ2hDYmRCyhjQkRS2hjQsQS2pgQsYROMRG5SERmisgTIrIi6HgaIyIXi8gyN9aLg46nIRE5041troh8N+h4GiMig0TkSRGZm4rfs4RugSZqhTdZfzyWqi5T1ZuBV3CqlqZlnDh1ED4B8nBKRaVVfKq60f13/CbgSy8tj+IsUdVpfsTX2A/aq5kvYAxwLrA+Zlo28C9gENAOeBenDvlZOEkb+8qPWW8O0Cld4wSy3PV6As+kW3zuOpcDK4Br0vXfMWa9uan4P2oFDlpA49cKj1t/XFV/CXwt3nZE5BSc+mmH0jlO136gfTrGp6oLgAUi8irwZy9j9DLOVLJT7uQ1Vlu8KdOAP/gWUXwtilNErhSRR4E/AQ/5HBu0PL6L3WGUHsWpTZcqLY2zm4jMBEaKyH/5HZwdoZPX4triqnqXT7E0pUVxqup8YL5/4XxGS+NbijNqSqq1NM69wM3+hXM8O0InrzW1xYOQ7nGme3xRaR2nJXTyjtUfF5F2OPXHFwQcUzzpHme6xxeV3nGmouUtLC/gWZxxt2px/lJPc6dPAD7Eaf38mcWZ2fFlWpyxL3vaypgQsVNuY0LEEtqYELGENiZELKGNCRFLaGNCxBLamBCxhDYmRCyhjQkRS2jjCRE5S0S2pXP1kLbAEtp4QlXX4fRrvi7oWNoyS2jjpd3AsKCDaMssoY2XpgPtRWRA0IG0VZbQxhMiMh44EXgVO0oHxhLaJE1E8oD7gVuAdcDwYCNquyyhjRd+DjylqluxhA6UJbRJioicDowDfutOsoQOkBU4MCZE7AhtTIhYQhsTIpbQxoSIJbQxIWIJbUyIWEIbEyKW0MaEyP8D4mK/qYfF0eUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5, 2.5))\n",
    "\n",
    "gs = model[1]\n",
    "ax.semilogx(np.array(gs.cv_results_['param_alpha']), -1*gs.cv_results_['mean_test_score'])\n",
    "\n",
    "ax.set_xlabel('$\\lambda$')\n",
    "ax.set_ylabel('MAE (Ha)')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the error can change by about a factor of 2 when we adjust $\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12min 17s, sys: 180 ms, total: 12min 18s\n",
      "Wall time: 34.8 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "data['u0_atom_pred'] = model.predict(data['fchl_rep'].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE = 0.0089 Ha = 5.61 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "mae = mean_absolute_error(data['u0_atom'], data['u0_atom_pred'])\n",
    "print(f'MAE = {mae:.4f} Ha = {mae*627.50947406:.2f} kcal/mol')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By tuning the hyperparameter, we reduced the error to ~5 kcal/mol!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benefit of Scalable Kernel\n",
    "Here, we repeat the \"split on molecule size\" test from the [Coulomb Matrix example](./fitting-coulomb-matrix.ipynb) to see whether the validation results are as poor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Made a training set of 100 entries and a test set of 881 by splitting based on molecule size.\n"
     ]
    }
   ],
   "source": [
    "size = np.percentile(data['n_atom'], 50)\n",
    "train_data = data.query(f'n_atom<={size}').sample(100)\n",
    "test_data = data.query(f'n_atom>{size}') \n",
    "print(f'Made a training set of {len(train_data)} entries and a test set of {len(test_data)} by splitting based on molecule size.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train on the smaller size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 17.2 s, sys: 100 ms, total: 17.3 s\n",
      "Wall time: 826 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Pipeline(steps=[('kernel', FCHLKernel()),\n",
       "                ('krr',\n",
       "                 GridSearchCV(cv=ShuffleSplit(n_splits=10, random_state=1, test_size=0.1, train_size=None),\n",
       "                              estimator=KernelRidge(kernel='precomputed'),\n",
       "                              param_grid={'alpha': array([1.00000000e-08, 1.81160919e-08, 3.28192787e-08, 5.94557071e-08,\n",
       "       1.07710506e-07, 1.95129342e-07, 3.53498111e-07, 6.40400427e-07,\n",
       "       1.16015530e-06, 2.10174801...78538e-06,\n",
       "       1.24960914e-05, 2.26380341e-05, 4.10112707e-05, 7.42963951e-05,\n",
       "       1.34596032e-04, 2.43835410e-04, 4.41734470e-04, 8.00250228e-04,\n",
       "       1.44974067e-03, 2.62636353e-03, 4.75794431e-03, 8.61953566e-03,\n",
       "       1.56152301e-02, 2.82886943e-02, 5.12480588e-02, 9.28414545e-02,\n",
       "       1.68192432e-01, 3.04698957e-01, 5.51995432e-01, 1.00000000e+00])},\n",
       "                              scoring='neg_mean_absolute_error'))])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "model.fit(train_data['fchl_rep'].tolist(), train_data['u0_atom']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate on the larger size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_99/260439986.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  test_data['u0_atom_pred'] = model.predict(test_data['fchl_rep'].tolist())\n"
     ]
    }
   ],
   "source": [
    "test_data['u0_atom_pred'] = model.predict(test_data['fchl_rep'].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAAD0CAYAAACsLwv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAohElEQVR4nO3deXzU1dX48c+Z7AkEAoRV1rIoiI9AWhdApQgKosBTweLyPEoFV1CrrWtxq0/7Ky51qVrc6gIq7qIoDIgrimwqMBpBBFmHsAdISJg5vz9mEjLJJJkkk8yS8369eDUzc+c7Z5oc7/d7v/eeK6qKMSY+OCIdgDEmfCyhjYkjltDGxBFLaGPiiCW0MXEkMdIB1LdWrVpply5dIh2GMWGxZcsWCgsL2bt3705VzS7/etwndJcuXVi2bFmkwzCmzh588EFuvPFGJk2axFNPPbUxWBs75TYmBjz33HPceOONjBs3jieeeKLSdpbQxkQ5VWXWrFkMHz6cF198kYSEhErbxv0ptzGxTkSYM2cOHo+HlJSUKttaD21MlFq6dClnn302e/bsITU1lYyMjGrfYz20MVHE6XLz2do8Ojv2cPNlY2natCkFBQVkZWWF9H5LaGOixPR5uTz5yU8c3uPGPfPPNE124HQ6ad++fcjHsIQ2Jgo4XW6e/HgdRQf24n71DjxFBYy793m6d+9eo+PYNbQxUeCztXl4FLxFhxBHAm3H3cn5wwbV+DiW0MZEgZM7NyPJISRltafdxMdI79inVsexhDYmwoqLi3ly2jU4vnwaAHEkcMSrfLY2r8bHsoQ2JoK8Xi9/+MMfmDNnDmcNPom0JN+kkbSkBAb3qDBVu1o2KGZMA5o+L5cFru2c2bstNw3vyQ033MCLL77Ivffeyx133Fx622pwj2yG9W5T4+NLvNcUy8nJUVucYaLB9Hm5/GvRutLHx25zMu+Fh7n++ut58MEHEZGQjyUiy1U1p/zzdsptTANZ4Noe8DgvpT1XXnklDzzwQI2SuSqW0MY0kDN7twXgSP5OALRjP7qNuR6HI3xpaAltTD1wutxMe2c1Tpe79LkTOzan9d41bP33JA6tXcLug0X8a9E6ps/LDdvnRsWgmIhMB84FioCfgMtUdW+QdjcAlwMKrPK3K2zAUI2pltPlZurLKyko9vDass1MHNSVT3/cwfKvFrN99jSSsjuT2qlvafsFru386axeYfnsaOmhncDxqnoC8CNwa/kGItIBmArkqOrxQALw+waN0pgQfLY2j4JiDwAFxR6eWLSO5Su+Yfvr95CQ2ZrW4+7GkZJe2r7kVDwcoiKhVXW+qh7xP/wKOKaSpolAmogkAunA1oaIz5iaGNwju/R+MkDxoX24Z0/DkZJBmwvuISG9WelrI/u2C1vvDFFyyl3ORODV8k+q6hYRuR/4BSgA5qvq/GAHEJHJwGSATp061WOoxlQ0rHcb2jdP5ae8gwAkpDej2SnjSOvan8TM1gzonEWf9pm1vtdclQbroUVkgYisDvJvdJk2twNHgJlB3p8FjAa6Au2BDBG5ONhnqeoMVc1R1Zzs7JrPtjGmLq6euYKf8g7iKdhPUd4GADJzRpPUsiMArq376yWZoQF7aFU9s6rXReR/gVHAUA0+2+VM4GdVzfO3fxM4FXgp3LEaU1tXz1zB3FXb8BYVsOO1uziyfwcdrngaR1JqaZuCYg+frc2L7YSuioicDdwMnK6qhypp9gtwsoik4zvlHgrYFDATNX73xGKWb9yDHikm782/UrR9HdljbwtIZqj9PO1QREVCA48BKYDTP2PmK1W9UkTaA0+r6khVXSIirwMr8J2WrwRmRCxiY8q4euYKXzJ7PeycM53Cjd/ScuQNpPc4GYDOLdIZ9V/tyS8srrfTbYiShFbVoGUZVHUrMLLM4zuBOxsqLmNCMX1eLnNXbQMgf8X7HPpxMVm/nUSTvkNL26QmOcI6ml2ZqEhoY2JVyTVziaYnjiAhI4uM4wYHtAvnveaqWEIbUwtOl5uHF65l9ZZ9ABxYvZC0bjkkpDcLSGaHwFVndG+Q3hmiZGKJMbHE6XJzzcwVpcmc/82H7Hr/IfZ//WaFtv++JKfBkhksoY2pEafLzV3vrKbI4wXg4A+fs3v+46R2HUDzwYHTIq4Z0r3eBr8qY6fcxoRo+rxcnli0Dq//ccHPK9k5535S2h9L9thbkYSk0rbhntIZKuuhjQmB0+XmiY+PJrOqlz0fP0tSy2Noff60gHvNnVuk8fhF/SMSp/XQxlTD6XJz8+vf4i0zf1HEQZtx9wDgSG0S0P6OUbUrwRsOltDGVKL8SDZA8d7t5C+fQ9YZl5HQpOJ+UyP7tmvw6+ayLKGNCcLpcnPVS8s5UqZb9hzYw45X/4K3MJ+m/UeRlNUu4D1tM1Midqpdwq6hjQni4YVrA5LZW3gA9+y/4Dm4h9bn31UhmQHyDhQFlByKBEtoY8pxutysKXOa7S0uZMfr91C8azPZY28npcOxQd/nqeVuF+FkCW1MObOWbKTs+t3iXZsp2rmRVufeRFrXfgFtB3TOIsHhK8Fbn6uoQmXX0MaUMX1eLp/8GNjLprTtTocrniYhrWnA802SE3jjqlPrvNtFOFlCG+NXdmcLVWXPgidJzGpPZs7oCskM8L8DuwK+kkORTuQSdsptjN/LX28s/XnfZy+Rv+J9PPm7KrRrnpbINUMabsFFTVgPbRq9kvvNuw8WA7B/6Tvs+/JVmpwwnOZnXBbQNtEhTB93YtT0yOVFRQ8tIveKyHci8o2IzPdXKgnW7mwRyRWRdSJyS0PHaeLP9Hm5TH5xWcAyyD0fPUV6z1NpcdY1AXtOpSY6eOLiAVGbzBAlCQ1MV9UTVPVE4D1gWvkGIpIA/AsYAfQGJohI7waN0sQVp8vNE4vWUbYkpffwIVK79KPVuX9CHAkB7R+9sH9UJzNEySm3qu4v8zADCFb18zfAOlVdDyAir+Ar6+uq/whNPJq1ZOPRxRZHipDEZDIHnEvTfiMrJHOkp3SGKlp6aETkPhHZBFxEkB4a6ABsKvN4s/+5YMeaLCLLRGRZXl5kb/Sb6ON0ubnsua/5JNf3t3F4+zq2/HsShZvXAFRI5uZpiRGf0hmqqCm0r6q3q2pHfEX2rw12iCDPBd2t3grtm8qUzNFelJuHFyjetYkds6eBI4HEZsF74PNODNpvRKWoKbRfxizgfSpW99wMdCzz+BhsbytTA06XmymzVpTO0T6yfwfuV6eBOGhzwb0kNm1V4T3RMPurJqLilFtEepR5eB7wQ5BmS4EeItJVRJLx7Tz5bkPEZ2JfSc9ceMR31ewpyMf96jS8hw/SZvw9JLWo2Asf36EZj0zoFxPXziWiYlAM+LuI9AK8wEbgSoByhfaPiMi1wDx8W8k+q6prIhaxiSmPLPwxYPWUIzmN1I59yOgzhOQ23Sq0H9m3XcxcN5clwbeRih85OTm6bJntmNOYlWxRA77RbO/hQyRkNK+0fd8OmcyZMrjS16OBiCxX1Zzyz0fFKbcx9WX6vNyjyez1kPfuP9g+62a8xYeDtncAU4f2bMAIw8sS2sQtp8vNf774GfAV9dv1wSMUrP2Kpv1H4UhKqdA+wSFcFYHSu+EULdfQxoSV0+Vm6ssrKSj2+FZOffQMB1cvpNnAC8kccG5A2wSB03pmc+FJnWM6mcES2sSpWUs2UlDsAeDANx+Qv+wdmg44l2YDJwS0G9IrPhK5hCW0iRtOl5tZSzaSd6AooIRQeq+BeA7to9mpFwQstojWJZB1YQlt4kLJflMlW9QAFG5eQ0q7niSkN6N5uZ55QOesuEtmsEExEyc+W5sXkMwF65fjfvk29n4+q0Lbzi3SeOOqUxsyvAZjPbSJOcFqeDVNPbqvVOFmF3lv/R/J2V1odvL5Ae9NTnBEdGeL+mYJbWJK2dHr15ZtZuKgruQXFrPyF9+95qIdP5P3+t0kZLai9bi7caRklL63Q/NU7jrv+LgZAAvGEtrElM/W5pWOXhcUewJ2g1Svh7x3/h+SnEab8fcGzAZLdEjcJzNYQpsYM7hHNq98van0etlb5jVxJJB93p8hIZHEZq1Ln09NdMREtZFwsEExE1OG9W7DwO4tA57zFOST/+18AJLbdCO5VaeA10/5VctGkcxgCW1iUO/2zUp/9hYVsuP1u9jtfJziPRWXxyc6hAtP6tyQ4UVUjU+5RSQDKFRVTz3EY0yVrp65grmrtgGgR4rJe+s+iratJXvMLSRlBRaL7dshk6lDezaa3hlCSGgRceArJnAR8GvgMJAiInnAXGCGqq6t1yhNo+d0ufnrey427j4E+AbAdr73AIUbVtJyxHWk9wy8r9ypRXrUL4GsD6Gcci8CfgXcCrRV1Y6q2hoYDHyFrzjBxXUJIpS63CLSUUQWicj3IrJGRK6ry2ea2FFSbaQkmQEOb3ZxKPcLsoZMpMkJwwLaC/CXUY2zwnO1BQ5EJElVi+vappr3Z5aU8hWRqUBvVb2yXJt2QDtVXSEiTYHlwBhVrbKMrxU4iH2XPfc1i3IrVm8tyttAcnaXgOcay4h2rQschJKodUlm//urrcutqttUdYX/53zgeyop42vi1/6l71CwfjlAhWSGxjWiHUyNBsVEJAvoAaSWPKeqn4YjEBG5D/gfYB8wpJq2XYB+wJJwfLaJTiVTPNOSfX+m+d/OZ89HT5HRZwhp3QZUaC/QqEa0gwk5oUXkcuA6fOVzvwFOBr4Efhvi+xcAbYO8dLuqvqOqtwO3i8it+Opyly/jW3KcJsAbwPXlevaybSYDkwE6deoUrImJcuVXTx3KXczueY+R2rU/LUdMrdA+MzWRB8ZH7yZyDaUmPfR1+Ea5v1LVISJyLHB3qG8OQ11uRCQJXzLPVNU3q/isGcAM8F1DhxqjiR6PLPyxNJkLNnxD3px/kNKuJ9ljbkMSkgLaDuicFberp2qqJhNLClW1EEBEUlT1ByAsC0pDqcstvpXpzwDfq+qD4fhcE52cLjdrth49+SpYt4SkrA5kn38njuTUCu37tM9syPCiWk166M0i0hx4G3CKyB7Ct3NFtXW5gYHAJcAqEfnG/77bVHVumGIwUWLWko14FVQVESFr6GT08EEcqU0qtE1OcMTUzhb1LeSEVtWx/h/vEpFFQDPgg3AEoaq/q+T5rcBI/8+fE3x/KxNHnC43i3LzOLLPzc45D9By5HUkteiAlEtmh/h65sY2E6w6ocwUe5RKNoUDzgQqjlAYU0Ml9cA+/TEPz8E9uF+9A++h/ain4h3R7CYp/N9/97VEDiKUHrrsrIy7qWT02ZjaKjui7T18EPfsO/Ec2E2bC/5a4V5zZmoiS+8IdXy18ak2oVX1+ZKfReT6so+NCYdZSzb6krm4kB2v30Pxzl9o/bu/kNLhuAptHxh/YsMHGENqunzSbgGZsPJNHtnpe+A5AiitRv2xwsSRJIfw1P/k2Gl2NaxiiYmoe99zUezxgMeDI7UJbS78O74FfkcJ8PjFAyyZQxDKoFg+R3vmdBEpuUEogKqq3QQ0tXL1zBVs3HWQPQufoihvI23G3Y0kBk4asRlgNRPKNXTThgjENB4lI9qLcvPY98XL5C+fQ+avx0JC4J/jkF7ZPHfZbyIUZWwKpYcWrWaNZShtjIHAEe39y95l3xezyOh7Js2HTAzYpiY5wdHoF1rURkgFDkRkiogErHIQkWQR+a2IPA/8b/2EZ+JNyRztg65P2LNwBmk9Tqbl2VMCktnUXigJfTbgAV4Wka0i4hKR9cBaYALwkKr+px5jNHHC6XKzaotvCCapdVcyep9B9nl/RhwJgG9b1xJFHi+fra1Y1MBULZRr6ELgceBx/2qnVkCBqu6t59hMnJg+L5e3V25m295Cjuxzk5DZmuRWnWh17k0B7U7rmc1X63dTUOwhLSnB5mjXQo1uW/krk2yrp1hMHJo+L5d/LVoHQJH7J7bPupVmp4wPuufUhSd15sKTOlfYt8qEzu5Dm3r19srNABTv3oJ79jQcqRlk9D49oE35cruWyLVnCW3qRcmtqV0HijiyfyfuV+8AoM0FfyUx03cqfXyHZlw3tIclcBjVKaH9c7v/GaZYTJwoe2tKvR52vH4X3sIDtJnwN5Ja+Oo6juzbjscv6h/hSONPXXvoPwL/DEMcJo6U3XxdHAk0H3QRjtQmpLTtDvhKBlky14+67m0VlpuHoRTaL9M2QURWish74fhsE347DxShR4oo3LQagPSep5DaqW/p61YyqP7UtYcO1+yw6ar6FygttD8NfxmiIK7DV5Pb/iqiSMk1M8A3G3eRN2c6Beu+pv2kf5PU/GixVysZVL9qsjgjWG+cHo4gQim074/lGOAc4D58p/smCgRcM6uy64OHKfjxS7KGTgpI5sa4eVxDi5rFGSEW2v8n8GegypisLnfDKilQoKrsXfQsB1ctoNmpE8jMGR3Qbt2OgxGKsPEIpYd+t6rXVfW8UD6oroX2RWQUsENVl4vIGdXEZHW5G8j0ebml+04VbljJ/qVv0bT/KJoNurBC24JiD5+tzbMeuh6Fcg19CrAZXwH8JdRyICwMhfYHAueJyEh8W/FkishLqlqnnS9N7Tldbp74eF3p49Qu/Wg15lbSe54SdLGFTeesf6GMcrfFt5Xs8cDDwDBgp6p+oqqfhCOIUArtq+qtqnqMqnbBt1/1R5bMkfXZ2jy8Cod+/JKinb8gImT0Glih4kiiQxjSK5tHJvSz3rmehXIN7QE+BD4UkRR8K6w+FpF7VPXRMMURSqF9E0Wmz8vlve+2UrB+OXnv/D/SfpVD6/++o/R1h8BVZ3Qnv7DY5mU3oJBuW/kT+Rx8ydwFeASodG+pmgql0H655z8GPg7X55uauXrmCuau2sbhLd+T9/b/kdSqI61GXh/QpmNWGn86Kyw7JZkaCGVQ7Hl8p9sfAHer6up6j8pELafLzdxV2yjK28CO1+8moUkL2oy/p8I2Nd2yK25bY+pfKD30JcBBoCcwtcxghxUJbIRKJo/s+3I2kphMmwv+SkJGVkAbh9g+zZESyjV0XaeHmjjhdLn59EffLapWI6/nSP4uEptVvDZOdNifTKTY//MmZPNXrmPH3EfwFh5AEpNJymoXtJ2VD4ocS2gTkvveWcm/b7+CA2s+oihvQ8BriQ64Zkh30pJ8tcHsfnPkWIEDU6np83JZ4NrOnvxDrPrPXyj85Xtajb6Z1I7HB7Qb3qcdfzqrFyd2bG7lgyLMEtoEVVILTL0eds65n8Kfl9Pi7Clk9BpYoW2rJsmAr3SQJXJk1fqUW0Ta+e9Pmzi0wLUdAM+hfRzemkvzMy6l6X+dFbRt09SkoM+bhleXHvpF4Fci8oaq3lRtaxNTmqQmoaokNmlB+4mP4kjJqLRtfmHFTdlNZNS6h/YvtugGPBe+cEw0cLrcfPX2f9g97zHU66mQzJ1bpJGc4PvTsQGw6FKrHlpEegOdgVWquia8IZlIuXrmCha6trNr5Xx2zX+K9F6DKrS5Zkh3/nRWL/++zjYAFm1qe8p9N/AaMFlEOquq7W0V40rmZx/68Ut2ffiobynkqBtLt6mBo8kMNgAWrWqb0E5VnQ3MDmcwJnIWr8ujcON35L37D5Lb9SB77G0BezUnOoQTOzaPXIAmJLW9hj5VRN4UkadExGp7xYFTu2ejnmKSW3eh9fl34UhOC3j9iFdt9lcMqG0PvVpV7xeRRKBPOAMyDcvpcvPR6l8Y268T7v1nsqxrv4ACBQ4Br9rgV6yobUKPEpHDwDxV/bauQYjIvcBofAUOdgCX+tdCl2/XHHga33JOBSaq6pd1/fzGyulyc9WT89jw/J947rSLGfHfF3Dtb3sy59ut7DpwmNN7tWZsvw42+BVDapvQFwD9gLEi0l1VJ9UxjlDrcj8MfKiq54tIMmEqI9xYfbj0BzbOvA3v4YM4sruwKDePr9bvrlAqyBI5dtT2Gvoa4A9Ab+D+ugYRSl1uEckETgOe8b+nyPaorr39+/fz7j+m4tm/k9bnTyO5dTfgaGVOE5tq20Nnqeo4fy/5EL4Er5MQ6nJ3A/KA50Tkv4DlwHWqasWea6i4uJjzzjuPDT+6uPfR/7AupQdfrNtFkcdr18oxrrYJXSQi/YFV+HrUatW1Lrc/1v7AFFVdIiIPA7cAfwnyWVZovwqJiYmcffbZXHHFFUyYMAHAJorECVGteR16EUkFpuA75Z6tqh+ELSCRzsD7qnp8uefbAl/5y/giIoOBW1T1nKqOl5OTo8uWLQtXeDGpJFkH/qolvZocpkuXLpEOydSRiCxX1Zzyz9f2GvpqVZ2uqpcB6+sWWsh1ubcDm/zlfgGGAq66fna8c7rcTH15Jc8v3sD4y67kuONPYNOmTZEOy9STGp1y+28bPQQcKyKFwLfA5cBldYwj1LrcU4CZ/mv39WH43Lj32do8Coo97Fv8CvuXvUtmzmi+359Ex0gHZupFjRLaP6p8mYichW+A6gTgjboGEWpdblX9BqhwmmEq1zQ1ifwV77Hv85lkHD+U5r/9Ay9//QvD+wQbzjCxrraDYqcAJwFbgRXhC8eEk9Pl5tEX3mC380nSup9EyxFTK2xTY+JLbRO6OfAV8AJwQ9iiMWHjdLm5f94PSLveNBt8MZm/Hlu6cqp3+2YRjs7Ul9om9G7/e3f4fzZRZPq8XP758lwczdqSkJZJ81N/H/C6a+u+CEVm6lu1518iUmHxhareAzyJb48r++uIIk6Xm4dnz2fby3ew64NHIh2OaWChXFC9WPKDiFxe8rN/wGqqqj5UH4GZ2nn7k2Vse2UajpR0Wpw5OWgb26YmfoWS0GV37r663GufhjEWU0dbtmzhjb9eCeqlzQX3kpjZGjj6CxR8VUdsJlj8CuUauuxUMin3mg2ZRpjT5S7dQC73xWnk79vLqdc8xJbEdqXrmCcO6mr7NDcSoSR0WxG5FN8kkvIJXfN5oyZsnC4318xcQZHH63vi+Itp0XEYmxLakexwMLB7Sy48qbMlcSMSSkLfhW8yx2XAMSKyBt/UzB+AVvUXmqnOZ2vzOHz4MPkr5tA0ZzSSnkVCum9r1yKPl44t0i2ZG5lQtpOdUfaxiByDb4ZYX+waOqJO7daCB2+5n4M/fEFSdlfSuvYj0SEc8aotg2ykanwfWlU3A5uBueEPx4Rq/prt3HbjVA7+8AVZv72ctK79AGjfPI0zemXb9XIjZZvVxSCny82ESVPY/eUbND/lAjJ/Pab0tc17DlkyN2I2Sh2D3lv8HXuWv0+TfueQOfjigNe8ipUQasSsh45Bo049gXcuexht1pakBAdeVbz++w3JCQ67dm7ELKFjyBtvvMHWrVs5duh4klseQ5HHi9ereMu0Gdi9pZ1uN2KW0DHC6XRy4YUXkpOTwxntB5Xee/YCCQ7B4x/ZtmmdjVtUJHQNCu3fgK9CiuIrUHiZqhY2ZKyRsGTJEsaOHUuvXr147733WLatiDe/cVNQ7LGZYCZArYoEhj0IkcyS2tz+Qvu9VfXKcm06AJ/7XysQkdnAXFX9T1XHjvUigWvWrGHw4MFkZWVx37Nv4dp79BrZqnQ2XpUVCYyKHjqUQvt+iUCaiBTj2zWjQi8ebxYvXkxaWhp3Pv4ydy/YSkGxh9eWbeaRCf24Z/Tx1R/ANCpRc9tKRO4TkU3ARfi2wgmgqlvw7dLxC7AN2Keq8ys51mQRWSYiy/LyYvMWTsmZ06RJk/j+++9ZV5hOQbEHsN0tTOUaLKFFZIGIrA7ybzSAqt6uqh2BmfgK7Zd/fxa+6+yuQHsgQ0QuLt/Of6wZqpqjqjnZ2bF3C2fv3r2cfvrpLFq0CIDMzEwG98gmLclXQsimdZrKNNgpt6qeGWLTWcD7VNw540zgZ1XNAxCRN4FTgZfCFmQUOHToEOeeey5LliyhuLi49PlhvdvwyIR+dt1sqhQV19Ai0kNV1/ofBi20j+9U+2QRSQcK8BXaj93RriCKi4sZN24cX3zxBa+88grDhw8PeH1Y7zaWyKZK0XIN/Xf/6fd3wHDgOvAV2heRuQCqugR4HV/Z4FX4Yp9RyfFijtfr5dJLL2Xu3Lk8+eSTjB8/PtIhmRgUFT10DQrt30nFU/GY53S5+eSHbbjzi/jb3/7G5MnBa4EZU52oSOjGzOlyc+3zX3JYkkg9YSL9z+0X6ZBMDIuWU+5G6/4HHmT9U1fjObCHwiNePl+3M9IhmRhmCR1BL7zwAvOf/QepbbrhSM+021GmzuyUO0LeffddJk6cyNChQ7nh/mdZsnG/3Y4ydWYJHQFffPEF48ePp3///rz11ls0bdqUc06MdFQmHtgpdwT07NmT888/nw8++ICmTZtGOhwTR6yHbkC//PILbdu2JTs7m5deiqsJbiZKWA/dQDZt2sSgQYO4/PLLq29sTC1ZQjeAnTt3Mnz4cPbt28cNN9h22qb+2Cl3mDld7oAFFPn5+YwYMYINGzYwf/58+vWziSOm/lhCh5HT5WbqyysDihDcfd1Elq9Yyd2PPsfgwYMjHaKJc3bKHUafrc0LKEIwa8lGtnQZSctRNzFrWzZOlzvCEZp4ZwkdRiVFCFS9eDYuB0BbdiHjuMFWZcQ0CEvoMBrWuw0P//5E2v3wGptfuZNeng1WZcQ0KLuGDrOlbz/DkndfZMqUKdz8h3EM+H6HVRkxDUdVo+YfcBO+ip+tKnn9bCAXWAfcEsoxBwwYoA3l8ccfV0Avvvhi9Xg8Dfa5pvEBlmmQv/eoOeUWkY7AMHylhoK9ngD8CxgB9AYmiEjvhouwauvXr2fKlCmMGjWKZ599Focjav6vNY1INP3VPQT8mcprcv8GWKeq61W1CHgFXxXQqNCtWzfmzp3L7NmzSUpKinQ4ppGKioQWkfOALar6bRXNOgCbyjze7H8uohYvXsy8efMAGD58OGlpaRGOyDRmDTYoJiILgLZBXroduA1fccAqDxHkuaC9uYhMBiYDdOrUqQZR1sx3333HOeecQ4cOHRg6dCiJiTbGaCIr4nW5RaQvvuL534oIwDHAChH5japuL9N0M9CxzONjqGQrHFWdgb8iaE5OTr1s3vXTTz8xfPhwMjIyeP/99y2ZTVSI+F+hqq4CWpc8FpENQI6qli+utRToISJdgS3A74ELGyrOsrZu3cqwYcM4cuQIixYtonNn28LVRIeouIauTLm63EfwbZEzD/gemK2qayIR17PPPkteXh4ffPABxx13XCRCMCaoqNhOtj7Vx3ayqsq6devo0aNHWI9rTKgq2042qnvoaHL48GEuv/xy1q5di4hYMpuoZAkdAo/HwyWXXMIzzzzDkiVLIh2OMZWyhK6GqnL11Vfz2muvcf/993PxxUF3sDUmKlhCV+P2229nxowZ3Hrrrdx4442RDseYKllCV+Hw4cN8/PHHXHHFFdx3332RDseYakX8PnS0UlVSUlJYsGABKSkp+Ce9GBPVrIcO4s033+Sss84iPz+f9PR0EhISIh2SMSGxhC5n4cKFTJgwgfz8fFsCaWKO/cWW8fXXXzN69Gh69uzJ+++/T0ZGRqRDMqZGLKH9XC4XI0aMoHXr1syfP58WLVpEOiRjaswS2k9V6dq1K06nk3bt2kU6HGNqpdGPch84cICMjAz69OnD0qVLbTTbxLRG3UPv27ePwYMHc8sttwBYMpuY12gTuqCggEFDz+a7Vatp1u3ESIdjTFg0yoQuLi5myIjRrF6+hBYjb+DFzVm2TY2JC40yoSdNmsSST5y0GH4VGb1Pt21qTNyIqkExEbkJmA5kly9B5K/b/QK+QoNeYIaqPlybzxkzZgyJLY7hi4xBFBR7bJsaEzeiJqGrK7QPHAFuVNUVItIUWC4iTlV1hfoZa9eupUePHowZM4YxY8ZU2MvZmFgXTafcVRbaV9VtqrrC/3M+vrpiIdflfuyxxzjuuOP4/PPPS58b1rsN94w+3pLZxI2oSOgQC+2Xbd8F6AcELR8iIpNFZJmILMvLy2PmzJlMmTKFc889l5NPPjl8gRsTZRqsSGAohfZVdV8VZXxLjtME+AS4T1XfrO5ze/TooT///DOnnXYac+fOJTU1tfZfwpgoUVmRwIhX/fQX2l8IHPI/VVJAv3yhfUQkCXgPmKeqD4ZyfIfDof379+ejjz4iMzMzjJEbEzlRm9DlVdZDi28a1/PAblW9vgbHywM2hjPGSrQCgp5VxJF4/46x9P06q2qFWzNRM8odjIi0B55W1ZHAQOASYJWIfONvcpuqzq3qGMG+dH0QkWXB/osZT+L9O8bD94u6hFbVLmV+3gqM9P/8OcE3rDPG+EXFKLcxJjwsocNnRqQDaADx/h1j/vtF3aCYMab2rIc2Jo5YQhsTRyyh60hEbhIRFZFWlbx+tojkisg6EbmloeOrLRG5V0S+E5FvRGS+/xZisHYbRGSVv1149+2tZzX4jjHzO7Rr6DrwrxB7GjgWGBBkMkwC8CO+VWSbgaXAhJqsEIsUEclU1f3+n6cCvVX1yiDtNlDFVN1oFsp3jLXfofXQdVPlCjHgN8A6VV2vqkXAK8DohgquLkr+0P0yqPw7xqwQv2NM/Q6jbmJJrCi7QqyK4oIdgE1lHm8GTqrv2MJFRO4D/gfYBwyppJkC80VEgX+rakzd+gnhO8bU79B66CqIyAIRWR3k32h8q8SmVXeIIM9FTU9XzfdDVW9X1Y7ATODaSg4zUFX7AyOAa0TktAYKPyRh+I5R/Tssz3roKqjqmcGe968Q6wqU9M7HACtEpPwKsc1AxzKPS1aSRYXKvl8Qs4D3gTuDHGOr/393iMhb+E5RPw1bkHUUhu8Y1b/D8qyHrgVVXaWqrVW1i3/u+Wagf/nlnvgGUHqISFcRSQZ+D7zbwOHWioj0KPPwPOCHIG0y/OWgEJEMYDiwumEirLtQviMx9ju0HjrMyq4QU9UjInItMA9IAJ5V1TWRjTBkfxeRXvgKMm4EroQKK+DaAG/5z1ISgVmq+mGE4q2Nar9jrP0O7baVMXHETrmNiSOW0MbEEUtoY+KIJbQxccQS2pg4YgltTByxhDYmjlhCxzARuUJEtvnX85b869tAn+3xf94aEflWRP4oIo5yr5X8G1Dm5+0isqXM4+Ryx00TkU/8yxYRkQPlXr9URB6rJrZkEflURBrdxKlG94XjzAnAHar6TAQ+u0BVTwQQkdb45kI3wzcXuvS1Mkra3gUcUNX7KznuROBNVfXUNjBVLRKRhcAF+BZdNBrWQ8e2vsA3kQ5CVXcAk4FrpYq1pCG6CHgn1MYi8raILPefKUwu89Lb/mM1KtZDx7Y+wHMi4vU/fjzYemQRyVLVPfUZiKqu959ytwbS5OjuJj+r6thQjuE//e6mqhvKPF32WAAtCFwcMVFVd4tIGrBURN5Q1V34Fon8unbfJnZZQscof/mjHap6QgjNHwIurd+IgKNrh4OdcoeiFbC33HMBxxKRS4Gy29VMFZGS/2B0BHoAu1TVIyJFItLUv594o2AJHbtOIPhyP0TkTnw92V7gS+BYEbkJ+BdwH5AC7AEeA17Ftw64D7AYX+2su1S1RssgRaQb4AF21OK7lCgAQt7vV0TOAM4ETlHVQyLycbn3pwCFdYgn5tg1dOzqS/A1yh2AJHzJfDK+3RRf8g9CTcG3xPEafIUN++EbgPoHvgGtp4DXgM7+Yy30H69KIpINPAk8pnVYvue/LEgQkVCTuhmwx5/Mx+L7viUxtQTyVLW4tvHEIuuhY1df4HQRGeF/rMBg4F7gOiAb3ynoCcC3/jZ9gIf916qH8I08vyW+fbd3qapXRI4HnvJfD3cHdlfy+SXXtknAEeBFIKQ9u6sxHxgELAih7YfAlSLyHZALfFXmtSFAlTuTxiNbDx1nRORGIBNoie+PfCPwO+DvQDdgPL5kfgS4CZiEL+nPVNX7ReRFVb3En9gTVfWPDRx/P+CPqnpJHY/zJnCrquaGJ7LYYAltoo6ITASer+296JJSQar6Qngji36W0MbEERsUMyaOWEIbE0csoY2JI5bQxsQRS2hj4ogltDFxxBLamDjy/wFdOcR2w/N3FAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 252x252 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.5, 3.5))\n",
    "\n",
    "# Plot the actual vs predicted\n",
    "ax.scatter(test_data['u0_atom'], test_data['u0_atom_pred'], s=10)\n",
    "\n",
    "# Place a y=x line on the graph. We fix\n",
    "#  the x and y limits first so that they do not\n",
    "#  re-adjust to create more whitespace after \n",
    "#  adding the new line\n",
    "ax.set_xlim(ax.get_xlim())\n",
    "ax.set_ylim(ax.get_xlim())\n",
    "ax.plot(ax.get_xlim(), ax.get_xlim(), 'k--')\n",
    "\n",
    "ax.set_xlabel('$E_{atom}$, DFT (Ha)')\n",
    "ax.set_ylabel('$E_{atom}$, ML (Ha)')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE = 0.0110 Ha = 6.90 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "mae = mean_absolute_error(test_data['u0_atom'], test_data['u0_atom_pred'])\n",
    "print(f'MAE = {mae:.4f} Ha = {mae*627.50947406:.2f} kcal/mol')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still achieve reasonable accuracy on these larger molecules, thanks to the scalable kernel. It is higher than what we observed in the \"random split test,\" but the difference is not as large as what we observed for the Coulomb Matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing a Learning Curve\n",
    "Learning curves show the improvement of a model as you add more training data.\n",
    "Here, we demonstrate how to make a learning curve to illustrate a few additional points about FCHL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the test sizes we will evaluate for the learning curve\n",
    "train_sizes = [10, 30, 100, 300, 1000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare a standard test set\n",
    "train_data, test_data = train_test_split(data, test_size=500, random_state=1) \n",
    "\n",
    "# Avoids warnings about setting values of a copy\n",
    "test_data = test_data.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [03:41<00:00, 44.27s/it]\n"
     ]
    }
   ],
   "source": [
    "# Loop over the test sizes\n",
    "results = []\n",
    "for ts in tqdm(train_sizes):\n",
    "    # Split off some training data\n",
    "    train_subset = train_data.sample(ts, random_state=1)\n",
    "    \n",
    "    # Fit the model\n",
    "    train_time = perf_counter()\n",
    "    model.fit(train_subset['fchl_rep'].tolist(), train_subset['u0_atom'])\n",
    "    train_time = perf_counter() - train_time\n",
    "    \n",
    "    # Evaluate the model on the test set\n",
    "    test_time = perf_counter()\n",
    "    test_data['u0_atom_pred'] = model.predict(test_data['fchl_rep'].tolist())\n",
    "    test_time = perf_counter() - test_time\n",
    "    \n",
    "    # Store performance information\n",
    "    results.append({\n",
    "        'train_size': ts,\n",
    "        'train_time': train_time,\n",
    "        'test_time': test_time,\n",
    "        'alpha': model[1].best_params_['alpha'],\n",
    "        'mae': mean_absolute_error(test_data['u0_atom'], test_data['u0_atom_pred'])\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = pd.DataFrame(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing performance to the original JCP paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "jcp_results = pd.DataFrame([\n",
    "    [100, 0.2769171852851224], [1005.3012238214807, 0.06426615177922405],\n",
    "    [2502.606983490659, 0.03831336782897019], [5029.126239917856, 0.02719291345844958],\n",
    "    [10026.471083195112, 0.01905330215897317], [20148.744549550407, 0.013094302135827789]\n",
    "], columns=['size', 'mae'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'MAE (eV)')"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAC5CAYAAADnCem7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfGUlEQVR4nO3deXhU5fXA8e/JBoFA0IAgRNkNu4DsUUGLyg4CVnFDBdRWuiG06M8q2tJaaItaUcHKUlFB2RRBqLRFlKDsq4hQQAm7IGtAQji/P+6QhjCTTJLJ3FnO53nmcebOXc6LOXPvfe+7iKpijIk+MW4HYIxxhyW/MVHKkt+YKGXJb0yUsuQ3JkpZ8hsTpeLcDqA0VK5cWWvVquV2GMa4bvXq1d+pahVv30Vk8teqVYtVq1a5HYYxrhORb3x9Z5f9xkSpiDzz+zS2Ppw6eOny8lfAiG3Bj8cYF0XXmd9b4he03JgIFl1nfhMysrOzyczM5MyZM26HEhHKli1Lamoq8fHxfm9jyW9ckZmZSYUKFahVqxYi4nY4YU1VOXz4MJmZmdSuXdvv7UL+sl9EyovIVBF5XUTuKbUDWe/GoDpz5gwpKSmW+AEgIqSkpBT5KsqV5BeRSSJyUEQ25VveRUS2ish2ERnpWdwXmKmqQ4BepRbUP3rBgc2ltntzKUv8wCnOv6VbZ/4pQJe8C0QkFhgPdAUaAQNEpBGQCuz2rJZToqOWv8L78oQk2L8RXrsePhwGpw6X6DAmPMTGxtK8efPc165duwBYsWIFN954I2lpaTRo0IDBgweTlZXFlClTGDp06EX76NSpU26bklq1avHdd98FuxjF5so9v6ouFZFa+Ra3Abar6g4AEZkO9AYycX4A1lHSH6uCHudlHYElz8PKv8OmmdDpSWg9CGL9r0Ax4SUxMZF169ZdtOzAgQPccccdTJ8+nfbt26OqzJo1ixMnTrgTZCkKpXv+GvzvDA9O0tcAZgP9RORVYJ6vjUXkYRFZJSKrDh06VPSjl7scuo2Bn2RA9Zaw8DfwajpsX1z0fZmAm7t2D+nP/5vaI+eT/vy/mbt2T6kcZ/z48QwcOJD27dsDzuV0//79qVq1aqkcz02hVNvv7aZFVfUU8GBhG6vqRGAiQKtWrYpfe3dFA7hvDny9EBY9CdP6wTVd4NbRULlesXdrim/u2j08MXsjp7Odu749R0/zxOyNAPRpUaPY+z19+jTNmzcHoHbt2syZM4dNmzYxcOBAn9vMmDGDzz77LPfz9u3bi318t4VS8mcCV+X5nArsdSUSEUjrCnVvhi9eg0/GwivtoN2jcOMIKJvsSliR7M4Jyy9Z1qPZldzXvhZjFn6Vm/gXnM7OYdS8zfRpUYMjp87yk2mrL/p+xiPtCz2mt8v+QuO8805efvnl3M+dOnUq0vahJJQu+1cC9UWktogkAHcBHxRlByLSU0QmHjt2LDARxZWB9F/Az1bDtXdCxsvwt+tgzT/gfMnqHo3/9h3z/gjraFZ2wI/VuHFjVq9eXfiKEcCVM7+IvAN0AiqLSCbwjKq+ISJDgUVALDBJVYv07E1V5wHzWrVqNSSgAVeoCr3HQ+vB8NFI+OBnsOJ16PonqNkhoIeKVgWdqatXSmTP0dOXLK9RKRGAy8sn+HWm98fQoUNp06YN3bt3p23btgBMmzaNzp07B2T/ocSVM7+qDlDVK1U1XlVTVfUNz/IFqnqNqtZV1dFuxFag6i3goYXQ7w3IOgyTu8J7D8LR3YVva4ptxG1pJMbHXrQsMT6WEbelBfxYVatWZfr06QwfPpy0tDQaNmzIp59+SsWKFf3avlmzZqSmppKamsqwYcMCHl8gSSSN2y8iPYGe9erVG7JtWyn30jubBctehGUvOJ/Tf+ncIiSUK93jRogtW7bQsGFDv9efu3YPYxdtZe/R01SvlMiI29JKVNkXibz9m4rIalVt5W39iEr+C1q1aqVBG8zj6G5Y/AxsmgUVa8APJ+CH45euZ92GL1LU5DeFK2ryh1KFX3iqdBX0nwQPfgTlUrwnPli3YRNyIir5A17bXxQ1O8DDS4J/XGOKKaKSX1XnqerDyckuPYePiS18HWNCREQlf8jLOuJ2BMbksuQPpvFt4MsitVsyptRY8gear27DiZdDhSvh3fvgvQfgVPh0/YxEhw8fzu3KW61aNWrUqEHz5s2pVKkSjRo18rrN008/zeLFhXf0UlUqV67M999/D8C+ffsQkYv6BFSpUoXDh/3rOr5r1y6aNGni17pFEUpt+0ssz3N+94Io6HFeTrbTLuCTMbBzKXQbC437On0JTFClpKTktusfNWoUSUlJDB8+nF27dtGjRw+v2zz33HN+7VtEaNu2LcuXL6dbt25kZGTQokULMjIyuP7669m6dSuVK1cmJSWl0H3l5JReM/KIOvO7XuFXmNh4p2PQI0uhUk2Y+RDMuBdOHHA7stA2tj6MSr70NbZ+qRwuJyeHIUOG0LhxY2699VZOn3aaFj/wwAPMnDkTgJEjR9KoUSOaNWvG8OHDL9lHeno6GRkZAGRkZDBs2DCWL1+e+7lDhw6oKiNGjKBJkyY0bdqUGTNmALBkyRJuuukm7r77bpo2bXrRfnfs2EGLFi1YuXJlicsZUckfNq5oCIM+hs7PwraP4ZW2sH6GjSPoS5CHXN+2bRuPPfYYmzdvplKlSsyaNeui748cOcKcOXPYvHkzGzZs4KmnnrpkHx06dMhN/hUrVtCnTx9273aagWdkZJCens7s2bNZt24d69evZ/HixYwYMYJ9+/blbjN69Gi+/PLL3H1u3bqVfv36MXnyZFq3bl3ickbUZX9YiY2D638Jad3g/cdgzsOweTb0GAcVq7sdXXB9NNIZRq04Jnf3vrxaU+j6fLF2Wbt27dx+/tddd13u8F4XVKxYkbJlyzJ48GC6d+/u9TahTZs2rF27llOnTpGdnU1SUhJ16tRh+/btZGRk8PjjjzNhwgQGDBhAbGwsVatWpWPHjqxcuZKKFSvSpk2bi0biPXToEL1792bWrFk0bty4WOXKz878bqtyjdNZ6LY/wI5PYHw7WDvNrgJcVKZMmdz3sbGxnDt37qLv4+LiWLFiBf369WPu3Ll06dIl/y4oV64c9erVY9KkSbRs2RKAdu3asWDBAg4ePEhaWhoFNa0vX778RZ+Tk5O56qqrWLZsWUmKdnE5AranEBASFX7FERML7R9zRgx6f6hzJbB5DvR8EZJT3Y6u9BV2hh5VQB3Og/MDG4sfTp48SVZWFt26daNdu3b4+ntLT0/nhRdeYNSoUQC0b9+ee++9l3bt2iEi3HjjjUyYMIGBAwdy5MgRli5dytixY/nqq68u2VdCQgJz587ltttuIykpibvvvrvE5YioM3/IV/gVJqUuPDAfuo6FbzKcq4BVk+0qIMScOHGCHj160KxZMzp27Mi4ceO8rpeens6OHTtyxwNs2bIlmZmZdOjgjAFx++2306xZM6699lpuvvlmxowZQ7Vq1Xwet3z58nz44YeMGzeO999/v8TlsF59oerITmfQkF2fQu2O0OtvcFlNt6MKmCL16rMJVv1S1F59EXXZH1Eurw33fwCrJ8PHT8Mr7eGWZ6HVIIiJqAu2wlmCl4oo+ysKMzExztwBP10OV7WBBcNhak84ssPtyEwEsOQPB5WudoYT7/U32L/BmU/g81fh/Hm3IzNhLKKS39X+/KVNBFreDz/9HGpdDwtHOmMIfhe+48ZHYn2TW4rzb2kVfuFIFdZPd2YVOuPjhy7EK8N27txJhQoVbKbeALgwRfeJEycumaLbKvwijQg0HwB1b4K/+BjBNsSHDUtNTSUzM5NiTa1mLlG2bFlSU4vWJsSSP5xV8P1MONTFx8dfcpYywRVR9/zGGP9Z8hsTpSz5I9mqyW5HYEKYJX+48zVsWGwCfPhLWD4+qOGY8BFRFX5h26uvBOZ2XnLR3PXgzGP3fM80ev/3GVj0JJw95YwgZI/UTB4RdeYP+159xTB20Vavc9eP+Xgn9JsE1w6A/4x2phSLwDYdpvgi6swfjfZ6mbo6d3lsHPR+BeITnUlFz2ZB1zHR1zHIeFVo8otIKnAXcANQHTgNbALmAx+pqjUwd5GvuesTE2I5fTaHxIRY6P5XSCgPGX+D7Cynj4DNLhT1CjwFiMhkYBJwFvgTMAD4KbAY6AJ8JiI3lnaQxjdvc9fHxQhZZ3PoM34ZWWfPOff6t/wOOj0B696CWYPg3FmXIjahorAz/19UdZOX5ZuA2SKSAFwd+LCMvy7MUZ9/7vqUpARW7fqecgme/8Ui0GkkxJeDj38L2afhjqkQX9bF6I2bCuzYIyLDgemqmhm8kEou4jv2FNH63UeZ9vk3jOrVmPJl4mDl32H+484IQQPecW4JTEQqqGNPYTU/NYDlIrJURH4iIpUDH54pbet2H2Xmmkx6vfwZW/Ydh9aDoc9rzhBhb/b13TPQRLQCk19Vf4VzWf9boBmwQUQ+EpH7RaRCMAI0JTewQy3eGtSW42fO0Wf8Mt7+4lv02rug/2TYsxqm9oJT/s0bZyJHoc981PGJqv4EuAp4AfgVYHNMhZEO9Sqz4Oc30Kb25Tw5ZyMLNu6Hxn3grrfh4BaY0t2mDYsyfj/wFZGmwHPAeJza/ydLKyhTOqpUKMPUB9vw5zuupUsTpzvwD3V+BPe8B0e/hcld4Ohul6M0wVLYo776IvJbEfkSeBvIAm5V1baq+kIwAiyKiB7GK0BiYoT+16USGyMcPHGGm//8CW8erIXeN8e59J/cFQ7/1+0wTRAUduZfBJQF7lTVpqo6WlVDdujYaGzeWxJxMTHUr5rEb+duYuin8ZwcMMdpBDS5q3MrYCJaYRV+dVT1/1R1o4jUFJHOACKSaBV+4e/y8glMGtiaJ7o2YOHm/XR79zhbu04HBCZ3g73r3A7RlCK/7vlFZAgwE5jgWZQKzC2lmEwQxcQIj3Ssy7uPtOdcznle2hAHDy5wnv1P7QXffuF2iKaU+Fvh9xiQDhwHUNVtgI+O5CYcXVfzMhb84gb+cHtTSKnLgf5zyCmXAm/e7swebCKOv8n/g6rmNgYXkTjA+odGmErlEkguF4+q8tiHh+iT9RSnk1LhrTvg60Vuh2cCzN8uvZ+IyJNAoojcgtO5Z17phWXcJCL8X/eGDH37DDcceJyPUsZR+e0f43UokBCfH8D45u+ZfyRwCNgIPAIsAJ4qraCM+1pcfRkLfn4DLRvU5eZDj3tPfAj5+QGMb36d+T199l/3vEyUSC4Xz4T7rmNqRgp87HY0JtAKa+Qzz9NwJt7Ld3VE5DkReaj0wjNuExEeSLfJNSJRYWf+IcAw4AUROYJz6V8WqA1sB15W1fdLN0RjTGkoMPlVdT/wa+DXIlILuBJnGK+vVTWr9MMzYWHbx1D/FrejMEXkd8ceVd2lqstVdZ0lfhTyMT/AOY1B37oDPhkD5204x3Bio/ca/3h5nHfyh3OMmbeG32S/Rvn/jIY9a+D21yCxUvDjM0VmYzibYksqE8dz/dtQ/q430K5jyPn6n5x+5UY4sNnt0IwfCqvtr1jAd0EZuNPzVOENEZkZjOOZYhDh+yYP8qvE0Rw/fozsCTdzfoP97wp1hZ35l1x4IyL/yvfd3MJ2LiKTROSgiGzKt7yLiGwVke0iMrKgfajqDlUdVNixjLsuL5/A878awqvXTGLtuZrEzB5E1gcjICfb7dCMD4Ulf96GXZcX8J0vU3DG9//fRiKxOKMBdQUaAQNEpJGINBWRD/O9rPNQGCmXEMeoe37E7p4z+Mf5rpRbMxGd2suGBwtRhSW/+njv7fOlG6suBY7kW9wG2O45o58FpgO9VXWjqvbI97K2o2GoX+vatH/sdf57wzhk71p0YkfOffO522GZfApL/itEZJiIPJ7n/YXPVYp5zBpA3oHiMj3LvBKRFBF5DWghIk8UsN7DIrJKRFYdOnSomKGZQKlftQJ1f/QQDF7MsexYmNydY0tfsclCQ0hhyf86UAFIyvP+wue/F/OY3m4XfP5FqOphVX1UVeuq6h8LWG+iqrZS1VZVqhT3d8kEXLUmZNz8Hsu0Kcn/foJ9Ux90Zgsyriushd+zvr4TkdbFPGYmzhDgF6QCe4u5LxMGurVpxI5aH/DmpN9wz8532D9uMykPvUt8Zesz4KYiPef3VMw9JyLbgFeLecyVQH0Rqe2Z6+8u4INi7it/fDZ6b4iqc0VF7nj8Zd6sM4bEU3uQ1zvC9sVuhxXVCk1+z8CdI0VkPfAmzkAet/ia/yvftu8Ay4E0EckUkUGqeg4YijMy8BbgXVUNSKsQG703tJWNj2XgwIc5fv/HxFW6Cqb157sFv7dmwS4pbKLODCAZp0Z+uqpuE5GdqhrS12s2UWcYOJvFwbcf4YpdH7C10g3UHvwmCUmXuR1VxCnJRJ2HcCr4qvK/2v2Qra61y/4wklCO5HsmsyD1l9T5PoND4zqwf9sat6OKKoWN298baAqsAZ4VkZ3AZSLSJhjBFZVd9oeXMvFxdBv8LKs7/YOEnCwqvtWF9QsnuR1W1PBnos5jqjpJVW8B2gHP4AzuYZO6mYBod1MPfnhoCbvi6nLt57+ChU9as+AgKPCev8ANRWqq6jcBjicg7J4/PJ07e4bYxb9FVkzkvMQRo+cuXclGCy6Sgu75C3zOLyKFPYLrVeyoSoGI9AR61qtXz+1QTDHEJZSFbmPJvrIl8e8/6n0lGy04YAobzKM9TlPcd4Av8K8zj2tUdR4wr1WrVkPcjsUUX3yLAeAr+U3AFJb81YBbgAHA3cB84J1APZc3xrinsNr+HFVdqKoDcSr7tgNLRORnQYnOGC9yzofs0+aw4k8LvzIi0heYhjNh50vA7NIOrDjsOX906PzXT1j+38NuhxH2ChvGayqQAbQEnlXV1qr6O1XdE5Toisie80cQH6MFZydWRlUZ8PrnPPP+JrLOenkiYPxSWPPe88Apz8e8KwqgqupzjD832aO+yHb6bA5jFn3FlIxdXHVZOcb0b0a7OiluhxWSit28V1VjVLWC51Uxz6tCqCa+iXyJCbE807Mx04e0QwRWf/O92yGFJRu334SttnVS+OgXN5AQ65zDln59iDJxMbS1qwC/RNS4/VbhF33KJcQRFxuDqvLSv7Zx58TPGfXBZqsL8ENEJb9V+EUvEeEfg9owsH1NpmTsouuLn7JiZ/6xY01eEZX8JrqVS4jj2d5NeGdIO86rcufE5azffdTtsEKW3fObiNO+bgoLf3Ejc9buoVmqcxV4+OQPpCSVcTmy0GLJbyJS+TJx3NuuJgDfHs6i64tLubP11TSoVoEX/7WNvUdPU71SIiNuS6NPC58jx0c0S34T8VKSEujbMpVJy3Y6DVQ8y/ccPc0TszcCROUPQETd81ttv/GmfJk4ftenCZWTEi4Zg+50dg5jF211JS63RVTyW22/Kcjhk2e9Lt97NDonEYmo5DemINUrJXpdHhcrrP02+loJWvKbqDHitjQS42MvWpYQG0PZuBhufyWDYe+u48DxMy5FF3yW/CZq9GlRgz/2bUqNSokIUKNSImP6N2P5k535aae6fLh+Hzf9eQmTPtvpdqhBYbX9Jqr0aVHDa83+r7s04K7WV/P7+V/ywzlnBqELPV5FQnr0umKz5DfG4+qUcky8v1Vu0s/bsI/pK77l6Z6NaFAt8jqx2mW/MflcONPnnD/P5r3H6fbipzz9/ia+P+X9aUG4Kva4/aEoz9DdQ7Zts7HdTcl9f+osLyz+mmlffEtSmThG9WqEIIxdtDUsWgkWNJhHRCX/BTaSjwm0rftP8NyHm6lXJYl3V2VyOjsn97vE+Fj+2LdpSP4AlGSiTmMMkFatAtMGtWXxlgMXJT6EbytBS35j/CQi7D3qvR1AOLYStOQ3pgh8tRKMjQm/VoKW/MYUga9WgsmJcWHXHsCe8xtTBBcq9fLX9vdodiVxnoFE/7BgC9WTy3Jvu5q5y0KRJb8xReSrlSDAuZzzbN1/golLd/De6kx+16cJLa++LMgR+id0f5aMCUNxsTFMebA14+9uyeGTZ+n7SgYjZ20IyQZClvzGBJiI0L3ZlSx+vCNDbqjN/A37yMr3eDAUWCMfY0rZsdPZJCfGo6r8YcEWksvF884Xu4PSQrCgRj4Rdc+fp3mv26EYkys5MR6AvcfO8PYX33Lq7P+uAtwcRzCiLvttGC8TympUSqSi54cgL7daCEZU8hsT6vYfC50Wgpb8xgSRrxaCVyaXZdbqTHLOB68OzpLfmCDy1kIwMT6WmxpW4fH31tP/tQy+2n88KLFY8hsTRN7GEfxj36b8vndTXrizOd8czqLHS58xdtFXnCnlx4P2qM+YEHLk1FlGz9/CrDWZdGlcjdfuu65E+4uaR33GhLvLyyfwlx9fS9+WNXIfER4/k01OjnJZ+YSAHssu+40JQen1KtOkhvPI+k8ffcWP/voJc9fuIZBX6nbmNybE3duuJpv3HueXM9Yxa00mna6pwqRlu0rcQtDu+Y0JAznnlWmff8Po+V9yNufinC1oDEEbw8+YMBcbIwzsUMvrfX9xWwha8hsTRg4e/8Hr8uK0ELTkNyaM+Goh6Gt5QSz5jQkjvloIjrgtrcj7stp+Y8KIrzEEi1Pbb8lvTJgpaAzBorDLfmOiVEQ+5xeRY8CFmTqTgWNe3lcGvivBYfLuqzjrePsu/7KCPnt7H4wyFbZeJJbL3+W+ypH3c7D/BmuqahWv36hqxL2AiX68XxWoYxRnHW/f5V9W0Gdv74NRpmgsl7/LfZUjX1mC+jdY0CtSL/vn+fE+kMcozjrevsu/rKDPpVEuf/cTbeXyd3lBsc/zsbwkSrSviLzs94eIrFIfzR7DVSSWCaxcpSVSz/z+mOh2AKUgEssEVq5SEbVnfmOiXTSf+Y2Japb8xkQpS35jopQlPyAidUTkDRGZ6XYsgSQifUTkdRF5X0RudTueQBGRhiLymojMFJGfuB1PoIhIeRFZLSI9gnG8iE1+EZkkIgdFZFO+5V1EZKuIbBeRkQCqukNVB7kTadEUsVxzVXUI8ABwpwvh+q2I5dqiqo8CPwZC9hFgUcrk8Rvg3WDFF7HJD0wBuuRdICKxwHigK9AIGCAijYIfWolMoejlesrzfSibQhHKJSK9gM+AfwU3zCKZgp9lEpHOwJfAgWAFF7HJr6pLgSP5FrcBtnvO9GeB6UDvoAdXAkUplzj+BHykqmuCHWtRFPX/l6p+oKodgHuCG6n/ilimm4B2wN3AEBEp9dyMti69NYDdeT5nAm1FJAUYDbQQkSdU9Y+uRFd8XssF/AzoDCSLSD1Vfc2N4ErA1/+vTkBfoAywIPhhlYjXMqnqUAAReQD4TlXPl3Yg0Zb84mWZquph4NFgBxNAvsr1EvBSsIMJIF/lWgIsCW4oAeO1TLlvVKcEK5CIvez3IRO4Ks/nVGCvS7EEkpUrfIRMmaIt+VcC9UWktogkAHcBH7gcUyBYucJHyJQpYpNfRN4BlgNpIpIpIoNU9RwwFFgEbAHeVdXNbsZZVFau8ClXqJfJOvYYE6Ui9sxvjCmYJb8xUcqS35goZclvTJSy5DcmSlnyGxOlLPnDkIikiMg6z2u/iOzJ8/nSCdwv3raViBTa5FdEMgIUazkReUtENorIJhH5TESSCtnmyQK+e8izrw2e/fX2LH/O0zPO+Mme84c5ERkFnFTVP+dZFudpTOI6EXkCqKKqwzyf04Bdqup9onlnnZOqeskPhIikAp8ALVX1mOdHpIqq7iyl8COanfkjhIhMEZG/ish/gD+JSBsRyRCRtZ7/pnnW6yQiH3rej/IMOLFERHaIyM/z7O9knvWXiDNqzlees7h4vuvmWfaZiLx0Yb/5XAnsufBBVbdeSHwRuVdEVniuWCaISKyIPA8kepa9lW9fVwAngJOefZ28kPie8vf3XNlcuAraKCLq+b6uiCwUZ6ScT0WkQQD+2cNatPXqi3TXAJ1VNUdEKgI3quo5z+XwH4B+XrZpgNOXvAKwVUReVdXsfOu0ABrjdEBZBqSLyCpggucYOz1NWb2ZBPxTRPrjDLwxVVW3iUhDnNGF0lU1W0ReAe5R1ZEiMlRVm3vZ13qcwS52isi/gNmqetGsNaq6CmgOICJjgYWeryYCj3qO3RZ4BbjZR8xRwZI/srynqjme98nAVBGpj9NlNN7HNvM9Z+IfROQgUBWn51leK1Q1E0BE1gG1cM6+O/Jccr8DPJx/56q6TkTqALfijC2wUkTaAz8CrvN8BkgEDhZUOM+PWhegtWf7cSJynaqOyr+uiPwYaAnc6rk96AC85zkWOGMBRDVL/shyKs/73wH/UdXbRaQWvvu/5733zsH734S3dbz1S/dKVU8Cs4HZInIe6AacxbkKeMLf/Xj2pcAKYIWIfAxMBkblXUdEGgPP4lyV5IgzKs5RH1cTUcvu+SNXMv+7136gFPb/FVDH88MCPgYIFZF0EbnM8z4BZ9y6b3BuAfqLyBWe7y4XkZqezbJF5JIrFRGpLiIt8yxq7tlX3nWScYbGul9VDwGo6nGcW4U7POuIiFxb9CJHFjvzR64xOJf9w4B/B3rnqnpaRH4KLBSR73DOxt7UBV71VBLGAPOBWaqqIvIUTn1ADJANPIaTzBOBDSKyRlXzjtEXD/xZRKoDZ4BDXDoCUx+gJvD6hUt8zxn/Hk8cT3n2Mx2nDiFq2aM+U2wikqSqJz2JPR7Ypqrj3I7L+Mcu+01JDPFUAG7Guc2Y4G44pijszG9MlLIzvzFRypLfmChlyW9MlLLkNyZKWfIbE6Us+Y2JUv8PRQcOtiqo4usAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.loglog(jcp_results['size'], jcp_results['mae'], 'o--', label='FCHL')\n",
    "ax.loglog(results['train_size'], results['mae'] * 27.2114, 's-', label='This Work')\n",
    "\n",
    "fig.set_size_inches(3.5, 2.5)\n",
    "ax.legend()\n",
    "ax.set_xlabel('Training Set Size')\n",
    "ax.set_ylabel('MAE (eV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our implementation agrees well with the original paper!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAADQCAYAAADS8b86AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJcklEQVR4nO3dd5wU9fnA8c9zjTvaHVW4o/cuTRQxCoqBqCgqMfYuMYlRf9aQaDSJxiSYBBMrNow1QhRRiRhRxE7viKAgcEc/7mh3XHt+f8wcLHe7V3d3tjzv12tf3M7Mzj47zHf3me98i6gqxhhjjDEmdBK8DsAYY4wxJtZZwmWMMcYYE2KWcBljjDHGhJglXMYYY4wxIWYJlzHGGGNMiFnCZYwxxhgTYpZwGWNCRkT+KyJXBXtbL4nI+SKyRUQOiMggr+Mx3hOR40RkvojsF5G/eh2P18TxvIjsFZEFXscTKSzhqgUR2SQioyssu1pEPvUqpmgkIioi3byOw/jnJhLljzIRKfB5fllt9qWqP1LVF4K9bW2IyEj3cxxwfxDXicg19djlw8BNqtpYVZcGK04TWfx931dhIrAbaKqqt4cwrJARkWkiUuSWk1wR+Z+I9Krj7k4BzgTaqeqwIIYZ1SzhijLulUNQ/99EJDGY+wslEUnyOoZY5yYSjVW1MbAZGOez7OXy7aLs/yLH/TxNgbuBp0WkT2124PN5OwKr6xJENJU1UysdgTVah5HEI6wc/cUtJ+2AncC02u7A/TwdgU2qerCOr49JlnAFiYjcKSL/qbDsnyIyxf17nog8JCILRCRfRN4SkeY+254kIp+LSJ6ILBeRkT7r5onIgyLyGXAI6FKD/U0Xke3uuvki0tdn3TQReUJEZovIQWCUiJwtIktFZJ97u+R+n+07ubVS17jr9orIjSJygoiscGN+tMJnv1ZE1rrbzhGRju7y+e4my90rqZ+4y88RkWXuvj4XkQE++9okIneLyArgoIgkuc+zfWoszqjjf52pIbemaKt77LcDz4tIMxF5R0R2uf/X74hIO5/XzBOR692/rxaRT0XkYXfbjSLyozpu21mO3sL5QEQeE5GXqvsM6pgJ7AX6iEiCiPxKRL4VkT0i8np5OfI5768Tkc3AJyJyAEjEOX+/dbfr7caeJyKrReRcnzj9lbVN4nxfrBCRgyLyrDi3pP7r83ma+eyjurL8mIi86772KxHp6rO+rzg1FbkiskNEfu0uD/i5TWVVnY8iMg24CrjL/U4bXcvz6kN3ud/vTHedivOdu95d/5iIiM/6G9zX7heRNSIy2F2eKSL/ccvnRhG5uSafV1UPAa8A/arbj4jcLyIzROQlEdkHXAc8Awx3j8fvfGLc4J6Ls0Qks8Ln+4WIrAfWy9HvmrtEZKeIbBOR8SJyloh84+7j1z6vHyYiX7hlcJuIPCoiKV4dv6oOrD1q+AA2AaMrLLsa+BRoCxwEMtzlSThXCEPc5/OAbPcEbgT8B3jJXZcF7AHOwkmCz3Sft/J57Wagr7vf5Kr2577mWqAJ0ACYAizzWTcNyAdGuO+XCowE+rvPBwA7gPHu9p0ABZ50t/0hUAjMBFq78e8ETnO3Hw9sAHq78d4DfO7z/gp083k+2H39iTg/Zle5x7qBz3FfBrQH0oCewBYg0ye+rl6fH7H48D3n3XOkBPize16lAS2AC4GG7vk2HZjp8/p5wPU+ZaUYuMH9f/4ZkANIHbb9AufWXgrO7Yt9vud/hc8wEtjq/p0AnO/uuydwK/AlzhV9A+Ap4NUK5/2/cMpYWsXzF6csbgB+7cZyOrAf6FlFWdvkvudxHC07S4BBbgwfAvfVoiznAsNwytrLwGvuuibANuB2932bACe66wJ+bnv4PferOx+nAQ/4vLZW5xU1+858B8gAOgC7gLHuuh/j/BacAAjQDaeGKQFYDPzWPTe7AN8BYwJ83iOfAWiMk3B9Ut1+gPvdYzPe3TbNPV6f+uz7dJxbroPd4/FPYH6Fz/c/oLn7+pE43zW/xSljN7if+RX3PO6L8xvUxX39EOAk99h1AtYCt4bz+NXonPL6pI6mB04BPADk+TwOlZ9YwH+BG9y/z8GpYi5/7TzgTz7P+wBFOIX3buDFCu81B7jK57W/r7A+4P78xJ3hnnDpPgXrX9V81inA392/O7mvz/JZvwf4ic/z/5Sf4O5xuM5nXYJ7nDrq0ZPfN+F6AvhDhfdfx9EEbhNwrc+6bjg/UqOBZK/Pi1h+UDnhKgJSq9h+ILDX5/k8jk2iNvisa+ieC21qsy3OF2YJ0NBn/UtUnXCV4ZTXXJzk/WJ33VrgDJ9t2+L8eJR/cSvul7rPNr4J1w+A7UCCz/pXgfvdvyuVNfeYXubz/D/AEz7Pf4lP0lrhtRlULsvP+Kw/C/ja/fsSYGmA/QT83F6fc5HyoHLCVdW5O41jE65anVfU7DvzFJ/1rwO/cv+eA9ziJ/4Tgc0Vlk0Cng/weafhJDF57jk9C+ha3X5wEq75FdZfzbEJ17M4tyvLnzd2j0cnn893us/6kUAB7u8ZTpKluBcM7rLFuJUCfj7LrcCbPs9Dfvxq8ojZe6UhNF5VPyh/IiJXA9e7T1/AufJ5GrgceLHCa7f4/P09TubeEieb/rGIjPNZnwx8FOC1Ve5PRHYDD+Jk7q1wfmxw3yvf3/5E5ETgTzg1Zik4VyHTK7zfDp+/C/w8b+z+3RF4RI7trSM4V/Pf+/kcHYGrROSXPstSgEyf50fiVdUNInIrTkHvKyJzgNtUNcfPvk1w7VLVwvInItIQ+DswFii/DdZERBJVtdTP67eX/6Gqh9xa/cZ+tqtq25ZArjq3PcptwakBDSRHVdv5Wd4ReFNEynyWleLUPvnuO5BMYIuq+r7+e5xzvarX16gsidPmq7qyvN3ntYc4ejzbA98GiLuqz50d4DXxrjbnbm3Pq5p8Z9b2/7kjkCkieT7LEnFqrQJ5WFXv8V0gIkNqsJ+qygg45WRJ+RNVPSAie3A+36YA+9jj8x1S4P4bqJz0AP4GDMVJhpNwEjJf4Th+VbI2XME1ExggIv1warherrDe9wehA06GvxvnRHtRVTN8Ho1U9U8+26uf9wu0v0uB83BqgNJxrqjAKcCB9vcKzhVNe1VNx7l9KNTNFuCnFT5Pmqp+XsX2D1bYvqGqvhooXlV9RVVPwSkUinOby4RexfPmdpxbcyeqalPgVHd5Xc+dmtgGNHeTvXJVJVtV2QL8qMK5l6qqvkmHv7JXLgdoL8d2ZOnAsUlLVa+vTk3KciBbcGooAq2r7nObuqvteVXb78yK7+Xv/3kLsLHCPpuo6ll1+CzV7ae6czwH57saABFphNMcIVjl5Anga6C7+z30a2r+HRTq43eEJVxB5F75z8BJXhao6uYKm1wuIn3cH4rfAzPcDP4lYJyIjBGRRBFJdRsN+rsir8n+mgCHcW77NQT+WIPwm+DUGhSKyDCcL/q6ehKYJG7jXhFJF5Ef+6zfgXM/vNzTwI0icqI4GonTiL+Jv52LSE8ROV1EGuBUgRfgXD2a8GuCc/zzxGkUfF+o31BVvwcWAfeLSIqIDAfGVfOyQJ4EHpSjnTpaich5tXj9VzhtN+8SkWRxOruMA16rYzwV1aUsl3sHaCMit4pIAxFp4tZkQ/0/t6labY9vdd+ZVXkGuENEhrjfn93c910A7BOnk0ua+9vST0ROqOVnCcZ+XgGuEZGB7vf2H4GvVHVTLWMJpAlOO84D4gxl8bNavDbUx+8IS7iC7wWcxucVbyfiLpuGU7WZCtwMoKpbcK5if43TmG8LcCfV///43R9OY8zvca4e1uA03qzOz4Hfi8h+nEaCr9fgNX6p6ps4NU6vidNrZRXwI59N7gdeEKdHyUWquginUeSjOL3HNuC0AQikAc7tz904n701zrEz4TcFp5Hrbpzz7L0wve9lwHCcROQB4N84iUltPYJTs/u+e+5/idN2o0ZUtQg4F+f83g08Dlypql/XIRZ/6lKWy2Pbj9MBZxxOOVkPjHJX1+tzm2rV6vjW4DszIFWdjnPb+RWcDhszgebuxfc4nHaVG3HOz2dwakprLBj7UdW5wL047RW34dQoXVybOKpxB04lwX6cC/h/1yK2kB4/X+U9LEyQiEgHnKrNNqq6z2f5PJxGvc8E6X2Cuj9jopmI/BunsXjIa9iMMaYurIYriNx2HLfhdMveV932xpi6EWcMuK7ijHc0FqeGeKbHYRljTEDWSzFI3EaAO3Cq/8d6HI4xsa4N8AZOw9utwM/UptkxxkQwu6VojDEmpojIeOBsnPadj6nq+95GZIzdUjTGGBMFROQ5caZ5WVVh+VhxpvfaICK/AlDVmap6A07nm594EK4xlcR0DVfLli21U6dOXodhotTixYt3q2orr+MINisXpj68KhcicirOTB//UtXyOf4SgW9wemNuBRYCl6jqGnf9X4GXVXWJ/706rEyY+qpJuYjJNlzijNg+rlu3bixatMjrcEyUEhF/o+JHlLrcOunUqZOVC1NnXpULVZ0vIp0qLB6GM+XOdwAi8hpwnoisxRk65r/VJVtgZcLUX03KRUzeUlTVt1V1Ynp6nYfLMMYzduvEmBrL4tgpYba6y36JMzr/BBG50d8LRWSiiCwSkUW7du0KfaQm7sVkwmVMlJtGhZ6u7q2Tx3AGQ+wDXCIifXw2ucddb0w88Td9i6rqP1R1iKreqKpP+nuhqk5V1aGqOrRVq5hrOWAiUEzeUoxUM5dmM3nOOnLyCsjMSOPOMT0ZPyir+hfGaBzGv1DcOhGRicBEgA4dOoQkbhP9ovC7YSvHzqPZDmfePmOCJljlwmq4wmTm0mwmvbGS7LwCFMjOK2DSGyuZuTS8c8VGShym1up86wTsat5UL0q/GxYC3UWks4ik4EwXM6umLxaRcSIyNT8/P2QBmugWzHJhNVxhMnnOOgqKj51fuaC4lDumL+eJed8C8MTlg+nSqjFvLcvm8Y++rbSPadeeQNv0NP69cDPPfbqp0vrXfzqc9IbJPP/ZRl5bsKXS+lm/HBEwjslz1kX6lWy8C3jrBPhHuIMxsSfSvxtE5FVgJNBSRLYC96nqsyJyEzAHSASeU9XVNd2nqr4NvD106NAbQhGziX7BLBeWcIVJTl6B3+UlZUrnlo0ASE50KhybpiUfWeYrKcFZnx5gvbua5o1S/K4XJGAcgZabiFHvWye+vXeNqSjSvxtU9ZIAy2cDs8McjokTwSwXUZVwudPnzMe5snnH63hqIzMjjWw//0FZGWk8ecWQY5aN6tmaUT1bB9zX2H5tGduvbcD15w3M4ryB/jPvQHFkZqQF3J+JCEdunQDZOLdOLq3NDuxq3lTFvhuMqaxF4xR2HyiqtLwu5SIsbbhq0829GncDr4cmytC6c0xPEircFEpLTuTOMT3DHkdacuIxyxKEsMdhAnNvnXwB9BSRrSJynaqWAOW3TtYCr9fm1om7X2uvYgK6dXT3SvetvfiOCicrE6Y615/SOWi/3WEZab42IwTj3Id/qMIurgUGAC2BVGB3TWq4hg4dqpEymJ2q0u++OZQpFBaXRkwvxZaNG/CzkV249pQuYY8j0onIYlUd6nUcwRZJ5cJElhmLtvDX/33D9vzCgN9RsVgurEyYikpKy0hMEESkRr0Ua1IuwnJLsTbd3FX1IeCcivsQkVFAI5wxiApEZLaqlvnZLiK7v3+3+yAHi0r50wX9uXiYt3GNH5RV6WQpK1PWbNtHvywbLNaYeLNsSx49jmvMhKHtmTC0ffUvMCaGqSp3zlhBg6QEHrqgv9/fzLrwcliIQN3c/VLV36jqrcArwNP+ki13u6nA74AlKSkpwYu2nhZszAVgWOfmHkfi3+PzNnDB45+zfEue16EYY8Io92ARVz+/gLv/s9LrUIyJCM9+upE3l2aTlZGGiL8O4nXjZcLlt5t7dS9S1WnV3U6MxKl9VOH49hl+ew9GgstO7EirJg34+ctLyDtUuYGgiX7WXsX4M3nO1+wvLOGmUdZ71ZhP1+/mj7PXMrZvG34R5DLhZcIVshGCI/GH5dITO/DWL0YENVsOpmaNUnjsssHs3F/I//17GWVloW/bZ8IrEi9EjLeWbt7Lawu3cO2ITvRs08TrcMIuEn8rjHe25B7ipleX0K11Yx6+6HgSKraWrycvE656jRBclUj7YSktU8LROaG+BrbP4N5z+vDRul088XHlgVeNMbGjtEy5961VtG7SgFtG9/A6HE9E2m+F8dbm3EM0Skli6hVDadwg+E3cwzUsREi6uVfxfhF11fLOihyG/XEuW3IPeR1Kta44qSNXn9yJoR2beR2KCbJIKxd1de6559KsWTMmTJgQcJuHH36Yvn370q9fP1566aUwRhc98guKaZSSxD1n9wnJj4sJr5qUi9mzZ9OzZ0+6d+/O448/HsboosOIbi356I6RdApR05+wDAvhlUjp6vubN1cya1kOy+77IYlBrqIMteLSsiMj4MebWOz+DpFTLurqo48+4sCBA7zwwgvMmDGj0vqVK1dy1VVX8fnnn6OqjBw5kjlz5pCRkRH+YCNc+fd/bZo6xGK5iPYyAdWXi5KSEnr37s2HH35IixYtGDp0KHPnzqVt28CDaMeLqfO/JSkhgWtGdKpzs5+alIuY/CWNtCv5BRtzGdyxWdQlW/+cu57LnvmK4lK/HUKN8cSoUaNo0iRwe6O1a9cyfPhwUlNTSUtLY+DAgbz33nthjDDyvfLVZrbnFyIiEduu1NROdeViwYIF9OnTh/bt29OwYUPOP/983nknqiZsCYl563by0H+/ZsnmvSF/r5hMuCLpvvzeg0Ws33kgYoeDqEqHFg1ZsDGXh+es8zoUE2HuuOMOxo4dG7T9nX766YgIzz///DHLX3nlFVJTUykuLq7xvvr378+8efPIy8tj7969fPTRR2RnZwct1mi3aFMuv35zJS99+b3XoXgu2BfnkVwucnJyaN/+aD+1du3axX252Lj7IDe/upRebZrylwkDQn7xEZMJVyRZuCmyx9+qynkDs7jipI48Nf875qze7nU4pp6C+eOycOFChg0bFoSoHEuXLiUzM5Pp06cfs3zJkiX079+f5OTkGu+rd+/e3HrrrZx++ulMmDCB4cOHk5RkbZTAGT373rdWk5meys9HdfU6HM8F++I8ksuFv+ZD8Vy7eeBwCRP/tYjEBGHqFUNomBL674iYTLgi6ZZiu2YNuWZEJwa08762rS7uOac3A9qlc8f05Xy/56DX4Zh6CMaPS3FxMSkpKcyfP58//OEPiAh9+/atV1zffvsteXl53HPPPXzwwQfk5eUdWbd48WKGDBkS+MUB3HDDDSxZsoS5c+eSnJxMt242xhTAi19+z9pt+7j3nD5h+YGJF9FQLrKystiy5ehY41u3biUzM7NeMUazT9fvYtOegzx66WDaN28YlveMyYQrkm4p9slsyn3j+tIgKbH6jSNQg6REHrt0MKnJiazK3ud1OMZjiYmJfPHFFwB89dVXbNu2jU8//bTSdvfffz89e/ZkwoQJ5Obm8v777zN48GD+9re/Vdp28eLFJCUlcc0119CuXTveeustwLkiX7p0aZ0Srp07dwKwbt06FixYwJgxY2q9j1izc38hf3v/G07t0Yqx/dp4HU5MiYZyMWzYMFavXs2WLVsoKCjgzTff5Oyzz67Dp40NY/u1Zd6doxjRrWXY3tMucUKosLiUddv30zezKUlR3NOvffOGzL9zFGkp0Zk0muBJSEhg27ZtNGnShBNOOMHvLYkPPviAbdu2sXjxYqZOncr555/P3r17eeWVV+jXr1+l7RcvXkzfvn1JTU3lwgsvZMaMGVx11VVs2LCB/Pz8Sj8so0ePZvny5Rw8eJB27doxffp0hg8fzsCBA1m2bBkA48ePJy8vj0aNGvH888/bLUUgKSGBs/q35caRXeP6VlIoREO5SEpK4u9//ztnnHEGZWVl3HLLLXFZwzV37Q4aJCVySveWZGWkhfW9Y/JbSETGAeO8vo2waNNeLn/2K164dhin9WjlaSz1VZ5svbMih+LSMs4f1M7jiIxXli5dyvHHHx/wR3vp0qVcfvnlNG7cmNtuu41Zs2bx85//3O+PChx7e+TCCy/ktNNOY9++fSxevJiUlJRKr/vggw/87qc82QL4/PPP6/DJYlvzRin8ecIAr8OIWdFQLsaNG8e4cePq8Oliw/od+7n51aX0aNOEEd1ahP3CI3qrXaoQKbcUF2zcQ4LA4A4ZnsYRLKrKK19t5lf/WcmaHLu9GG2C1bZx2bJlDBo0KOD6Xr168e677wIwd+5cCgsLmTJlCps2bfK7ve/tkRNPPJHWrVsza9YslixZQr9+/YikSeijUXFpGbe9vozVOd63aY00wWzva+UisuUXFDPxxcWkpSTy+GWDPanljcmEK1Is2JRLn8ymNEmteU+SSCYiPHLxINLTkvn5y4vZV1jzLsnGe8G6EFm+fDkDBgSuKRk3bhyFhYV06NCB3/72t7z++uv85S9/YdSoUTz88MPHbLtx40Zyc3OP/LCICBdccAEzZsyoc4N5c6xpn23ijSXZZO8t8DqUiBPMi3MrF5GrtEy55bWlbMk9xBOXD6FtenhvJZaLyVuKkaCopIylm/O47MSOXocSVK2aNODRSwdzydNfcveMFZ5dKRjvlJSU8PXXX5OTk0PDhg39juA+ZcoUpkyZcuR5hw4dOPfccyttV94w+Pjjjz+ybMKECfzwhz8kJSWFiy66KBQfIW5szy9kygffcHqv1pzZ5zivw4lpVi4i17srtzFv3S4eGN+PEzp5N0RTTNZwRcKwECuz8zhcUsawzrE3J+Gwzs25e2xP/rtqO59t2ON1OCbMHnzwQV577TXatWvHpEmT6rWvxYsX06dPH1JTU48sGzFiBBkZGezbt4/BgwfXN9y49uDstRSXKfeN62MXRiFm5SJyjRvQluevOYHLTuzgaRw2l2KIFBaXsvj7vfTLSic9LTZuKfpSVT7/dk9Yu9SGWyzOGQexMW+cqd6X3+3h4qlfcuvo7tw6ukfQ9huL5cLKRGz6Zsd+UhITQjYZta+4nUsxEqQmJzKiW8uYTLbAaVNQnmytydnHngOHPY7IGONrYPsM7jm7NzeeZiPKm/iz92AR172wkIkvLqKsLDIqlizhCoHSMmXKB9+wbvt+r0MJuf2FxVw89Qtu/fcySiPkpDYm3qkqqcmJXP+DLqQm2/h5Jr6UlJbxy1eXsiP/MH+6cAAJCZFxO90SrhBYt30/Uz5YHxfdsJukJjPprN58sn43//xwvdfhmCpEQttGE3o5eQWMmTKfxd/neh1KxLMyEZv+/N7XfLphNw+M78fgDpHTjjomEy6vC1H5hNVe9oYIp4tPaM+Fg9vxyNz1fPzNLq/DMQFEyvh0JrQeeHcNm3MP0bpJavUbxzkrE7Hnf2t28PQnG7lyeEcuOqG91+EcIyYTLq8L0YJNubRNT6VdM2/G+gg3EeGB8f3oeVwTbn1tKdvzC70OyZi4NP+bXcxeuZ2bRnUL24S8xkSSU7q15M4xPbn3nD5eh1KJjcMVZKrKgo25nNw1/NMGeKl89N4Zi7fSorGNgGxMuB0uKeW+Wavp3LIRN5zaxetwjAmr3INFJCcKTVKT+cUob6f1CyQma7i8tGv/YQqKSuPmdqKvLq0ac9fYXiQnJlBQVOp1OMbElbeW5bBx90HuP7cvDZKsobyJH8WlZdz40mIuefrLiO68ZTVcQda6aSrLfnsmJRH8nx5qm3Yf5JKnv+Ses/tw9oC2XodjTFz48ZB2tGuWxsldY3dsPGP8eeCdNSzYmMuUnwwkMUJ6JPpjNVwhkJSYENddsTMz0mibnspdM5bz7a4DXodjTMzLO1SEiFiyZeLO6wu38MIX33PDDzozflCW1+FUKWoSLhEZKSKfiMiTIjLS63gCuW7aQl5dsNnrMDyVkpTAo5cOpkFyIj97aTGHikq8DsmYmPXh1zs45c8fsWJrntehGBNWSzfv5Z6ZqzilW0vuHtvL63CqFZaES0SeE5GdIrKqwvKxIrJORDaIyK+q2Y0CB4BUYGuoYq2PHfsKmfv1Tg4UWoKRmZHGlJ8MZP3OA9zz5ipieQopY7xSWOw0lD+uaQN6tWnqdTjGhNVxTVMZ3ac1/7xkEEmJkV9/FK42XNOAR4F/lS8QkUTgMeBMnARqoYjMAhKBhyq8/lrgE1X9WESOA/4GXBaGuGtlwUZn/K1hneOvwbw/p/ZoxS1ndGfx93s5XFIWN7dZRaQdcDHwAyATKABWAe8C/1XVMo/iGgeM69YtMnvwmNp7Yt63bMkt4JXrTyQlKfJ/cIwJhqKSMpIShMyMNB6/bIjX4dRYWBIuVZ0vIp0qLB4GbFDV7wBE5DXgPFV9CDinit3tBRqEJNB6Wrgpl4YpifTNtCvNcr88vTtARDdkDCYReR7IAt4B/gzsxKmV7QGMBX4jIr9S1fnhjk1V3wbeHjp06A3hfm8TfN/vOcgTH3/LuOMzOTmGJ5EPJbsIiT6qym/fWsWu/YeZeuXQqPpt8bKXYhawxef5VuDEQBuLyAXAGCADp7Ys0HYTgYkAHTp0CEacNbZgYy6DOzSLiqrNcCkvDNvzC3nov2v5/bn9SG8YmxN6u/6qqqv8LF8FvCEiKUB4T0wTk+au3UlKYgL3nN3b61BqTEQSgOM5WvO7WlV3eBWPXYREn5e/2sxrC7fw85FdoyrZAm8TLn9HKmBDH1V9A3ijup2q6lQR2QaMS0lJCVtdY1mZ0rNNk4iatymS5OQX8O6KbRw8XMrTVw6J2UFh/SVbItIMaK+qK1S1CNgQ/shMrLn2lM6MOz6TVk0issL/GCLSFbgbGA2sB3bh1vyKyCHgKeAFr263m+iwYGMu989azaierbj9hz29DqfWvEy4tgK+Ex21A3KCsWMvrloSEoRHLh4UrreLOoM7NOPXZ/Xm9++s4fjfvc/+whIyM9K4c0xPT7ryzlyazeQ568jJKwhJHCIyDzgXp4wtA3aJyMeqelvQ3sTEpYKiUjbuPkifzKZRkWy5HgCeAH6qFXrQiEhr4FLgCuAFD2IzUWBbfgE/f3kx7Zs3ZMrFg6Kudgu8TbgWAt1FpDOQjdPI+NJg7NiL+/IHD5fQqIGNI1uVZg2TSRDY5/bizM4rYNIbKygpLWPcwMxjtk0UISkxAVWlqLTyRW/5+rIypbis8vqkhAQSE8Tv+reX5XDvW6soKC7ziWMlQDCTrnRV3Sci1wPPq+p9IrIiWDs38evxeRt4fN63fHT7SDq0iI75ElX1kirW7QSmhC8aEy18L4xbNmmAAE9fOYT0tOhslhKWDEFEXgVGAi1FZCtwn6o+KyI3AXNweiY+p6qrg/F+XtRwXfHsV7Rs3ICpVw4N11tGnYff/4aKA/AXFJdxx4wV3DHj2Fzk/0b34JbR3dm+r5DhD31YaV/3nN2b63/Qhe92H2T03z6utP5PF/Tn4mEdWJmdz3mPfVZtbAXFpUyesy6YCVeSiLQFLgJ+E6ydmvi2cfdBnvr4O849PjNqki1fIvJj4D1V3S8i9wCDgQdUdYnHoZkIM3NpNpPeWElBsTNN3K79h0lLTmBV9j66tW7icXR1E65ein6vblR1NjA72O8X7hqugqJSVmzN5/of2ISxVcnJKwi47s4xx96PLx9ao3GDpErrAIa6c1W2aJTid33/dukAtElPrbR+8px1tY6vDn6PczHxqaouFJEuOG1XjKmT8t5ZDZISmHRW5A/yGMC9qjpdRE7B6QT1MM6txoAdpkx8mjxn3ZFkq1xBcVmwL4zDKibvgYW7hmvplr2UlCnDOluD+apkZqSR7SepycpICzi7e3UzvzdrlFLl+uOaplZa/8pXm/3GkZmRFnA/taWq04HpPs+/Ay4M2huYuPPequ18sn43943rQ+smqV6HU1flv6BnA0+o6lsicr+H8ZgIFegCOMgXxmEVk+MXiMg4EZman58flvdbuHEvIjCkow14WpU7x/QkrcLgp2nJiX5rqKI1DhG5R0QCnggicrqIVDXOnDF+7dhXyMD2GVxxUkevQ6mPbBF5CudW+2wRaUCM/g6Z+mnZ2H+HkGBeGIeb1XAFwcJNufRq0zRqG/KFS3k1cCh7B0ZAHCuBt0WkEFjC0e7v3YGBwAfAH4PxRia+XD2iM1cM7xSVvbN8XIQzAPDDqprntnO80+OYTITZtPsgh4pKEI4dK8qLC/RgismEK9xtuC4/qQPFpTZXYE2MH5QVEfffQxWHqr4FvCUi3YERQFtgH/ASMFFVo7c+3Hji210H2LjrIKP7HBe1yZaINFbVA6p6CJ/xFFV1G7DNdxuvYjSR45sd+2nUIIlbRnfnhc+/9/QCPZhiMuEKdw3X2H5tw/E2Joqo6nqskbypp/KG8iu35vPJ3adHcy36WyKyDHgLWKyqBwHcziSjcGq+ngZmeBah8ZyqIiL8sG8bftC9FWkpiUw8tavXYQWN3Tuvp1XZ+azJ2ed1GCZOiUgXEXlWROyHKga9s2Ibn23Yw51je0VzsoWqngHMBX4KrBaRfBHZg1Pz2wa4SlXtHI5jJaVl/OKVJby1LBuAtJTEal4RfWIy4Qpno/lH5q7n5y8vDvn7mPghIs+JyE4RWVVh+VgRWSciG0TkV+D0flTV67yJ1ITSgcMlPPDuGvpnpXPpsOifflNVZ6vqZaraSVXTVbWFqp6sqg+q6nav4zPeUVV+O2s1s1duJ7+g2OtwQiYmEy5VfVtVJ6anp4f0fcrKlEWbcjmhk/VONEE1Dadh8REikgg8BvwI6ANcIiJ9wh+aCZdHPviGnfsP84fx/aK27VakC3ePduPf4/O+5ZWvNnPjaV25cngnr8MJmZhMuMJlw64D7D1UfGSQTmPKiUgPEZlbXkslIgPckbWrparzgdwKi4cBG9warSLgNeC8WsQzUUQWiciiXbt21fRlxkO92zblxtO6MrB9htehxKxwXZybwN5YstUZzHRgJndFcQ/EmrCEqx4WbHR+Ey3hMn48DUwCigFUdQXOfKF1lQVs8Xm+FcgSkRYi8iQwSEQmBXqxqk5V1aGqOrRVq1b1CMOEywWD23H32KgdUd6YGtm4+yAnd23BXyYcT0KM1+TGZC/FcA0LsXBTLq2bNKBD8+ib08yEXENVXSByzBdIST325++bSFV1D3BjjXbgwaTupvbeXbGNnfsLuTL6x9zyy53Wp7uqPi8irYDGqrrR67hMeJX3SLz9hz0pKikjJSn2639i8hOGq5r4D+P78fw1J1DhR9UYgN0i0hV33D4RmYA73lAdbQXa+zxvB+TUZgd2+yTy7Sss5r5Zq5m5LMdvhh3tROQ+4G6c2l+AZJyeiiaOZOcVMO7RT1mV7bSdi4dkC2K0hitcmqYm0zfTfryMX78ApgK9RCQb2AhcXo/9LQS6i0hnIBvn9uSl9Y7SRISZS7OZPGfdkTk+rxzeIVZvr5wPDMKZhQFVzRGRJt6GZMIp/1AxVz+3gO37CklOjI9Eq1x8fdog+mzDbv4xdz0FRaXVb2zijtu4fTTQCuilqqeo6qaavFZEXgW+AHqKyFYRuU5VS4CbgDnAWuB1VV1dm5isR1Zkmrk0m0lvrDxmQvUn5n3HzKXZHkYVMkWqqhyt+W3kcTwmjAqLS7nhxUV8v+cQU68YSs828ZVrWw1XHb2zIod3V2zjplHWHsZUJiIZwJVAJyCp/Lazqt5c3WtV9ZIAy2cDs+saU7hnYDA1M3nOOgqKj71wKygudXpuRfE0JgG87k5enSEiNwDX4nQwMTGurEy5ffpyFmzM5R+XDGJ41xZehxR2lnDV0VcbnfG3YrTa39TfbOBLnMmsyzyOxUSwnDz/02sGWh7NVPVhETkTZ37RnsBvVfV/HodlwqCotIyDh0v49Vm9OPf4TK/D8URMJlyh7o21+8Bhvtt1kIuGtq9+YxOvUlX1Nq+D8GW9FCNT2/RUcvILKy3PzEjzIJrQU9X/ichXuL8/ItJcVSuOO2diSGmZkpqcyLNXnUA811HEZBuuUPfGWrTJ+W6wEeZNFV4UkRtEpK2INC9/eBmQ9VKMPAcPl1BcVkZShV+htORE7ozBQSBF5KcisgNYASwCFrv/mhj17optXPD4Z+w5cJjEBInrXv0xmXCF2o59h2nWMJn+WfbDZQIqAibjNH5fjP2wGD8eeHctuw8UceNpXcjKSEOArIw0Hrqgfyy23wK4A+jrzqfYRVU7q2oXr4MyobFgYy7/9/oykhMTaNQgJm+o1YodgTq46uROXH5Sx5gclNAEzW1AN1Xd7XUgJjJ9sGYHry7YzE9P68IdY3pxx5i4GFX+W+CQ10GY0Fu/Yz/Xv7CQ9s3SeOaqoaQmJ3odkucs4aojS7ZMNVYTYT8s1oYrcuzaf5i7/7OC3m2bctuZPbwOJ5wmAZ+7bbgOly+sSe9dEz127Cvk6ucX0iA5kWnXDCOjYYrXIUUES7hq6ZP1u5g8Zx3/uHgQnVraEDImoFJgmYh8RIT8sNiwEJFj2ZY8ikrLeOTigTRIiqsr/6eAD7HeuzGtqKSMFo1T+OP5/WlvU98dETUJl4gkAH8AmgKLVPUFL+L48rs9rMnZR+umDbx4exM9ZroPYyo5s89xfPar02mamux1KOFWEmm9d03wFJeWkShC++YNeesXI+K6gbw/tUq4RKQ90BfoB/THafw4tAavew44B9ipqv18lo8FHgESgWdU9U9V7OY8IAvIxZlXzhMLN+6lb1Y6DVOiJlc1HvDqgsBEto27D7ImZx9nD2gbj8kWwEciMhF4m2Nrfm1YiCinqtz9nxWUlil/v2igjVHpR7VZg4j8FLgK6AM0AN4FVgGzgAdr+D7TgEeBf/nsNxF4DDgTJ4FaKCKzcJKvhyq8/lqcQfK+UNWnRGQGMLeG7x00hcWlLNuSx1Undwz3W5soISKvq+pFIrISd/oSX6o6wIOwAGvD5bXi0jJu/fcyNu0+yCndWpLeMC4TrvL5Pyf5LFPAeipGub++/w1vLMnm/0b3sGQrgJpU00wCfgLsBv4EpAHPqermmr6Jqs4XkU4VFg8DNqjqdwAi8hpwnqo+hFMbdgwR2YrT1R6c9jF+uVdPEwE6dOhQ0xBrZMXWfIpKyxjWOf6mJDA19nf330rnsNesDZe3/vnhBpZvyePxywbHa7KFqnb2OgYTfC9/9T2PfrSBi09oz81n2AVdIDVJuM5R1VXu3z92bwO+LSLTgEdUta4NH7OALT7PtwInVrH9G8A/ReQHwPxAG6nqVGAqwNChQyvVMNRHg6QExvQ9jqEdmwVztya2PAYMVtXvvQ7ERI7F3+/l0Q/Xc8HgLM7q39brcMJORE5X1Q9F5AJ/61X1jXDHZILjgzU7uHfmKkb1bMUD4/tZu60qVJtw+SRb5c/fc3te3QN8Bgyv43v7+18JmCCp6iHguhrtOES3To5vn8FTV1TbZM3EN/u2MccoLC7ltteXkZmRxu/O7et1OF45Fad34jg/6xTngjpoRKQL8BsgXVUnBHPf5liNGiRxcteWPHrpYJISbSz1qtSp5beqHgbuFZEX6/HeWwHfyQjbATn12F9IlZYpew4cpnXTVK9DMZEtS0T+EWiljTcUf1KTE7ntzB5kZqTRJD4bygOkAKjqNXXdQW06X7lNVa5z2/uaECgoKiUtJZHhXVtwUpfmVrNVA/VKR1X1m3q8fCHQXUQ6i0gKcDFOQ/x6C8WccWu37WPYH+cyZ/X2oO3TxKQCjk7l4+/hGREZJyJT8/PzvQwjrhwucZqbnjcwK97nXh0bhH1Mq7gfn85XP8Lp2HWJiPQJwnuZKuw+cJgfPTKfFz7fBGDJVg2FZWwDEXkVGAm0dBu/36eqz4rITcAcnCuT51R1dZDeL+i3FBdsdHotD2hn8yeaKu2J1CEhrNF8eO3cX8j4Rz/j12f35pwBmV6H47VEEWlGgFvuNRkWojadr4A11e0vlB2sYtmhohKum7aQbfmF9Lffw1oJS8KlqpcEWD4bmB2C9wv6D8uCjbm0a5ZG2/S0YO3SxKai6jcxsU5VuWvGCvYcLKJXmyZehxMJeuHU8AZqu1vXYSH8dr4SkRY4wxYNEpFJbu/3Y980hB2sYlVJaRm/fGUpK7PzefLyIQzuYB3IaiMmR+8Mdg2XqrJwUy6n9WgVlP2Z2KWqJ3kdg/HeS19+z7x1u/jduX3p1toSLmCNqg4KwX79JnCquge4MQTvF7dUlXvfWs3cr3fyh/H9+GHfNl6HFHViMuEKdg3Xd7sPsudgESd0jus2GMaYGtiw8wAPzl7LaT1aceVwGyQ5xKKq81U0mrk0m8lz1pGTV0DTtGTO7N2aK06y87ourA9nDbRs1ICHf3w8I3taDZcxpmqfbdhNo5QkJk8YYI2Jj3okRPutV+cr60hStZlLs5n0xkqy8wpQIL+gmE837GHm0myvQ4tKMZlwBbsQpTdMZsKQdtZ+y9SKiCSKSKaIdCh/eB2TCb2rTu7Eh7ePtCFkfKjqtPruw+189QXQU0S2ish1qloClHe+Wgu8XpvOV6Ho0R5LJs9ZR0HxsRO7FBSXMnnOOo8iim52S7EG/rN4K8M6N6d984bB2J2JAyLyS+A+YAdQPhuDAjaXYoxasnkvpWXKCZ2ax+3UPaEU7s5XBrLzCvwuzwmw3FQtJmu4giknr4Dbpy/n/TU7vA7FRJdbgJ6q2ldV+7sPz5ItsKv5UNpfWMwtry3lrhkrKCmt62xnsU9ERtRkmfHeM598F3BdZobd7amLmEy4gnlLceEmZ3iYE63BvKmdLYA1DIkTv3t7Ddl7C3j4xwNsepOq/bOGy8LC2nD598GaHTzw7loGtGtKavKx53NaciJ3junpUWTRzW4pVmPBxlwaN0iid9umQYjMxJHvgHki8i5wuHyhqv7Nu5BMKPx35TZmLN7KL0/vxpCOdmHmj4gMB04GWonIbT6rmuIMfO0JGwzYv1G9WvPQBf358ZB2vLNi25FeipkZadw5pifjB2V5HWJUismEK5gWbsplcMdmJCZYbyNTK5vdR4r7MDFo1/7DTHpzJQPapXPzGd29DieSpQCNcX5zfAcm2wfY5NIRoKS0jMlz1nH1iE60TU/jkmFOH5/xg7IswQoSS7iqkF9QzPqdBzhvoJ1spnZU9Xdex2BCr0WjFH4xshun925Nst1KDEhVPwY+FpFpqvo9gIgkAI1VdZ+30ZnC4lJufnUp76/ZQbvmDW2crRCJyYQrWL2x0tOSWXLPmcEJysQFEZmiqreKyNs4vRKPoarnehCWCYGikjJSkhK44dS6zkoTlx4SkRuBUpypftJF5G+qOtmLYKznLhw4XMLEfy3i82/3cP+4PpZshVBMXpIFszdWs0YpNGtkd4RMjb3o/vsw8Fc/DxMD1u/Yz2mTPzoyqb2psT5ujdZ4nKEcOgBXeBVMvPfczTtUxGXPfMVXG3P564+P5+oRnb0OKabFZA1XsPz+7TUc3z7dbimaGlPVxe6/H3sdiwmNopIybnltGYdLyujU0sbmq6VkEUnGSbgeVdViEbGJoz0iIgjwxGWDbW7EMLCEK4CDh0t44YtN3JhitwtM7YlId+AhoA9wZMhxVfXshLLbJ8Hx9w++Yc22fUy9Ygitm9ho8rX0FLAJWA7MF5GOOA3nTRhl5xXQolEK6WnJvPGzk0mwTmFhEZO3FINh6ea8I6NGG1MHzwNPACXAKOBfHL3d6Il4v30SDAs25vLkx99y8QntrUagDlT1H6qapapnqeN7nPJhwuSbHfs5/7HP+PUbKwEs2QqjmEy4gjGY3YKNe0gQGNKxWRAjM3EkTVXnAqKq36vq/cDpHsdk6mnO6u10aN6Qe8/p43UoUUlEjhORZ0Xkv+7zPsBVHsYTVwOfLt+Sx0VPfQHAT0/r6nE08ScmE65gXMkv2JRLn8ymNEm1OdFMnRS63d7Xi8hNInI+0NrroEz93HN2b978+QgaNbDWGHU0DWei6Uz3+TfArV4FE0+1vp9/u5tLn/6SJqlJTL9xOD3bNKn+RSaoYjLhqi9VJTkxgRHdWnodioletwINgZuBIcDleHglb+rnk/W72LDzACJCc+u1XGsiUp6htlTV13EndFfVEpwhIkwIFRaXctu/l5OZkcaMG0+mY4tGXocUl+wyzQ8R4cXrTvQ6DBOlRCQRuEhV7wQOANd4HJKph+35hdz0ylL6tG3KqxNP8jqcaLUAGAwcFJEWuGPUichJ2JyjIZeanMhzV59A2/RUG+bIQ5Zw+aGqiFhDQlN7IpKkqiUiMkRERFWty3sUKytT7pi+nKKSMv54QX+vw4lm5V+otwGzgK4i8hnQCpvaJ2Re/GITuw8U8X9n9qBPps0H7DVLuPy48aXFpCUnMuXiQV6HYqJP+ZX8UuAtEZkOHCxfqapveBWYqb3nP9/Epxt288fz+9O5pd2GqQffSavfxBn0VHAmdh8NrPAqsFikqjz20QYefv8bRvc+jtIytfmAI4AlXBWUlSlffLuHs/q39ToUE92aA3tweiYqzo+LApZwRYn1O/bz5/e+ZnTv1lwyrL3X4US7RJzJqyv+6ns6cmwsjk2nqjz036+ZOv87zh+UxV8mDLBkK0JETcIlIj8ALsOJuY+qnhyK91m3Yz/7CksY1tnG3zJ10tq9kl/F0USrnN1ejCLtmzfk+lM6c+0pna2JQf1tU9Xfex1ERar6NvD20KFDb/A6lmC5Z+YqXv5qM1cO78j94/raOFsRJCwJl4g8B5wD7FTVfj7LxwKP4Fz9PKOqfwq0D1X9BPhERMYDC0MVa/ncaDbgqamjQFfyYAlX1CgpLSM1OZG7xvbyOpRYYb/6YTKkYzOaN0rhtjN72IVChAlXDdc04FGc0baBIz25HgPOBLYCC0VkFs4P1kMVXn+tqu50/74UuD5UgS7YlEvb9FTaNUsL1VuY2BaRV/KmejOXZjN5zjqy8wpITBDuGtuTn55qg0MGyRleBxDLCopKWZmdz7DOzblgcDuvwzEBhCXhUtX5ItKpwuJhwAZV/Q5ARF4DzlPVh3BqwyoRkQ5AvjvbvF8iMhGYCNChQ4dax3pKt5YMbJdhVwamruzEiUIzl2Yz6Y2VFBQ7Q0KVlil//983HNcklfGDbPL6+lLVXK9jiFX5BcVcN20hq3P2Mf+uUbRq0sDrkEwAXrbhygK2+DzfClQ3+NV1OHPUBaSqU0VkGzAuJSVlSG2DumRY7ZM0Y3yE9UpeRBoBjwNFwDxVfTmc7x8rJs9ZdyTZKldYXMbkOess4TIRa9f+w1z13ALW79zPIxcPsmQrwnk50nyt27io6n2q+nl1O67rdA3ZeQXsOXC4Vq8xxlcwruRF5DkR2SkiqyosHysi60Rkg4j8yl18ATBDVW8Azq3ve8ernLyCWi03xmvZeQVc9NQXbNx9kGeuOsF61kcBLxOurYBvX+t2QE4wdlzXCUkf+eAbRv/tY2ysSuOxacBY3wU+bR5/BPQBLnEn/m3H0ZpimyKljtpmpPpdnplhbTlNZJq+aAu7DxzmpeuHcVqPVl6HY2rAy4RrIdBdRDqLSApwMc4IxPVW1xquBRtzGdqpubXfMp5S1flAxZqyI20eVbUIeA04D+fCpbyVbMDyLCITRWSRiCzatWtXKMKOOnsOHOb215ezc38hd43pRVpy4jHr05ITuXNMT4+iM+FQ14tzL5WWORUCN5/endk3/4AhHa1HfbQIS8IlIq8CXwA9RWSriFznTlp6E87M8WuB11V1dZDer9aFaOf+QjbtOcQwGw7CRCZ/bR6zcAZSvVBEngDeDvRiVZ2qqkNVdWirVnY1/OHXOxgz5RPeXp7Dss15jB+UxUMX9CcrIw0BsjLSeOiC/tZ+K8bV9eLcKws25jJ2yny25B4iIUFo39zTcWNNLYWrl+IlAZbPxpniIdjvV+vB7BZu3AvACTbgqYlMfts8qupBajg5diyOql1bh4pKePDdtbz81WZ6tWnCS9cPo1cbZ4658YOyLMEyEeujdTv52UuLycxIs5Hjo5SXtxRDpi41XAs27qFhSiJ9bYJPE5nq3eYx2q7mQ+Ev763jlQWbmXhqF966acSRZMuYSPb28hxueGER3Vo3ZvpPh1vbwigVNVP71EZdariuPaUzp/ZoRXJiTOagJvodafMIZOO0ebzU25CiQ0lpGXkFxbRs3ICbz+jOmL5tGN61hddhGRNQ+SC8OXkFNGuUQu7BIoZ1as4zVw+laWqy1+GZOorJ7KIuNVwdWzTijN7HhTAqY2omVG0eo7GBcH1t2n2QHz/1Bde9sIjSMqV5oxRLtkxEKx+ENzuvAAVyDxaRlCBMGNLOkq0oF5MJV21vnazOyef1hVsoKLJe9cZ7qnqJqrZV1WRVbaeqz7rLZ6tqD1XtqqoP1mG/cXNLUVV5dcFmzvrHJ3y78wDXn9LZ2r2YqOBvEN6SMuWRues9isgES0zeUqytt5dv49lPv+PcgZleh2KMqae8Q0XcMX0FH6zdwcldW/DXi46nbbq1eTHRIdsG4Y1ZMZlw1bY31sJNufTPSie1wjg8xsSSeOml2CApkZy8Au45uzfXjuhMgtVsmSiQf6iYP7y7JuB6aygf/eL+lmJhcSkrtuYxrLO16zCxLZZvKR4qKuGv76/jUFEJaSmJzLppBNf/oIslWyZq/HrmSt5cms2ZvVuTmnzsT7MNwhsbYrKGqzaWbs6juFQZ1rmZ16EYY+pg6ea9/N+/l/F97iH6ZqYztl8bkqy3sakBr2t98w4VUVKmtGzcgLvG9OSnp3ZhQLuMY3opZmakceeYnjZGXAyI+4Tr6+37EMGmRzAmypSUlvHPDzfw6EcbaNM0lVdvOImTulhNtam5ugwhFCxzVm/nN2+uYmjHZjx5xRA6tmh0ZJ0NwhubYjLhqs1VyzUjOnPewCzS06y7rYltXl/NB9vv3l7Di19+z/mDsvjdeX2ty7yJCrkHi7h/1mpmLc+hd9um3HR6bJRHU72YTLhqe9XSvFFKiCMyxnteXs0Hi6pSWFxGWkoiN/ygC8M6N2fc8da72ESHJZv3MvFfi8gvKOb/Rvfg56O62mDbcSQmE66aWp2TzyMfrOdXP+pFl1aNvQ7HGFOFXfsP86v/rEBEePrKIXRo0ZAOLWzyXhM9OrdoRP+sdO7+US+bVioOxXVq/fmGPby/ZgeNG8R13mlMxPvfmh2MnTKfTzbs5mQbKd5ECVXl7eU5XP38AkpKy2jWKIXnrxlmyVaciutMY8GmXDq2aEjrpqleh2JMyEVjG66Dh0t44N01vLpgC33aNuXViwfS47gmXodlTLV27i/k3pmrmLN6B8e3z2DvoWJaNWngdVjGQ3GbcJWVKYs25dr8iSZuRGMbrsLiUuau3cmNp3Xl/87sToMkG5zYRDZV5a1lOdz/9moOFZUy6Ue9uO6UzjZUiYnNhKsmV/Lf7jrA3kPFDOtsw0EYE0mKS8uYvmgrFw1tR4vGDZh7+2k0sR6IJkoUlyqPfrSBLi0b8ZcJx9OttbUPNo6YTLhqciW/r7CE49tnMKyTJVzGeKXiAI9Xn9yJd1ZuY/mWPJo3SmFsvzaWbJmIp6rMWp7DGb2Po3GDJF667kRaNWlgE6abY8RkwlUTQzo2461fjPA6DGPi1syl2Ux6YyUFxaWAM2nvg7PXkpacwGOXDmZsvzYeR2hM9bblFzDpjZXMW7eL35zVmxtO7UKbdGsXbCqL24SrqKSMlCS7p27iR6Q1mv/ze18fSbZ8paelcPaAth5EZEzNqSqvL9rCA++spaRMuX9cH64c3snrsEwEi8uEa+veQ5zx14/5+08GclZ/+2I38cHLRvO79h9m/je7+GbnftbvOMA3O/azLb/Q77Y79vlfbkwkmTxnHY/P+5aTujTnzxcOOGZqHmP8icuEa8HGXA6XlNHJCogxQVNYXMqGnQdYv3M/3+w4wPod+7n65M6c0r0l63fu5/bpy0lJTKBLq0YM6tCM/IJi9heWVNpPZkaaB9GbeFTbWl9V5VBRKY0aJPGTE9rTNj2Vy07sSIK11TI1EDUJl4h0AB4FdgPfqOqfaruP8ga62XkFCPD1tn30ybQB6Ex8q9hw/c4xPaucOLewuJTvdh1k/c79tG/ekMEdmrEl9xCnTv4IVWebpAShS6tG7C8sBmBQ+2Z8cNtpdGrR8Ej3+IptuADSkhO5c0zP0H1YY3xUVetbsVxce0on5q7dSVpyIs9cNZSOLRpxxXC7aDc1F5aES0SeA84BdqpqP5/lY4FHgETgmWqSqB7Au6r6lIj8q7YxVPxyV+A3M1eRkCA2K7uJW/4ark96YyUAP+rfhv2FJbRs3IDSMuWmV5awbvt+Nu05SJmbWF12YgcGd2hGZkYat5zRne6tm9DjuMZ0atnomDni0lISK3WPLy93tUn2jAkHf+XiD++spUGicN+5/ap5tTH+hauGaxpO7dSRRElEEoHHgDOBrcBCEZmFk3w9VOH11wJLgd+IyE+AF2sbwOQ56yo10C0oLmXynHX2BW/iVqByccf05dw+fTkje7Ti2atPIDFB2F9YQo/jmnDOgLZ0P64JPY5rQueWzhV+YoJw6+getX7/8YOyrPyZiOOvXABkNErh0hM7eBCRiQVhSbhUdb6IdKqweBiwQVW/AxCR14DzVPUhnNqwY4jIHcB97r5mAM/7ey8RmQhMBOjQ4WjByMkr8BtboOXGxBp/7VUCnf8lZcpNo7oxqEPGkWUvXX9iiCM0JjIEKhc79x0OcyQmlng5LkIWsMXn+VZ3WSDvATeLyJPApkAbqepUVR2qqkNbtWp1ZHmghrjWQNfEC1V9W1UnpqenH1kW6PzPykjjjjE9beorE5fs98KEgpcJl79uHRpoY1VdpaoTVPVGVb2jyh2LjBORqfn5+UeW3TmmJ2nJx87DZg10TbyzcmFMZVYuTCh4mXBtBdr7PG8H5ARjx/6u5McPyuKhC/qTlZGG4FzBP3RBf2s/YuKalQtjKrNyYULBy2EhFgLdRaQzkA1cDFwajB0HGlvFGugaU5mVC2Mqs3Jhgi0sNVwi8irwBdBTRLaKyHWqWgLcBMwB1gKvq+rqcMRjjDHGGBNO4eqleEmA5bOB2SF4P8+mMDHGGGOMqUhUA7ZTj3oikg+sd5+mA+Wt6FvijFgfLL77Dsb2gdb7W16TZb7PI+U41PUYBFpX1Weu+Lymx6CjqrYKsC5qRWi5qMm2wSwXVZ0fwTwO4fpuCLTOykUNiMguIA/7noy1chHO38zqy4WqxuwDmBrg70Whep9gbB9ovb/lNVkWicehrsegLp85nMcgGh7ReD5UtU1dykU150fQjkO4vhtqehysXNTs2EZLuQjl+RDK4xCvv5le9lIMh7cD/B3K9wnG9oHW+1tek2WReBzqegwCravqM1d8HspjEA2i8Xyoapu6lIvqzpdgCdd3Q6B1Vi5qLhrLRSjPh5q8f13F5W9mTN9SDEREFqnqUK/j8JodBzsGvuxYOOw42DHwZcfCYceh/scg1mu4ApnqdQARwo6DHQNfdiwcdhzsGPiyY+Gw41DPYxCXNVzGGGOMMeEUrzVcxhhjjDFhYwmXMcYYY0yIWcJljDHGGBNilnAZY4wxxoSYJVyAiHQRkWdFZIbXsXhJRMaLyNMi8paI/NDreLwgIr1F5EkRmSEiP/M6Hi9ZubAyUc7KxVFWLqxclKt1uQjWyLGR9gCeA3YCqyosHwusAzYAv6qwbobXcUfIcWgGPOt17B4fg4RYOgb1PBYxVS6sTNTrOFi5OLrOyoWVi/J1NSoXnn+4EB60U4HBvgcNSAS+BboAKcByoI/P+pgqQPU4Dn8FBnsdu1fHADgX+By41OvYvT4W7vqYKhdWJup2HKxcWLmwclG/chGztxRVdT6QW2HxMGCDqn6nqkXAa8B5YQ8ujGpzHMTxZ+C/qrok3LGGSm3PBVWdpaonA5eFN9LQs3JhZaKclYujrFxYuSgXynIRswlXAFnAFp/nW4EsEWkhIk8Cg0RkkjehhZXf4wD8EhgNTBCRG70ILIwCnQsjReQfIvIUMNub0MLOyoWViXJWLo6ycmHlolxQykVSqKKLUOJnmarqHiAeTppygY7DP4B/hDsYjwQ6BvOAeeENxXNWLqxMlLNycZSVCysX5YJSLuKthmsr0N7neTsgx6NYvGTHwY6BLzsWdgzK2XE4yo6FHYNyQTkO8ZZwLQS6i0hnEUkBLgZmeRyTF+w42DHwZcfCjkE5Ow5H2bGwY1AuKMchZhMuEXkV+ALoKSJbReQ6VS0BbgLmAGuB11V1tZdxhpodBzsGvuxY2DEoZ8fhKDsWdgzKhfI4iNut0RhjjDHGhEjM1nAZY4wxxkQKS7iMMcYYY0LMEi5jjDHGmBCzhMsYY4wxJsQs4TLGGGOMCTFLuIwxxhhjQswSrgjmztm1zH1sF5Fsn+cp1bx2qIjUauoFEblWRFaKyAoRWSUi57nLfy8io+vzWYwJFisXxlRm5SLy2ThcUUJE7gcOqOrDPsuS3AHZgrH/dsDHwGBVzReRxkArVd0YjP0bEwpWLoypzMpFZLIarigjItNE5G8i8hHwZxEZJiKfi8hS99+e7nYjReQd9+/7ReQ5EZknIt+JyM1+dt0a2A8cAFDVA+WFx33PCe5VUPkV00oRUXd9VxF5T0QWi8gnItIrHMfCmHJWLoypzMpFZEnyOgBTJz2A0apaKiJNgVNVtcStxv0jcKGf1/QCRgFNgHUi8oSqFvusXw7sADaKyFzgDVV923cHqroIGAggIpOB99xVU4EbVXW9iJwIPA6cHqTPakxNWbkwpjIrFxHCEq7oNF1VS92/04EXRKQ7oEBygNe8q6qHgcMishM4DmcGdADcwjgWOAE4A/i7iAxR1fsr7khELgIGAz90q5JPBqaLSPkmDer7AY2pAysXxlRm5SJCWMIVnQ76/P0H4CNVPV9EOgHzArzmsM/fpfj5v1enQd8CYIGI/A94HrjfdxsR6Qv8DucqqVREEoA8VR1Yp09iTPBYuTCmMisXEcLacEW/dCDb/fvquu5ERDJFZLDPooHA9xW2SQdeA65U1V0AqroPp1r5x+42IiLH1zUOY4LEyoUxlVm58JAlXNHvL8BDIvIZkFiP/SQDD4vI1yKyDPgJcEuFbcYDHYGnyxtDussvA64TkeXAauC8esRhTDBYuTCmMisXHrJhIYwxxhhjQsxquIwxxhhjQswSLmOMMcaYELOEyxhjjDEmxCzhMsYYY4wJMUu4jDHGGGNCzBIuY4wxxpgQs4TLGGOMMSbE/h+yWH2Ixaw5TgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 612x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(8.5, 3.))\n",
    "\n",
    "ax = axs[0]\n",
    "ax.loglog(results['train_size'], results['alpha'], '--o')\n",
    "ax.set_ylabel('$\\lambda$')\n",
    "ax.set_title('Hyperparameters')\n",
    "\n",
    "ax = axs[1]\n",
    "ax.loglog(results['train_size'], results['train_time'], '--o')\n",
    "p = linregress(np.log(results['train_size'].iloc[-3:]), np.log(results['train_time'].iloc[-3:]))\n",
    "ax.text(0.1, 0.8, f'$t \\propto N^{{{p[0]:.1f}}}$', transform=ax.transAxes, fontsize=14)\n",
    "ax.set_ylabel('Train Time (s)')\n",
    "ax.set_title('Training Performance')\n",
    "\n",
    "ax = axs[2]\n",
    "ax.loglog(results['train_size'], results['test_time'], '--o')\n",
    "p = linregress(np.log(results['train_size'].iloc[-3:]), np.log(results['test_time'].iloc[-3:]))\n",
    "ax.text(0.1, 0.8, f'$t \\propto N^{{{p[0]:.1f}}}$', transform=ax.transAxes, fontsize=14)\n",
    "ax.set_ylabel('Test Time (s)')\n",
    "ax.set_title('Inference Performance')\n",
    "\n",
    "for ax in axs:\n",
    "    ax.set_xlabel('Train Size')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('figures/fchl-scaling.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can note a few things:\n",
    "\n",
    "1. The hyperparameters do change as a function of training set size. Be careful when fitting this model to always re-tune $\\lambda$\n",
    "2. The training time scales with $O(N^2)$, *i.e.* with the number of training entries. The quadratic scaling is a result of the time it takes to compute the kernel matrix\n",
    "3. The evaluation/test time scales linearly with the training set size, which is driven by the cost of comparing a point to each training entry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "Going back to our learning objectives, let's review that FHCL is...\n",
    "\n",
    "- [x] kernel method (*uses features based on similiarity to training set*)\n",
    "- [x] created with atomistic, scalable kernels (*similarity functions that express atomic contributions to molecular properties*)\n",
    "- [x] that use \"alchemical\" many-body similarity functions (*functions that capture similarity between atoms' environments*)\n",
    "- [x] to regress atomic properties with ridge regression (*using similiarty scores as inputs to a penalized linear regression*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}