{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recreating Ling _IMMI_ (2017)\n",
    "In this notebook, we will recreate some key results from [Ling et al. _IMMI_ (2017)](https://link.springer.com/article/10.1007/s40192-017-0098-z), which studied the application of random-forest-based uncertainties to materials design. We will show that the errors produced from the Random Forest implemented in lolo (the code used by Ling et al.) are well-calibrated and that the uncertainties can be used with Sequential Learning to quickly find optimal materials within a search space.\n",
    "\n",
    "Note: This notebook will require you to install [lolopy](https://pypi.org/project/lolopy/) and establish an account with Citrination to get an an API key (see [Quickstart](https://citrineinformatics.github.io/api-documentation/quickstart/index.html)), and set it as an environment variable named CITRINE_KEY. Also, the uncertainity calculations do not currently function on Windows. \n",
    "\n",
    "Last used with matminer version 0.4.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "from matminer.data_retrieval.retrieve_Citrine import CitrineDataRetrieval\n",
    "from matminer.featurizers.conversions import StrToComposition\n",
    "from matminer.featurizers.base import MultipleFeaturizer\n",
    "from matminer.featurizers import composition as cf\n",
    "from lolopy.learners import RandomForestRegressor\n",
    "from sklearn.model_selection import KFold\n",
    "from pymatgen import Composition\n",
    "from scipy.stats import norm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the random seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the Datasets\n",
    "The Ling Paper used 4 different datasets to test the uncertainty estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdr = CitrineDataRetrieval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 165/165 [00:03<00:00, 46.07it/s]\n"
     ]
    }
   ],
   "source": [
    "data = cdr.get_dataframe(criteria={'data_set_id': 150888}, print_properties_options=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the composition and class variable from strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f898fcd82d484dd78c460573a779f103",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>HBox</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='StrToComposition', max=165), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "data = StrToComposition(target_col_id='composition').featurize_dataframe(data, \"chemicalFormula\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['ZT'] = pd.to_numeric(data['ZT'], errors='coerce')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.reset_index(drop=True, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Features\n",
    "Every dataset except the steel fatigue dataset uses the composition-based features of [Ward et al.](https://www.nature.com/articles/npjcompumats201628)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "f =  MultipleFeaturizer([cf.Stoichiometry(), cf.ElementProperty.from_preset(\"magpie\"),\n",
    "                         cf.ValenceOrbital(props=['avg']), cf.IonProperty(fast=True)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e7cbc101a00949beb07d4537a8fd8c20",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>HBox</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='MultipleFeaturizer', max=165), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "X = np.array(f.featurize_many(data['composition']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the Residuals and RF Uncertainty\n",
    "As described in the Ling paper, ideally-calibrated uncertainty estimaes should have a particular relationship with the errors of a machine learning model. Specifically, the distribution of $r(x)/\\sigma(x)$ where $r(x)$ is the residual of the prediction and $\\sigma(x)$ is the uncertainty of the prediction for x should have a Gaussian distribution with zero mean and unit standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = RandomForestRegressor()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the errors from 8-fold cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = data['ZT'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_resid = []\n",
    "y_uncer = []\n",
    "for train_id, test_id in KFold(8, shuffle=True).split(X):\n",
    "    model.fit(X[train_id], y[train_id])\n",
    "    yf_pred, yf_std = model.predict(X[test_id], return_std=True)\n",
    "    y_resid.extend(yf_pred - y[test_id])\n",
    "    y_uncer.extend(yf_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the normalized residuals ($r(x)/\\sigma(x)$) against the normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcYAAACICAYAAAB9R0uBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcFMX9//HXm2WXBTEIgooRRePx88ALifFKTLxvDd6a4EmIGolXQOOBGqOIxiMeXzyIiKKoCGo8QKICigqogGBCQI6AAiIgCCrs8fn9UbU4jDOzs8vOzO7yeT4e+9jp6uvTPV1d3dU9VTIznHPOORc0KXQAzjnnXH3iBaNzzjmXwAtG55xzLoEXjM4551wCLxidc865BF4wOueccwkaZcEoaaWk7TKMnyPp0BzH0EfSE7lcR30k6f8kXVfoOGpC0kGSphc6Dpcfkq6R9Eih43D1V70vGCVdLemVpLQZadJOBzCzlmY2K6Y/Jukv67H+cyS9Xdv5a7G+jpJMUtOk9PXajvWIp0YXEWbWw8xuznLZ6/vdzJH0bbwQqvq7L4v5TNL2CTGPNbOdahtHNesqyPdWH0h6S9IySc0KHUsiM/urmV2Qi2UruFTSVEmrJM2X9KykTrlYX1xnnZ6jsrmor23eayjqfcEIjAEOkFQEIGkLoBjYOylt+zitqwPJBXM9dly8EKr6u6TQAblwgQccBBhwfEGDya97gJ7ApUAbYEdgOHBMIYPKkazyXqpzSU3PL3k/H5lZvf4DSoBvgM5x+FTgH8DopLSZCfMYoaDsDpQBa4CVwEtx/BzgSmAKsBwYApSmWf85wNtpxm0JvAgsBWYCFyaM6wM8kTB8PDAN+Ap4C9g5zTI7xvibJqU/BvwlMSbgDmAZMBs4KmHaNnEffR7HD08YdywwKcYxDtg9YdwcoFfcL6uBp4BK4Nu4//4Up3sWWBj33Rhg1zRxHgzMB64AvgAWAOfGcT/4boCrgKFJ2/134O40+2oOcGiacdvHY2Q58CUwJKaPift3VVzvaVVxJi33qrgfVgGPApsDrwJfA6OA1gnTp9wfqbYx4bgZCiyO392lCcv6KTARWAEsAv5W6DxYy3x7PfAO8Dfgn0njjgY+ifvyM+DKmN4W+Gc8NpcCY4EmiXm6muPsTwnH2YlxPf+Ny7omVd7k+/zWDfhfPFb+nDBtc2AgIR/9O65jfppt3gGoAH6aYb+0Ah6P3/1c4NqEbTyHzPn6HGBW3G+zgbOAnYHv4npXAl/FaY8BPorH0TygT4pzzA+2GTgyHq9lcXmTa5H3zonf/V1x3/8lTVqTuP1z4/f2ONAqKcbzY4xj8nr8FjoDZZnJ3gQui5/vA84DbklKG5Aw/dpMREIGSvpSxxNOUG3iAd8jw5ecrmAcDTwAlAJ7xoP9kBSZb0fCCfYwwt3unwgFaUmKZVYdENUVjGXAhUAR8HtCIag4/mVCYd86ru8XMX3veADuG+frFvdFs4T9MgnoADRPlwHi/t8YaAbcDUxKE+fBQDlwU4zjaMJFTutU3w3QPu6nTeJw0xhv51pkzqeAPxMyXylwYKrjIyHO5ILxPUJh+OMYw4fAXnGb3wBuqOn+iMNNgA8IBUcJsB3hZHdEHP8u8Jv4uSXws0Lnv1rm2ZnARUDneKxunjBuAXBQ/Nwa2Dt+vhX4v3isFBPuOKuO6eoKxvK4T4sJ+WIxMDh+L7sSCo/tUuTNjnHZDxMKwT0IF4U7x/G3EfJ5a2ArwsVSuoKxBzC3mv3yOPBCjKsjoeA+v7p8DWxEKOR2SsgruybM93bSeg4GOsXjbXfCRdaJWW7z2v2TYTvmkLlgLAf+QMjDzdOknRePk+0Ix/rzwKCkGB+P2948n8dvQ6hKhXBg/jx+PohwJTk2KW10DZd5r5l9bmZLCXcre9ZkZkkdgAOBXmb2nZlNAh4BfpNi8tOAl83sdTMrI1wRNgf2r2HMieaa2cNmVkG4om0PbC6pPXAUoaBfZmZlZla1by4E+pvZ+2ZWYWYDCRniZwnLvdfM5pnZt+lWbGYDzOxrM1tNyER7SGqVZvIy4KYYxyuEq9CUz/PMbAHhjuuUmHQk8KWZfZBhPwyX9FXC34UJ690G2DJ+PzV9BvN3M1tkZp8RjrX3zeyjuM3DCIVkVdw12R9dgHZmdpOZrbHwLPxh4PSEuLeX1NbMVprZezWMu+AkHUjY98/E7+5T4MyEScqAXST9KB6jHyaktwe2icfLWItnySyUAbfE/PU04e7znvi9TCPU1uyeYf4bzexbM5sMTCYUFhBqo/4a45wP3JthGZsSCv2U4qOf04CrY1xzgDtZ95yRMl/HcZXAbpKam9mCuF0pmdlbZvaxmVWa2RTCheIvstzmbKXLewCfm9nfzaw84VySnHYWoUZklpmtBK4GTk+qNu1jZqsynY9yoaEUjGOAAyW1JpxUZhCqAfePabtR8+eLCxM+f0O4YqmJLYGlZvZ1Qtpcwh1GqmnnVg2YWSWheiPVtOXxf3FSejEh81dZG7+ZfRM/tiTc7S01s2Uplr0NcEXiwRyn3zJhmnkp5ltLUpGk2yR9KmkF4coRwokolSVmVp4wXN2+HgicHT+fDQzKFA/hKniThL+HY/qfCFfa4yVNk3ReNctJtijh87cphltCrfbHNsCWSd/BNXx/8jufUMPwH0kTJB1bw7jrg27ASDP7Mg4PjmlVuhJqD+ZKGi1pv5jej3AHMVLSLEm9a7DOJbEwgfD9QJrvLI1054MtWTdPZMofSwgFWTptCbUEcxPSks8ZKfO1ma0iFKo9gAWSXpb0/9KtSNK+kt6UtFjS8jhf8jG5vufAdHkPUu+n5LR1zovxc1O+zwvplpNz1RaMVS+4FNi7hLr57oR6asxsBaGaoTvhSmR2mnmzveKsqc+BNpI2TkjbmvDMJNW021QNSBKhQEo17QJCAdgxKX1b1j2I0pkX49okzbhbkg7mFmb2VMI0yfsrefhM4ATgUMJ3UhWnsogtWarvZjiwu6TdCM9Dn6zFcjGzhWZ2oZltCfwOeCDxTdQ6VN3+SN7GecDspO9gYzM7OsY9w8zOADYD+gLPSdooB3HnhKTmhLusX0haKGkhcBnhLnoPADObYGYnELZxOPBMTP/azK4ws+2A44DLJR0SF/0N0CJhVVvkZ4tYQKhCrdIhw7T/AraStE+a8V/yfU1GlXTnjB8wsxFmdhih8P0PoaYBUuejwYT3HzqYWStCFXW2ebQuzpmplpGcts55kbAvyln3giZX5++MsrljnCmpn6Rdch5NGvE2eiJwOaFaq8rbMS3T3eIiQh32+pCk0sQ/M5tHuGu9NabtTrjaT3UifwY4RtIhkooJL6OsjvOvI171DgVukbSppGJJZwC7EF7+yChWR75KKAhax/mrqpwfBnrEq0lJ2kjSMUmFe7Lk/bdxjH0J4UT11+piqsGyMbPvgOcIGXu8mf2vNguWdIqkqhPaMkIGq7qjqItjokp1+yN5XeOBFZJ6SWoe7zh3k9Qlxn22pHaxVuGrOE8FDceJhHh3ITye2JPwgshY4LeSSiSdJalVrPZcEadH0rGSto8XjlXpVds+CTgz7q8j+WG1YK48A1wd89KPgbRvPcearAeApyQdHLe1VNLpknrHvP0MIW9vLGkbwvmr2t87S9pc0vHxImk14ZFE4vG8laSShFk2JtQcfSfpp6xblV2dRUBHSbmuUXwKuEzStpJaEvLOkKQapoLIZsN3JzwgfkTSe5K6S/pRjuNKZTThCjPxWdHYmJapYHyU8DzjK0nDa7nu/QlVMWv/Yj34GYQ7hM8Jz51uMLPXk2c2s+mEasG/E64ajyO86rwmzfouIry5NYXw4sclwDFmtijN9Ml+Q7gy/U+c/48xjomE54z3EQqLmYSH4pncClwb99+VhIfhcwlXuZ8QXlKprXTfzUDCiwPVVaMCvKR1f0s1LKZ3Ad6XtJJw5dwzoVahDzAwrvfU9Ygfqt8f62xjPDkeRygwZhOOh0cId5sQnqtOi3HfA5weLxYaim7AP8zsf/GufaGZLSQcc2fFaX4DzIlVzz34vup8B8IbvysJtUQPmNlbcVxPwn77Ki6ntnm5pm4ivPE6O8b2HKFgSudSwrbeT4j1U+AkwnsMEF4+WUV44eptwgXggCziaEK4oP6ccG74BeE8AeFlsGnAQklV1dcXATdJ+prwUtIzWayjyrPx/xJJH2aYLl3ey9YAQh4fQ9i/3xH2T8FVvfGV3cThzuMpYBPCAXKzmc3MUWxuAyVpa0KhvkWsMneuXpD0e8LFSr7uWF0BZPWMMd7CDyNcwd5JqBp6CXgl48zO1VCsvrkceNoLRVdoktpLOkBSE0k7Ee7aanpn5BqYbFoTmEH4HWE/M0t8JvZcwrMr59ZbfH6yiFA1eWSBw3EOwluk/Qkvv31F+CnIAwWNyOVctVWpkg5M/g2YpAPM7J2cRuacc84VQDYF44dmtnd1ac4551xjkLYqVeFHt/sD7SRdnjDqR4Tmigqibdu21rFjx0Kt3rm8+OCDD740s3b5Wp/nK7chyDZfZXrGWEJoCaEp4TcxVVYAJ69feLXXsWNHJk6cWKjVO5cXkrJpzKHOeL5yG4Js81XagtFC+5qjJT1mZnnNpM4551yhZKpKvdvM/gjcJ+kHDyLNbEPqY80559wGIlNValWrI3fkIxDnnNsQdOz98jrDc25rjH0YN2yZqlI/iP/Xduek0JNFh9iNiXPOOdfoZNPyzVuSfiSpDaHPrn9I+lvuQ3POOefyL5tGxFvFprl+TWgcuDOhix3nnHOu0cmmYGyq0Cv8qcA/cxyPc845V1DZFIw3ASOAmWY2QdJ2hPZTM5I0QNIXkqamGS9J90qaKWmKJG9JxznnXMFVWzCa2bNmtruZXRSHZ5lZ1yyW/RiZG4I+itD/2g5Ad+DBLJbpnHPO5VS1vWtIakfo3LZj4vRmdl6m+cxsjKSOGSY5AXjcQmOt70naRFL72AO9c845VxDZdDv1AjCW0Ht1RR2u+8fAvITh+THtBwWjpO6Eu0q23nrrOgzBuQ2X5yvnUsumYGxhZr1ysG6lSEvZ1YeZPQQ8BLDPPvtk7g7EOZcVz1fOpZbNyzf/lHR0DtY9H+iQMLwV8HkO1uOcc85lLZuCsSehcPxO0gpJX0taUQfrfhH4bXw79WfAcn++6JxzrtCqrUo1s42rmyYVSU8BBwNtJc0HbgCK4zL/D3gFOBqYCXwDnFub9TjnnHN1KZu3UgWcBWxrZjdL6gC0N7PxmeYzszOqGW/AxTUJ1jnnnMu1bKpSHwD2A86MwyuB+3MWkXPOOVdA2byVuq+Z7S3pIwAzWyapJMdxOeeccwWRzR1jmaQi4k8p4g/+K3MalXPOOVcg2RSM9wLDgM0k3QK8Dfw1p1E555xzBZLNW6lPSvoAOITwo/wTzezfOY/MOeecK4Bs3krtBPw/4Avg314oOueca8zSFoySWhHaSe0ATCHcLXaS9D/ghNh5sXPOOdeoZHrGeDMwEdjBzE4ysxOBHYEJwC35CM4555zLt0xVqYcCu5vZ2jdQzaxC0jXAxzmPzDnnnCuATHeMa8ysPDkxpq3OXUjOOedc4WS6YyyVtBc/7B5KQLPcheScc84VTqaCcQHwtzTjFuYgFuecc67g0haMZvbLfAbinHPO1QfZtHzjnHPObTC8YHTOOecSeMHonHPOJai2YJQ0VNIxkrwQdc451+hl0x/jg8C5wL2SngUeM7P/5DYsVx917P1yyvQ5tx2T50iccy53qr0LNLNRZnYWsDcwB3hd0jhJ50oqznWAzjnnXD5lVT0qaVPgHOAC4CPgHkJB+XrOInPOOecKIJtup54ndDs1CDjOzBbEUUMkTcxlcM4551y+ZfOM8REzeyUxQVIzM1ttZvvkKC7nnNsgJD679+f19UM2Val/SZH2bl0H4pxzztUHaQtGSVtI6gw0l7SXpL3j38FAi2wWLulISdMlzZTUO8X4cyQtljQp/l1Q6y1xzjnn6kCmqtQjCC/cbMW6jYl/DVxT3YIlFQH3A4cB84EJkl40s0+SJh1iZpfUJGjnnHMuVzI1Ij4QGCipq5kNrcWyfwrMNLNZAJKeBk4AkgtG55xzrt5IWzBKOtvMngA6Sro8ebyZpeuSqsqPgXkJw/OBfVNM11XSz4H/ApeZ2bzkCSR1B7oDbL311tWs1jmXDc9XzqWW6eWbjeL/lsDGKf6qk9zBMYAlDb8EdDSz3YFRwMBUCzKzh8xsHzPbp127dlms2jlXHc9XzqWWqSq1f/x/Yy2XPR/okDC8FfB50jqWJAw+DPSt5bqcc865OpGpKvXeTDOa2aXVLHsCsIOkbYHPgNOBM5PW0T6hwYDjgX9XG7FzzjmXQ5neSv1gfRZsZuWSLgFGAEXAADObJukmYKKZvQhcKul4oBxYSngL1jnnnCuY6t5KXS+xxZxXktKuT/h8NXD1+q7HOeecqyuZfuB/d/z/kqQXk//yF6LLt8rKSvr27ctBBx3EwoULuf766znssMNYvXBmoUNzrkEbOnQoC5/sxXfzprJy6hssfLIX30wfV+iwXJJMVamD4v878hGIqz9GjBiBmfG//f/Ez+7+ANiX8h1/wjcTX6Rk8+0AIaV66dg5l87cuXMZMWIEm51yI01KSqEDtNhxf74aO4jSbfdCxc0KHaKL0t4xmtkH8f9oQtuoywjPAd+Naa6RGT16NBdccAFHHXUUvXv3Rvr+8Gi6cVta//I8Vn/2HxYP/yuV360sYKTONRwzZ87k+OOPp3379jz00EOhUIyalJTS5pALsfI1LHrqGmbO9FqZ+qDaRsQlHQN8CtwL3AfMlHRUrgNz+TVlyhRuuOEGbr/99ozTlW61CxvvcSRfvnQHZsk/S3XOJVqyZAlnn302/fr1o6SkJO10RS1asekRF3P22WezfPnyPEboUsmm26k7gV+a2UwAST8BXgZezWVgLn/MjE022YRBgwbRpk2baqdvvl1nijfdCitf49U/zqVhZpSUlNC/f3922mmnaqcv3rQD8/f6Hbvf/CY0aYKaFHk3VAWSTbdTX1QVitEs4IscxePyrKKigq5du1JcXEyHDh2qnyFq2mpzlo99glXT38lhdM41XNdeey3jx49njz32yHqepq02Y9X0t/nq7cE5jMxVJ9Nbqb+W9GtgmqRXYhdR3QjNuE3IW4Qup2688UYOOOAA2rdvX+N5Wx10Fl9PfJFPPvF24Z1LNGzYMObOncuvfvWrGs+70S4HU758Id/817u9LZRMVanHJXxeBPwifl4MtM5ZRC6vdt11V0499dRazdukuJS2x15Oq1at6jgq5xq2oqIi+vfvX6u3tyWx6RF/oHzZ59VP7HIi0w/8z81nIC6/li1bxm233UbfvuvXPG3TVpvT5frhrPr4X7T+1fkA/lzEbbAqKiq49NJLueeee2jaNJtXOFJrUlJK09btOeWUUxgyZAhNmmTz1MvVlWq/OUmlwPnArsDa94zN7LwcxuVy7A9/+APdunWjY++X10mvTaFW3HYbylcs5ru5UyjdZve6CtG5Bufuu+9m6623Xq9CsUqTkuZ07tyZe+65h8suu6wOonPZyuYyZBCwBXAEMJrQS8bXuQzK5dbixYvZYYcdOOyww+pkeZJoc/jvKV++sE6W51xDVFlZyfz587nyyit/MK5j75fX/tXEVVdd5T+LKoBsCsbtzew6YFVsP/UYoFNuw3K5smDBAqZOncoNN9xQp8statGKjXb5JSsmvFCny3WuIVizZg3PP/88d911F0VFRXW23KKiIi6//HLuvPNOysrK6my5LrNsCsaqb+MrSbsBrYCOOYvI5UxlZSXdu3entLS0+olrQU2LKftyLs8//3xOlu9cfdWnTx+WLFlS/YS11KJFC/r06ZOz5bt1ZVMwPiSpNXAd8CLwCd6hcIM0dOhQ9ttvP/bbb7+craP1Id157rnnqKyszNk6nKtPZsyYwaxZs+jevXvO1tGjRw8WLFjA4sWLc7YO971qnxCb2SPx42hgu9yG43Klqr3G4uLijNPV9BlIsiYlpQwePJg33niDAw88MGMzWM41dIsWLaJNmzY8+eSTOW1YXxIDBgxg8uTJlJeX1+p3xy572bSVuqmkv0v6UNIHku6WtGk+gnN1Y+nSpWvbYMzXa9+zZ8/m6qu9q03XeFVUVNCtWzdmzZpVp88VEyW/tPPdd99x7rnnUlFRkZP1uSCbs+TThCbgugInA18CQ3IZlKtbF110EbfeeiubbbZZ3tZ53nnnsWjRIsaPH5+3dTqXT/369ePII4+kS5cueVvnvvvuy+GHH87DDz+ct3VuiFTdq8CSPjCzzklpE81sn5xGlsY+++xjEydOLMSqGxwz47PPPqOkpCRtobi+VaepVP0WcvXq1UB4BrPbbrvV+Xoas5jv8pbHPF/VzLx58ygtLaVt27Ypq1Bzma/MjI5/eoE1i+fSbIvt1xnnMss2X2Vzx/impNMlNYl/pxJ613D1mJnRq1cvBg0alNc7xUTNmjVj5cqV9OjRgw8//LAgMThX1wYPHsyVV16ZtlDMNUmgJix/dwirPvGucXMhUyPiX0taAfwOGAysiX9PA94MQz03YMAAKisr6d27d0Hj6NzvPebs+TsO6HoB21w1vKCxOLe+JkyYwFNPPcXjjz9ekEKxipoU0e64q1j1yVuUL/fOjupaprZSN85nIK5ulJeX8/rrr3PmmWdSWlqat2qeTJq2bMNmp95I2eI5vPPOOxxwwAF5Xb9z68vMePXVVznssMN49tlnadbsh/2Q5jtfqWkJ7bpeR+U3Kxg6dChdu3bN6/obs6wa9JN0PPDzOPiWmf0zdyG52poxYwbdu3ena9euHHXUUQWLI9UJQmpC0Uatuf322+nUqRM33HBDtT8dca4+WLp0KRdffDFbbLEFhxxyyDoNZOS7MEwmNUElzRkzZjRDhw7l/vvvp3Vr7/xofWXTiPhtQBfgyZjUU9KBZlbYOjq3VllZGStWrGDKlCk88MAD7LzzzgXPsKkUbdSa4cOHM2TIECorK5kxYwY77LBDocNyLiUz44svvmD27Nmcf/75HHrooQWNJ12eblLcjHvuvIeRI0fy7bffsluvZ2naesu1tUX+Yk7NZXPHeDSwp5lVAkgaCHwEVFswSjoSuAcoAh4xs9uSxjcDHgc6A0uA08xsTk02oKFKPMiTD9xM4xJtfcVQVrw/lG9nfchjd97AGWecUfeB1rFtr34F2Jir3hvOgYuGs3TpUq6++moOPvjgQofm3FoPPvggTzzxBEcccQTXX399ocPJyuGHH46Z8fVHr7BmwQxa7nkkLXereUfJLsuqVGATYGn8nFWvtJKKgPuBw4D5wARJL5pZYnfv5wPLzGx7SacTmpo7LcuYNihlZWUUFxczbNgwxo0bx7Jly1Cb4ynZYgda7X86Z5xxfKFDrJGi0pY89dRTLFq0iKVLl/LSSy8xcOBA9t13X0477TS22GILiouLC/qCg2v8KioqkMS4ceMYNWoUEydO5Nlnn2Xbbbdl5MiRbLTRRoUOsUYk0eaQC6lc8y1lS+ZRtvQzjjnmGLp06cIJJ5xAp06dkJSzBgkai2wKxluBjyS9CYjwrDGbJk1+Csw0s1kAkp4GTiC0tVrlBKBP/PwccJ8k2Xr0s/Lqq69y2rUPANByr6Ppf8JWjBw5Egh9EE6bNo033ngDgOfL9mTNwpmsnjeVbvt35KqrrmLcuHG88847APTu3Zs333yT999/H4Brr72W1157jarfe/Xp04fhw4czadIkAG6++WaGDBnC1KlTAXih+BesmvovypbM47f7daRfv3489NBDTJ8+nSXvzaXNId1ZMXE4G+9xL2aVbHpUT74aO4iyL/8HFeVsdkofWnY6lPKvFtKktCWbdb2Om3ddyUknncRee+3Fzje+QYvtf1rbXVVw61YNNaF88+MZ894M7p45kjsO3oj+/ftjZtxxxx1MnjyZYcOG0aJFC2666SYmTZrEqFGjKCoqonfv3owfP54xY8Ygid69ezN27FjGjRsHwJ///Gdef/11JkyYAMB1113Hyy+/vPYnJDfeeCNDhw5lypQpANxyyy0MHjyYadOmAdC3b18GDBjA9OnTAbjjjjt48MEH+fTTT4HQB9/dd9/NnDlzkMS9995L3759mT9/Pk2bNuWuu+7illtuYeHChZSWltKvX7987N46NXXqVPr37w/Asccey+abb86jjz4KwIknnkirVq0YOHAgACeffDLNmjXjySfD05fTTjsNSTz99NMAnHXWWaxevZrnnnsOgG7durF8+XKGDw9vLZ933nksXryYl156CYDu3bszb948Xn31VSA0WDFjxoy1+bpnz558/PHHjBo1CoArr7yS8ePH88Ybb1BRUcF1113HqFGj6NnvH1SWrabN4b/n4q0X8dprr1FUVMTgwYNZvnw5BxxwAH/84x/Z+cZwfuCtt3K5S3OqSUlzmrXfEYAnb3+SiRMnUlZWxkcffUSvXr2oqKige/fubLXVVtx22220aNGC888/n0022YRHHnmEpk2bcuaZZ1JSUsITTzwBhO+tsrKSIUNC+y6//e1vWbVq1dpOA84991yWLl3KCy+EHnYuvPBCPv/8c15+OeTzHj16MHv2bEaMGAHAxRdfzPTp09d+bz179mTy5Mm8Fff7FVdcwfvvv8/bb78NQK9evRgzZgzvvvsuANdcc806+fr666+nXbt2dbL/Mv7AX+FyfSugnPCcUcD7ZlZtx3uSTgaONLML4vBvgH3N7JKEaabGaebH4U/jNF8mLas7UNVC707A9Ky3MP/aEloHakwa2zY1hO3ZxszqJpen4fmq4Hyb8i+rfJXxjtHMTNLw2PLNizUMIFUdWHIpnM00mNlDwEM1XH9BFLJVoFxpbNvU2LantjxfFZZvU/2VTcs370mqTWOA84EOCcNbAZ+nm0ZSU8Lzy6U455xzBZJNwfhLQuH4qaQpkj6WNCWL+SYAO0jaVlIJcDo/vOt8EegWP58MvLE+zxedc8659ZXNyze1+qW4mZVLugQYQfi5xgAzmybpJmCimb0IPAoMkjSTcKd4em3WVc80iKqpGmps29TYtme6QCyZAAAGHElEQVRD0Bi/M9+meirtyzeSSoEewPbAx8CjZlaex9icc865vMtUMA4ByoCxhLvGuWbWM4+xOeecc3mXqWD82Mw6xc9NgfFmtnc+g3POOefyLdPLN2VVH7wKtWYk9ZH0maRJ8e/oQsdUG5KOlDRd0kxJjaJtXElz4gtkkyR5z7wNiOer+qux5atMd4wVwKqqQaA58E38bGb2o7xE2ABJ6gOsNLM7Ch1LbcUm/f5LQpN+wBlJTfo1OJLmAPskNyLh6j/PV/VXY8tXmfpj9Mb0NmzZNOnnnKsZz1cNQDa/Y3S1c0n83ecASQ2xg7QfA/MShufHtIbOgJGSPohNormGxfNV/dSo8pUXjLUkaZSkqSn+TgAeBH4C7AksAO4saLC1k1VzfQ3QAfElsqOAiyX9vLoZXP54vmqwGlW+yrbbKZfEzLLqtVTSw8A/cxxOLmTTpF+DY2afx/9fSBpGqNoaU9ioXBXPVw1TY8tXfseYA5LaJwyeBEwtVCzrIZsm/RoUSRtJ2rjqM3A4DfO72SB5vqqfGmO+8jvG3Lhd0p6EKpI5wO8KG07NpWvSr8Bhra/NgWGx8+OmwGAze62wIbka8HxVPzW6fJWxP0bnnHNuQ+NVqc4551wCLxidc865BF4wOueccwm8YHTOOecSeMHonHPOJfCCMUckmaQ7E4avjI0g5zOGxySdHD8/ImmX9VxeR0k/+H1STP82tqz/iaTHJRXXch2vSNokRXofSVfWZdyuYfK85Xkr17xgzJ3VwK8lta3NzAp9YNYZM7sgxy34f2pmewKdCK15nFqbhZjZ0Wb2VZ1G5hobz1u14Hkre14w5k458BBwWfIISdtI+ldsDPlfkraO6Y9J+pukN4G+8UpuoKSRsb+zX0u6PfZ79lrVlaOk6yVNiG1KPqT4S9ukdb4laR9Jx+v7/uymS5odx3eWNDo2AjyiqpWRmD5Z0rvAxdVttJlVAOOJDSNLKpLUL8Y3RdLvYnp7SWNiHFMlHRTT51Sd8CT9OcY4CtgpeVvi57YKXd5UXb2OlfRh/Ns/u6/KNTCet/C8lUteMObW/cBZklolpd8HPG5muwNPAvcmjNsRONTMrojDPwGOIXRN8wTwppl1Ar6N6QD3mVkXM9uN0G/msekCMrMXzWzPeAU6GbgjngT+DpxsZp2BAcAtcZZ/AJea2X7ZbLCkUmBfoKrli/OB5WbWBegCXChpW+BMYESMYw9gUtJyOhOay9oL+HWctzpfAIfFxoxPY9396hoXz1uet3LGm4TLITNbIelx4FJCZquyH+GABBgE3J4w7tl4ZVjlVTMrk/QxoQmpqkzxMdAxfv6lpD8BLYA2wDTgpUyxxem/NbP7Je0G7Aa8Hi+Ii4AF8aSziZmNToj1qDSL/ImkScAOwHNmNiWmHw7srvg8BmgVp5kADIgnjuFmNilpeQcBw8zsmxhvNu1JFgP3KTQbVkE4EbpGyPMW4HkrZ7xgzL27gQ8JV4fpJLbLtypp3GoAM6uUVGbft+FXCTSNV5EPEHrPnqfwEkJppoAkHQKcAlR1DSNgWvKVq8KD+mzbDPzUzPaM1URvSTrezF6My/6DmY1IEcfPCVfmgyT1M7PHkyZJt+5yvq/tSNzWy4BFhKvkJsB3WcbuGibPW563csKrUnPMzJYCzxCqPaqMI1RlAJwFvL0eq6g6eL+U1BI4OdPEkrYhZPZTzazqSns60E7SfnGaYkm7xgf1yyUdmBBrRma2AOgNXB2TRgC/T3hms6NCa/zbAF+Y2cPAo8DeSYsaA5wkqblCy/3HJYybA3SOnxO3txWwwMwqgd8Qrs5dI+V5y/NWrnjBmB93Aolv0F0KnCtpCuEg61nbBccM9jCh+mc4oRolk3OATQmt4U+S9IqZrSFkgr6SJhOeSVQ9XD8XuD++IPBtqgWmMBxoER/6PwJ8Anyo8Fp3f0JNxcHAJEkfAV2Be5K260NgSIxlKDA2YfQdhBPCONbdrw8A3SS9R6jqSb5DcI2P5y3PW3XOe9dwzjnnEvgdo3POOZfAC0bnnHMugReMzjnnXAIvGJ1zzrkEXjA655xzCbxgdM455xJ4weicc84l+P9mVFqtlHWxTwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 468x144 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2, sharey=True)\n",
    "\n",
    "x = np.linspace(-8, 8, 50)\n",
    "\n",
    "# Plot the RF uncertainty\n",
    "resid = np.divide(y_resid, y_uncer)\n",
    "axs[0].hist(resid, x, density=True)\n",
    "axs[0].set_title('With Lolo Uncertainty Estimates')\n",
    "\n",
    "# Plot assuming constant errors\n",
    "resid = np.divide(y_resid, np.sqrt(np.power(y_resid, 2).mean()))\n",
    "axs[1].hist(resid, x, density=True)\n",
    "axs[1].set_title('Assuming Constant Error')\n",
    "\n",
    "for ax in axs:\n",
    "    ax.plot(x, norm.pdf(x), 'k--', lw=0.75)\n",
    "    ax.set_xlabel('Normalized Residual')\n",
    "\n",
    "axs[0].set_ylabel('Probability Density')\n",
    "\n",
    "fig.set_size_inches(6.5, 2)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we compare the error distribution using the Lolo uncertainty estimates (_left_) and the assumption that all entries have the same error (_right_). The normalized residuals for the uncertainty estimates have a distribution closer to the unit normal distribution, which means - as expected - that it better captures which predictions will have a higher error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequential Learning\n",
    "One important use of model uncertainties is to employ them to guide which experiments to pick to find optimal materials with minimal experiments/computations. As described in the Ling paper (and [other](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.054303) [nice](https://www.nature.com/articles/srep19660) [articles](https://link.springer.com/article/10.1007/s10822-015-9832-9)), it is not always best to pick the experiment that the model predicts to have the best properties if you can perform more than one experiment sequentially. Rather, it can be better to pick entries with large uncertainities that, when tested and added to the training set, can improve the models predictions for the next experiments. \n",
    "\n",
    "Here, we demonstrate one approach for picking experiments: Maximum Liklihood of Improvement (MLI). In contrast to picking the material with the best predicted properties (an approach we refer to Maximum Expected Improvment (MEI)), the MLI approach pickes the material with with the highest liklihood of being better than the best material in the training set - a measure that uses both the predicted value and the uncertainty. The MLI  method is equivalent to the [Expected Improvement metric common in Bayesian Optimization](https://sigopt.com/blog/expected-improvement-vs-knowledge-gradient), and balances a tradeoff between picking entries with favorable predictions (*exploitation*) and those with high uncertainties (*exploration*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Pick an initial training set\n",
    "We'll start with a small set of entries from the training set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Picked 10 training entries\n"
     ]
    }
   ],
   "source": [
    "in_train = np.zeros(len(data), dtype=np.bool)\n",
    "in_train[np.random.choice(len(data), 10, replace=False)] = True\n",
    "print('Picked {} training entries'.format(in_train.sum()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert not np.isclose(max(y), max(y[in_train]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Demonstrate picking the entries based on MLI and MEI\n",
    "Just to give a visual of how the selection process works"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.fit(X[in_train], y[in_train])\n",
    "y_pred, y_std = model.predict(X[~in_train], return_std=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For MEI, we pick the highest predicted value. For MLI, we pick the material that has the highest probability of being better than any material in the training set. As we assume the predictions to be normally distributed, the probability of materials can be computed from the [Z-score](https://en.wikipedia.org/wiki/Standard_score) $Z = (y - y^*)/\\sigma$ where $y^*$ is the maximum of the $y$ of the training set. Formally, the probability can be computed from the Z-score using the cumulative distribution function of the normal distribution. For our purposes, we can use the Z-score becuase the probability is a monotonic function of the Z-score (stated simply: the material with the highest probability will have the highest Z-score)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "mei_selection = np.argmax(y_pred)\n",
    "mli_selection = np.argmax(np.divide(y_pred - np.max(y[in_train]), y_std))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted ZT of material #13 selected based on MEI: 0.14 +/- 0.08\n",
      "Predicted ZT of material #1 selected based on MLI: 0.12 +/- 0.12\n"
     ]
    }
   ],
   "source": [
    "print('Predicted ZT of material #{} selected based on MEI: {:.2f} +/- {:.2f}'.format(mei_selection, y_pred[mei_selection], y_std[mei_selection]))\n",
    "print('Predicted ZT of material #{} selected based on MLI: {:.2f} +/- {:.2f}'.format(mli_selection, y_pred[mli_selection], y_std[mli_selection]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this particular iteration, the MEI and MLI strategies pick the same material. Depending on the random seed of this notebook and that used by lolo, you may see that the material picked by MLI has a lower predicted $ZT$ but a higher variance. According to the logic behind MLI, picking that entry will (1) yield a higher liklihood of finding a well-performing material and (2) lead to an improved model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Run an iterative search\n",
    "Starting with the same 32 materials in the training set, we will iteratively pick materials, add them to the training set, and retrain the model using 3 different strategies for picking entries: MEI, MLI, and randomly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_steps = 32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_inds = set(range(len(y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random Selection\n",
    "Just pick an entry at random, no need to train a model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_train = [set(np.where(in_train)[0].tolist())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_steps):\n",
    "    # Get the current train set and search space\n",
    "    train_inds = set(random_train[-1])  # Last iteration\n",
    "    search_inds = sorted(all_inds.difference(train_inds))\n",
    "    \n",
    "    # Pick an entry at random\n",
    "    train_inds.add(np.random.choice(search_inds))\n",
    "    \n",
    "    # Add it to the list of training sets\n",
    "    random_train.append(train_inds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Maximum Expected Improvement\n",
    "Pick the entry with the largest predicted value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "mei_train = [set(np.where(in_train)[0].tolist())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_steps):\n",
    "    # Get the current train set and search space\n",
    "    train_inds = sorted(set(mei_train[-1]))  # Last iteration\n",
    "    search_inds = sorted(all_inds.difference(train_inds))\n",
    "    \n",
    "    # Pick entry with the largest maximum value\n",
    "    model.fit(X[train_inds], y[train_inds])\n",
    "    y_pred = model.predict(X[search_inds])\n",
    "    train_inds.append(search_inds[np.argmax(y_pred)])\n",
    "    \n",
    "    # Add it to the list of training sets\n",
    "    mei_train.append(set(train_inds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Maximum Likelihood of Improvement\n",
    "Pick the entry with the largest probability of improvement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "mli_train = [set(np.where(in_train)[0].tolist())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_steps):\n",
    "    # Get the current train set and search space\n",
    "    train_inds = sorted(set(mei_train[-1]))  # Last iteration\n",
    "    search_inds = sorted(all_inds.difference(train_inds))\n",
    "    \n",
    "    # Pick entry with the largest maximum value\n",
    "    model.fit(X[train_inds], y[train_inds])\n",
    "    y_pred, y_std = model.predict(X[search_inds], return_std=True)\n",
    "    train_inds.append(search_inds[np.argmax(np.divide(y_pred - np.max(y[train_inds]), y_std))])\n",
    "    \n",
    "    # Add it to the list of training sets\n",
    "    mli_train.append(set(train_inds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPQAAACICAYAAADQz8ZEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGi9JREFUeJztnXt8VdW177+/BEgCxCAQqRIxQLFQRWOJr5YqeusDK9Dbq1VqOSC9x7fHnoO9tdpa5NrTeo61LVfbHo+t1F7k4bWCPaJoESpWBQlSiIRHxQCx1QSERJS8x/1jrR02YSd77+z3zvx+Puuz15przrnGXHuPPeeaa4wxZWY4HI7sICfVAjgcjvjhFNrhyCKcQjscWYRTaIcji3AK7XBkEU6hHY4swim0w5FFOIV2OLIIp9AORxbRJ9UC9JShQ4daaWlpqsVwOJJCRUXFPjMrDpcvYxW6tLSUDRs2pFoMhyMpSNodST435HY4soiM7aF7wvyN81lZvTLVYjh6StNHcPhAqqVICN8/8zbOm3BjzPX0KoV+ec/LtLS3cNYJZ6VaFEdP2P1naG6Bwk+lWpK4Uzgg7ONxRPQqhW5obuCCkguY+/m5qRal19DS0kJNTQ2NjY2xV5Z/OSgXBsbnx59WGFRVVZGfn09JSQl9+/btUTW9TqGP63dcqsXoVdTU1FBYWEhpaSmSYqvs/VbIL4JBI+IjXJphZuzfv5+amhpGjhzZozp6zaRYY2sjTW1NHJfnFDqZNDY2MmTIkNiV2dqhvRVye9ZzZQKSGDJkSEyjmV6j0A3NDQCuh04BMSszQFuL95nTL/a60phY71XYIbekf+nuvJk9FEEdlwM/B3KBx8zsxyHyfA2YCxjwFzP7erh6o6GhySl0RhNQ6CzuoeNBJD10ob+VAzcDw/3tJuCz4QpLygUeASb7+adL+mynPGOA7wJfMLPTgG9F0YaIcD10htPW7H32QKFzc3MpKyvj9NNPZ8qUKRw8eDAuIlVXV3P66afHpa54EVahzew+M7sPGAp8zszmmNkcYAJQEsE1zgH+ama7zKwZWAxM65TnH4FHzOyAf83aaBoRCR0K7Z6hM5P2QA8d/ZC7oKCATZs2UVlZyeDBg3nkkUfiLFz6EM0s9wigOei4GSiNoNxwYG/QcQ1wbqc8pwJI+jPesHyumb3QuSJJNwA3AIwYEd1Mp+uhU899f3ibrX9r6FnhtiZoa4V+649K/uxJx/GDKadFXM3555/P5s2bATh06BDTpk3jwIEDtLS0cP/99zNt2jSqq6uZPHkyEydO5LXXXmP48OEsX76cgoICKioqmD17Nv3792fixIkd9TY2NnLzzTezYcMG+vTpw0MPPcRFF13EggULWLZsGW1tbVRWVjJnzhyam5v53e9+R15eHitWrGDw4ME9uychiGZS7HfAeklzJf0AWAf8NoJyoZ7yO8cO7gOMASYB04HHJA06ppDZo2ZWbmblxcXRvYsMPEMX5RVFVc6RJphBjBNGbW1trFq1iqlTpwKQn5/PM888w8aNG1m9ejVz5swhENZ6586d3Hrrrbz99tsMGjSIp59+GoDrr7+e+fPn8/rrrx9Vd6DX37JlC4sWLWLmzJkds9WVlZU8+eSTrF+/nnvuuYf+/fvz1ltvcf755/PEE0/E1KbORNxDm9kPJT0PfNFPut7M3oqgaA1wctBxCfC3EHneMLMW4F1J2/EU/M1I5QtHoIce2HdgvKp0REk0Pekx1G2HnFwY8umoix4+fJiysjKqq6uZMGECl1xyCeC997377rt55ZVXyMnJ4b333uODDz4AYOTIkZSVlQEwYcIEqqurqa+v5+DBg1x44YUAzJgxg+effx6AV199ldtvvx2AsWPHcsopp7Bjxw4ALrroIgoLCyksLKSoqIgpU6YAMH78+I7RQryI6rWVmW00s5/7WyTKDJ5SjpE0UlI/4Frg2U55lgEXAUgaijcE3xWNbOFoaG6gsG8huTm58azWkSzaWno8wx14ht69ezfNzc0dvenChQupq6ujoqKCTZs2MWzYsI5eNS8vr6N8bm4ura2tmFmXr5W6W7AiuK6cnJyO45ycHFpbW3vUpq6IWKEl5Un6uqS7Jd0b2MKVM7NW4DZgJVAFLDWztyXNkzTVz7YS2C9pK7Aa+LaZ7Y++OV3T0NTgJsQyFWv3JsVifAddVFTE/PnzefDBB2lpaaG+vp4TTjiBvn37snr1anbv7t5DcdCgQRQVFfHqq68C3h9CgAsuuKDjeMeOHezZs4fPfOYzMcnbE6KZFFsO1AMVQFM0FzGzFcCKTmn3Bu0b8C/+lhDqm+vdhFim0ub3YnF4B33WWWdx5plnsnjxYq677jqmTJlCeXk5ZWVljB07Nmz5xx9/vGNS7LLLLutIv+WWW7jpppsYP348ffr0YcGCBUf1zMlCka5tJanSzNLmpVt5eblFE+BgxooZ5OXm8dhljyVQKkdnqqqqGDduXGyVNB2C/Tth8GjIz/4/5VD3TFKFmZWHKxvNM/RrksZHK1y60NDshtwZS7uzEouUaIbcE4HrJe3CG3ILb7R8RkIkizPO0yqDicFKrLcRjUJPTpgUSaChySl0xtLWAsrxfKEd3RKNQs/sIn1ePARJJI2tjTS3N7shd6YSeGUVD6+tLCcahf44aD8fuBLvNVTa48w+M5y25qx3m4wX0ViK/ST4WNKDHGsgkpZ0uE66HjozaWuBvPxUS5ERxBLgoD8wKl6CJBLXQ2cwZt4sdwwTYpKYMWNGx3FrayvFxcVceeWVACxYsIDi4mLKyso6tq1bt6ale2Q4Iu6hJW3hiFNFLlBMBjw/A9Q31QNQ1M85ZmQcbT13mwwwYMAAKisrOXz4MAUFBbz00ksMHz78qDzXXHMNDz/88FFp1dXVPb5mqojmGfrKoP1W4APfrDPtcT10mvD8XfD+lujKWBu0fAJ9CiAnxM/1U+Nh8jEBcI5h8uTJPPfcc1x11VUsWrSI6dOns3bt2uhkyQAiHnKb2W5gEDAF+O9EEK0kXXDBDTIYa/c+Y5zhvvbaa1m8eDGNjY1s3ryZc8892iV/yZIlRw25Dx8+HNP1UkU0Q+478CKL/N5PWijpUTP7PwmRLI4418k0IYKe9BgO1ULDezBsPOT2POr0GWecQXV1NYsWLeKKK6445nyoIXcmEs0d+iZwrpl9DCDpAeB1IP0Vusm5TmYsbS1AjucLHSNTp07lzjvvZM2aNezfH1dnvrQhGoUW0BZ03EboaCRph7PjzmDamuNmVDJ79myKiooYP348a9asiV22NCQahX4cWCfpGf/4K8Cv4y9S/HF23BlMDIENOlNSUsIdd9wR8tySJUs6/JwBfvGLX3DSSSfF5brJJBrDkock/Qn4Al7PHGkIopTjghtkMO0t0C+2uY9Dhw4dkzZp0iQmTZoEwKxZs5g1a1bIspWVlTFdO9lEEmh/hJntATCzCrwABxlFQ3MDo/uPTrUYjmgxOzLkdkREJK+tlgV2JD2dQFkSRn2Ti1aSkTg/6KiJRKGDZyMywtQzGDNzk2KZShysxHobkSi0dbGfETS2NdLS3uJ66EykY4E610NHSiSTYmdKasDrqQv8fTgSsSStNcUtUpfBdEQqcT10pIRVaDPLaGsMZ/aZwbS1AIqLUUlvIevXhw4otPO0ykDiFKkkEvfJ22677ZhypaWl7Nu3L6ZrJ5vsV2gX3CBzaW+Oy3A72H0SCOk+mS1E8h76QTO7MxnCJALnOpk+PLD+AbZ9uC3yAs0fe8PtPl1HKxk7eCzfOec7Yaty7pNHuDjhUiQQp9CZjHnRPuNAOPfJbKHn/mgZQn1TPUIU9itMtSi9nkh60g7aWuCDSigqgQHRLR0cinDuk9lCpK+t3gW2AJVBn9v85V/TmobmBgb2G0hOnP7pHUki8Moqju+gnfukx2bgy8DpwHjgUmAO3hKxe9NpvatQOE+rDCUBVmLOfdLHzP6Gt0j7i4E0eQvlRr/6dpJxK2ZkKG3xt+Puzn1ywYIFLFvW4bbAG2+8EbfrJpNIFDpkXBZ/Cdid8RUn/jg77gylvRnPqCT2aZ6euk9mYtTPSB4sH5D05YRLkiDckDtDccvf9IhIFLoF+KWk/9n5hKRF8RcpvjQ0NVCU56zEMo44RirpTUQynvk78FXgBUklZjY36NxnEiJVnOhwneyuh25sgB0roT0jQoxnHhqDfbwPRdvTtjZCL3zV6D3J9pxIJ8WqJX0B+IOkEuAGs0DA5PBIuhz4Od6KG4+ZWch4rpKuAp4CzjazDZHW3xURuU5u+DX8cW6sl3J0Qf55P2J/QTtDBvSJXqn79q71rMyM/fv3k5/f83ZHotDyL7Zf0sXAUuBZSV8jgqifknKBR4BLgBrgTUnPmtnWTvkKgX8C1kXXhK4JLIHT7aTY+5VwXAnM+q94XdYRRElbOzUHGqk7GPH//xHqDwAH4i5TOpOfn09JSUmPy0ei0B2BAM2sUdJXgF8Ca4BIZpvOAf5qZrsAJC0GpgFbO+X738C/AXGzG4/I7LNuGww7DQaPjNdlHUH0BUbGbujliJCwk2JmNrvTcbuZ3Qg8B0SiBcOBvUHHNX5aB5LOAk42s267SUk3SNogaUNdXV3YC4cNbtDWCvt2wAljw9blcGQCPbaHNLP7iEyhQw3LO578JeUAP8WzPgt3zUfNrNzMyouLw//thw1u8OEuz8TwhIxZpsvh6JaYDJz9BezCUQOcHHRcgmd1FqAQz6x0jaRq4Dy8Z/TyWGSDCIbctf6ov9j10I7sIBkeC2/i2X2PlNQPuBZ4NnDSzOrNbKiZlZpZKfAGMDUes9xhh9y1VYCgOK3fvjkcEROxQvuL04VN64y/hvRtwEqgClhqZm9LmidpajTCRktDc0P3rpN1Vd5kWN+CRIrhcCSNaAxlLwE6O7RODpF2DGa2AljRKe3eLvJOikKmbmlobqCwX2HXrpO1Ve752ZFVhO2hJd0saQswVtLmoC3gI522dGsl1toE+99xz8+OrCKSHvpJ4HngR8BdQekfmdmHCZEqTnS7SN2+nWBtcMK45ArlcCSQSN5D15tZNfB74EN/ZnsG8Jj//jhtqW/uZk2rOj9YnVNoRxYRzSz3983sI0kTgcuA3wK/SoxY8aHb4Aa1Wz1f2yFjkiuUw5FAolHoNv/zy8AvzWw5kNZrlHQb3KB2GwweDX3SugkOR1REo9DvSfoP4GvACkl5UZZPKmFdJ2u3uuG2I+uIRiG/hvcu+XIzOwgMBr6dEKniwOHWw7S2t4ZW6OZP4EC1U2hH1hHxe2gz+wRvYixw/He84AdpSceaVqGilezbDphTaEfWEbFC+1E+rwNGmdk8SSOAT5nZ+oRJFwPd2nHXVnmfxU6hE4GZsaqqlpoDn6RalIzh4rHDGDGkf8z1RGMp9gugHW9pnHnAR8DTwNkxS5EAul2krrbKi/c8eFSSpcp+WtvamfuHt/m/b+xJtSgZxfDj+yddoc81s89JegvAzA74zhZpSdgeeuipkJv1KwEllYbGFm5duJG1O/dx44WjuPGC0eFD2jgAGJAXn99iNLW0+OGEDEBSMV6PnZZ0hB8KpdB122DEeUmWKLvZ++EnzF7wJu/u+5gH/sd4rjl7RKpF6pVEo9DzgWeAYZJ+CFwN3JMQqeJAl8ENGhugfi8Uz0q+UFlKxe4PueGJClrbjSe+eQ6fHz001SL1WqKZ5V4oqQL4EnA8MM3MqhImWYwEXCcH9h149Im67d5nL/Sy+vkfd/Lk+khiUkTH/kPNlBxfwG9mnc2o4oHhCzgSRiQLvp8H/Bj4EC+Q303AUOB2Sf9gZi8kVsSe0dDUhetkIEpJL4sjtmnvQX62agdnlw5m1NABca27ML8Pt0z6NMcPSNsplV5DpGtb3Q0UAS/jGZaskzQWWASkp0J3ZSVWtw36FMCg0qTLlCra240fLK9k6MA8fj2znMJ8tyJFthKJpVgfM3vRzJ4C3jezdQBmti2xosVGQ3MXS+DUbvVCDuWkrdVq3FmyYS9/qannnivGOWXOciL5VQfPZB/udC62dTsSSJc9dO22XvX8fPCTZv7thW2cUzqYaWUnpVocR4KJZMh9pqQGvHC8Bf4+/nHarlXS0NTAiQNOPDrxkw/h0Pu9yuTzwRe309DYyn3TTot+KRpHxhFWoc0sNxmCxJuQPXQvC2pQ+V49C9ftYeb5pYw70S2p2xvIygdJMwsd3KBjhjv7Fbq93bh3eSVDBvTjny85NdXiOJJEVir04dbDtFrrsUYltVWQdxwcNzx0wSzi6Y01bNxzkLsmj6OowE2E9Ray0ph5xzve+np5L8zj0FN3d6Tn00Qlo/nG3BdTJVrSONzSxoRTjuerZ2X/n5fjCFmp0EWFQ7i05SRyC4dSOeBoy6Xtg77INYUnd1Eye+iTK/7h/FJyctxEWG9Csa4YnyrKy8ttw4aYV8txODICSRVmFna9t6x8hnY4eitOoR2OLCJjh9yS6oBwrkNDgX1JECeZuDZlDvFs1ylmFnZR9IxV6EiQtCGS545MwrUpc0hFu9yQ2+HIIpxCOxxZRLYr9KOpFiABuDZlDklvV1Y/QzscvY1s76Edjl6FU2iHI4vISoWWdLmk7ZL+KumuVMvTUyT9RlKtpMqgtMGSXpK00/88PpUyRoukkyWtllQl6W1Jd/jpGdsuSfmS1kv6i9+m+/z0kZLW+W1akoyFKbJOof3FAB4BJgOfBaZLytSYQwuAyzul3QWsMrMxwCr/OJNoBeaY2TjgPOBW//vJ5HY1AReb2ZlAGXC5Hy33AeCnfpsOAN9MtCBZp9DAOcBfzWyXmTUDi4FpKZapR5jZK3jhk4OZBvzW3/8t8JWkChUjZvZ3M9vo738EVAHDyeB2mcch/7CvvxneOnD/z09PSpuyUaGHA3uDjmv8tGxhmL+Ub2BJ3xNSLE+PkVQKnAWsI8PbJSlX0iagFngJeAc4aGatfpak/A6zUaFDOQC7d3NphqSBeKuXfsvMGsLlT3fMrM3MyoASvFFiqDhXCf8dZqNC1wDBEQxKgL+lSJZE8IGkEwH8z9oUyxM1kvriKfNCM/u9n5zx7QIws4PAGrz5gUGSAkFEkvI7zEaFfhMY488w9gOuBZ5NsUzx5Flgpr8/E1ieQlmiRl4s4V8DVWb2UNCpjG2XpGJJg/z9Arz136qA1cBVfraktCkrLcUkXQH8DMgFfmNmP0yxSD1C0iJgEp4b3gfAD4BlwFJgBLAHuNrMOk+cpS2SJgJrgS0cWcThbrzn6Ixsl6Qz8Ca9cvE6yaVmNk/SKLxJ2cHAW8A3zKwpobJko0I7HL2VbBxyOxy9FqfQDkcW4RTa4cginEI7HFmEU2iHI4twCh0CSSbpJ0HHd0qaG6e6F0i6KnzOmK9zte/RtLpTeqnfvtuD0h6WNCuO154kqV7SpqDtS/Gqv4trzkvCNWZJSutFtp1Ch6YJ+KqkoakWJBjfkyxSvgncYmYXhThXC9yRYHe+tWZWFrT9MVEXkpRrZvcm8ho+swCn0BlIK148qH/ufKJzDyvpkP85SdKfJC2VtEPSjyVd5/vJbpE0OqiaL0la6+e70i+fK+nfJb0pabOkG4PqXS3pSTxjjM7yTPfrr5T0gJ92LzAR+JWkfw/Rvjo8F8WZnU9IGi3pBUkVvoxjfdl2yWOQpHZJF/j510r6dCQ3VdLZftvyJQ3wfYdP99v4iqRnJG2V9CtJOX6ZSyW9LmmjpKd8G3AkVUu6V9KrwNXB34t/7l/9chskfU7SSknvSLopSJ5vB93vgA9zqT+y+U9fvhclFfh1lwML/RFHgf8db/XLPxjJPUg4Zua2ThtwCDgOqAaKgDuBuf65BcBVwXn9z0nAQeBEIA94D7jPP3cH8LOg8i/g/ZmOwbM9zwduAL7n58kDNgAj/Xo/BkaGkPMkPKuqYryFB18GvuKfWwOUhyhTClT6dW/Ds256GJjln18FjPH3zwVe9vdfAE4DrsQzr73Hl/PdENeYBNQDm4K20f65+4EH8XzWvxuUvxEY5cvzEp7J5FDgFWCAn+87wL3+fjXwv4Ku2fG9+Odu9vd/CmwGCv37VOunX4r3py3/u/gv4AL//rQCZX6+pXgWXkfdUzzrr+0cMc4alOrfrZll5+qT8cDMGiQ9AfwTcDjCYm+a7wIo6R0gsG7tFiB46LvUzNqBnZJ2AWPxfmBnBPX+RXgK3wysN7N3Q1zvbGCNmdX511yI96NcFkH73pW0Hvh6IM3v/T4PPCV1OK3l+Z9r/bpHAj8C/hH4E55yh2KtmV0ZIn2eX6YR794GWG9mu3w5FuGNMBrxglT82ZenH/B6UJkl3TQxYL+/BRhonu/1R5Ia5dldX+pvb/n5BuLd7z14f1Kb/PQKPCXvTIMv32OSnsP7Q0g5TqG752fARuDxoLRW/EcVeb+y4OfQYDvd9qDjdo6+153tbQ2vp7jdzFYGn5A0Ca+HDkWsa8X+K54D/iv+cQ6eD29ZiLxrgZvwRgX3At/G61lfCZG3OwbjKU9fvJFJoG1d3ZOXzGx6F3V1dV/g6Hvf+Xvp49f9IzP7j+BC8ny0g/O3AQWdKzezVknnAP8NzwHoNryABinFPUN3g3nOAUs5OnRMNTDB35+G98OMlqsl5fjP1aPwhm4rgZvluRYi6VRJA8LUsw64UNJQf8JsOl6vGRFmtg3YijeMxjy/5HclXe3LIElnBl3r80C7mTXiDaNvxFP0aHgU+D6wEC9ET4Bz5HnI5QDXAK8CbwBfCDyjS+ov6dQor9cVK4HZQc/kwyWFC6rwEd7QPTCaKTKzFcC38EIPpRzXQ4fnJ3j/vgH+E1juD1dX0X0v0RXb8RRvGHCTmTVKegxvaLfR7/nrCBOyxsz+Lum7eG56AlaYWbQuej/kyLAT4Drgl5K+h/dntRj4i5k1SdqLp2TgKfJ0QkzU+XxRXgSPAPcD/YFWM3vS/wN6TdLFeL3m68CPgfF4vf4zZtYu73XaIkmBof/3gB1RtvEYzOxFSeOA1/3h/CHgG3g9clcswJtoPIwXs265pHy8e3/MBGoqcN5WjpTjP1bc2cUztyMK3JDb4cgiXA/tcGQRrod2OLIIp9AORxbhFNrhyCKcQjscWYRTaIcji/j/y1aig/Ilg4AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x144 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for train_inds, label in zip([random_train, mei_train, mli_train], ['Random', 'MEI', 'MLI']):\n",
    "    ax.plot(np.arange(len(train_inds)), [max(y[list(t)]) for t in train_inds], label=label)\n",
    "\n",
    "ax.set_xlabel('Number of New Experiments')\n",
    "ax.set_ylabel('Best $ZT$ Found')\n",
    "    \n",
    "fig.set_size_inches(3.5, 2)\n",
    "ax.legend()\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this particular case, we find that the MLI strategy finds the best material more quickly than the Random or MEI approaches. In Ling 2017, they evaluate the performance of these strategies over many iterations and find that, on average, MLI finds the optimal materials as fast or better than any other approach. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}