AL_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7b46cd07",
   "metadata": {},
   "source": [
    "# Active learning of potential energy surfaces \n",
    "Active learning (AL) algorithms allows one to build machine learning algorithms with a reduced amount of data. In this tutorial we will learn some AL algorithms and how to use them to build potential-energy surfaces of molecules. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efb60503",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install a pip package in the current Jupyter kernel\n",
    "import sys\n",
    "!{sys.executable} -m pip install numpy\n",
    "!{sys.executable} -m pip install tensorflow\n",
    "!{sys.executable} -m pip install Cython\n",
    "!{sys.executable} -m pip install scipy\n",
    "!{sys.executable} -m pip install h5py\n",
    "!pip install git+https://github.com/ntucllab/libact.git\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48895f51",
   "metadata": {},
   "source": [
    "## Reading data \n",
    "We will consider the example of building the PES of a pyrrole water dimer. For this we have apriori computed data from electronic strcutre theory. Let's start by reading and visualizing the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "153612be",
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py as hp \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "from sklearn.ensemble import RandomForestRegressor as RFR\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from libact.base.dataset import Dataset\n",
    "from libact.base.interfaces import QueryStrategy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b55fbc3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'dataset.h5'\n",
    "items = {}\n",
    "def print_attrs(name, obj):\n",
    "    for index, elem in enumerate(name):\n",
    "        if elem == '/' in name:\n",
    "            key = name[0:index]\n",
    "            x = f.get(name)\n",
    "            x = np.asarray(x)\n",
    "            items[key] = x\n",
    "f = hp.File(filename, 'r')\n",
    "f.visititems(print_attrs)\n",
    "X = (items['ExpDistance'])\n",
    "Y = (items['Energies'])\n",
    "Y += np.abs(np.amin(Y))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00a4d9e3",
   "metadata": {},
   "source": [
    "We will now plot the distribution of energies "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d5dc3f03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAESCAYAAADXMlMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAd20lEQVR4nO3de7gdVZnn8e+vw02BlkCOyBDCTRqJCAGOQZqrjYZAC7SP9GPSXFuYjAS0sZtWMozgoCPgrcVHAgQ6AtNNALGBjBMIICAtckswQAIGQoySCCQQkOsACe/8Uesklc257J3ae1ftc36f56nnVK26vfukst+zqmqtpYjAzMxsff1Z2QGYmVlncyIxM7NCnEjMzKwQJxIzMyvEicTMzApxIjEzs0I2KDuAdhgxYkTssMMOZYdhZtZR5s6d+0JEdA203ZBIJDvssANz5swpOwwzs44i6ff1bOdbW2ZmVogTiZmZFeJEYmZmhTiRmJlZIU4kZmZWiBOJmZkV4kRiZmaFOJGYmVkhTiTWcpMnT2by5Mllh2FmLeJEYmZmhTiRmJlZIaX0tSVpOvAZYHlE7N7L+n8Gjk2LGwC7AV0RsVLSEuBVYDWwKiK62xO1mZn1pqwayZXA+L5WRsR3I2JMRIwBpgC/jIiVuU0+mdY7iZiZlayURBIR9wArB9wwMxGY0cJwzMysgEo/I5H0frKay89yxQHcJmmupEnlRGZmZj2qPh7JkcC9Nbe1DoiIZZI+CNwu6bephrOOlGQmAYwaNao90ZqZDUGVrpEAE6i5rRURy9LP5cCNwNjedoyIaRHRHRHdXV0DDvBlZmbrqbKJRNIHgIOBm3Nlm0ravGceGAfMLydCMzOD8l7/nQEcAoyQtBQ4F9gQICIuTZt9FrgtIl7P7bo1cKMkyGK/JiJubVfcZmb2XqUkkoiYWMc2V5K9JpwvWwzs2ZqozMxsfVT21paZmXUGJxIzMyvEicTMzApxIjEzs0KcSMzMrBAnEjMzK8SJxMzMCnEiMTOzQpxIrG08brvZ4OREYmZmhTiRmJlZIU4kZmZWiBOJmZkV4kRiZmaFOJGYmVkhTiRmZlaIE4mZmRXiRGJmZoU4kZiZWSFOJGZmVkgpiUTSdEnLJc3vY/0hkv4kaV6azsmtGy9poaRFks5qX9RmZtabsmokVwLjB9jmPyNiTJrOA5A0DLgYOBwYDUyUNLqlkZqZWb9KSSQRcQ+wcj12HQssiojFEfE2cC1wdFODMzOzhlT5Gcl+kh6RdIukj6aybYFnctssTWVmZlaSDcoOoA8PA9tHxGuSjgBuAnZp5ACSJgGTAEaNGtX0AM3MLFPJGklEvBIRr6X5WcCGkkYAy4DtcpuOTGW9HWNaRHRHRHdXV1fLYzYzG6oqmUgkfUiS0vxYsjhfBB4CdpG0o6SNgAnAzPIiNTOzUm5tSZoBHAKMkLQUOBfYECAiLgWOAU6VtAp4E5gQEQGsknQ6MBsYBkyPiAUlfAQzM0tKSSQRMXGA9T8GftzHulnArFbEZWZmjavkrS0zM+scTiRmZlaIE4mZmRXiRGJmZoU4kZiZWSFOJGZmVogTiZmZFeJEYmZmhTiRmJlZIXW1bJc0EjgV+EvgQ6n4OeBe4LKIeKavfc3MbHAbsEYi6QDgCeBvgQXANWla0FMmaf9WBmlmZtVVT43kh8BPIuLLva2UdFHa5uPNC8vMzDpFPc9IPko2TnpfLgF2b044ZmbWaepJJM8C/d262j9tY2ZmQ1A9t7a+B1yaBpi6HXg+lW8NfBo4CTijFcGZmVn1DZhIImKqpBeBrwAnkw0oBbAamAucEBHXty5EMzOrsrpe/42I64DrJG0IjEjFL0TEOy2LzMzMOkJDIySmxOHnIWZmtkbRBom/Bi51g0Qzs6GraIPEY3CDRDOzIa2UBomSpgOfAZZHxHvaoEg6FvgaIOBV4NSIeCStW5LKVgOrIqK73vOamVnzldUg8UpgfD/rfwccHBEfA74JTKtZ/8mIGOMkYmZWvlIaJEbEPcDKftb/OiJeSov3AyMbOb6ZmbVPJzRIPBm4JbccwG2Sgqzn4draCgCSJgGTAEaNGtXC8MzMhrZKN0iU9Ml0zgNyxQdExDJJHwRul/TbVMOpjXsa6ZZYd3d3tCI+MzOrcINESXsAVwCHR8SLuViWpZ/LJd0IjAXek0jMzKw9GhohMSLeiYhn09TKJDIK+A/g+Ih4Mle+qaTNe+aBccD8VsVhZmYDK6VBoqQZwCHACElLgXOBDQEi4lLgHGArYKokWPua79bAjalsA+CaiLi1kXObmVlzDZhIUoPEW8jezLoNuDOt2pqsQeKXJB0eEffWe9KImDjA+lOAU3opXwzsWe95zMys9TxCopmZFeIREs3MrBCPkGhmZoV0QoNEMzOrsEo3SDQzs+qrbINEMzPrDB4h0czMCmmoZXtvJG2XxhcxM7MhqHAiAbYETmzCcczMrAPV07L9hAE2cR/tZmZDWD3PSK4E3iAbB6Q3zajVmJlZh6onCfyR7BXfzXub6L+xopmZDXL1JJK5wN79rA9AzQnHzMw6Tb0t2zfrZ/0i4JPNCcfMzDpNPS3b/3OA9a8Dv2xaRGZm1lH8oNxaavLkyWWHYGYt5kRiZmaFOJGYmVkhTiRmZlZIKYlE0nRJyyXN72O9JP1I0iJJj0raO7fuRElPpclds5iZlaysGsmVwPh+1h8O7JKmSWTD+SJpS+BcYF9gLHCupOEtjdTMzPpVSiKJiHuAlf1scjRwdWTuB7aQtA1wGHB7RKyMiJfIRmzsLyGZmVmLVfUZybbAM7nlpamsr3IzMytJVRNJYZImSZojac6KFSvKDsfMbNCqO5FI+lw/677WnHDWWAZsl1semcr6Kn+PiJgWEd0R0d3V1dXk8Gx9TZ482Y0UzQaZRmok/ybpCknv7ymQNFLSXcBXmhzXTOCE9PbWJ4A/RcSzwGxgnKTh6SH7uFRmZmYlaWTM9n2Ba4B5kv4O2JnsbaoHgD0bOamkGcAhwAhJS8nexNoQICIuBWYBR5B1CPkG8Pdp3UpJ3wQeSoc6LyL6e2hvZmYtVnciiYhHJXUDU4H7yLqPPzMiftToSSNi4gDrAzitj3XTAY8Rb2ZWEY0+bN8TOJispvA2MFbS5k2PygY9PycxGzwaedh+DnAPcDNZQtkH+AjwmKQDWxOemZlVXSPPSL4IHBkRt6XlhelB+LeAO4CNmx2cdS7XOMyGjkYSyR4R8UK+ICJWAWdJmtXcsMzMrFPUfWurNonUrLunOeGYmVmnGbQt283MrD2cSMzMrBAnEiuNu0sxGxycSMzMrBAnEjMzK8SJxMzMCnEiMTOzQtY7kUh6XzMDMTOzzlSkRvKr2gJJHylwPDMz60CNdJECgKQjgdHAZpK2i4j8GOrX0eDYJGawtm+uqVOnlhyJmTWq4UQCzCcb7nYEcLWkUWTD3T4LvNPE2MzMrAM0nEgi4nfAVEnze/rYkrQtsD1ZkjEzsyGkyDOS+ySdJOmfgTHA/RHxSnPCsk7QrFbpbt1u1tmKJJJrgW7gTeAzwMOSdm1KVGZm1jHW5xlJj50i4nM9C5LGAJcDBxUNyszMOkeRGsmrkj7csxAR84Dh9ewoabykhZIWSTqrl/X/Imlemp6U9HJu3ercupkF4jczsyYoUiM5Hbg5jY74OLAb8PuBdpI0DLgY+DSwFHhI0syIeLxnm4j4Sm77LwF75Q7xZkSMKRC3tcDkyZP96q7ZEFUkkTwF7A38DVkSeRr4Rh37jQUWRcRiAEnXAkeTJaPeTATOLRCntYkfmpsNTUUSya8iYh+yRojAmpbtvx1gv22BfCPGpcC+vW0oaXtgR+DOXPEmkuYAq4ALIuKmPvadBEwCGDVq1AAhmZnZ+mr4GYmkIyV9jdSyvWb1db3tU8AE4IaIWJ0r2z4iuoG/A34oaefedoyIaRHRHRHdXV1dTQ7LzMx6FG3ZflWqNTTSsn1Z2r/HyFTWmwnAafmCiFiWfi6WdDfZ85OnG4jfzMyaaMBEImlYvkbQhJbtDwG7SNqRLIFMIKtd1J73I2Rvgd2XKxsOvBERb0kaAewPfKeOc5qZWYvUUyN5TdKjwNw0PQw81pNEYE0toa9axToiYpWk04HZwDBgekQskHQeMCciel7pnQBcGxGR23034DJJ75Ldlrsg/7aXmZm1Xz2J5Atkb2ftQ1Zz+HPgLUnzySWXiJhb70kjYhYwq6bsnJrlb/Sy36+Bj9V7HjMza70BE0lEzABm9CxL2oUsqeyVfn4e2LyeY5mZ2eCzPr3/PiXpWbJbS38BbAwsb3ZgZmbWGep+/VfSn0s6XtLNwArgfLKW7OPI2oaYmdkQVM9bWycBxwCHAn8EfgZ8OyIeaG1oZmbWCeq5tTWd7I2sM4CfRMTbLY3IzMw6Sj23tu4CNgUuIevx92FJl0v6oqSPS9qotSGamVmV1fPW1qEAknYiG8hq7zR9FtgSeEfSgojYu5WBmplZNdX91lbqrXcxcH1PmaQdWJtczMxsCKrnYfsPgBuBeyPi3fy6iFgCLAFuaEVwVi093cR73BEzy6unRvI+svHZN5L0f4GbgNkR8WYrA7NqyY814nFHzCxvwIftEXFqRGwL/DXZ21vfAl6QNFPSFyS5j3YzsyGs7gaJEfFgRJwdEbsDewK/BE4Clkr6laQzUy/AZmY2hAyYSNIY6+uIiEUR8f2IOIhsPJHpwAFkw+KamdkQUqQb+VUAEbGCLJFMb1mUZmZWWaV0I2+Dhx+8m5m7kTczs0LcjbyZmRXibuStUnyrzKzzuBt5MzMrpJ4ayXSydiNnALtFxFeLJhFJ4yUtlLRI0lm9rD9J0gpJ89J0Sm7diZKeStOJReIwM7Pi6nlGchfZg/VLgB9JWsDat7XmAo80MkZJapdyMfBpYCnwkKSZEfF4zabXRcTpNftuCZxL1lFkAHPTvi/Ve34zM2uuerpIOTQitgQ+DBwP3AZsT9ZVygOkMUoaOOdYYFFELE4J6Frg6Dr3PQy4PSJWpuRxOzC+gXObmVmTldGN/LbAM7nlpcC+vWz3OUkHAU8CX4mIZ/rY1w/6zcxKVKjtRwu7kf8/wIyIeEvSfwOuAv6qkQNImgRMAhg1alSTwxs6/BaVmQ2k7td/m2gZsF1ueWQqWyMiXoyIt9LiFWQNH+vaN3eMaRHRHRHdXV3uoHh9OImYWT3KSCQPAbtI2jGN9z4BmJnfQNI2ucWjgCfS/GxgnKThkoaTtWGZ3YaYzcysD23v1iQiVkk6nSwBDAOmR8QCSecBcyJiJvBlSUcBq4CVZN3VExErJX2TLBkBnBcRK9v9GczMbK1S+seKiFnArJqyc3LzU4ApfezrnobNzCqkjFtbZmY2iDiRmJlZIU4kZmZWiBOJmZkV4kRiZmaFOJGYmVkhTiRmZlaIx1m39yi7a5Se80+dOrXUOMysPq6RGJB9eZedQMysMzmRWGU5uZl1BicSMzMrxInEzMwKcSIxM7NCnEjMzKwQJxJbhx9um1mjnEjMzKwQJxKrPNeSzKrNicTMzApxIhnC3ODPzJrBicTMzAopJZFIGi9poaRFks7qZf0/Snpc0qOSfiFp+9y61ZLmpWlmeyM3M7NabU8kkoYBFwOHA6OBiZJG12z2G6A7IvYAbgC+k1v3ZkSMSdNRbQl6kPBtLDNrhTJqJGOBRRGxOCLeBq4Fjs5vEBF3RcQbafF+YGSbYzQzszqVkUi2BZ7JLS9NZX05Gbglt7yJpDmS7pf0Ny2Iz8zMGlDpga0kHQd0AwfnirePiGWSdgLulPRYRDzdy76TgEkAo0aNaku8nWAwDBo1GD6D2WBSRo1kGbBdbnlkKluHpE8BZwNHRcRbPeURsSz9XAzcDezV20kiYlpEdEdEd1dXV/OiHyQ67XmJX1U2q64yEslDwC6SdpS0ETABWOftK0l7AZeRJZHlufLhkjZO8yOA/YHH2xb5IOUvaDMrou23tiJilaTTgdnAMGB6RCyQdB4wJyJmAt8FNgN+KgngD+kNrd2AyyS9S5YEL4gIJxIzsxKV8owkImYBs2rKzsnNf6qP/X4NfKy10ZmZWSPcst3MzApxIjEzs0KcSMzMrBAnkkFoqLyFNVQ+p1nVVbpBojWHv3DNrJVcIxkEhlJjvaHyOc06iWskg4i/ZM2sDK6RdLi+ksdQqqXUGqqf26wsTiRmZlaIE4kNSkO5RmbWbn5GYoOGE4dZOZxIOpS/NOuXH7+kt7FMPL6JWTFOJB1i8uTJ6/VFN9gTTiOfr3bbwf67MWsXJ5KK85ddMf79mbWeH7ZX0EAPiv3l2Br532vPv0G9v+veajuN7Fvv9v63typyjaSD+Eukedb3lliznqM0+lymv/ZCPdY3Nj8jsqKcSMzWUysSe38vA/S3jTXGv8PmciKpMNdA2q+VyaFH0VpIUfkXN8q6xvpKmI38bqqUBIZ6YnIiqYihfiEOBuvzpVzPLat69629hookjP6+qIu+Ql37LKqK13y7/j8Olv/3TiQtNNBFMlguIutbb1/gra4FrO/x+3s9ureEVE+top5rfKBk2sr/HwMlsr7WtzoBdtp3Q2mJRNJ44CJgGHBFRFxQs35j4GpgH+BF4PMRsSStmwKcDKwGvhwRs9sYelMU+UvUyjUY/o3Wp4bSqmPXe7xGbwmu7y3E2v0arW0V/fLvr6bZ18sVZTe0LSWRSBoGXAx8GlgKPCRpZkQ8ntvsZOCliPiwpAnAhcDnJY0GJgAfBf4LcIekv4iI1e39FPXr60KoXW/WbO24tuq9DdeqL9i8/mr/7X42VXtrsbZWV8/zofVpBlDG94kiov0nlfYDvhERh6XlKQARcX5um9lpm/skbQA8B3QBZ+W3zW/X1/m6u7tjzpw56xXrQH8B5Mtry5wgzAaPvv5PV/3/epEELmluRHQPuF1JieQYYHxEnJKWjwf2jYjTc9vMT9ssTctPA/sC3wDuj4h/S+X/CtwSETfUnGMSMCkt7gosLBDyCOCFAvu3W6fFC50Xs+NtrU6LFzov5nri3T4iugY60KB92B4R04BpzTiWpDn1ZOWq6LR4ofNidryt1WnxQufF3Mx4y+oiZRmwXW55ZCrrdZt0a+sDZA/d69nXzMzapKxE8hCwi6QdJW1E9vB8Zs02M4ET0/wxwJ2R3YebCUyQtLGkHYFdgAfbFLeZmdUo5dZWRKySdDowm+z13+kRsUDSecCciJgJ/CvwvyUtAlaSJRvSdtcDjwOrgNPa8MZWU26RtVGnxQudF7Pjba1Oixc6L+amxVvKw3YzMxs83I28mZkV4kRiZmaFOJH0Q9J4SQslLZJ0VsmxTJe0PLWv6SnbUtLtkp5KP4enckn6UYr7UUl75/Y5MW3/lKQTeztXk+LdTtJdkh6XtEDSP1Q5ZkmbSHpQ0iMp3v+ZyneU9ECK67r0cgjpZY/rUvkDknbIHWtKKl8o6bBWxJs71zBJv5H08w6Jd4mkxyTNkzQnlVXymkjn2ULSDZJ+K+kJSftVNV5Ju6bfa8/0iqQz2hJvRHjqZSJ7CeBpYCdgI+ARYHSJ8RwE7A3Mz5V9BzgrzZ8FXJjmjwBuAQR8AngglW8JLE4/h6f54S2Kdxtg7zS/OfAkMLqqMafzbpbmNwQeSHFcD0xI5ZcCp6b5ycClaX4CcF2aH52ulY2BHdM1NKyF18U/AtcAP0/LVY93CTCipqyS10Q611XAKWl+I2CLKsebi3sYWW8g27cj3pZ9kE6fgP2A2bnlKcCUkmPagXUTyUJgmzS/DbAwzV8GTKzdDpgIXJYrX2e7Fsd+M1nfapWPGXg/8DBZTwovABvUXhNkbxzul+Y3SNup9jrJb9eCOEcCvwD+Cvh5On9l403HX8J7E0klrwmytmu/I72UVPV4a2IcB9zbrnh9a6tv2wLP5JaXprIq2Toink3zzwFbp/m+Yi/lM6XbKHuR/ZVf2ZjTbaJ5wHLgdrK/zl+OiFW9nHtNXGn9n4Ct2hkv8EPgq8C7aXmriscLEMBtkuYq68YIqntN7AisAH6Sbh9eIWnTCsebNwGYkeZbHq8TySAR2Z8OlXuXW9JmwM+AMyLilfy6qsUcEasjYgzZX/pjgY+UG1HfJH0GWB4Rc8uOpUEHRMTewOHAaZIOyq+s2DWxAdnt5EsiYi/gdVKnsT0qFi8A6bnYUcBPa9e1Kl4nkr51Qlcsz0vaBiD9XJ7K+4q9rZ9J0oZkSeTfI+I/OiFmgIh4GbiL7NbQFsq66Kk9d9ld+OwPHCVpCXAt2e2tiyocLwARsSz9XA7cSJawq3pNLAWWRsQDafkGssRS1Xh7HA48HBHPp+WWx+tE0rd6unEpW74bmRPJnkP0lJ+Q3sr4BPCnVLWdDYyTNDy9uTEulTWdJJH1TvBERPyg6jFL6pK0RZp/H9nznCfIEsoxfcRbWhc+ETElIkZGxA5k1+adEXFsVeMFkLSppM175sn+LedT0WsiIp4DnpG0ayo6lKxHjUrGmzORtbe1euJqbbytfODT6RPZWw1Pkt0rP7vkWGYAzwLvkP2ldDLZPe5fAE8BdwBbpm1FNnDY08BjQHfuOF8AFqXp71sY7wFkVehHgXlpOqKqMQN7AL9J8c4HzknlO5F9sS4iu1WwcSrfJC0vSut3yh3r7PQ5FgKHt+HaOIS1b21VNt4U2yNpWtDzf6qq10Q6zxhgTroubiJ7i6nK8W5KVtP8QK6s5fG6ixQzMyvEt7bMzKwQJxIzMyvEicTMzApxIjEzs0KcSMzMrBAnEjMzK8SJxMz6JOlGSS9JuqHsWKy6nEjMrD8XASeUHYRVmxOJWQWl7imel7RzmXFExN3Aq72tk/RTSf/U3oisipxIrKNJulJS9DLdX3ZsBf13YFZEPF12IP04Dzhb0gfKDsTKtcHAm5hV3h3A8TVlb7fyhJI2ioiWnEPS+4FTgCNbcfyac82j9++BcRHxx/72jYjHJC0GjiPrs8mGKNdIbDB4KyKeq5lWAki6W9JUSd+W9IKyce+/J2nNtZ96P/2qpKclvalsTPHj8idIx7kk7bsCuDeVbyrpakmvpVtRUyT9PNWUTpD0oqSNa47175L660n6CLIOL++t2U+S/knZONpvSVoq6fya+L4vaaWkFZL+IfXqe7GklyX9QdI6CTcixkTE7r1M/SaRnJlkvc3aEOZEYkPBscAq4C+B04EzgM/n1n+LrDfl08jGMD8fuEzSX9cc5ziyHlMPZO0D6O8DBwOfJRsTZM+0HrLedv8MOLrnAOk20GfJutjvy4HA3Hhvj6rfBr6e4vso8LesO5LdsWTPM/YFLiAbQfEmsh6su8nGH79CaWyKJnkQGJu63rehqlXdGXvy1I4JuJIsSbxWM12Y1t8N3Fezz+3AFWl+U+BN4MCabX5I9oyC3HEerdlmM7JbaBNyZZsCLwFXpuUfA7fm1p9KNtzpBv18ppuAq3o51/8DvtjHPut8TrKEtwKYmSvbMMV7TAO/3zvScd4gG75gv5r1e5DVnnYu+1rwVN7kZyQ2GNwDTKopezk3/2jNuj8CH0zzo8nG6rhVUr4GsCGwpGa/2mFtd07brRkIKiJelzQ/t83lwMOSRkbEUrJxHq6KteOq9+Z9wPM1ZaOBjcnGlejLms8ZESFpOdk4Ez1l70h6ibWffUAR8akBNnkzF7MNUU4kNhi8ERGL+ln/Ts1ysPa2bs/PI4E/DLDf640GFhGPSHoYOEnSTWS3mI7rfy9eIBtAqVG9fc7+PnszbJl+rmjiMa3DOJHYUPc48BawfUTc2eC+T5N9UX8cWAxr3rjaPa3rcTnwVWAEcG9ELBzguL8BTqopeyLFeSjZSHdVsTuwLNaOD25DkBOJDQYbS/pQTdnqiBjwr+SIeFXS94DvSRLZbbLNgE8A70bEtH72fU3SdOBCSS+QDYX8P8j+4s/fJpsB/IDs+cgX6/g8s9Mxt4qIF3NxXgScL+mtFOdWwD4RcUkdx2yVA2nt+OPWAZxIbDD4FNmXeN4yYGSd+3+d7JnEmcAlwCtkY8x/p459zyR7wD6T7CH/vwBbkz0YB9YkgeuBY4DrBzpgZO0zHgQmsG77jClkD/K/TvbZngeuriPGlpC0CdkbaIeVFYNVg8dsN2ui1Gbk98B3I+L7ufJbgKUR8V/rPM54sn6uRkfE6pYEW5Ck04CjI2Jc2bFYuVwjMStA0l7AbmRvbm0OfC39vC6tH052+2ccWRuTukTErZIuJqt5/L7JYTfLO8CXyg7CyucaiVkBKZFcDuxK1p5lHnBmRMxN65eQvdn0vyLiwpLCNGspJxIzMyvEXaSYmVkhTiRmZlaIE4mZmRXiRGJmZoU4kZiZWSFOJGZmVogTiZmZFeJEYmZmhTiRmJlZIU4kZmZWyP8HRXACOD0FEo4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute the histogram of the energies\n",
    "bins = np.linspace(np.amin(Y), np.amax(Y), 200)    \n",
    "hists, edges = np.histogram(Y, bins=bins, density=True)\n",
    "\n",
    "#plot distribution of pool\n",
    "plt.bar(edges[:-1], 1000*hists, width=np.diff(edges), color='dimgray')\n",
    "plt.xlabel(\"Energy (cm$^{-1})$\", fontsize=14)\n",
    "plt.ylabel(\"$N/N_{tot}$  x  1000\", fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "bfeffe26",
   "metadata": {},
   "source": [
    "The plot shows that the majority of geometries correpond to electronic energies around 1800 $\\text{cm}^{-1}$. Let's now start to fit the data. Since active learning is an iterative procedure, one needs to be able to train models fast on the data. That is why we consider regression tree enesemble to fit the data.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b8bae4f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generalization error of the model when trained on 51750 training examples is: 27.195416083620533\n"
     ]
    }
   ],
   "source": [
    "X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.1, random_state=3)\n",
    "# define the Random forest model \n",
    "Trees = RFR(n_estimators=100, max_features=12, max_depth=None,\n",
    "            min_samples_leaf=1, min_samples_split=2, bootstrap=False, \n",
    "            n_jobs=-1)\n",
    "# fit the data \n",
    "Trees.fit(X, Y)\n",
    "# evaluate the test error on a test dataset \n",
    "Error_RFR = np.sqrt(mean_squared_error(Y_test, Trees.predict(X_test)))\n",
    "print(f\"Generalization error of the model when trained on {Y.shape[0]} training examples is: {Error_RFR}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cb2b8a0",
   "metadata": {},
   "source": [
    "## Active learning \n",
    "The question now is whether we can achieve the same accuracy with a less amount of data using active learning algorithms. Let's start by writing an object that implements a query-by-committee strategy using the trees of a random forest regressor as the committee. We start by defining a labeler, i.e., an object that returns the energies corresponding to certain molecular geometries. In real application, the labeler is an electronic-structure solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7877fa06",
   "metadata": {},
   "outputs": [],
   "source": [
    "Labeler = Y\n",
    "# define the size of the initially labeled data\n",
    "initial_size = 2000\n",
    "Indices = np.random.choice(np.arange(np.shape(Y)[0]), size=initial_size, replace=False) # to determine initial data\n",
    "Y = Y.tolist()\n",
    "for i in range(np.shape(Y)[0]):\n",
    "    if i not in Indices:\n",
    "        Y[i] = None\n",
    "Y = np.array(Y)\n",
    "# create a dataset object \n",
    "dataset = Dataset(X=X, y=Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "29ec357d",
   "metadata": {},
   "outputs": [],
   "source": [
    "class QBCR(QueryStrategy):\n",
    "    \"\"\"QBC query strategy object.\n",
    "    Args:\n",
    "    BatchSize: int. The number of elements to query in a single AL iteration.\n",
    "    Trees: an ensemble of trees sklearn instance composed of n trees.\n",
    "    dataset: an object containing the labeled and unlabeled data.\n",
    "    Attribs:\n",
    "    make_query: choose n=BatchSize data points to query from the pool.\n",
    "    \"\"\"\n",
    "    def __init__(self, *args, **kwargs):\n",
    "\n",
    "        super().__init__(*args, **kwargs)\n",
    "        self.BatchSize = kwargs.pop('BatchSize', 1)\n",
    "        self.Trees = kwargs.pop('Trees', None)\n",
    "        if self.Trees is None:\n",
    "            raise TypeError(\n",
    "                \"__init__() missing required keyword-only argument: 'Trees'\"\n",
    "            )\n",
    "\n",
    "    def _standard_deviation(self, predictions):\n",
    "        return np.std(predictions, axis=0)\n",
    "\n",
    "    def make_query(self):\n",
    "        # get the unlabeled entries, i.e. the pool data\n",
    "        unlabeled_entry_ids, X_pool = self.dataset.get_unlabeled_entries()\n",
    "        if len(unlabeled_entry_ids)>self.BatchSize:\n",
    "            # Predict the energies of the unlabeled data\n",
    "            predictions = np.array([estimator.predict(X_pool) for estimator in self.Trees.estimators_])\n",
    "            # compute the standard deviation of predictions for every unlabeled datapoint\n",
    "            stddev = self._standard_deviation(predictions)\n",
    "            ask_idx = np.flip(np.argsort(stddev))[0:self.BatchSize]\n",
    "        else:\n",
    "            #The number of remaining unlabeled data is less than batch size, query the remaining data\n",
    "            ask_idx = np.arange(len(unlabeled_entry_ids))\n",
    "\n",
    "        return unlabeled_entry_ids[ask_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "66906608",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0, size of labeled data: 2000, test error: 394.58520212764165 inverse centimeter\n",
      "Iteration: 1, size of labeled data: 3000, test error: 310.4995294371538 inverse centimeter\n",
      "Iteration: 2, size of labeled data: 4000, test error: 271.0478379481371 inverse centimeter\n",
      "Iteration: 3, size of labeled data: 5000, test error: 229.19607949561689 inverse centimeter\n",
      "Iteration: 4, size of labeled data: 6000, test error: 203.99898777265594 inverse centimeter\n",
      "Iteration: 5, size of labeled data: 7000, test error: 183.55387834992166 inverse centimeter\n",
      "Iteration: 6, size of labeled data: 8000, test error: 171.4468349584541 inverse centimeter\n",
      "Iteration: 7, size of labeled data: 9000, test error: 158.09904513347527 inverse centimeter\n",
      "Iteration: 8, size of labeled data: 10000, test error: 144.49418840535301 inverse centimeter\n",
      "Iteration: 9, size of labeled data: 11000, test error: 133.7115262938106 inverse centimeter\n",
      "Iteration: 10, size of labeled data: 12000, test error: 124.1293604726151 inverse centimeter\n",
      "Iteration: 11, size of labeled data: 13000, test error: 115.56733752452732 inverse centimeter\n",
      "Iteration: 12, size of labeled data: 14000, test error: 105.71522706840804 inverse centimeter\n",
      "Iteration: 13, size of labeled data: 15000, test error: 101.51092947576868 inverse centimeter\n",
      "Iteration: 14, size of labeled data: 16000, test error: 95.16085721268111 inverse centimeter\n",
      "Iteration: 15, size of labeled data: 17000, test error: 90.45814566906718 inverse centimeter\n",
      "Iteration: 16, size of labeled data: 18000, test error: 85.61012076202093 inverse centimeter\n",
      "Iteration: 17, size of labeled data: 19000, test error: 79.50964375529628 inverse centimeter\n",
      "Iteration: 18, size of labeled data: 20000, test error: 73.99404871813064 inverse centimeter\n",
      "Iteration: 19, size of labeled data: 21000, test error: 70.23851972629323 inverse centimeter\n",
      "Iteration: 20, size of labeled data: 22000, test error: 67.32901785109118 inverse centimeter\n",
      "Iteration: 21, size of labeled data: 23000, test error: 64.11562160070423 inverse centimeter\n",
      "Iteration: 22, size of labeled data: 24000, test error: 60.90575815464525 inverse centimeter\n",
      "Iteration: 23, size of labeled data: 25000, test error: 58.860086105908785 inverse centimeter\n",
      "Iteration: 24, size of labeled data: 26000, test error: 54.785659089662566 inverse centimeter\n",
      "Iteration: 25, size of labeled data: 27000, test error: 51.51711663053002 inverse centimeter\n",
      "Iteration: 26, size of labeled data: 28000, test error: 49.503245452302096 inverse centimeter\n",
      "Iteration: 27, size of labeled data: 29000, test error: 45.715142912983765 inverse centimeter\n",
      "Iteration: 28, size of labeled data: 30000, test error: 43.44226814004741 inverse centimeter\n",
      "Iteration: 29, size of labeled data: 31000, test error: 41.88621934219971 inverse centimeter\n",
      "Iteration: 30, size of labeled data: 32000, test error: 40.48607114057524 inverse centimeter\n",
      "Iteration: 31, size of labeled data: 33000, test error: 38.05916437599893 inverse centimeter\n",
      "Iteration: 32, size of labeled data: 34000, test error: 35.41858908177194 inverse centimeter\n",
      "Iteration: 33, size of labeled data: 35000, test error: 33.67757181479994 inverse centimeter\n",
      "Iteration: 34, size of labeled data: 36000, test error: 33.65896155897505 inverse centimeter\n",
      "Iteration: 35, size of labeled data: 37000, test error: 32.635704388535146 inverse centimeter\n",
      "Iteration: 36, size of labeled data: 38000, test error: 30.445362267248345 inverse centimeter\n",
      "Iteration: 37, size of labeled data: 39000, test error: 30.451610301656917 inverse centimeter\n",
      "Iteration: 38, size of labeled data: 40000, test error: 29.906711872830034 inverse centimeter\n",
      "Iteration: 39, size of labeled data: 41000, test error: 30.3366646528376 inverse centimeter\n",
      "Iteration: 40, size of labeled data: 42000, test error: 27.855053049858075 inverse centimeter\n",
      "Iteration: 41, size of labeled data: 43000, test error: 27.719765104626703 inverse centimeter\n",
      "Iteration: 42, size of labeled data: 44000, test error: 28.117241405142856 inverse centimeter\n",
      "Iteration: 43, size of labeled data: 45000, test error: 29.35087455654231 inverse centimeter\n",
      "Iteration: 44, size of labeled data: 46000, test error: 27.04933133988998 inverse centimeter\n",
      "Iteration: 45, size of labeled data: 47000, test error: 26.843780465450475 inverse centimeter\n",
      "Iteration: 46, size of labeled data: 48000, test error: 27.96906734388468 inverse centimeter\n",
      "Iteration: 47, size of labeled data: 49000, test error: 25.95221986463523 inverse centimeter\n",
      "Iteration: 48, size of labeled data: 50000, test error: 28.83917455117561 inverse centimeter\n",
      "Iteration: 49, size of labeled data: 51000, test error: 27.491416929738737 inverse centimeter\n"
     ]
    }
   ],
   "source": [
    "BatchSize = 1000 # number of elements to query at each active learning iteration\n",
    "iters = 0\n",
    "errors = []\n",
    "\n",
    "while dataset.len_unlabeled()>0: # While Pool isn't emtpy\n",
    "    # get labeled entries in the pool\n",
    "    X_L, y_L = dataset.get_labeled_entries()\n",
    "    # train a random forest regressor on the labeled data\n",
    "    Trees.fit(X_L, y_L)\n",
    "    # Initialize a QBC object\n",
    "    query_strategy = QBCR(Trees=Trees, BatchSize=BatchSize, dataset=dataset)\n",
    "    # Compute the root mean squared error of the RFR model on the OOP data\n",
    "    Error_RFR = np.sqrt(mean_squared_error(Y_test, Trees.predict(X_test)))\n",
    "    errors.append(Error_RFR)\n",
    "    print(f\"Iteration: {str(iters)}, size of labeled data: {X_L.shape[0]}, test error: {Error_RFR} inverse centimeter\")\n",
    "    # for analysis purposes we'll collect some data for AL Iteration 10 \n",
    "    if iters == 10: \n",
    "        y_sampled = y_L\n",
    "        yp = Trees.predict(X_test)\n",
    "    #Query new datapoints\n",
    "    ids = query_strategy.make_query()\n",
    "    # Label these datapoints and update the dataset\n",
    "    dataset.update(ids, Labeler[ids])\n",
    "    iters += 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28bd5c71",
   "metadata": {},
   "source": [
    "Now, let's plot these errors as a function of the number of training data. \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "566d5b83",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n_training_data = [initial_size]\n",
    "for i in range(1,len(errors)): n_training_data.append(n_training_data[i-1]+BatchSize)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3687b760",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArWElEQVR4nO3deXydZZ3//9cnJ3uzNVv3tA3dKAUKFKisbRkHqIzI/BwER2SUB4yKow46gs58Zxzm67iL8gVRFBW+XxcQVBgWtUJL2Qq00BW6pCtNl6RLli5Js3x+f5wr4bSmyelyck5y3s/H437kvq97OZ+7TfLJfV3XfV3m7oiIiABkJDsAERFJHUoKIiLSTUlBRES6KSmIiEg3JQUREemWmewATkR5ebmPGzcu2WGIiAwoS5Ys2eXuFT3tG9BJYdy4cSxevDjZYYiIDChmtvlo+1R9JCIi3ZQURESkW8KTgplFzOxNM3sybI83s1fNrMbMHjaz7FCeE7Zrwv5xiY5NREQO1x9PCp8F3o7Z/gZwl7tPAPYCN4Xym4C9ofyucJyIiPSjhCYFMxsNvA/4Sdg2YA7waDjkQeADYf3qsE3Yf1k4XkRE+kminxS+B3wR6AzbZUCDu7eH7a3AqLA+CngHIOxvDMcfxsxuMbPFZra4vr4+gaGLiKSfhCUFM7sKqHP3JSfzuu5+v7vPcPcZFRU9drMVEZHjlMgnhQuB95vZJuDXRKuNvg+UmFnX+xGjgdqwXguMAQj7i4HdiQhszY5mvvGH1TQebEvE5UVEBqyEJQV3/5K7j3b3ccB1wHPu/vfAfOCD4bAbgcfD+hNhm7D/OU/QZA9b9hzgvgXr2bhrfyIuLyIyYCXjPYXbgdvMrIZom8EDofwBoCyU3wbckagAxpblA7B5t5KCiEisfhnmwt0XAAvC+gbgvB6OaQH+rj/iqSqNJoUtuw/0x8eJiAwYaflGc25WhOFFuWzeo6QgIhIrLZMCQFVZvqqPRESOkLZJYWxpPptVfSQicpj0TQpl+dQ1t3LgUHvfB4uIpIk0TgpDgGj3VBERiUrjpNDVLVVJQUSkS/omhdLwpKCkICLSLW2TQnF+FsV5WWzeox5IIiJd0jYpAIwrUw8kEZFYaZ0UqsqGKCmIiMRI66QwtjSf2oaDtHV09n2wiEgaSOukUFWWT0ens63hYLJDERFJCWmdFMaFdxU2qQpJRARI86TQ9a7CFo2BJCICpHlSqCzMITcrQ43NIiJBWicFM6OqNF9DaIuIBGmdFCA6BpKG0BYRiVJSKM1ny54DJGg6aBGRASVhScHMcs3sNTNbZmarzOw/Q/nPzWyjmS0Ny/RQbmZ2t5nVmNlyMzs7UbHFGluWT0tbJ3XNrf3xcSIiKS2RczS3AnPcfZ+ZZQEvmtkzYd+/uPujRxx/JTAxLOcD94WvCdU1hPbm3QcYVpSb6I8TEUlpCXtS8Kh9YTMrLL3V0VwNPBTOWwSUmNmIRMXXpatb6ia1K4iIJLZNwcwiZrYUqAPmufurYddXQxXRXWaWE8pGAe/EnL41lB15zVvMbLGZLa6vrz/hGEeW5BHJMA2hLSJCgpOCu3e4+3RgNHCemU0DvgRMAc4FSoHbj/Ga97v7DHefUVFRccIxZkUyGFWSp26pIiL0U+8jd28A5gNXuPv2UEXUCvwMOC8cVguMiTltdChLuLFl+XqrWUSExPY+qjCzkrCeB7wXWN3VTmBmBnwAWBlOeQL4aOiFNBNodPftiYov1tiyfI1/JCJCYnsfjQAeNLMI0eTziLs/aWbPmVkFYMBS4BPh+KeBuUANcAD4WAJjO8zY0iE0Hmyj8UAbxflZ/fWxIiIpJ2FJwd2XA2f1UD7nKMc7cGui4ulNVeiBtHnPfs7IL0lGCCIiKSHt32iGd4fQ1sB4IpLulBSAqtLwpKDGZhFJc0oKQF52hMrCHD0piEjaU1IIxpZpCG0RESWFYGzZEL3VLCJpT0khGFuaz46mFlraOpIdiohI0igpBF3dUreoCklE0piSQjBW3VJFRJQUuowrU7dUERElhaAkP5ui3Ew9KYhIWlNSiDG2bIi6pYpIWlNSiFGlIbRFJM0pKcQYV5bP1r0Hae/oTHYoIiJJoaQQY2zpENo7nW0NLckORUQkKZQUYsQOoS0iko6UFGKcUlEAwKptTUmOREQkOZQUYlQU5nDqiCLmr65LdigiIkmhpHCEOVMqWLx5L40H2pIdiohIv0tYUjCzXDN7zcyWmdkqM/vPUD7ezF41sxoze9jMskN5TtiuCfvHJSq23syZUklHp7NwXX0yPl5EJKkS+aTQCsxx9zOB6cAVZjYT+AZwl7tPAPYCN4XjbwL2hvK7wnH9bvqYoQzNz1IVkoikpYQlBY/aFzazwuLAHODRUP4g8IGwfnXYJuy/zMwsUfEdTSTDmDW5kgVr6+no9P7+eBGRpEpom4KZRcxsKVAHzAPWAw3u3h4O2QqMCuujgHcAwv5GoKyHa95iZovNbHF9fWKqeGZPqWTP/kMs29qQkOuLiKSqhCYFd+9w9+nAaOA8YMpJuOb97j7D3WdUVFSc6OV6dOnECiIZxnNvqwpJRNJLv/Q+cvcGYD7wHqDEzDLDrtFAbVivBcYAhP3FwO7+iO9IxflZnFM1lOfUriAiaSaRvY8qzKwkrOcB7wXeJpocPhgOuxF4PKw/EbYJ+59z96RV6s+eUslb25vY0aghL0QkfSTySWEEMN/MlgOvA/Pc/UngduA2M6sh2mbwQDj+AaAslN8G3JHA2Pp02amVAMxfo6cFEUkfmX0fcnzcfTlwVg/lG4i2LxxZ3gL8XaLiOVYTKwsYVZLHc6vruP68qmSHIyLSL/RG81GYGXOmVPLiul20tHUkOxwRkX6hpNCLOVMqOdjWwasb9yQ7FBGRfqGk0Iv3nFJGblaG3m4WkbShpNCL3KwIF5xSznOr60hiRygRkX6jpNCHOVMq2bLnAOvrNfGOiAx+Sgp9mD0l2jX1udU7kxyJiEjiKSn0YVRJHlOGF+rtZhFJC32+pxDeSv4oMC72eHf/TMKiSjGzp1Ty44UbaGppoyg3K9nhiIgkTDxPCk8TTQgrgCUxS9qYM6WS9k7nhbW7kh2KiEhCxfNGc66735bwSFLYWWNKKMnPYv6aOt53xohkhyMikjDxPCn8XzO72cxGmFlp15LwyFJIZiSDCyeU88K6enVNFZFBLZ6kcAj4FvAK71YdLU5kUKno0okV7GxqZc3O5mSHIiKSMPFUH30emODuaV2hfvGkcgAWrq1nyvCiJEcjIpIY8Twp1AAHEh1IqhtRnMekYQUsVGOziAxi8Twp7AeWmtl8oLWrMJ26pHa5eGIF/3fRZg4e6iAvO5LscERETrp4nhR+D3wVeJk07ZLa5ZJJFRxq7+TVjUmZJVREJOHieVJ4FGhx9w4AM4sAOQmNKkWdP76UnMwMFq7dxazJlckOR0TkpIvnSeFZIC9mOw/4c2LCSW25WRHOG1/KwnX1yQ5FRCQh4kkKue6+r2sjrOf3dZKZjTGz+Wb2lpmtMrPPhvKvmFmtmS0Ny9yYc75kZjVmtsbMLj+eG0q0SydVUFO3j20NB5MdiojISRdPUthvZmd3bZjZOUA8vxHbgc+7+1RgJnCrmU0N++5y9+lheTpcdypwHXAacAXwg1BVlVIumVQBRLumiogMNvEkhc8BvzGzF8zsReBh4NN9neTu2939jbDeDLwNjOrllKuBX7t7q7tvJNoV9rw44utXEysLGF6UqyokERmUjpoUzGwkgLu/DkwBPgl8AjjV3Y+p95GZjQPOAl4NRZ82s+Vm9lMzGxrKRgHvxJy2lR6SiJndYmaLzWxxfX3//2I2My6eWM6L63bR0akhL0RkcOntSeEnZrbIzL4OXAisdveV7t52LB9gZgXAY8Dn3L0JuA84BZgObAe+cyzXc/f73X2Gu8+oqKg4llNPmksmVdDU0s6yrQ1J+XwRkUQ5alJw97nALGABcA2wyMx+G/5Sr4rn4maWRTQh/MLdfxuuu9PdO9y9E/gx71YR1QJjYk4fHcpSzkUTyjFTu4KIDD69tim4e4u7/8HdP+vuM4iOg5QJ3GNmr/V2rpkZ8ADwtrt/N6Y8duzpa4CVYf0J4DozyzGz8cBEoNfPSJahQ7I5Y3SJkoKIDDrxvLzWLTQA/4Boz6DsPg6/ELgBWGFmS0PZl4HrzWw64MAm4B/DtVeZ2SPAW0R7Lt3a9cJcKrp0Yjn3zK+h8UAbxfmajU1EBodjSgpdzGyFu5/e2zHu/iJgPex6updzvkp0SI2Ud/GkCu5+roaX1u9i7umaeEdEBoejJgUz+9uj7QKGJyacgWP6mBIKczJ5YV29koKIDBq9PSk8DPyCaDXPkXITE87AkRXJ4IIJZSxcuwt3J9qEIiIysPWWFJYD33b3lUfuMLO/SlxIA8clkyr446qdrK/fz4TKgmSHIyJywnrrffQ5oOko+645+aEMPJdM1JAXIjK49PaewgvuvuUo+9JujuaejCnNp7p8CPPX1CU7FBGRkyKesY+kF1edMYIXa3axoX5f3weLiKQ4JYUTdMN7xpEVyeDHL2xMdigiIies16RgZhlmdm1/BTMQVRTm8MFzRvPYG1upb27t+wQRkRTW1zAXncAX+ymWAevmi6tp6+jkwZc3JTsUEZETEk/10Z/N7AthJrXSriXhkQ0g48uHcPnU4Tz0yib2t7YnOxwRkeMWT1L4EHArsBBYEhb1PjrCP15aTVNLO79+/Z2+DxYRSVF9JgV3H9/DUt0fwQ0kZ1UN5bzxpfz0xY20dXQmOxwRkePSZ1Iwsywz+4yZPRqWT4d5EuQIn7i0mtqGgzy1fHuyQxEROS7xVB/dB5xDGDI7rN+XyKAGqlmTKplYWcAPn1+Pu6bqFJGBJ56kcK673+juz4XlY8C5iQ5sIMrIMG65pJrVO5pZuG5XssMRETlm8SSFDjM7pWvDzKqBlJ38Jtmunj6KYUU53L9wfbJDERE5ZvEkhS8A881sgZk9DzxHdFpO6UF2ZgYfv3A8L9XsZsXWxmSHIyJyTPp6ozkCnEl0vuTPAP8ETHb3+X1dOLzXMN/M3jKzVWb22VBeambzzGxd+Do0lJuZ3W1mNWa23MzOPuG7S5Lrz6+iMCeTH+lpQUQGmL7eaO4Arnf3VndfHpZ4x3JoBz7v7lOBmcCtZjYVuAN41t0nAs+GbYAriSaficAtDODG7KLcLD48s4qnV2xn9Y6jjT4uIpJ64qk+esnM7jGzi83s7K6lr5Pcfbu7vxHWm4G3gVHA1cCD4bAHgQ+E9auBhzxqEVBiZgN2nstPXnoKhblZfPWpt9UTSUQGjN5mXusyPXy9M6bMgTnxfoiZjQPOAl4Fhrl7V0f+HcCwsD4KiH0deGsoG5Cd/kvys/nsZRO588m3mL+mjjlThvV9kohIksXTpvCEu88+YjmWhFAAPAZ8zt0Pq0vx6J/Qx/RntJndYmaLzWxxfX1qz3h2w3vGUl0+hP/91Nt6y1lEBoS42hSO9+LhzefHgF+4+29D8c6uaqHwtWvaslpgTMzpo0PZkTHd7+4z3H1GRUXF8YbWL7IiGXx57qlsqN/PLxZtTnY4IiJ9SlibgpkZ8ADwtrt/N2bXE8CNYf1G4PGY8o+GXkgzgcaYaqYB67JTK7lwQhnfe3YdjQfakh2OiEiv4kkK04HTiLYpfCcs347jvAuBG4A5ZrY0LHOBrwPvNbN1wF+FbYCngQ1ADfBj4FPHcB8py8z4t/dNpelgG99/dl2ywxER6VWfDc3uPvt4LuzuLwJ2lN2X9XC8Ex2ie9A5dUQRHzp3DA+9somPzKyiuqIg2SGJiPQonlFSh5nZA2b2TNieamY3JT60weW2904mNyvCfz+9OtmhiIgcVTzVRz8H/giMDNtrgc8lKJ5Bq6Iwh0/NPoU/v72Tl2s0WJ6IpKZ4kkK5uz8CdAK4ezsaEO+4fPzC8YwemsedT75FR6deaBOR1BNPUthvZmWE9wm6egYlNKpBKjcrwpfnnsrqHc382+9X6E1nEUk58bzRfBvR7qKnmNlLQAXwwYRGNYjNPX0Et84+hXvnryc3K8K/XzWVaO9dEZHki6f30RtmdikwmWhvojXurg73J+ALfz2ZA4c6+NlLm8jPjvAvl09JdkgiIkB8Twpd7QirEhxL2jAz/v2qqbS0dXDv/PXkZ2dy6+wJyQ5LRCS+pCAnn5nxvz9wOgcPdfCtP64hNyvCTReNT3ZYIpLmlBSSKJJhfPvvzqSlrZP/evIt8rMjXH9eVbLDEpE0Fs/La8/GUybHJzOSwd3Xn8WsyRV8+Xcr+J9l25IdkoiksaMmBTPLNbNSoNzMhoZpNEvD3Aij+i3CNJCdmcEPP3IO544t5fOPLOP1TXuSHZKIpKnenhT+EVgCTAlfu5bHgXsSH1p6yc2K8KMbzmH00DxufmgxG+r3JTskEUlDR00K7v59dx8PfMHdq919fFjOdHclhQQYOiSbn33sXDLM+NjPX2f3vninwxYROTnieaN5h5kVApjZv5nZb+OZT0GOz9iyIfzkxhnsaGzh5ocW09KmEUVEpP/EkxT+l7s3m9lFROc/eAC4L7Fhpbezq4byvQ9N5813Gvjnh5fSqXGSRKSfxJMUuv5UfR9wv7s/BWQnLiQBuPL0Efzr3FN5ZuUOvv4HDbctIv0jnvcUas3sR8B7gW+YWQ7xJRM5QTddNJ4tew5w/8INjCnN54aZY5MdkogMcvH8cr+W6HwKl7t7A1AK/Esig5KoruEw5kyp5CtPrGLBmrpkhyQig1yfScHdDwB1wEWhqB3oc7JhM/upmdWZ2cqYsq+YWe0RczZ37fuSmdWY2Rozu/zYb2Vw6nq5bdKwQj79yzdZs6M52SGJyCAWzxvN/wHcDnwpFGUB/y+Oa/8cuKKH8rvcfXpYng6fMRW4DjgtnPMDM4vE8RlpoSAnkwdunEF+doSP//x16ppbkh2SiAxS8VQfXQO8H9gP4O7bgMK+TnL3hUC8r+ZeDfza3VvdfSNQA5wX57lpYWRJHg/ceC579h/iloeWqKuqiCREPEnhkEenCOuaeW3ICX7mp81seaheGhrKRgHvxByzlaMMpWFmt5jZYjNbXF9ff4KhDCynjy7me9dNZ9nWBj7/yDJ1VRWRky6epPBI6H1UYmY3A38GfnKcn3cfcAowHdgOfOdYL+Du97v7DHefUVFRcZxhDFyXnzacL105hadWbOe789YmOxwRGWTimXnt22b2XqCJ6Oxr/+7u847nw9x9Z9e6mf0YeDJs1gJjYg4dHcqkBzdfXM3GXfu5Z34NY8vy+bsZY/o+SUQkDn0mBTP7hrvfDszroeyYmNkId98eNq8BunomPQH80sy+C4wEJgKvHev104WZcefV03hnz0Fuf2w57nDtuUoMInLi4qk+em8PZVf2dZKZ/Qp4BZhsZlvN7Cbgm2a2wsyWA7OBfwZw91XAI8BbwB+AW91dLam9yIpkcP9Hz+GiiRV88bHl/HjhhmSHJCKDgEXbkHvYYfZJ4FNANbA+Zlch8JK7fyTx4fVuxowZvnjx4mSHkVSH2jv550eW8tTy7dw6+xS+8NeTMbNkhyUiKczMlrj7jJ729VZ99EvgGeBrwB0x5c3urllgUkR2ZgZ3X3cWRblZ3Dt/PY0H27jz/dPIyFBiEJFjd9Sk4O6NQCNwff+FI8cjkmH89zXTKM7L4ofPr6fpYDvfufZMsiIaokpEjk08A+LJAGBm3HHlFIrzsvjGH1bT3NLG//nw2RTk6L9YROKnPyUHmU/OOoWv/e3pLFy3iw/e9zLv7DmQ7JBEZABRUhiErj+vip9/7Fy2NRzk6ntf4rWNagISkfgoKQxSF0+s4Pe3XkhJXhZ//5NFPPz6lmSHJCIDgJLCIFZdUcDvPnUhM6vLuP2xFdz5P2/R3tGZ7LBEJIUpKQxyxflZ/OwfzuVjF47jpy9t5OMPLqbhwKFkhyUiKUpJIQ1kRjL4j785ja/97em8sn4X77v7RZZsVjuDiPwlJYU0cv15VTz6iQvIyIBrf7SIe+fXaPhtETmMkkKaOXNMCU995mKumDacb/1xDTf+7DXqm1uTHZaIpAglhTRUlJvFPdefxX9fczqvbdzD3Ltf4KWaXckOS0RSgJJCmjIzPnx+FY9/+kKK87L4yAOvcs9z6zjaAIkikh6UFNLclOFFPPHpC3n/mSP59p/W8vnfLKO1XaOWi6QrDYwj5Gdn8r0PTae6vIC7/ryWrXsO8qMbzmHokOxkhyYi/UxPCgJEq5M++1cT+f5101m6tYFrfvAS6+v3JTssEelnSgpymKunj+JXN59Pc0s719z7Ei+vVwO0SDpRUpC/cM7YUn5/64VUFuXy0Qde4+5n17GjsSXZYYlIP0hYUjCzn5pZnZmtjCkrNbN5ZrYufB0ays3M7jazGjNbbmZnJyouic+Y0nwe++QFzJpcwXfnreU9X3+Wv//JIh5dspV9re3JDk9EEuSoczSf8IXNLgH2AQ+5+7RQ9k1gj7t/3czuAIa6++1mNhf4J2AucD7wfXc/v6/P0BzN/WND/T5+v3Qbv3tzK+/sOUheVoS/Pm0YHzp3DBecUp7s8ETkGPU2R3PCkkL44HHAkzFJYQ0wy923m9kIYIG7TzazH4X1Xx15XG/XV1LoX+7Oks17+e2btTy1fDuNB9u4eGI5d1w5hdNGFic7PBGJU29Job/bFIbF/KLfAQwL66OAd2KO2xrK/oKZ3WJmi81scX19feIilb9gZswYV8p/X3M6r375Mv7XVVNZUdvIVf/nRW57ZCm1DQeTHaKInKCkNTR79BHlmB9T3P1+d5/h7jMqKioSEJnEIzcrwk0Xjef5f5nNP15yCk8u387sby/ga8+8TePBtmSHJyLHqb+Tws5QbUT4WhfKa4ExMceNDmWS4orzsrjjyinM/8IsrjpjBPcv3MDF33iO/3ryLb3nIDIA9XdSeAK4MazfCDweU/7R0AtpJtDYV3uCpJZRJXl899rpPPlPF3HxxAoefHkTl33neT7840U8tXw7h9o145vIQJDI3ke/AmYB5cBO4D+A3wOPAFXAZuBad99jZgbcA1wBHAA+5u59tiCroTl11TW38JvFW/nlq1uobThIeUEOHzp3NNefV8XoofnJDk8krSWt91GiKSmkvo5OZ+G6en6xaAvPrd4JwJwpw/joe8Zy0YRyMjIsyRGKpJ/ekoIGxJOEimQYsydXMntyJbUNB/nlq5v59Wvv8Oe3d1JdPoSPzBzL/3fOaIrzspIdqoigJwVJgtb2Dp5esZ2HXtnMm1sayMuK8OHzq/jUrFMoK8hJdngig56qjyRlraxt5KcvbeT3b9aSlxXhpourufni8RTm6slBJFGUFCTl1dQ18915a3l6xQ6G5mfxqVkTuOE9Y8nNiiQ7NJFBR0lBBowVWxv55h9X88K6XQwvyuXTcybwwXNGKzmInERKCjLgvLJ+N9/642re2NJAZWEON19czYfPr2JIjvpGiJwoJQUZkNydV9bv5t4FNbxUs5vivCz+4YJx/MMF4zRVqMgJUFKQAe/NLXv5wYL1zHtrJ/nZEa6dMYa/OXMkZ40p0bsOIsdISUEGjTU7mrlvQQ1PrdhOW4czvCiXy08bxhXTRnDe+FIiShAifVJSkEGnqaWN596u45mV21mwpp7W9k7KhmTz3qnDuHBCOTOry6go1DsPIj1RUpBB7cChdhasqeeZlTuYv7que7rQCZUFzKwuZWZ1GTOryyjXi3EigJKCpJH2jk5WbWvilQ27WbRhN69v3MP+Qx0ATBleyKWTKrhkUgUzxg0lJ1PdXCU9KSlI2mrv6GTltiZeXr+LF9buYvHmPbR1OHlZEWZWl3JJSBLV5UOIDtYrMvgpKYgE+1vbWbRhN8+vrWfh2no27T4AROeDuGRSBZdOKueCCeUUaZgNGcSUFESOYsvuAyxcF00QL6/fzb7WdiIZxlljSji/upQpw4uYMryQ8eVDyIwkbfZakZNKSUEkDm0dnby5pYGFa+tZuK6eVdua6OiM/nxkRzKYUFnAlOGFnD66mCumDWdEcV6SIxY5PkoKIsehtb2D9XX7Wb2jiTU7mlm9o5k1O5rZ0dSCGZw/vpSrp4/iymnDKcnXG9YycCgpiJxEm3bt5/Gl23h8WS0b6veTFTEunVTJ+6eP5NJJFZowSFJeyiUFM9sENAMdQLu7zzCzUuBhYBywiej8zXt7u46SgiSTu7NqWxOPL63liWXb2NnUSiTDOLuqhFmTK7l0UgVTRxRpGA5JOamaFGa4+66Ysm8Ce9z962Z2BzDU3W/v7TpKCpIqOjqdN7bsZcGaOp5fW8/K2iYAygtyuGRiOdOrSpg8rJDJwwtV1SRJN1CSwhpglrtvN7MRwAJ3n9zbdZQUJFXVN7fywrp6nl9bzwvrdrFn/6HufZWFOUweXsjkYYWcO76USydVaL4I6VepmBQ2AnsBB37k7vebWYO7l4T9Buzt2j7i3FuAWwCqqqrO2bx5c7/FLXI83J2dTa2s2dnMmh1NrNmxjzU7m1i3cx+t7Z3kZ0eYPaWSudNGMHtKBfnZmjNCEisVk8Iod681s0pgHvBPwBOxScDM9rr70N6uoycFGcjaOjp5dcMenl65nT+t2sGufYfIzcrg0kkVXDZlGKePLmZCZQFZej9CTrKUSwqHBWD2FWAfcDOqPpI01dHpvL5pD8+s2M4zK3dQ19wKQHZmBqeOKGLayCKmjSpm6ogixpUNoThfPZzk+KVUUjCzIUCGuzeH9XnAncBlwO6YhuZSd/9ib9dSUpDBqLPT2bBrP6u2NbKytpGVtU2s3NZIc0t79zHFeVlUleZHl7Lo11EleYwamseokjy1UUivUi0pVAO/C5uZwC/d/atmVgY8AlQBm4l2Sd3T27WUFCRduDtb9hxg9Y5mtuw+wJY9B9i85wDv7DnA1r0HaOs4/Oe4vCC7O0lMHlbEmWOKOXN0iaYxFaD3pNDvLVruvgE4s4fy3USfFkTkCGbG2LIhjC0b8hf7Ojqd7Y0Hqd17kNqGmK8NB3lrWxPPrNxB199+VaX5nDmmhDNHF3NWVQmnjSzWU4UcRt0cRAa4SIYxemg+o4fm97i/uaWNlbVNLNvawLJ3Gnhj817+Z9k2IDqm07RRRcwYV8rZVUM5Z+xQzViX5pLe0HwiVH0kcnzqmlt4c0s0QSzevJcVWxs51NEJRIcRr64YQnX5EKorChhfPoTqiiGMLM7T29mDREpVH4lI8lUW5nL5acO5/LThQHTwv5W1TSzZvIeVtU1s3LWfx96o7Z7aFCDDok8lAEZIDhYtL8nLpqwgm/KCHMoKsqkIX8eVDeG0UcWMLM7VJEYDhJKCiJCTGeGcsdHqoy7uTn1zKxt27WdD/X62NRyk0x3v3h/92unO3v2H2L3/ELv2tbJuZzO79h3qfvIAKMnPYuqIIk4bWcTUkUVHHXY8w4z87AgFOZnk50QYkp1JXlaEjAzD3Wlt76SlrYOWtk4OtnVwqL2TMaV5euHvJNK/pIj0yMyoLMqlsiiXmdVlx3Suu9PU0k5N3T7e2tbIW9ubeGtbEw+9spnW9s6+L3CEnMwMDnV00lNtdyTDmDSskOljSjhrTAnTq0qYUFFwzFVdre0d1DW1sqOphe2NLexsbGFHUwt7DxzilIoCpo0qZtrIIsoKEt/m4u5s3n2A5bWNFOdlMbI4l5EleQzJSfyvbLUpiEi/ae/oZMOu/exqboUefmd3dDoHDnVw4FA7+1vf/drS3kFOJIPc7Ah5WRFys6JfIxnGup3NvPlOtBG9KbzLUZCTSWVhDhkZRoZFn0AiGdGl053Wtk5a2qNPHC1tHbS2d3Koh2SVlxWhOC+LHU0t3WUji3OjCWJUMUPzs8iMZJCZYWRGjMyM6Pq+1nbq97VS19RKfXMrdc0t3S8kTqwsZPLwAiYNK2TK8CLGlw8hK2Js2n2ARRt2s2jDbl7dsOewz+xSnJfFyJI8Rhbn8v7pI7l6+qjj+n9Qm4KIpITMSAaThhUyaVjhSb92Z6ezcfd+lm5pYNnWBvYeaKOz0+nodDrco+vuGJAbEktuVgY5mdH1IdkRhhXlMrw4ugwryqUoNxMzo/Fg2+EvE9Y28qe3dvYZU2FOJhVFOVQW5nDm6BI63Fm7o5n5a+q6Z/XLihiFuVndgyaWF+Qws7qUmdVlnFVVwsFDHdQ2HGR7YwvbGg6yreEgtQ0t1Ickc7LpSUFE5Dh0PcW0d3bS3uG0dzodnZ20dTj52REqC3PJy+75HZDW9g421O9nzY5m1uxsZldzK9OrSphZXUZ1+ZCEN8rrSUFE5CTLz8487gbunMwIp44o4tQRRSc5qhOn4RdFRKSbkoKIiHRTUhARkW5KCiIi0k1JQUREuikpiIhINyUFERHppqQgIiLdBvQbzWZWT3TqToByYFcSw0kW3Xd60X2nn0Tc+1h3r+hpx4BOCrHMbPHRXtsezHTf6UX3nX76+95VfSQiIt2UFEREpNtgSgr3JzuAJNF9pxfdd/rp13sfNG0KIiJy4gbTk4KIiJwgJQUREek24JOCmV1hZmvMrMbM7kh2PMfLzH5qZnVmtjKmrNTM5pnZuvB1aCg3M7s73PNyMzs75pwbw/HrzOzGmPJzzGxFOOduS/TUTnEwszFmNt/M3jKzVWb22VA+2O8718xeM7Nl4b7/M5SPN7NXQ6wPm1l2KM8J2zVh/7iYa30plK8xs8tjylP258LMImb2ppk9GbbT5b43he/FpWa2OJSl3ve6uw/YBYgA64FqIBtYBkxNdlzHeS+XAGcDK2PKvgncEdbvAL4R1ucCzxCd+nwm8GooLwU2hK9Dw/rQsO+1cKyFc69MgXseAZwd1guBtcDUNLhvAwrCehbwaojxEeC6UP5D4JNh/VPAD8P6dcDDYX1q+J7PAcaHn4VIqv9cALcBvwSeDNvpct+bgPIjylLue32gPymcB9S4+wZ3PwT8Grg6yTEdF3dfCOw5ovhq4MGw/iDwgZjyhzxqEVBiZiOAy4F57r7H3fcC84Arwr4id1/k0e+eh2KulTTuvt3d3wjrzcDbwCgG/327u+8Lm1lhcWAO8GgoP/K+u/49HgUuC38FXg382t1b3X0jUEP0ZyJlfy7MbDTwPuAnYdtIg/vuRcp9rw/0pDAKeCdme2soGyyGufv2sL4DGBbWj3bfvZVv7aE8ZYSqgbOI/tU86O87VKEsBeqI/mCvBxrcvT0cEhtr9/2F/Y1AGcf+75EKvgd8EegM22Wkx31DNPH/ycyWmNktoSzlvtePb9Zp6Xfu7mY2KPsPm1kB8BjwOXdviq0KHaz37e4dwHQzKwF+B0xJbkSJZ2ZXAXXuvsTMZiU5nGS4yN1rzawSmGdmq2N3psr3+kB/UqgFxsRsjw5lg8XO8FhI+FoXyo92372Vj+6hPOnMLItoQviFu/82FA/6++7i7g3AfOA9RKsIuv5Qi421+/7C/mJgN8f+75FsFwLvN7NNRKt25gDfZ/DfNwDuXhu+1hH9Q+A8UvF7PdmNLyeyEH3S2UC0samrYem0ZMd1AvczjsMbmr/F4Y1Q3wzr7+PwRqjX/N1GqI1EG6CGhvVS77kRam4K3K8Rrfv83hHlg/2+K4CSsJ4HvABcBfyGwxtcPxXWb+XwBtdHwvppHN7guoFoY2vK/1wAs3i3oXnQ3zcwBCiMWX8ZuCIVv9eT/o91Ev6x5xLttbIe+Ndkx3MC9/ErYDvQRrQ+8Cai9afPAuuAP8f85xtwb7jnFcCMmOt8nGjDWw3wsZjyGcDKcM49hLfZk3zPFxGtZ10OLA3L3DS47zOAN8N9rwT+PZRXhx/sGqK/KHNCeW7Yrgn7q2Ou9a/h3tYQ09sk1X8uODwpDPr7Dve4LCyrumJLxe91DXMhIiLdBnqbgoiInERKCiIi0k1JQUREuikpiIhINyUFERHppqQgac3MZpnZBcdx3gwzuzuO414+vsj6vO4CM+t1Mncz+5yZ5Sfi82XwUlKQdDcL6DEpxLxl+xfcfbG7f6avi7v7MSeck+hzgJKCHBMlBRlwzGycmb1tZj8O8xH8yczywr7pZrYojEH/u67x6Y92HeATwD+HMe4vNrOfm9kPzexV4Jtmdp6ZvRLG/3/ZzCaHc2fFzAfwFYvOh7HAzDaY2WdiPmNfzPELzOxRM1ttZr/oGu/ezOaGsiVhHPwne4g1z8x+He77d0TfhO7ad5+ZLbbD52b4DDASmG9m8492nMhfSPabflq0HOtCdDiQdmB62H4E+EhYXw5cGtbv5IghNHq41leAL8Rs/xx4EoiE7SIgM6z/FfBYWJ/Fu2/kfoXosAU5QDnR8Xmywr59Mcc3Eh2TJgN4hegb3blER70cH477Vdd1j4jzNuCnYf2McP8zwnbXW7ARYAFwRtjeRMz4/Uc7TouW2EVPCjJQbXT3pWF9CTDOzIqJjin0fCh/kOjkRcfqNx4dxRSig7D9xqIz4t1FdNydnjzl0fH9dxEd1GxYD8e85u5b3b2T6JAe44iOjrrBo/MCQDQp9OQS4P8BuPtyosmvy7Vm9gbRoTNOIzoJTU/iPU7SmJKCDFStMesdnNxh4PfHrP8XMN/dpwF/Q/Qv++ON56THbGbjgS8Al7n7GcBTPcUY73EiSgoyaLh7I7DXzC4ORTcAz/dyCkAz0alAj6aYd4cg/ocTCrBna4Bqe3f+4Q8d5biFwIcBzGwa0SokiFZv7QcazWwYcGXMObH31ttxIt2UFGSwuRH4lpktB6YTbVfAzO40s/f3cPz/ANd0NTT3sP+bwNfM7E0SMCmVux8kOhfxH8xsCdFf5I09HHofUGBmbxO9pyXh/GVEq4NWE533+KWYc+4P153fx3Ei3TRKqkiSmVmBu+8LvZHuBda5+13JjkvSk54URJLv5jBf8yqi1VU/Sm44ks70pCAiIt30pCAiIt2UFEREpJuSgoiIdFNSEBGRbkoKIiLS7f8HAeLlNgBtlioAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(n_training_data, errors)\n",
    "plt.xlabel(\"no. training data\")\n",
    "plt.ylabel(\"test error 1/cm\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f23449f1",
   "metadata": {},
   "source": [
    "Now we plot the distribution of the points collected by QBC strategy\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fd483061",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAESCAYAAADuVeJ5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAb4klEQVR4nO3df7QcZZ3n8ffHQBAB5UcyDockJGBAI6sQrgFHUXcQJuiQjEccw4jAipMD1+zIKOOEYeTsMqiLA4zOmQQJLIK7IwGZFe+y0fgLx5VVSILhR8JELjFKIj8CooI6QPC7f9Rzk07Tuffp3O6u6u7P65w+t+qpp6u+fW/d+vZTVc9TigjMzMzG8pKyAzAzs+7ghGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaWZY+yA2iVSZMmxfTp08sOw8ysq6xZs+aJiJicU7dnEsb06dNZvXp12WGYmXUVST/JretTUmZmlsUJw8zMsjhhmJlZFicMMzPL4oRhZmZZnDDMzCyLE4aZmWVxwjAzsyxOGDZug4ODDA4Olh2GmbWZE4aZmWVxwjAzsyxOGGZmlsUJw8zMsjhhmJlZFicMMzPL4oRhZmZZnDDMzCyLE4aZmWVxwjAzsyxOGGZmlsUJw8zMsjhhmJlZllIShqS5kjZIGpa0eBd1/lTSeknrJH2x0zGamdnO9uj0BiVNAJYAJwGbgVWShiJifU2dmcCFwJsi4ilJv9fpOM3MbGdltDDmAMMRsTEingOWA/Pr6vw5sCQingKIiMc7HKOZmdUpI2EcAjxcM785ldU6AjhC0h2SfiBpbqMVSVooabWk1Vu3bm1TuGZmBtW96L0HMBN4G3A6cI2k/esrRcSyiBiIiIHJkyd3NkIzsz5TRsLYAkytmZ+SymptBoYi4vmI+DHwI4oEYmZmJSkjYawCZkqaIWkisAAYqqtzK0XrAkmTKE5RbexgjGZmVqfjCSMitgGLgJXAA8DNEbFO0iWS5qVqK4EnJa0Hbgf+KiKe7HSsZma2Q8dvqwWIiBXAirqyi2umA/hIepmZWQVU9aK3mZlVjBOGmZllccIwM7MsThhmZpbFCcPMzLI4YZiZWRYnDDMzy+KEYWZmWZwwLMvg4GDZIZhZyZwwzMwsixOGmZllccIwM7MsThhmZpbFCcN2y+DgoC+Em/UZJwwzM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaWxQnDzMyyOGGYmVmWUhKGpLmSNkgalrS4wfKzJW2VtDa9PlhGnGZmtsMend6gpAnAEuAkYDOwStJQRKyvq3pTRCzqdHxmZtZYGS2MOcBwRGyMiOeA5cD8EuIwM7MmlJEwDgEerpnfnMrqvVvSvZJukTS1M6GZmdmuVPWi9/8GpkfE64BvADc0qiRpoaTVklZv3bq1owGamfWbrIQhaYqkT0i6XdID6XW7pEt349v/FqD2PVNS2XYR8WREPJtmrwWObbSiiFgWEQMRMTB58uQmwzAzs2aMmTAkvRl4AHgPsA74YnqtGymT9KYmtrkKmClphqSJwAJgqG6bB9fMzkvbNzOzEuXcJfUZ4PMR8ReNFkr6bKrzhpwNRsQ2SYuAlcAE4LqIWCfpEmB1RAwBfyFpHrAN+Dlwds66zcysfXISxmuB942y/CpgYTMbjYgVwIq6sotrpi8ELmxmnWZm1l451zAeAUY75fSmVMfMzHpYTgvjcuBzkuZQ3LH0WCp/JUXnu7OB89sRnJmZVceYCSMilkp6EvhL4ByK6w4ALwBrgDMj4ub2hWhmZlWQNTRIRNwE3CRpT2BSKn4iIp5vW2RmZlYpTY0llRKEr1eYmfWhrIQhaQpwHvAHwO+n4keB/wd8LiIe3tV7zcysN4y3495pNN9xz8zMulDHO+6ZmVl3yumH8VqK51fsylXAUa0Jx8oyODhYdghmVnHuuGdmZlnccc/MzLK4456ZmWVxxz0zM8vijntmZpbFHffMzCyLO+6ZmVkWd9yzbCN9NZYuXVpyJGZWBnfcMzOzLO64Z+PiHuJm/cMd98zMLIs77pmZWRZ33DMzsyzuuGdmZllyLnqPStJUSdc1+Z65kjZIGpa0eJR675YUkgbGG6eZmY3PuBMGcCBwVm5lSRMobtM9BZgFnC5pVoN6+wEfBu5sQYxmZjZOY56SknTmGFWmNbnNOcBwRGxM618OzAfW19X7O+Ay4K+aXL+1yODgoDvpmdl2Odcwrgd+A8QuljfbSjkEqB17ajNwXG0FSbOBqRHxfyQ5YZiZVUDOwf5nFLfO7tfoxeid+pom6SXAlcBHM+oulLRa0uqtW7e2MgwzM6uTkzDWALNHWR6AmtjmFmBqzfyUVDZiP4qhRr4jaRNwPDDU6MJ3RCyLiIGIGJg8eXITIdh4uHe3WX/KSRiXA3eMsnwY+I9NbHMVMFPSDEkTgQXA0MjCiPhlREyKiOkRMR34ATAvIlY3sQ0rWU5SGRwcdPIx6yI5Pb3/7xjLfw38a+4GI2KbpEXASope49dFxDpJlwCrI2Jo9DWYmVkZWnFbbdMiYkVEHBERh0fEJ1LZxY2SRUS8za2L7uDWgllvKyVhmJlZ93HCMDOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLEt2wpD07lGW/XVrwrEydWKoDnfuM+tezbQw/qekayW9bKRA0hRJtwN/2frQzMysSppJGMdRjBy7VtKApPcC9wL/Dry+HcGZmVl15DxACYCIuDcNMb4U+D7FsOYXRMQ/tis4MzOrjmYver8eeCvFkObPAXPSs7fNzKzHNXPR+2Lgu8BXKBLHscCrgfskndCe8MzMrCqyT0kB5wKnRsTX0/wGSccDlwLfBPZqdXDWO3x3lFn3ayZhvC4inqgtiIhtwGJJK1oblpmZVU32Kan6ZFG37LutCcfMzKrKPb3NMviUmpkThnUhH7zNyuGEYaVzAjDrDk4YZmaWpZm7pKwP+du/mY1wC8PMzLKUkjAkzZW0QdKwpMUNlp8r6T5JayV9T9KsMuLsV50Y5tya47+HVcFuJwxJe+/m+yYAS4BTgFnA6Q0Swhcj4j9ExNHAp4ErdzdO25kPPGa2u8bTwvhefYGkV2e8bw4wHBEbI+I5YDkwv7ZCRPyqZnYfipFxrYe5VWNWfU1f9JZ0KkXLYF9JUyPi4ZrFNzH2szEOAWrfs5niWRv12/kQ8BFgIvCHu4hlIbAQYNq0abkfwTrMicCsN+zOXVL3A1OBScAXJE0DtgCPAM+3KrCIWAIskfRnwN8CZzWoswxYBjAwMOBWCDsOzkuXLm2qfjtj6aRmP7+Z5Ws6YUTEj4Glku4fGUNK0iHAoRTJZCxbKBLOiCmpbFeWA1c1G6eZmbXWePphfF/S2cBkYD3w1Yj4Xcb7VgEzJc2gSBQLgD+rrSBpZkQ8mGbfCTyImZmVajwXvZcDA8BvgT8G7pZ05FhvSkOiLwJWAg8AN0fEOkmXSJqXqi2StE7SWorrGC86HWVmZp01nhbGYRHx7pEZSUcD1wBvGeuNEbECWFFXdnHN9IfHEZeVyBe4zXrXeFoYT0t61chMRKwFDhh3RGZmVknjaWEsAr6Snra3HngN8JOWRGVmZpUznhbGg8BsYDXFHVIPAX/aiqDMzKx6xtPC+F5EHEvRWQ/Y3tP738YdlXWErzeYWTPK6OltZmZdaLw9vW+QdCht6Olt1iru/W3WGmMmDEkTIuKFkfkW9PQ2M7MulNPCeEbSvcCa9LobuG8kWQBExBZGH97DzMy6XE7C+ADF3VDHUgzh8XLgWUn3U5NEImJN26K0vuJTSGbVNGbCiIgbgRtH5iXNpEgex6Sf7wX2y1mXmZl1r6b7YaRBAW8D7gGeBvYCHm9xXGbZuuXhS90Qo9loshOGpJdLer+krwBbgU9R9Ow+meKhSGZm1sNy7pI6GzgNOBH4GfAvwCcj4s72hmZmZlWSc93hOoo7oM4HPp+ew21mZn0m55TU7cA+FE+9e1rS3ZKukXSupDdImtjeEM3MrArGTBgRcWJEHAi8Cng/8HWKTnqXAneSkkhbo7S+5IvEZtWSfStsRGwENgI3j5RJmk7x1L3ZLY/MzMwqZcwWhqQrJZ0g6UV1I2JTRNwSEX/TnvDM8vVCi6QXPoP1rpxrGHtTPL/7MUnXS/oTSXu3OS6zvtItfUk6qRW/E/9OWyvnGsZ5EXEI8E6Ku6UuBZ6QNCTpA5ImtztIMzMrX3bHvYi4KyIuioijKJ558a/A2cBmSd+TdEEatdas8jr5zdOtB+sVTQ9vDhARw8AVwBWphXEqMC8tvrzlUdpOBgcH+35gvt09ADfzPv+ezXaW08J4RtKdkpZKOkfSMZK2J5qI2BoR10XEn0REVrKQNFfSBknDkhY3WP4RSesl3SvpW+khTWZmVqKOD28uaQKwBDgJ2AyskjQUEetrqv0QGIiI30g6D/g0xai4Zl3NQ7dbNytjePM5wHDq14Gk5cB8YHvCiIjba+r/ADgjc91m2/m6gVlrNf0Mi4h4UNIjFKezjqD54c0PAR6umd8MHDdK/XOArzYbp1VX1Q7kvlZhlic7YUh6OUVL4DSKIc0fB76cpu9oR3CSzqDoSf7WXSxfCCwEmDZtWjtCsA7r9N1LZpavjOHNtwBTa+an0OB54JLeDlwEvDUinm20oohYBiwDGBgYiN2Mx8zMMuTcJXUdRb+L84HXRMTHxvksjFXATEkz0ki3C4Ch2gqSjgGuBuZFhJ/mZy3RqhaFWybWrzo+vHlEbAMWASuBB4CbI2KdpEskjfTl+HtgX+BLktZKGtrF6sxsDO44aK2Sc5fUiQCSDmPHyLSzgXcBBwLPS1oXEdkj1kbECmBFXdnFNdNvz12XmVnV9cqNFR7e3Cqtqv0WWhmXv/1bt2j6ttpaEbEJ2ATc0opgzKy6SdIse/BBM3sxtw6snzhh9DBf7DT//a2VxnVKysxGN9oB2wdz6zZuYfQBH5jMrBWcMMzw6TuzHE4YfcIHxP7Riudge1+xRnwNo0f53LmZtZpbGBXng7uZVYUThlkbjHVap91fBBptv9NDx/fil51e/EzNcMIwM7MsThg9oFe/zXUb/x16i/+eL+aEUSHeQa0s4933vN/2ByeMCvI/X3dpZ6Jv17rL2sf8pai7OWGYmVkWJwwz266Zb/85dWvruGXR/dxxz7qCDzZja/XvqFefy1H/9LsyPme3/m7dwjBrsW5Jbr6eYM1ywjDrUj7YW6f5lFQX6Nbmay+o0kG51/aDXvs8/cAtDDOrhColZ2vMLQwz6zlOPu1RSgtD0lxJGyQNS1rcYPlbJN0taZuk08qI0czMdtbxFoakCcAS4CRgM7BK0lBErK+p9lPgbOCCTsdXFf6GZLaDr3dUQxmnpOYAwxGxEUDScmA+sD1hRMSmtOx3JcTXVZxYrBN8wDYo55TUIcDDNfObU5mNk5OHmbVTV1/0lrQQWAgwbdq0kqOxXtHLibfRZ+vlz9sP6nuut1MZLYwtwNSa+SmprGkRsSwiBiJiYPLkyS0JzsysbFVN4mUkjFXATEkzJE0EFgBDJcRh1jOqeoCx3tLxhBER24BFwErgAeDmiFgn6RJJ8wAkvUHSZuA9wNWS1nU6znZp9Jxl/7O3Xi/8TnvhM1i+bvh7l9IPIyJWRMQREXF4RHwilV0cEUNpelVETImIfSLioIh4bRlxmllnlfkFqhsO2GXz0CBmZpalq++S6jd+GI1Zde1OX5Vu+z92wiiBO0HZePiLQ//q5C20jfiUlJmZZXELwyyTv813t0Yts935tp773la2BqrS4dItDDPrSlW4Rb3fbot3wqiIftrprHvtzn7abwfVVqni782npMysVFU7KNquOWGUyP8oZuXqxv/BMmN2wuiQbtwxzcrSzf8v3Rz7WJwwzKyn9fIBvNOcMNrIHfSsF5VxJ1I/q9Lnd8Jokyr9kS2f/27NqcrvqypxjKYXvkD6tlozG1U3HIytM5wwzMwsi09JmVnXqnrrp+rxNcstDDPrKVU/SFc9vtG4hWFmVsdP/WvMLYwWqOKYL2ad5v+D3ueE0QH+JzJrjar+L1U1rlZzwjAzsyy+htFi/fJNw6yf+P+64ITRQt6pzPL+D8bzv1KV/7OqxNFJpZySkjRX0gZJw5IWN1i+l6Sb0vI7JU0vIUwzM6vR8YQhaQKwBDgFmAWcLmlWXbVzgKci4lXAPwCXdTZKMzOrV8YpqTnAcERsBJC0HJgPrK+pMx/4L2n6FuCfJCkiol1B1T6wvbapWT9QWD82Q83MANTGY3DjDUqnAXMj4oNp/v3AcRGxqKbO/anO5jT/UKrzRN26FgIL0+yRwIbdCGkS8MSYtaql22LutnjBMXdCt8UL3RdzTryHRsTknJV19UXviFgGLBvPOiStjoiBFoXUEd0Wc7fFC465E7otXui+mFsdbxkXvbcAU2vmp6SyhnUk7QG8AniyI9GZmVlDZSSMVcBMSTMkTQQWAEN1dYaAs9L0acC323n9wszMxtbxU1IRsU3SImAlMAG4LiLWSboEWB0RQ8B/B/6HpGHg5xRJpV3GdUqrJN0Wc7fFC465E7otXui+mFsab8cvepuZWXfyWFJmZpbFCcPMzLL0dcIYa4iSDsZxnaTHU/+TkbIDJX1D0oPp5wGpXJL+McV8r6TZNe85K9V/UNJZjbbVwpinSrpd0npJ6yR9uMpxS3qppLsk3ZPi/a+pfEYafmY4DUczMZXvcngaSRem8g2S/qgd8dbFPkHSDyXdVvWYJW2SdJ+ktZJWp7JK7hM129pf0i2S/k3SA5LeWOWYJR2Zfr8jr19JOr8jMUdEX74oLrg/BBwGTATuAWaVFMtbgNnA/TVlnwYWp+nFwGVp+h3AVwEBxwN3pvIDgY3p5wFp+oA2xnwwMDtN7wf8iGKol0rGnba7b5reE7gzxXEzsCCVfw44L00PAp9L0wuAm9L0rLSv7AXMSPvQhDbvHx8BvgjcluYrGzOwCZhUV1bJfaImvhuAD6bpicD+VY+5JvYJwKPAoZ2Iua0fpsov4I3Aypr5C4ELS4xnOjsnjA3AwWn6YGBDmr4aOL2+HnA6cHVN+U71OhD/V4CTuiFu4GXA3cBxFL1g96jfJyju4ntjmt4j1VP9flJbr02xTgG+BfwhcFuKobIx0zhhVHafoOjj9WPSDUDdEHNdnCcDd3Qq5n4+JXUI8HDN/OZUVhWvjIhH0vSjwCvT9K7iLu3zpFMfx1B8a69s3OnUzlrgceAbFN+0fxER2xpse3tcafkvgYM6GW/yGeBjwO/S/EEVjzmAr0tao2LoHqjwPkHR4toKfD6d9rtW0j4Vj7nWAuDGNN32mPs5YXSNKNJ/Je9/lrQv8C/A+RHxq9plVYs7Il6IiKMpvrXPAV5dbkSjk/THwOMRsabsWJrw5oiYTTEa9YckvaV2YdX2CYqW2Gzgqog4Bvg1xemc7SoYMwDp2tU84Ev1y9oVcz8njJwhSsr0mKSDAdLPx1P5ruLu+OeRtCdFsvjniPhf3RJ3RPwCuJ3idM7+Koafqd/2roan6WS8bwLmSdoELKc4LfXZKsccEVvSz8eBL1Mk5irvE5uBzRFxZ5q/hSKBVDnmEacAd0fEY2m+7TH3c8LIGaKkTLXDo5xFcY1gpPzMdOfD8cAvUzN0JXCypAPS3REnp7K2kCSKHvkPRMSVVY9b0mRJ+6fpvSmutzxAkThO20W8jYanGQIWpDuSZgAzgbtaHS9ARFwYEVMiYjrF/vntiHhfVWOWtI+k/UamKf6W91PRfQIgIh4FHpZ0ZCo6keJRC5WNucbp7DgdNRJbe2Nu90WZKr8o7h74EcW57ItKjONG4BHgeYpvPOdQnHv+FvAg8E3gwFRXFA+gegi4DxioWc8HgOH0+k9tjvnNFE3ee4G16fWOqsYNvA74YYr3fuDiVH4YxcFzmKJpv1cqf2maH07LD6tZ10Xpc2wATunQPvI2dtwlVcmYU1z3pNe6kf+pqu4TNds6Glid9o1bKe4YqnrM+1C0Hl9RU9b2mD00iJmZZennU1JmZtYEJwwzM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZn1O0pclPSXplrJjsWpzwjCzzwJnlh2EVZ8ThllJ0pAMj0k6vMw4IuI7wNONlkn6kqSPdjYiqyonDKs8SddLigavH5Qd2zj9DbAiIh4qO5BRXAJcJOkVZQdi5dtj7CpmlfBN4P11Zc+1c4OSJkZEW7Yh6WXAB4FT27H+um2tpfH/+skR8bPR3hsR90naCJxBMR6R9TG3MKxbPBsRj9a9fg4g6TuSlkr6pKQnVDwf/XJJ2/fvNFLnxyQ9JOm3Kp47fUbtBtJ6rkrv3Qrckcr3kfQFSc+kU0gXSrottXzOlPSkpL3q1vXPkkYb/fgdFIM33lH3Pkn6qIpnLD8rabOkT9XFd4Wkn0vaKunDaRTaJZJ+IemnknZKrBFxdEQc1eA1arKoMUQxMqr1OScM6xXvA7YBfwAsAs4H3luz/FKKUYA/RPGM608BV0t6Z916zqAY3fMEdlwIvgJ4K/AuimdSvD4th2J02JcA80dWkE7fvIti+PddOQFYEy8e/fOTwMdTfK8F3sPOT0V7H8X1huOA/0bxRL5bKUZdHqB4PvW1Ss9FaJG7gDlpWHjrZ+0cgtcvv1rxAq6nSAbP1L1GHnL/HeD7de/5BnBtmt4H+C1wQl2dz1BcQ6BmPffW1dmX4tTXgpqyfYCngOvT/D8BX6tZfh7FIzL3GOUz3Qrc0GBb/w6cu4v37PQ5KRLbVmCopmzPFO9pTfx+v5nW8xuK4fXfWLf8dRStocPL3hf8KvflaxjWLb4LLKwr+0XN9L11y34G/F6ankXxrIivSar9Rr8nsKnuffWPQz081dv+wKGI+LWk+2vqXAPcLWlKRGymeMbADbHjuduN7A08Vlc2C9iL4pkGu7L9c0ZESHqc4hkHI2XPS3qKHZ99TBHx9jGq/LYmZutjThjWLX4TEcOjLH++bj7Yccp15OepwE/HeN+vmw0sIu6RdDdwtqRbKU4NnTH6u3iC4kE9zWr0OUf77K1wYPq5tYXrtC7khGH9YD3wLHBoRHy7yfc+RHFAfgOwEbbf4XRUWjbiGuBjwCTgjojYMMZ6fwicXVf2QIrzRIqnplXFUcCW2PHsaOtTThjWLfaS9Pt1ZS9ExJjfeiPiaUmXA5dLEsXprX2B44HfRcSyUd77jKTrgMskPUHxKN2/pfgGX3t660bgSorrF+dmfJ6VaZ0HRcSTNXF+FviUpGdTnAcBx0bEVRnrbJcTaP/zqa0LOGFYt3g7xcG61hZgSub7P05xzeAC4CrgVxTPIf90xnsvoLjQPURxsf0fgFdSXKAGth/sbwZOA24ea4VR9G+4C1jAzv0bLqS4oP5xis/2GPCFjBjbQtJLKe74+qOyYrDq8DO9zZqU+lz8BPj7iLiipvyrwOaI+PPM9cylGMdpVkS80JZgx0nSh4D5EXFy2bFY+dzCMBuDpGOA11DcKbUf8Nfp501p+QEUp21OpuijkSUiviZpCUVL4ictDrtVngf+c9lBWDW4hWE2hpQwrgGOpOgPsha4ICLWpOWbKO4k+kREXFZSmGZt54RhZmZZPDSImZllccIwM7MsThhmZpbFCcPMzLI4YZiZWRYnDDMzy+KEYWZmWZwwzMwsixOGmZllccIwM7Ms/x/Tw13zQRxXggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute the histogram of the energies\n",
    "bins = np.linspace(np.amin(y_sampled), np.amax(y_sampled), 200)    \n",
    "hists, edges = np.histogram(y_sampled, bins=bins, density=True)\n",
    "\n",
    "#plot distribution of pool\n",
    "plt.bar(edges[:-1], 1000*hists, width=np.diff(edges), color='dimgray')\n",
    "plt.xlabel(\"Energy (cm$^{-1})$\", fontsize=14)\n",
    "plt.ylabel(\"$N/N_{tot}$  x  1000\", fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "18cf641c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '|E_p-E_test|')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzxUlEQVR4nO2dfXQc5Xnof4/Wa5AdiuzgUBA2Ni7HnFAnNlGxXff2AGniAAFU8kEc09A0hXPb3jSEXDd2wgmQ0otz3SQkOfekhSQtKQ5x+IgCwb0OF0x76sRuZGTjmOBgMBiLLxVbDsECy9Jz/9gZeTSa2Z2ZndmZlZ7fOTrafXd29tnd2fd53+dTVBXDMAzDiENL3gIYhmEYzYcpD8MwDCM2pjwMwzCM2JjyMAzDMGJjysMwDMOIzaS8BciKk046SWfPnp23GIZhGE3Ftm3b/ktVZ9Q6btwqj9mzZ9Pd3Z23GIZhGE2FiDwX5TgzWxmGYRixMeVhGIZhxMaUh2EYhhEbUx6GYRhGbEx5GIZhGLEZt9FWRm26enpZu3E3L/QPcGpbKyuXzaNzYXveYhmG0QSY8pigdPX0svq+nQwMDgHQ2z/A6vt2ApgCMQyjJma2mqCs3bh7RHG4DAwOsXbj7pwkMgyjmTDlMUF5oX8g1rhhGIYXUx4TlFPbWmONG4ZheDHlMUFZuWwereXSqLHWcomVy+blJJFhGM2EOcwnKK5T3KKtDMNIgimPCUznwnZTFoZhJMLMVoZhGEZsTHkYhmEYsTHlYRiGYcTGlIdhGIYRG1MehmEYRmxMeRiGYRixyVR5iMh3ROQVEfmFZ2y6iDwkIk85/6c54yIiXxeRPSLyuIic43nOVc7xT4nIVVnKbBiGYdQm653HPwPv842tAh5W1TOBh537ABcCZzp/1wDfhIqyAW4AFgHnAje4CscwDMPIh0yVh6r+O3DAN3wZcIdz+w6g0zP+Xa2wBWgTkVOAZcBDqnpAVQ8CDzFWIRmGYRgNJA+fx8mq+qJz+yXgZOd2O/C857j9zljY+BhE5BoR6RaR7r6+vnSlNgzDMEbI1WGuqgpoiue7TVU7VLVjxowZaZ3WMAzD8JGH8njZMUfh/H/FGe8FZnqOO80ZCxs3DMMwciIP5XE/4EZMXQX8yDP+MSfqajFwyDFvbQTeKyLTHEf5e50xwzAMIycyraorIncB5wEnich+KlFTa4AfiMgngOeADzuHbwAuAvYAh4GPA6jqARH5W+DnznFfVFW/E94wDMNoIFJxO4w/Ojo6tLu7O28xDMMwmgoR2aaqHbWOs34eRup09fRakynDGOeY8jBSpaunl9X37WRgcAiA3v4BVt+3E8AUiGGMI6y2lZEqazfuHlEcLgODQ6zduDsniQzDyAJTHkaqvNA/EGvcMIzmxJSHkSqntrXGGjcMozkx5WGkyspl82gtl0aNtZZLrFw2LyeJDMPIAnOYG6niOsUt2sowxjemPIzU6VzYbsrCMMY5ZrYyDMMwYmPKwzAMw4iNKQ/DMAwjNubzMAxjwmMldeJjysOoyfVdO7lr6/MMqVISYfmimdzcOT9vsQwjFaykTjLMbGVU5fqundy5ZR9DTvXlIVXu3LKP67t25iyZYaSDldRJhikPoyp3bX0+1rhhNBtWUicZpjyMqgyF9HsJGzeMZsNK6iTDlIdRlZJIrHHDaDaspE4yTHkYVVm+aGasccNoNjoXtnPL5fNpb2tFgPa2Vm65fL45y2tg0VZGVdyoKou2am4sYq46VlInPtbD3DDGOW7EnJ8rF88yBWKMIWoPczNbGcY4xyLmjCww5WEY4xyLmDOywJSHYYxzLGLOyAJTHoYxzrGIOSMLLNrKMMY5FjFnZIFFWxmGYRgjWLSVYRiGkRmmPAzDMIzY5KY8ROTTIrJLRH4hIneJyPEiMkdEtorIHhFZLyKTnWOPc+7vcR6fnZfchmEYRk7KQ0Tagb8GOlT1d4ES8BHgS8BXVfV3gIPAJ5ynfAI46Ix/1TnOMAzDyIk8zVaTgFYRmQRMAV4ELgDucR6/A+h0bl/m3Md5/N0iFqRuGIaRF7koD1XtBf4e2EdFaRwCtgH9qnrUOWw/4FYqaweed5571Dn+rf7zisg1ItItIt19fX3ZvgnDMIwJTF5mq2lUdhNzgFOBqcD76j2vqt6mqh2q2jFjxox6T2cYhmGEkFeS4B8Be1W1D0BE7gOWAm0iMsnZXZwG9DrH9wIzgf2OmetE4NXGi20YhpEfXT29rN24mxf6Bzi1rZWVy+blVko+L5/HPmCxiExxfBfvBp4ANgEfdI65CviRc/t+5z7O44/oeM1uNAzDCKCrp5fV9+2kt38ABXr7B1h93066enprPjcL8vJ5bKXi+H4M2OnIcRvwWeA6EdlDxafxbecp3wbe6oxfB6xquNCGYRg5snbjbgYGh0aNDQwOsXbj7lzkya22lareANzgG34GODfg2DeADzVCLsMwjCLyQv9ArPGssQxzwzCMJuDUttZY41ljysMwDKMJWLlsHq3l0qix1nKJlcvm5SKPlWQ3DMNoAtyoqqJEW5nyMAzDaBI6F7bnpiz8mNnKMAzDiI0pD8MwDCM2pjwMwzCM2JjyMAzDMGJjysMwDMOIjSkPwzAMIzamPAzDMIzYmPIwDMMwYmPKwzAMw4iNKQ/DMAwjNlaexDCMEYrUqc4oNqY8DMMAjnWqcxsOuZ3qAFMgxhjMbGUYBlC8TnVGsTHlYRgGULxOdUaxqWm2EpG9gALi/B/1sPNfgVtV9evpimcYRqM4ta2V3gBFkVenOqPY1FQeqjqnEYIYhpEvK5fNG+XzgHw71RnFJrLZSkS+FGXMMIzmpHNhO7dcPp/2tlYEaG9r5ZbL56fqLO/q6WXpmkeYs+pBlq55hK6e3tTObTQWUfVbokIOFHlMVc/xjT2uqu/IRLI66ejo0O7u7rzFMAzDwR/NBZWdTdoKyqgPEdmmqh21jqu58xCRvxCRncA8EXnc87cXeDwNYQ3DGP9YNNf4Ikqex/eAfwVuAVZ5xl9T1QOZSGUYxrjDornGFzV3Hqp6SFWfBa4HXlLV54A5wJUi0pateIZhjBfCorYsmqs5iZPncS8wJCK/A9wGzKSyKzEMw6jJymXzaC2XRo1ZNFfzEqc8ybCqHhWRy4FvqOo3RKQnK8GMaFzftZO7tj7PkColEZYvmsnNnfPzFsswxuA6xa121vggjvIYFJHlwMeAS5yxcvoiGVG5vmsnd27ZN3J/SHXkvikQo4h0Lmw3ZTFOiGO2+jiwBPg7Vd0rInOAf0n6wiLSJiL3iMiTIvJLEVkiItNF5CERecr5P805VkTk6yKyx4n0OqfW+ScCd219Pta4YRhGWkTeeajqEyLyWWCWc38vUE+S4NeA/6uqHxSRycAU4HPAw6q6RkRWUYnu+ixwIXCm87cI+Kbzf0IzFJKjEzaelBW3/4zNTx8LrFs6dzrrrl6S6msYhtFcRFYeInIJ8PfAZGCOiCwAvqiql8Z9URE5EfhD4E8BVPUIcERELgPOcw67A3iUivK4DPiuVjIatzi7llNU9cW4r10E0piMr+/aGfpYSST0sbj4ZQXY/PQBVtz+s1gyW58IwxhfxDFb3QicC/QDqOp24IyErzsH6AP+SUR6RORbIjIVONmjEF4CTnZutwNeW8x+Z2wUInKNiHSLSHdfX19C0bKl2mQcFb+vw8/yRTMTy+fHL2ut8SDczOLe/gGUY30irDSFYTQvsRzmqnpIRq9qh+t43XOAT6rqVhH5GqMTEFFVFZFY9hdVvY1KGDEdHR3p2m4S4l9xB1UthXiTcTWfxpWLZxXOWV4ts9h2H4bRnMTZeewSkY8CJRE5U0S+Afw04evuB/ar6lbn/j1UlMnLInIKgPP/FefxXip5JS6nOWOFJmjFnQbVfBpFUxxgmcWGMR6Jozw+CZwNvEklOfAQ8KkkL6qqLwHPi4ibHfRu4AngfuAqZ+wq4EfO7fuBjzlRV4uBQ83g7whacUfl+q6dzF29gdmrHmTu6g2jfBxhPo00fR0uS+dOjzUehGUWG8b4I47yuFhVP6+qv+f8XQ/EdpZ7+CSwTkQeBxYA/wtYA7xHRJ4C/si5D7ABeAbYA9wO/GUdr9sw4qysvZOx69Nwdxhu/oarQMJ8Gmn6OlzWXb1kjKKI6+C3zOLxg5VUN1zqLck+ZqwoFKEk+9I1jwSaqsotMOjxFvkn47mrNwSapkoiPH3LRUDzZZZbtFXzYyXVJwZRS7JHaUN7IXAR0C4i3jazvwUcTS7i+CesM1utH1uU/I2bO+fnoiySKgHLLM6WRuTiWOCD4SWK2eoFoBt4A9jm+bsfWJadaM1P0s5sjfRpxMFCbotJGuHfUbDAB8NLlB7mO4AdIvI9VR0MO05E7lXVD6Qq3TggyYp7+aKZgXkcWfg04lDPyrNZzFbNIqeXNHJxohAWam6BDxOTOOVJQhWHQ9KEQcOHa44qmk8j6crTbyt3dyxA7hOzV1m0TSnzmzeOMjhcMQ8WSc4iEGaGtcCHiUmcJMFaFCIpb7yQl0+jGklXnkW1lfuV2sHDY9dHRZCzKFhJdcNLmsrDGOckXXkW1VYeNQ8nbzlrsXTu9EATVZxcnKhY4IPhEifPoxb5enONzEkaAFDUJMGoSiFvOWuRRi6OYcQlSqjub6nqr0Mem6Wqrmf3s6lKZhSSJCvPotrKq9UacymCnFEwRWE0mihmq0ep1J1CRB5W1Xd7HutyH1PVn6QtnBEd1/Hb2z9ASYQhVdoLYpMuqq08SKmVS8LUyZM4NDBYGDkNo4hEUR5ec5TfiGqmqhzxKgzhWMSCm0xYpGihotrKj5vUMqI8pk0pc8MlZxdSTsMoGlF8HhpyO+i+0SC8CXsQ/kW40ULGaNzPr3/gWITVG4NJOwwYxsQjys7jbSJyHZVdhnsb5/6MzCQzqhKnYm/Ro4XyoKjhw0a+NGOSaF5EUR63AycE3Ab4VuoSGTXp6umN1Ruk6NFCeVDU8GEjP4qczFpEopQnuakRghjRcC/wOJx/lm0Q/VipDcOP7UbjkSjPQ0QeS1sQIxpJGkxterKY/dyjkkUPCesxYvix3Wg8kiYJWpRVTiS5kP3PaaaGPkGVfK9dv72uirGuXXtgcGikUnHUhEdj/FLUZNaiklR5PJiqFBOcOJN52IXc3tZKe4SLv9nKqofttDY/fWBUa96o+KPUhlRHdhymOCY2ae9Gm2mRloRYykNEfltELgW2ishvZyTThCLuZH7+WTPGbPtayyXOP2sGh4+M7c3lv/ir2XXTJK0fTrWd1l1bn499vka9f6P5SFp+J4hmW6QlIXJhRBH5c+ALwCNUzFbfEJEvqup3shJuIhDFSReWDAiVL+KcWSdy77beMedpay1z46Wjk94aYddNM2qlWgmRIVWWrnkk1q7B7NpGNdJKZp0Izvc4VXVXAgtV9VUAEXkr8FPAlEcd1JrM/BNxUJbmlmcOBraunXrcpDEXaiOijNL84axcNo9r128PfTyuYrIoq+Rc37Uz1x4zcXMw0szZiPveJ8IiJY7Z6lXgNc/915wxIwLXd+1k7uoNzF71IHNXbxix19dy0kWJrgrreR50oTYiyijNH07nwvaapcWjmp26enp5/c3apj1jLNd37eTOLftGrrUhVe7csi+R3ykJcc1AaZqNkrz3ieB8j6M89lDxddwoIjcAW4Bfich1nqxzI4Cwi2/F7T+rOZnXs1IJulDTtOvGed1q47VYd/USrlw8q2oP96jdDL3lSKBSz8qirGoT5l9K4ndKQlxfVZq+rSTvfSKEgscxWz3t/Ln8yPl/QsCxExr/dvnFQ8ET2+anD/DTZw7w+2dM59lXBwK311HKhgchMOZC9cv11SsWZDJpZlGC3e2suHTNI6l1MwSYMnmsac8YS9juNmw8beLuZtPc/SZ570WtJJ0mcXqYV800F5FvqOon6xepuQlyFldDtaJErlw8K9CGGjQRR0EZ7QPo6ull5T07GBw6VnH32vXb6X7uwMjrdvX08rn7HuewUyBQBFYsmkXH6dNj/Qiy/OEEfR7ivJ9qzvOJYIPOErfMf9B4I4jrq0rTt5X0vRe1knRapNmGdmmK52pakmSAQ2ULHKQ8vBNxnB1ISYQ5qx4cmbhvemDXiOLwcueWfXScXvEpXPeD7Qx7DlGtPH7Xfz7P0HC8Mu9Z/XD8n4c3+qyabOYor4/li2Zy55Z9geONIO5uNs3db97vvaik2YbWIPlKttYWOMiGWut8XkfhwcODoceu3bibtRt3j1Ico87leyDvvIjOhe1sXnUB7W2tY6LPwmSbCDboLLm5c/4ov1NJJHS3nAVxfXVp+vbyfu9FRTQlm6WIPKaq56RyshTo6OjQ7u7uhr9umE2+3ALV2kWURHj6lotin9ePPw8kKnGfJ8DeNRePGmt0KOecVQ8GyhwkG1i5bcOIgohsU9WOWselufOIbfwUkZKI9IjIj537c0Rkq4jsEZH1IjLZGT/Oub/HeXx2inKnStgKd+2HFnDrFQuYUg7+yGttgaPuaJIojpJIbPONwqjM8TxCOeNGdXU/d4CXDr2BAi8deoPu5w5kJpthjHfqUh4i4vWZfC3BKT4F/NJz/0vAV1X1d4CDwCec8U8AB53xrzrHFZJq2+XOhe088bcX1twCB5X2yNI2P6TKymXzaImp/r2x89/bOtYmDNFDOZOUM4ljiso7T8Ewxhs1zVYi8h+q+gfO7X9R1T/xPJbYVCUipwF3AH8HXAdcAvQBv62qR0VkCXCjqi4TkY3O7Z85CuslYIZWET4vs1W9dPX0svLuHQx6/AzlFmH2SVN46pXXM3nN9rZWNq+6gK6e3qrZ3GG0tZbH5E94eTbAhOTFH6EGFSXgKt1q5qaopqi5qzeERsxUMxcmwcxj44+J9J1GNVtFibaa6rl9tv91Ykk1mluBv+FYnshbgX5VdVOA9wPut9MOPA/gKJZDzvH/NUoYkWuAawBmzZpVh2iNxXthwljT0+CwZqY44FizqM6F7XzmBztCJ9lhxwnvp5riiBLKWSuhK6xOlvvcKD/oOLH69UwU1o1u/GHfaTBRlEe1rUkib7uIvB94RVW3ich5Sc4RKIzqbcBtUNl5pHXeLAnaadSDACfW2An4uXdbLx2nT6f7uQMMh0yyyxfNZNOTfbETFqOEM1bLwQhTLJ//4U6GlUg/6GomML9yq3eiKGpBvIm0ck6bon6neRPF59EmIn8sIh9wbl/u/H0AODHh6y4FLhWRZ4HvAxdQ8Zm0efwopwHur74XmAkjfpYTGSd1tW68f1dqigMq2vzGS/0bxOoMDA7xufse584t+wJXA65PJszHMG1KOfC8reWWmtFWXT29tITsTk5taw1VLK8fGQr8QX96/fZRyqJW216/cqu3rEURkxEnQnnwLCnid1oEoiiPfwMuBd7v3L7E+Xs/8O9JXlRVV6vqaao6G/gI8IiqrgA2AR90DruKYyVQ7nfu4zz+SDV/RzMRZ4cQhfa21kSrocMhccQlkREF4AYDtHoixt48OsTbTzkhUKnccvk7qr6mO6kFmY5cx3eSKLCVd+8YmRirJW0GxerXO1EUsSCe9TCpjyJ+p0UgitnqhmoPiojrXOhX1V/XKc9nge+LyM1AD/BtZ/zbwL+IyB7gABWFM2EIy8FoAbxTfhZJb/6Jvfu5Awx4FM2wU15l6dzw+lxhhE3sJZFRCV1xnfiDwzpiUgib9AUCd0X1ZqJnUderXmzlXB9F/E6LQBTlcQfVfRvu3PbPwHfjCqCqjwKPOrefAc4NOOYN4ENxz503UezM06aUA7O/W6RSHsR9XvdzB8Yk4IXVnErbHOGtGRUWervlmYOxo5bCJq9h1ZHPqXNhOzfevyv2Ds09d1xlUO9EUcSCeFaapT6K+J0WgZrKQ1XPb4Qg442ojteL33FKYN2cjy4abVLpXNhetfaV/3XTxCt7mtVVo05qN156duzikO45zj9rBut8vpxqhRSjTBS1FgVFK4hnK+f6Kdp3WgTSLIxoeIgSoeEmrgWx6cm+mq8RNIklLcxYC1f2NKurxpnUjpvUEvl9lVuElcvm0dXTy73begO7L0K4Qq82UTRj2GZRVs4W8TW+MOWREbXszNUUR9jz3/OVR0PzPdxJLAvF4X2NKxfPSq3CaNikBhVT2Qv9A7RNKfObN46OSZocUg0t5DipJCPnrfV5xA25bNawzTxWzl5l4f8em0HpGtUx5ZERtUwytcp2+E031RSHy8DgUOjOIA28kVe1CiAGrTIhePXrzSL3O8eD/EGDw1o1O3VgcDiWIo2Tu2LO52j4d2hB32MzKF0jHFMeGVHLJFNrgj/4+puj+nFEzTAfUqW1XMpkB+LK7Hb1CyPItLPy7h0gjGpG5fXNxN011VKPcRVpV09vpEnMnM/RiGo+zVvpmiktOdbPIyPcnIi21mMJdMd78iNq+QgODw6PSuiKSkmED7wru4s/StHCoIljcFjHNKNyV55Z+WmGVClHrPYYNefB+oJEo968mEZgyZP1YTuPjHnz6LGciIOHB0cUQVh3siDiTKxDqty7rTdxX49auO1rb3pgFzdcUslk96/c4qwm3W6AWdDWWub1I0drH0j0ya4ozudqFGE1HbZD85K30g3zX930wK7cP79mILVmUEWjCFV1wxo4uVVsaznNi06LMMZp3Voucdyklsh5GSURfqt1UtVOh0kJy6EJwv1Omp1aFYrzlKNcEqZOnsShgcFCTMphzcT85PH55UmaVXWNhNRyrnacPp0HH38x8cR55tumsueV1zPZYUQhKNppYHCI48stkf0uQ6qZKA6A/ojndUN7xwNFiQZrhh1alN0RFMOxX4TdpB/zeWRItZo4XT29rLxnx5iJs7XcQrlU25DT3tbKQ9edV0hHbf/hwTENscKKJ2ZJ2GfjdYO0tZZZ+6F35v5DTItmiQZL0vwrbYL8V2FYYcux2M4jQ6pFXN30wK4xDmSA48slbrjk7Jr1nHr7ByoKaNm8RA2csuRUpzijPys76zwUP2Gff6NMEHmsFosSDVYtmRLCe7TkvTt6/c2jgSbXoha2zHPRYzuPDKnWkjbMVHPw8CA3PbAr0vnTLkOSFm5zKS/uZ5EkEz0J06aUq37+WZPXarEo0WDVJrwiVfntXNjO5lUXsHfNxWxedQE3Xnp2IT4/L0XdTdrOI2OSZPZG9QEMDA5x4/3RFE0j+fGOF0PrcDVql+RGguVVkyiv1WLevgZ3txXmS6g24eU9GUL+n18QRdlN+jHl0eSk3Q8kDfoHBgOT7rp6sgsh9tP93IFYfc7rbT0bNVy5ERNkXgozimnSnfCKOBm6FK0IYlELW5rZKici5q6N4G3A1AwEmSDWbtzdsMgwt/xLFPNRPSamsOee2BocIJD1BJmnI7pWsqcb1VYU01qzkKf5tRq282gAQSvTjy4KLjAYxkBIp7+iErTCjtv/vB7csiRh5qPP3ff4yHfSElDGJKqJKez8QeHKWU+QeVf8rbmrchZMRTQNQTHDYV2KthsCUx6ZE/aDvuXyik+gmZMEq6HA3NUbRoomNjqssCTC9V07QxXW4cFhDjuPhdW/imJiCjum//AgX71iQUMno7yjcmrlTQwOHevwWITJ0KssTnSqEQTVXstbzqJiyiMmcVcn1X7Qm1ddwKYn+xq6Im8kQ6rcuWUfe/t+w5ZnDjb0tU96S7luxRzFxFTNmdnoCTLvqJwg23xestTCv6gL8h0WIRy2yJjyiEESs0CYYnDHi/JjypLNTx9o+Gu+/NqRus8RFHIMlV4sbkl6AUotwpAn3T4v+33eUTlec1TYdV8Up3izVP0tMs3lhc2ZJPHpYXkN7nitH1Ncx7qRHkHdHN16ZK6pS4GhYWXq5FLuzswiOKLdvIlbr1iQuyzVaIaqv0EUITPfxXYeMUhiFqjV83vlsnmsvHvHmE55bsmMrp5ePr1+e271qyYyQd9rWBOvNwaH2bvm4qxFqkqRHNFFkiWIZqj66yfvgAg/pjxikMQs0B7ynHbvc/y7C8/9u7v3meLIiaDvtdpiYPaqB0M7KzaKIjiiXYoki58g/0y5RXjL8ZPoP1yMqr9+8g6I8GNmqxgkMQvUes7ajbvH1Lhyo1IgH39B2pRysr0tnTs98XOFYJ9HrfIqbpDA9V3FLB1jVAjKnVj7oXfS84X3jpQqKZLigPwDIvyY8ohBkmSdWs8p2gWRBcNBtdszRoB1Vy+peoz3O1k6d/qoDaAC927rHWNTXr5oZqTXr9Wj3sgff12roikLP9WqdOeBma1ikmQrXu051UxheZdcTos8zG4rFs8CwhtCTZtSHtX8aemaR8bIGWQScM1RbrRVGFF7pxtGVIpWpsR2HjlTzaxVxKKHSWhQId1RuJP8DZecPaY/SrkkI4UTXeLsAG/unM/Tt1zEs2surhlNVy9Fiq4pMhPhcypamRLbeeRMtaiUovXpSEJruYSgHG5geRWhspNwP88rfm8mm57sqxr1kzRHYvEZ0wL9UovPmFbXe4B0o2uKXHqjXooWhZQlRQpCMOVRAIp0QaRJuzNJfToHJegqgt7+gZFM86Vzp4f6QaZMHrsJj2ISePbV4B1Lz75+5q7ewJBq4gistKJrxvvkWrQopIlCLmYrEZkpIptE5AkR2SUin3LGp4vIQyLylPN/mjMuIvJ1EdkjIo+LyDl5yN1o8mjdmiYv9A9UGls10Gw1uSShPpbNTx9gxe0/GzO+4vaf8dQrr48ZP23a8TUnnzBz1+HB4RG/R9IIrLSCKYrUfCkJtUxSEyHopIjk5fM4CnxGVd8OLAb+SkTeDqwCHlbVM4GHnfsAFwJnOn/XAN9svMiN54ZLzs4tzDUNlEpjq0b4jpfOnc6zay7mSEBrXy9BJqawcOggheInTqRL3AistKJrwibR3v4B5q7ekGpYcdq+hyjl8osWhTRRyEV5qOqLqvqYc/s14JdAO3AZcIdz2B1Ap3P7MuC7WmEL0CYipzRW6sbTubCdL3/onU2/A8mScotw6xULaoblZkVQwEMYcSOw0io3Um0STTMvJYvWu1F2TUUoyzIRyT3aSkRmAwuBrcDJqvqi89BLwMnO7XbAu2zb74z5z3WNiHSLSHdf39i6RM1I58J2er7w3rzFKCQCXHHuzFGmpbQVba1JNSgCJizQKm4EVtTomlqr/SgKLo28lCzMY1FMUkWLQpoo5OowF5G3APcC16rqr8Xz41JVFZFYSzVVvQ24DaCjo6PpA+29ETLGWCqJfPvpOH36qL4M/iq3XoKyzpfOnR5qunKd7dWc3f6AB7d4op+oCYbVzu0nijM8SrXbNPJSsvA9RI2CG49BJ97qzXmXvQkit52HiJSpKI51qnqfM/yya45y/r/ijPcC3l/eac7YuKWrp5eV9+wYMQEYwQwMDnPt+u0jn1P/wCAtBO9AwqKt1l29pGopkzu37At0tIdxc+d8rlw8a2SnURLhysWzMvnhR13tu9nUWealhJnHlNo7uDCSmqSi+F6KnBvir95cxLI3ojlkwkpli3EHcEBVr/WMrwVeVdU1IrIKmK6qfyMiFwP/A7gIWAR8XVXPrfYaHR0d2t3dndl7yJqFX/xJYGa0EY32ttZRGeQu1fIdZq96MNK5i7QKnLPqwcDFhUBgld+wXVEays2/C0rrNeLmqATJ0VoujTJlRTkmT9wwbz8lEZ6+5aJMX1tEtqlqR63j8jJbLQX+BNgpItudsc8Ba4AfiMgngOeADzuPbaCiOPYAh4GPN1TaHDDFUR9BppKunl6uW78dN12xt3+A65wclM6F7ZQCepkH4a4Cobo5K+qkV08CX9zkRn95lTQVoStzWHLrXVufT/Q61UxSQZ9dlLyPoueG1GrlUARyUR6q+h+ER/+/O+B4Bf4qU6GMcUXQ5Ln6vsfx57kPO+OdC9tZvmhmrNa11SbDMF9E93MHRrLd26aUeXNwaFT2fdwEviT1jm7unJ/ZrqlaZYS0J76wzzhs5+NdUBQ9NyRsIZNW2Zs0yD3aygimrTU4aqi13EKpONdPIQmbPAdCSqS4466vIirVJsOwle26LftG/DMHDw8Glm2JE6FUxEijrOt9uYR9xmGv411QFD03JCy4IknQRVZYeZIUSbN+0I2Xnj2mw2AL4ROgUWHalDKq8On127nx/l2IMNLcpxreWlhnvm1qpATBapNhWFRT1LV3nBVw0SKNwnZwaU98YZ/RkCqt5VLV3VjRKtT6ydK8mBamPFIi7fpBQQUTDx85ar6QKly5eBb3busd+Q76B459VrVajnprYbWWS5EUSJarwKgr4CIWPGzUxBfm72n3+D7CPpeit8mFbM2LaZBLtFUjaHS01dI1j4ReyEFRP7UImhTGQ5XdrHDNfF6FUe/5ph43aeTzn/3WVrY8czDyZBg1ciuIqFE/RYgYylN5FeH9j0eKHm017kjTARe0i1l5z4665BvPtJZLzDhhciRTU1T6BwZHFFFv/wAHXj/Clz/8zswnpbbWMjdeenak6KKWAKdqIyOG3Fwkt42y9zqNovjqVTpxdw9FT7prNkx5pETSfhAw+qJ2rej+/aC/z7lxjIHBoVQVR9hr3PTArsgTXLWs9SDaE+QvhDnsGxUxdNMDu8Zcl4NDyud/uLPqhJ6miTeqv8ef3xI13NoIx6KtYhKWlRrUD6LauIs/k1TJp22rUZuDhwcjZyHXylr3E6WHdlB0URCNihgK87+9fmSoanHEPErEh9Xusl7zyTHlEYNqVUPDVr61VsR28ebPcZNaKEf8JcSZ4NZdvYRn11xMe4TJfOEXf1JTMUXZUbSWS5x/1ozMy27EKZPhVwx55Fg0Q9Jds2HKIwZZrJjs4s2fN48OEzUCOskEF+U5Bw8P1ixfHrajKImM5Hh84F3t3LutN9Wy6H7CSpxUw/sZ5JFj0ajck4mEKY8YZLFisos3e6L224hCkgku6nNqLUTCigR++cPvZO+ai9m86gI2PdmXuUkoyW7Z+xnk0X+jGZLumg1THjGotmIKM3vUMofYxZstJRFuuTwdh2jSCS5Ow6iwmlxL1zzCp9dv5/hyC22t5dBs8kaYhKrtlq9cPGvMey23CIePHB0xowENz4pvZKXjiYJFW3l4z1ceHeWjOPNtU3nouvOASp/roGgq74QSlIex9kMLqr7mzZ3zY5sAjOgsX1RpFlWtl0UYInDqia115zBE6afhcqKvLI0/Mung4UFayyW+esWCQFnqifqLSrW6Szd3zh/TX+XQwOCIc723f4DrfrCdr3x4QaL8p3ooetJds2E7Dwe/4oCKs/s9X3mUFbf/LDDs8rhJLSMrps6F7dx6xYJRq6lbQ37gfuJE5RjREEaXAE+yY1ixaBabV10wYhKqZ2Xs9tO49YoFVGtL77dixvWzNcIkVMsE5L7XvWsu5sjRoTHRg8NaWWil3T/daCy283CoFi0V9tibR4fHlDxIMsGsu3pJqIIyknF8uUTH6ceUsvu9eEuyeznzbVN5pu9wQxLISi3CcEjeTr8v/DWuGaoRZTfilB8JKvzoYrkWzY0pjwZSLat2zoy3mPJIkaCkvmrmq8NHhlNtsuP/rs8/awabnuyrabbym5eSmKEaUSgxTRNQ0j4fRr6Y8mgQXT29fObuHSO9tXv7B/jM3cdKOVi+R/q4SX3e7nFhk3dv/wCzVz0YadfhVwz+uleLz5jGY/sOjcqgjuLXCjIvNbL6axblO0SgVjS6has3J+bzcDjzbVNDx8N8EnF8FZ//4c4RxeEyNFwp5QDVf0AWzJucG+/fBRxzPNeiVq/ooETRzU8fGNVrevPTByJlgntxo8L8O4ZG9evIqmf2ikW1+6NYuHpzYsrD4aHrzhujQNxoq6BSE0vnTmfd1Usin//1I8GTiTte7Qdk67Lk9A8MjuwU4kzoYTvBuOeJihsVFoTXAe2WGk87ezyr8h3+ENkgzpgxZdT9sBJARrEws5UHNyw3iDiKIglxW6Aa0XFNTHFodNHBO7fs48c7XhzVvCrLgoJ+sizf4fpH5q7eEHi+Z/oOj9zO8j0a6WI7jwYRtu5yx2/unF/VPGUb++S4vok4RGllmjb9Tj6Eaw67dv12Ftx0rOZVlgUFG1G+I4qCyqNoopEM23lkQFBU1YrFswJ3Fis8PbPDjim1yBh/iQFTyi1Mm3ocLzj+hzDappRjN9NafMa0wO8xyIEdxNK503n21YHY0VZ++gcGR1beaWaP+0PDTz5hMi+/dmTMcVErIERxtldLLnTJo2iikQzbeaRMWOXdjtOn1yyPcHPnfKZOHlvGYmhYbecRwBtHh1m5bB57a1SuPeQ0dQoLighi1wuvBX6PMLa0xtK508d8r+uuXjIqwfDmzvkjSYLlalmCAbgr77Bdjz8rvRZBOUUvv3aEk0+YnKh8R1Rne5T6UnkUTTSSYW1o68S/4po8SRgISIyK2o52zqoHzUEeA5FK971avd2nTSkzZfKk2Ct/P0nbCrv4u+/FYdqUMocOD45JciyXhLUfjN7lsFqL3GfXXBxbrjBfRklkTO5MrR2KtZbNH2tD2wCCupMNDNbnaA1LCjOCUQ1vSuTl4OHBMdnbSajXfLJ24+7EXSHD3ufgkDas9WwQcZzttZILG5Ehb6SDKY86iBPGGHXbHZYUdny5JdIkaYQTppiF0eHQ5ZJwdFgDk9uSmk9c/0lWC4M8fQK1fBlxkw8bkSGfN2n0cM8bUx51EDWMsUVqF+bz9zGfOrnE4SNDI87WH+94MQWJJzbnnzWDe7f1jlHMH3hXO5ue7OOF/gHappT5zRtHAxVHUGZ3lIkxyBQThlAJnIjrXI+j1ML6qyct0BkWZr580UzrHe6jq6eXmx7YNWoh2KzhyOYwr4OoYYzDCt3PhdetCupj/vqRIVYsnsXKZfO4d1sv/QO266iXTU/2BWZru87svWsuZsrkSQyGRLadM+vEUT/uqI7iOImFCuzt+w2bV10QGgTgv+rilitJI+nVS7VeGetCcpfCxscz7iIiyILQjOHItvOogziJfdWKv1XL7g3qDGck44X+gTE2dfcH645XM/9seebgqPvVvjfvdx3XpOTuCsJMmN6dUlKTR9pJr2G+jLC9+UQMCqm1iGi2cOSmirYSkfcBXwNKwLdUdU3YsfVGW0WpiloS4YwZU9jT93rN4m9Q6S/xva37sJQNw2gc06aUmVySwDyW8UxJ4MsfjtZTyEvUaKumUR4iUgJ+BbwH2A/8HFiuqk8EHV+P8ohjowZL4jMMo7hEbUrnElV5NJPP41xgj6o+o6pHgO8Dl2XxQnGL30VRHJbkZxhGHmTlS2km5dEOeI3M+52xEUTkGhHpFpHuvr6+xC+Uhe3R9iWGYeRBVr6UZlIeNVHV21S1Q1U7ZsyYkfg8WZRCsJ4FhmHkQValXZpJefQC3uI4pzljqbNy2Txay2NrTIVRqlGraOnc6ZELzBmGYaRJFl0noblCdX8OnCkic6gojY8AH83ihYJKJIRFWy1fNJOO06eHZg/74+ct2sowGotFW2WTeNg00VYAInIRcCuVUN3vqOrfhR3bqMKIhmEY44lxWRhRVTcAG/KWwzAMY6LTTD4PwzAMoyCY8jAMwzBiY8rDMAzDiI0pD8MwDCM2TRVtFQcR6QOeS/j0k4D/SlGcRtBsMjebvNB8Mpu82dNsMkeR93RVrZllPW6VRz2ISHeUULUi0WwyN5u80Hwym7zZ02wypymvma0MwzCM2JjyMAzDMGJjyiOY2/IWIAHNJnOzyQvNJ7PJmz3NJnNq8prPwzAMw4iN7TwMwzCM2JjyMAzDMGJjysOHiLxPRHaLyB4RWZWjHN8RkVdE5Beeseki8pCIPOX8n+aMi4h83ZH5cRE5x/Ocq5zjnxKRqzKUd6aIbBKRJ0Rkl4h8qglkPl5E/lNEdjgy3+SMzxGRrY5s60VksjN+nHN/j/P4bM+5Vjvju0VkWVYyO69VEpEeEflxk8j7rIjsFJHtItLtjBX5umgTkXtE5EkR+aWILCmqvCIyz/lc3b9fi8i1DZFXVe3P+aNS6v1p4AxgMrADeHtOsvwhcA7wC8/Y/wZWObdXAV9ybl8E/CuVVumLga3O+HTgGef/NOf2tIzkPQU4x7l9AvAr4O0Fl1mAtzi3y8BWR5YfAB9xxv8B+Avn9l8C/+Dc/giw3rn9dudaOQ6Y41xDpQyvjeuA7wE/du4XXd5ngZN8Y0W+Lu4A/ty5PRloK7K8HrlLwEvA6Y2QN7M30ox/wBJgo+f+amB1jvLMZrTy2A2c4tw+Bdjt3P5HYLn/OGA58I+e8VHHZSz7j4D3NIvMwBTgMWARlQzcSf5rAtgILHFuT3KOE/914j0uAzlPAx4GLgB+7Lx+YeV1zv8sY5VHIa8L4ERgL04wUdHl9cn4XmBzo+Q1s9Vo2oHnPff3O2NF4WRVfdG5/RJwsnM7TO5c3o9jHllIZSVfaJkdE9B24BXgISqr8H5VPRrw+iOyOY8fAt7aYJlvBf4GGHbuv7Xg8gIo8BMR2SYi1zhjRb0u5gB9wD85psFvicjUAsvr5SPAXc7tzOU15dGkaGV5ULg4axF5C3AvcK2q/tr7WBFlVtUhVV1AZUV/LnBWvhKFIyLvB15R1W15yxKTP1DVc4ALgb8SkT/0Pliw62ISFXPxN1V1IfA6FbPPCAWTFwDHz3UpcLf/sazkNeUxml5gpuf+ac5YUXhZRE4BcP6/4oyHyd3Q9yMiZSqKY52q3tcMMruoaj+wiYrZp01E3C6b3tcfkc15/ETg1QbKvBS4VESeBb5PxXT1tQLLC4Cq9jr/XwF+SEVJF/W62A/sV9Wtzv17qCiTosrrciHwmKq+7NzPXF5THqP5OXCmE70ymco28P6cZfJyP+BGQVxFxa/gjn/MiaRYDBxytqwbgfeKyDQn2uK9zljqiIgA3wZ+qapfaRKZZ4hIm3O7lYqP5pdUlMgHQ2R238sHgUecVd39wEec6KY5wJnAf6Ytr6quVtXTVHU2lWvzEVVdUVR5AURkqoic4N6m8n3+goJeF6r6EvC8iMxzht4NPFFUeT0s55jJypUrW3mzdOA04x+VaIRfUbF9fz5HOe4CXgQGqayGPkHFXv0w8BTw/4DpzrEC/B9H5p1Ah+c8fwbscf4+nqG8f0Bla/w4sN35u6jgMr8D6HFk/gXwBWf8DCqT6R4qZoDjnPHjnft7nMfP8Jzr88572Q1c2IDr4zyORVsVVl5Hth3O3y73N1Xw62IB0O1cF11Uoo+KLO9UKjvKEz1jmctr5UkMwzCM2JjZyjAMw4iNKQ/DMAwjNqY8DMMwjNiY8jAMwzBiY8rDMAzDiI0pD8OoAxEZ8lU1HVOJWUT+VEROTXj+80Tk9+uX1DDSZVLtQwzDqMKAVsqbVONPqeSRvJDg/OcBvwF+muC5hpEZtvMwjAwRkQ8CHcA6Z2fSKiLvEpF/cwoFbvSUkfhrqfRDeVxEvu8UmPzvwKed5/63HN+KYYzCkgQNow5EZIhKpq7LLaq63nfMo8D/VNVup/7XvwGXqWqfiFwBLFPVPxORF4A5qvqmiLSpar+I3Aj8RlX/vjHvyDCiYWYrw6iPKGYrL/OA3wUeqpQDo0SlDA1UymGsE5EuKmUxDKOwmPIwjMYiwC5VXRLw2MVUOkheAnxeROY3VDLDiIH5PAwje16j0poXKoUIZ4jIEqiUsReRs0WkBZipqpuAz1Ipn/4W33MNozCY8jCM+mj1hequCTjmn4F/cDoWlqiUR/+SiOygUn34953xO0VkJ5VKv1/XSo+RB4A/Noe5UTTMYW4YhmHExnYehmEYRmxMeRiGYRixMeVhGIZhxMaUh2EYhhEbUx6GYRhGbEx5GIZhGLEx5WEYhmHE5v8DNDV0dVVNYBYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(Y_test, np.abs(yp-Y_test))\n",
    "plt.xlabel(\"E test\")\n",
    "plt.ylabel(\"|E_p-E_test|\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a63fa143",
   "metadata": {},
   "source": [
    "## Summary and conclusion\n",
    "We have seen that the query by commmittee algorithm does indeed reduce the number of electronic structure calculations needed for building a PES. However, it collects a lot of outliers, which makes it perform poorly in low-energy regions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06911595",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}