Energy_Prediction_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import multiprocessing\n",
    "\n",
    "import scipy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import RandomizedSearchCV as CV\n",
    "\n",
    "from pymatgen.io.vasp import Poscar\n",
    "from pymatgen.io.ase import AseAtomsAdaptor\n",
    "\n",
    "from ase.visualize import view\n",
    "from sklearn.svm import SVR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Requirements:\n",
    "numpy, scipy, pandas, pymatgen, sklearn, ase\n",
    "\n",
    "\n",
    "Optional: tqdm for progress bars, thundersvm for GPU support"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from thundersvm import SVR # much faster with GPU support\n",
    "\n",
    "def tw(iterable, **kwargs):\n",
    "    \"Add a progress bar using tqdm if available\"\n",
    "    try:\n",
    "        import tqdm\n",
    "        tqdm.tqdm._instances.clear()\n",
    "        # clear existing instances\n",
    "    except AttributeError:\n",
    "        pass\n",
    "    try:\n",
    "        return tqdm.tqdm_notebook(iterable, **kwargs)\n",
    "    except NameError:\n",
    "        return iterable  # could not import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up (basic settings here)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CdTe dataset is available at /ufrc/hennig/sxiexie/CdTe_Archive\n",
    "\n",
    "\n",
    "Big picture:\n",
    "1) Parse structure and energy files\n",
    "\n",
    "2) Generate representations\n",
    "\n",
    "3) Optimize and fit models with support vector regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_cores = 1  # number of cores to use for multiprocessing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### uncomment the following line to use all available cores\n",
    "n_cores = multiprocessing.cpu_count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# directory of POSCARs and energy files\n",
    "working_directory = '/home/jason/Downloads'\n",
    "directory_path = working_directory + \"/LiGe_Archive/both/\"\n",
    "\n",
    "structure_ext = 'poscar'  #file extension names\n",
    "target_ext = 'energy'\n",
    "\n",
    "element_list = ['Li', 'Ge']  # list of elements to consider"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# File Parsing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = (os.listdir('/home/jason/Downloads/LiGe_Archive'))\n",
    "pos = []\n",
    "ene = []\n",
    "for b in a:\n",
    "    b = b.split('.')\n",
    "    if b[-1]=='energy':\n",
    "        ene.append(b[0])#+'.'+b[1])\n",
    "    elif b[-1]=='poscar':\n",
    "        pos.append(b[0])#+'.'+b[1])\n",
    "same = (set(ene) & set(pos))\n",
    "for b in same:\n",
    "    os.rename('/home/jason/Downloads/LiGe_Archive/'+b+'.energy','/home/jason/Downloads/LiGe_Archive/both/'+b+'.energy')\n",
    "    os.rename('/home/jason/Downloads/LiGe_Archive/'+b+'.poscar','/home/jason/Downloads/LiGe_Archive/both/'+b+'.poscar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e01333d2c58745af9d0dc1c4e9fe75ec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=2515), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "def process_entry(args):\n",
    "    \"\"\"parse files\"\"\"\n",
    "    entry_name, element_list, directory_path, structure_ext, target_ext = args\n",
    "    target_file = '{}/{}.{}'.format(directory_path,\n",
    "                                    entry_name,\n",
    "                                    target_ext)\n",
    "    with open(target_file) as f:\n",
    "        lines = f.read().splitlines()\n",
    "    local_energies = [float(line) for line in lines]\n",
    "    # This dataset has per-atom energies but we won't use them\n",
    "    energy = np.sum(local_energies)\n",
    "    \n",
    "    structure_file = '{}/{}.{}'.format(directory_path,\n",
    "                                       entry_name,\n",
    "                                       structure_ext)\n",
    "    structure = Poscar.from_file(structure_file).structure\n",
    "    entry_data = [len(structure), energy]\n",
    "    for element in element_list:\n",
    "        element_percent = structure.composition.get_atomic_fraction(element)\n",
    "        entry_data.append(element_percent)\n",
    "    return structure, entry_data\n",
    "\n",
    "file_stems = [os.path.splitext(f)[0]\n",
    "              for f in os.listdir(directory_path) \n",
    "              if structure_ext in f]\n",
    "# check the working directory for files matching the ext\n",
    "\n",
    "columns = ['Size', 'Total Energy'] + element_list  # pandas header\n",
    "\n",
    "# prepare zipped arguments for parallel parsing\n",
    "zipped_args = [(f, element_list, directory_path, structure_ext, target_ext)\n",
    "               for f in file_stems]\n",
    "parsed_data = list(tw(multiprocessing.Pool(n_cores).imap(process_entry,\n",
    "                                                         zipped_args),\n",
    "                      total=len(zipped_args)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unpack parallel parsing results\n",
    "structure_list, table_data = zip(*parsed_data)\n",
    "structures = {k: v for k, v in zip(file_stems, structure_list)}\n",
    "df = pd.DataFrame(columns=columns, data=table_data, index=file_stems)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get Reference Energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While learning on total energies is desirable in some scenarios, we often want to predict formation energy per atom to compare relative stabilities. Here, we normalize the total energy of each structure by the reference energies. These are the lowest-energy elemental structures in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference_energies = {}\n",
    "for element in element_list:\n",
    "    pure_entries = df[df[element] == 1.0]\n",
    "    pure_entry_energies = pure_entries['Total Energy'].values\n",
    "    pure_entry_energies /= pure_entries['Size'].values\n",
    "    minimum_energy = np.min(pure_entry_energies)\n",
    "    reference_energies[element] = minimum_energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Li': -0.6352893700000001, 'Ge': -3.95538245}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reference_energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute Formation Energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_formation_energy(entry_data, element_list, reference_energies):\n",
    "    element_fractions = {element: entry_data[element] \n",
    "                         for element in element_list}\n",
    "    reference_contributions = [(reference_energies[element] \n",
    "                               * element_fractions[element])\n",
    "                               for element in element_list]\n",
    "    formation_energy = (entry_data['Total Energy']/entry_data['Size'] \n",
    "                        - np.sum(reference_contributions))\n",
    "    return formation_energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Formation Energy'] = df.apply(get_formation_energy, axis=1, \n",
    "                                  args=(element_list, reference_energies))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up DSCRIBE's SOAP Handler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Smooth Overlap of Atomic Positions representation is a well-established method for capturing atomic environments for machine learning. It has few hyperparameters, and the default choices here are reasonable for many systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dscribe.descriptors import SOAP\n",
    "rcut = 7  # maximum pairwise distance to consider. \n",
    "# Stick to the second or third nearest-neighbor distance.\n",
    "\n",
    "# maximum degree and orders to calculation in spherical harmonics power-series solution\n",
    "nmax = 6\n",
    "lmax = 8\n",
    "\n",
    "# Setting up the SOAP descriptor\n",
    "soap = SOAP(species=element_list,\n",
    "            periodic=True,\n",
    "            rcut=rcut,\n",
    "            nmax=nmax,\n",
    "            lmax=lmax,\n",
    "            rbf='gto',\n",
    "            sigma=0.125,\n",
    "            average='off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "adaptor = AseAtomsAdaptor()\n",
    "entry_names = sorted(structures.keys())\n",
    "flattened_entry_list = [adaptor.get_atoms(structures[entry_name])\n",
    "                        for entry_name in entry_names]\n",
    "soap_raw = soap.create(flattened_entry_list, n_jobs=n_cores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2515"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#soap_raw.shape\n",
    "len(entry_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Shape of passed values is (21544, 702), indices imply (2515, 702)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/managers.py\u001b[0m in \u001b[0;36mcreate_block_manager_from_blocks\u001b[0;34m(blocks, axes)\u001b[0m\n\u001b[1;32m   1680\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1681\u001b[0;31m         \u001b[0mmgr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mBlockManager\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mblocks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1682\u001b[0m         \u001b[0mmgr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_consolidate_inplace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/managers.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, blocks, axes, do_integrity_check)\u001b[0m\n\u001b[1;32m    142\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdo_integrity_check\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 143\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_verify_integrity\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    144\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/managers.py\u001b[0m in \u001b[0;36m_verify_integrity\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    344\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_verify_integrity\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mmgr_shape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                 \u001b[0mconstruction_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtot_items\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maxes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mtot_items\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/managers.py\u001b[0m in \u001b[0;36mconstruction_error\u001b[0;34m(tot_items, block_shape, axes, e)\u001b[0m\n\u001b[1;32m   1718\u001b[0m     raise ValueError(\n\u001b[0;32m-> 1719\u001b[0;31m         \u001b[0;34m\"Shape of passed values is {0}, indices imply {1}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpassed\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mimplied\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1720\u001b[0m     )\n",
      "\u001b[0;31mValueError\u001b[0m: Shape of passed values is (21544, 702), indices imply (2515, 702)",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-6da3cec469d9>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msoap_raw\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mentry_names\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'SOAP_{}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msoap_raw\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, data, index, columns, dtype, copy)\u001b[0m\n\u001b[1;32m    438\u001b[0m                 \u001b[0mmgr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minit_dict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    439\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 440\u001b[0;31m                 \u001b[0mmgr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minit_ndarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    441\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    442\u001b[0m         \u001b[0;31m# For data is list-like, or Iterable (will consume into list)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/construction.py\u001b[0m in \u001b[0;36minit_ndarray\u001b[0;34m(values, index, columns, dtype, copy)\u001b[0m\n\u001b[1;32m    211\u001b[0m         \u001b[0mblock_values\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    212\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 213\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mcreate_block_manager_from_blocks\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mblock_values\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    214\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    215\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/managers.py\u001b[0m in \u001b[0;36mcreate_block_manager_from_blocks\u001b[0;34m(blocks, axes)\u001b[0m\n\u001b[1;32m   1686\u001b[0m         \u001b[0mblocks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"values\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mb\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mblocks\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1687\u001b[0m         \u001b[0mtot_items\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mb\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mblocks\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1688\u001b[0;31m         \u001b[0mconstruction_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtot_items\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblocks\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxes\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1689\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1690\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/pandas/core/internals/managers.py\u001b[0m in \u001b[0;36mconstruction_error\u001b[0;34m(tot_items, block_shape, axes, e)\u001b[0m\n\u001b[1;32m   1717\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Empty data passed with indices specified.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1718\u001b[0m     raise ValueError(\n\u001b[0;32m-> 1719\u001b[0;31m         \u001b[0;34m\"Shape of passed values is {0}, indices imply {1}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpassed\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mimplied\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1720\u001b[0m     )\n\u001b[1;32m   1721\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Shape of passed values is (21544, 702), indices imply (2515, 702)"
     ]
    }
   ],
   "source": [
    "data = pd.DataFrame(data=soap_raw, index=entry_names,columns=['SOAP_{}'.format(i) for i in range(len(soap_raw[0]))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each entry has its own 567-dimensional vector representing the chemical environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing SOAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hb5dn48e+t4Z3YGc5OSCChEDaEhL1HgJZRoEAXfQulCwqltA19C2X8+ha6B5QCpUApe7QNmwAhbLJDEgJZZC8ncbyHxv374xzbsiPJx0ayZOv+XJcvS0ePdZ5jW7r1rPsRVcUYY4zpKl+mK2CMMaZ3sgBijDGmWyyAGGOM6RYLIMYYY7rFAogxxphuCWS6Aj1p8ODBOnbs2ExXwxhjepV58+ZtV9XyjsdzKoCMHTuWuXPnZroaxhjTq4jI2njHrQvLGGNMt1gAMcYY0y0WQIwxxnSLBRBjjDHdYgHEGGNMt1gAMcYY0y0WQIwxxnSLBRAPbrjrNG6/+8JMV8MYY7JKTi0k7K65gU0MDO/MdDWMMSarWAvEGGNMt1gA8cAX9YNKpqthjDFZxbqwPGgY91NKmqoyXQ1jjMkqFkA8qPSXURAIZboaxhiTVawLywNRRbEuLGOMiWUBxBhjTLdYAPFAUDTTlTDGmCxjAcQDJ4DYr8oYY2Jl9F1RRKaKyCcislJEpsV5PF9EHncf/0BExnZ4fIyI1IrIdWmtJzaL1xhjOspYABERP3AncAYwEbhERCZ2KHYZUKmq44E/ALd3ePwPwIvprivWgWWMMbvJZAtkMrBSVVerajPwGHBOhzLnAA+6t58CThYRARCRc4HVwNJ0V1TAZmEZY0wHmQwgI4H1Mfc3uMfillHVMFAFDBKRYuCnwM2dnURErhCRuSIyt6KioptVtRaIMV68/697ePJHN2S6GqaHZDKAxPtI3/GdOlGZm4E/qGptZydR1XtUdZKqTiovL+9GNUHUWiDGeDHv7fFsqzsx09UwPSSTK9E3AKNj7o8CNiUos0FEAkApsBOYAlwgIr8GyoCoiDSq6h3pqOi3/rUTpAKmpuPZjTGmd8pkAJkDTBCRccBG4GLgyx3KTAcuBd4DLgBeV1UFjm0pICI3AbXpCh7OSfIQtawvxhgTK2PviqoaFpErgZcBP/APVV0qIrcAc1V1OnAf8JCIrMRpeVycodpm5rTGGJPFMvqxWlVfAF7ocOzGmNuNQNKtAFX1prRUrv1ZiD8cY4wxucuWVxtjjOkWCyCeWAvEGGM6sgDihQ2BGGPMbiyAeGItEGOM6cgCiCfWBDHGmI4sgHhmLRBjjIllAcQT68IyxpiOLIB4otaLZYwxHVgA8UqsBWKMMbEsgHhiXVjGGNORBRBPrP/KGGM6sgDimbVAjPEqEg5nugqmB1gA8UBsEN2YLtFIJNNVMD3AAohXNohujGeRUHOmq2B6gO2S5MHTF+UxqLGZyzJdEWN6CVVrgeQCCyAebMkbbEMgxnSBRqKZroLpAdaF5ZFaBDHGM2uB5AYLIB4IagHEmC6I2CB6TrAAYoxJOY1aF1YusADigc9aIMZ0iY2B5AYLIJ4oavHDGM+iWBdWLrAA4oFgg+jGdEnYAkgusADiiS1DN6YrVK0LKxdYAPFA1MZAjOmKqA2i5wQLIB4IoJbKxBjPNGpdWLnAAogHzjoQY4xXNo03N1gAMcaknAWQ3GABxANbB2JM11gAyQ1JkymKyCjgYuBYYATQACwBngde1JyZamEBxJiuyJm3hhyXMICIyP3ASOA54HZgG1AA7A1MBf5XRKap6ps9UdHMUsvGa0wXqK0DyQnJWiC/U9UlcY4vAZ4RkTxgzGc5uYhMBf4E+IG/q+ptHR7PB/4JHAbsAC5S1TUiMhm4p6UYcJOq/vuz1CVpPdUWEhrTFZGobWmbCxKOgcQLHiIyQEQOdB9vVtWV3T2xiPiBO4EzgInAJSIysUOxy4BKVR0P/AGnJQROEJukqgfjtIbuFpG07W1is7CM6RqN2ismF3Q6iC4ib4hIfxEZCCwC7heR36fg3JOBlaq6WlWbgceAczqUOQd40L39FHCyiIiq1qtqy0ecAtK8VNxSmRjTRTaInhO8zMIqVdVq4IvA/ap6GHBKCs49Elgfc3+DeyxuGTdgVAGDAERkiogsBRYD34kJKO2IyBUiMldE5lZUVHSzqooNghjjna1Ezw1eAkhARIYDX8IZUE+VeO/IHVsSCcuo6gequh9wOHC9iBTEO4mq3qOqk1R1Unl5ebcrag1yY7yzaby5wUsAuRl4Gae7aY6I7AmsSMG5NwCjY+6PAjYlKuOOcZQCO2MLqOoyoA7YPwV1iktULZWJMV1g03hzg5cAsllVD1TV7wGo6mogFWMgc4AJIjLOndF1MTC9Q5npwKXu7QuA11VV3Z8JAIjIHsDngDUpqFNcYl1YxnSJtUByg5cA8hePx7rEHbO4Eqd1swx4QlWXisgtInK2W+w+YJCIrASuBaa5x48BFonIQuDfwPdUdftnrVPS+qbzyY3pY1TtFZMLki0kPBI4CigXkWtjHuqPs27jM1PVF4AXOhy7MeZ2I3BhnJ97CHgoFXXwQmwlujFdohFbSJgLkq2dyANK3DL9Yo5X43QnGWNMXNYCyQ0JA4iqzgJmicgDqrq2B+uUdXy2oZQxXWKD6LkhWRfWH1X1GuAOEdnt44Sqnh3nx/osCyDGeGeD6LkhWRdWyxjDb3uiItlMUNTihzGe2ULC3JCsC2ue+31Wz1Unm1kEMcYzGwPJCZ0mIBSRo4GbgD3c8gKoqu6Z3qplD0umaEzXWDLF3OAlg+19wA+BeUBOzs2zdO7GdI0NoucGLwGkSlVfTHtNspitAzGmayK2DiQneAkgM0XkN8AzQFPLQVWdn7ZaZR1rjhvTFTYLKzd4CSBT3O+TYo4pcFLqq5OdbD8QY7rIPnPlhE4DiKqe2BMVyWbONF4LIMZ4FrUurFzgZRbWjfGOq+otqa9OlrJPU8Z0iaUyyQ1eurDqYm4XAJ/HyZ6bM6wLy5iusWm8ucFLF9bvYu+LyG/Zfd+OPs1mYRnTNYoNoucCL/uBdFQE5MwiQnB3JEQIR6Oc//R8Vu6qz3SVjMlqUWuB5AQvYyCLaRsF8APlQO6Mf8S4a8F65s3ZzCVba5nzveMyXR1jspeNgeQEL2Mgn4+5HQa2ursJ5oyWMZDKynUMCGykoW54pqtkTFazMZDc4GUMJKf3AoG2XFi1n75CeMJ/OGTXOOALma6WMdnLUpnkhO6MgeQkRfCH1zO6Yl/Cvo2Zro4xWc3iR26wAOKBqAJCU6g/8w/7GdGyr2a6SsZkNUummBs6DSAicqWIDOiJymSrlmm8hXmDAdhRPjTDNTImu9lCwtzgpQUyDJgjIk+IyFSR3Mvp4aRzB1Hn12VpTYxJzgbRc0OnAURVfw5MwNkX5BvAChH5PxHZK811yyoWNIzxzlogucHTGIg6/w1b3K8wMAB4SkR+nca6ZQ1x2h8xR+zFYUwyUdsPJCd4WUj4A+BSYDvwd+DHqhoSER+wAvhJequYeS3rQIQo/Zs2ELVPV8YkZS2Q3OBlIeEg4Isd14OoalREPp/gZ/oWdwyk345NXL3oGBqLX4QzvpLpWhmTtdTSueeEpF1Ybivj/ESLCVU1J7LytnRh9atw7pdt3yej9TEm29k03tyQNICo81+wSETG9FB9slSHbLzWOjcmKZuFlRu8dGENB5aKyGxi9gZR1bPTVqss07ofiNinKmO8sDGQ3OAlgNycrpOLyFTgTzhZfv+uqrd1eDwf+CdwGLADuEhV14jIqcBtQB7QjDOw/3ra6qktjQ6bymuMF9aFlRu8JFOclY4Ti4gfuBM4FdiAs1hxuqp+FFPsMqBSVceLyMXA7cBFODPCvqCqm0Rkf+BlYGQ66glxpvHamhBjkrIurNzgJZXJESIyR0RqRaRZRCIiUp2Cc08GVqrqalVtBh4DzulQ5hzgQff2U8DJIiKqukBVN7nHlwIFbmslPbSlC6vlgAUQY5LRqLVAcoGXhYR3AJfgrPkoBC53j31WI4H1Mfc3sHsrorWMuwdJFc604ljnAwtUtSneSUTkChGZKyJzKyoqulXRti1tLXAY44WNgeQGryvRVwJ+VY2o6v3ACSk4d7x3447/dUnLiMh+ON1a3050ElW9R1Unqeqk8vLybldU41bPGBOPLbbNDV4G0etFJA9Y6KYu2QwUp+DcG4DRMfdHAZsSlNkgIgGgFNgJICKjgH8DX1fVVSmoT0It6dxV7EVhjCc2BpITvLRAvuaWuxJnGu9onG6jz2oOMEFExrkB6mJgeocy03HSqABcALyuqioiZcDzwPWq+k4K6tKpaLsuLOvKMiYZm4WVG5K2QNyZUr9U1a8CjaRwSq+qhkXkSpwZVH7gH6q6VERuAeaq6nScDMAPichKnJbHxe6PXwmMB24QkRvcY6ep6rZU1W93AmqBwxhPrAsrJyQNIKoaEZFyEclzZ0qllKq+ALzQ4diNMbcbgQvj/Nz/A/5fquuTiKitRDemK2wMJDd4GQNZA7wjItNpvxL99+mqVLZpCR0xo/cZqokxvYRN480JXgLIJvfLB/Rzj+XUxwsBovhiZmMZY5KJ2iB6TvASQD5S1SdjD4jIbt1KfZq2JTKxl4UxHtgLJSd4mYV1vcdjfVZrMkVjjCe2kDA3JGyBiMgZwJnASBH5c8xD/XG2tc0ZziB6uyMZqokxvYWNgeSCZF1Ym4C5wNnAvJjjNcAP01mpbKTijwkiFkCMScbG0HNDwgCiqotwNpN62M1DlbNilw9aw9yYxNYFIqz3R/mSdWHlhGRdWE+o6peABSK75/BQ1QPTWrOs4lx+WNVb8jBjctTjJc5ysQtVeeAbv6FwQBUX/aHHlmyZHpasC+tq9/vne6Ii2cxSYBnTRarUFRxGXUOmK2LSKeEHalXd7H5fCzQBBwEHAk3usZwhLS2Q/sXcetFANhxSE7dcYzjCwX9/m6dWpjGjijG9gH3myg1eNpS6HJgNfBEnoeH7IvLNdFcsm7S0QCpGDgDgrX1GxC339sZd7FpZxbTHFvZU1YzJThZBcoKXhYQ/Bg5R1R0AIjIIeBf4Rzorll3avxo0QSqTaN0WxtavZUt+/ABjTK6wLW1zg5cAsgFn6m6LGtrvJNjnSZxb8Xz0wZt8IbyKdcGPgHPTXCtjss+4/B2M6FcDuyyA5IJks7CudW9uBD4Qkf/ifBQ/B6dLK2e0dGG1TG1P+NJYMZeC6FfYf2VODREZ02rbccNZ5juQc1+110AuSNYCaUmcuMr9avHf9FUnW7XlwkomGG35deantTbGZKsdPmfbaOvCyg3JFhLe7G4odZuq/rgH65R1WqfxWhp3YzxR23wtJySdhaWqEeDQHqpLr5E4sWJrpOmpqhiTlSL2GsgJXgbRF7qbST1J+w2lnklbrbKMeEzLEEHxg8UPk/M6ph81fZOXADIQ2AGcFHNMgdwJIO53bf2eKEJY5DAGIGLJFHNCpwFEVf+nJyrSq3QaJyyQmNwWUosguaDTACIiBcBlwH5AQctxVc2Z1egtXVidbSrV1tNlAcTktqhY2tFc4OWv/BAwDDgdmAWMov3Cwj6vZRZWW1dWggBhrxljAIjYLKyc4OUtb7yq3gDUqeqDwFnAAemtVrbRJPdijlur3Rggcbof07d4CSAh9/suEdkfKAXGpq1GWai15SHS4UiHcpb33RgAIhY/coKXWVj3iMgA4OfAdKAEuCGttco2mvRuK7ExEGMAJ+2PP9OVMGnnZRbW392bbwJ7prc62al1EN3igjGeqA2i54SEf2UR+apI4v8CEdlLRI5JT7WyS1sXVscjHcs5x60jy+S6sMWPnJCsBTIIZz/0ecA8oAJnGu944HhgOzAt7TXMBq0tkOQBQtVSmRgDELVB9JyQLJnin0TkDpwV6EfjbGfbACwDvqaq63qmipnXcRC98xkm9uIxuS1qL4GckHQMxE2mOMP9ylmt+4F0EjiaIj7ePWstR81v6IFaGZO9rAWSGzLaUykiU0XkExFZKSK7dYeJSL6IPO4+/oGIjHWPDxKRmSJS67aS0qx9F1YizeVlvFdyCA8fMy79VTImi3WWtcH0DRkLIO5eI3cCZwATgUtEZGKHYpcBlao6HvgDcLt7vBFnKvF1PVJX93ung+huS8VePCbXRX32GsgFmWyBTAZWqupqVW0GHsPZLjfWOcCD7u2ngJNFRFS1TlXfxgkk6dehCyvhLCuPad+N6fMsgOSETgOIiAwVkftE5EX3/kQRuSwF5x4JrI+5v8E9FreMqoaBKpzZYZ6JyBUiMldE5lZUVHSrotLSheVrCSAJXhzWAjEGgKi9BnKClxbIA8DLwAj3/nLgmhScO95/WMeP8F7KJKWq96jqJFWdVF5e3pUf3e2MnQ0M2tR3YxxqL4ac4OXPPFhVn8DJTtDSEoik4NwbgNEx90cBmxKVEZEATh6unSk4d5e0hI1OZ5ZYF5YxAKjfIkgu8PJXrhORQbifw0XkCJyupM9qDjBBRMaJSB5wMU6urVjTgUvd2xcAr6tm4F26QyqTznYktC4sk+vso1Ru8JJM8VqcN/K9ROQdoBznzfwzUdWwiFyJ0z3mB/6hqktF5BZgrqpOB+4DHhKRlTgtj4tbfl5E1gD9gTwRORc4TVU/+qz1iqdjCyTRi6PTQXZjcoSXLqxnb72FwrIiTrmqRyZTmjTwkkxxvogcD3wO5730E1UNdfJjnqjqC8ALHY7dGHO7Ebgwwc+OTUUdvGidnttJF5a0zjyxFkhPmvHq9ag+zWmnLs90VYzLyzrCdRuPgY1wSvqrY9LEy5a2X+9w6FARQVX/maY6ZR/tuJAw/qujrYvL9CSf7wkAIpEwfr+XRrVJO5+9CnKBl1fb4TG3C4CTgflAzgSQ3bqwEg6B2MBhJjSE86lsHEAkErIAkiWi9lLICV66sK6KvS8ipTj7pOeO3bqwkrdArAurZ/15wbdZXjmeC05tBgozXZ2cFQmH2+7YSyAndOdzQj0wIdUVyWYtCwk7HSTvbKGhSYvlleMBaA41Z7gmuU0jMbP7bSV6TvAyBvIsbe+ZPpy8VU+ks1LZqrPVtZ3tF2LSqylsASRb2ELC3OClw/i3MbfDwFpV3ZCm+mSltllYrUfilrMU1plxyoEreXXYZJqbmzJdlZwWjba1QDrfM8f0BZ1+TlDVWTFf7+Ra8HB425GwNa7Yi4e56ys562/v0hRORdKC5J4bfjyNUsiuXZVpP5dJrF0Xlr0EckKyPc9rRKQ6zleNiFT3ZCUzrjUXls+9m2gQ3bqwWlz++AKWrqnkpdXbe+ycm7dt6bFzmd05+8+5t60LKyck29K2X09WJLu1DKJ3UqqTWVq5pCHk7Mq4eeta2Htoj5xz2w5rgWSSRqKtt20/kNzg+XOCiAwRkTEtX+msVLbxulGUdtJCySVHVs/m+s3zWLVkTlqef+mOWhZvr2l3bFd1/BRta9Ys45lnfkk0Go37uEmN2BaIdePmBi/7gZwtIiuAT4FZwBrgxTTXK7u409vDPj+QrAvL/d4Tdcpyx24cSbjwGMpXpie9yFm/mcUXfvtmu2O1TfFnYb0372q2RN9kx441aamLcUQisYPoMcdj14eYPsVLC+RW4AhguaqOw1mJ/k5aa5Vlgg3OJ9cGf17Scmq5sGI4bybhzvr9uqk4UktJuBYAUffv0xw/gDyy+jR+PfcHbK+tS0tdjENjWnixs7DaDa6bPsVLAAmp6g7AJyI+VZ0JHJzmemWVlgBS7y8AEk9RtBZILCffpj+cnt/G0IOiDDqsdbd6ACIJFhJuqB0OQFPI3sjSKdEYSKK/i+n9vKwD2SUiJcCbwMMiso3WTp3c4Ktz3qDqfQVJy7V1bVkLpCWAFKcmcfNuPhq5b+ttwQkhicY4Ig3KoPw6tm7bDHvn1GefHhWOtv2tY2dhRcPp+R8wmeelBXIO0AD8EHgJWAV8IZ2VyjYF9W4AkaLWY3fO/ohDfvMsTeGYZrvfAkcLFedNwxfxp/1c0trmi9/aOeCIajaetDdzFi1Le11ymcaMdcS20kPNjZmojukBCVsgInIH8Iiqvhtz+MH0Vyn7FDSF8GmkNYAowpbX7+Gbawp4aWaEc0491z1O6+O5Tt0WCNH0ZMed9tQjSFTgxJgWhcQPIAtKPwdAYw8sasxtGucWRCK7t0BsYL1vSNYCWQH8TkTWiMjtIpKzbf9if5Qi6mkUJ9OrIoz5ZAo+TmPj26+3lrMdCduoOG8QvjQFkGBkKgE9HWhrgSSaOSo4rcSI2jTedIoNClFf21uLhuIEEBsX6RMSBhBV/ZOqHgkcj7Od7P0iskxEbhSRvXushlkgIEqhNrQ7ppIPQCimf9cWEkI0EuHG338DxPm074sE037OlgCiCXrLWh6PRO1NK52iGtOFBSzOC7PBHyEcJ0dZpMm6tfoCL7mw1qrq7ap6CPBl4DwgpzqTfQF/uwDidFE5n2Y1NoGcz1og85fM4Z5DrqGxv5+VR35INJL+DZ5ax0ASxG2f+3hexAJIOkVCsYPowktFIR7t10w0TndV7LhI2AbZey0vCwmDIvIFEXkYZwHhcuD8tNcsi/iDAQqi7T9Ftaw98EVjwoVtKMXHSxcA8NLpe/LomBPYtU8+73zwfJrP6rZAfPH/nVtbKIH0D+jnsrqGtnU2sZmp483Cim2BhBtsfU46zZ17By++dDL19TtT/tzJkimeKiL/ADYAVwAvAHup6kWq+p+U1ySL+X0+8qLtP0XNOX0Lt140kHDMfEVL5w47Kp3kiU3uosu6fkHenTczpefYuLMtZUk0Gm1tYWiC/EutLRTLz5RW1VW7Wm+rCKfus4ITh68mHGcMpDmmBdLcUN8j9ctVVVUfEwisJT+/f8qfO1kL5GfAe8C+qvoFVX1YVXPyo4L4IC9mPYMivDlkfwAihUXtjsd+z0V1jY388oFlFIWdN4iGwgChpvyUnuPN19sy6VTW1NGa7NLfWQskd/8uPaGqprb1torw7B4n8OKBx7Luw4Xc+Z3Xeeu+O1sfDzXW8UJRM28WhAg35nYAUVXef/9JGtLUEguFthMOF+H3p747Odkg+omqeq+qpr7d08v4/EIwdqwDCLiJ46KBtkFiS+cOGqqnsXAohe7UzYYCP+FIaqdsrlm2oPX2ps1rWsN1NFEXlmrSx3PJyy//khdf/GVanrtqZ1s25NjW+Lo5mwBY+VbbOGJTQx1L8yJ8UBAm1NB+gkquWbtuDnX103j11e+l5fkj0Uqi0fQkV7dXlAfiF/I6vAn63VxPoYK2/FhtvVm5+0nXH3J+T/nuwGltQYDi6ppkP9Jl4dq2T6wb1qxu27PeWiCdmr5mK8+uq0jLc9fUtW0TtK2k7Q1LfLE7YjtiB9EbG1L7/9HbaNT5Xfj9Sz2Vf/bmm3njrj95KltbW4XIJiD13VdgAcQTn9/XrgUCgt9tgYTyg9z+txNpaqjpdMOpXOALuVNm3U/7q/oNhSZvLZC3XnyJHZs3d36OmOR8G9d+QmsXVoIxjpYxkkiCAJNLnv/0dF749NS0pLaviwnsc4oOanugdT+EtkkM4cZ6Sk/1s9eJ1TTW5XYACYUb2fbKgUQbvG3JvG7zsSxddICnsjPfuJL8/HrQ9MyGtFeUBz6/kBfp0IXlpgNrLgiwYtgX+PsN30Ql2Pp4b/HqG39gzdr5qXvCJqgt30Kj3/ldrAyOp3HgoE5/bPXKT7iwYBjXzO58fkZDqLD1du32bW3rPHw+Hvj3z3nmpd93+Ann8bAFkFY7K9en/Dkb6p2uqCHRrYSlrWs3GnFbgLQFkO1r1rDVN4ylefvS3Ng7hlaXLPkvCxY8ltLnDIXqWfnK6+zcdRUVz30xpc8NQNN6KtaVUxMuT/1zYwHEk8bqEMF2c9kFvzuNtykY4LnSk4hEhvPJoM7fKLNJKNTE5S/tzc/+8wxz5977mZ/vmWcfp55i/nDSRN4tOaz1eFNR8iSUAC/OeQOAJSV7dFo2EFIaSnbSWLKT5qam1vZexCfMXD2URxYUtysvrQEkt6fxhsNt62CWrFqS8udvSVlSFmm/43W4yfnwpfh58kc3snzWiyyZ3dZCaa6vJRMi4TAPfesW3n3wb57Kb912LTsr/5cZr34/ZXV4Y9YhbN2ynN+UNbAsOI6q6uQt8HCTt1ZKi4rp5/DigpuZ9t75VNalfh2UBRAPhu07jGBsqmqkdRC9Md9pGlaOGciHBfu7JXpHF9amHZuYtnU2I1aEqaq+7TM/34Z1Kyko2L2p3FCUeB+VOx+7m//3yJ+5ddgUAAaG4+8qCNDYWM2atTPo16+Q3541nt+cNR7CkXYtkJcPPZa3j5jS7udaBtFzuQXy3nu3M+vNfflpzSKmVS/jyXdWpfwcUff3XBpqHxDCUefvr9KfbXUnMOuBamI7hJszNAtr57oVVPuPYeE7e3oq37C5iPr1/fD5Xuq87K6dzHv6X52W8/nCbA4N4ugpa9kwoolFC/+atPyuzWtbbz/9kxuTll396Qu8X7ofswrDlOcrA4qT72fUHbn7iuqCY775PYIxifiapACfuxK9MRjkqjnTqS5r61ZJ9RjIhx8+ziszvM3QaG5s4vFp04hGkvdxh8PNzJj+MJH84xkXPo3t9QO6Xb+G6gqa6ipprNpOU9mw1uP91AkGDQWJ05ncOnQKdww/rvW+JMhXFY1GOefWx/jmI+v4zdQLW49vCw6ibQwk/r+zj/bjMqlQUfEJoVDvScfx+sfz+cGrt3PveQdx04WnUFGRurrv3PkpL750Ii17h/Vvbj+rqp5+rA1EaPY7A+vN+eU0+drOX+ex79+L7dtX8fwLF1BTs63Tsls+dgat1df5+EBDQw33zrmNny/5P1ZWjqW5OXnQe+JHd/H+jBGsmfNWwjJRd1y1fuAoXis7irVHl1Nbn7wFUrnuU5qD9YT9jWypPiFhueXLn+XTT6/i7QKnVdjU4GP6f1Pb/QYZDiAiMlVEPhGRlSIyLc7j+SLyuPv4ByIyNuax693jn4jI6emuayBmDKReilkfGAnAzEGT+MvhZ7NgROq3ia+t3cqLL53J3I/vYd66JnbuXJO0/MGfho8AACAASURBVPwlc7n7lp/wymGjuPvHV8YvM+sh5r/5EE89+2XWrVzPC2evZ8WUxdz2wY86rc/lT/2Zw159jlsf/0u746fNnslZ7zzP+vxSqsvafk+l0WpEo9Tlx18H8snSRbsdq/UXxykJby5aTGCA8vmV7bPo/OuE81vDdSimi2rLxx/y6exZ/OPZJ2iSvN0e/yyam+v5cPGZzHj1ws4LA+FIlMaQt4kEK1fO5NnnTmH+gn9+liru5rllp/OVivfY5Hf+bz83PHUBZNGiO8jLW9eazLJfY/vn3jhiCOtOWcWi0W2zvyKFbQPnm3YlbnV2xauvfpO33/8a72wu5v3Zv+m0fMWa9Wyb+A7h0jlJy82efQfvvncwiwqgyq/8as61VFauTfozzdHxAMy65yVU44+Kzpv3J6Z/cjKNxU4AWxXYkzc+Sdwaaqyv58PFH3LPeUX86oIRbD3gw8TXtv1DttaVM7loC2ft9TFfHD6fEf2HJK1zd6Q/UVECIuIH7gROxVntPkdEpqvqRzHFLgMqVXW8iFwM3A5cJCITgYuB/YARwKsisreqpi1fd7DDbnb14rzRNbgp3j/Ka9vgKFkLZN3SRcx+9GHOv/V2pJOV62+/93te33gUH44cy6Ly/Zn4wX8554yr45atrq7izIoAnPo/ALx31r4E7/glX7v8GnZV7uTpmf9l5cdbCZf4CUaV0ui+bD9oIPMKD2LeWDiifF6cui7mg0f+SfXYkbwxOMiMQU5L4c4ho7h0zWr+/tw/EfWzh6+McMCPv7yIh0ef0Prz+RqiiDrqg21N58VLFvDq8rns3FHBjLETITC29bHB0W3U+pzfZ21zmAK/j807dnDXvx4nL7SFo2QwL124e1LolpdnTX7bWMsf3n6O+sIg60oHscnv5P4Mp6gFsn3nGn7/zo/YZ8g8zvRQ/pLn72Vp8XiuXf8Rl3/jqqRlX57/Kn/96LscsHoFjx3S+XPvqqqkpKiEQDBxK+/mX/yQC4Pj+ejUsa3H1oxMPF5XW1vN9a88xF55JbwU8DG6qYY7p36TvPz2Y1nhcBObN88jpCU8+cm57Nvo/CX8HVqR6w+t543Soyg5roaJjzcQkABrD9kC7APApqrkn+ZVleZIlPwkqWiqqjbwxqelVHEOT44/ncYlz3Lqye3LLH7hKfKLSxh90GSmTb+XtWMG837xiXAAbLnnVn51xQ2tZSPhMKGGWgr6lfHqsjk8suxXfH7wYt4/eD/KQ5XM+aSZzw/dl3hm/O7XVA0r49kTYWNgKrtuvJqrb/1z6+PbVn5E2cix7Kz6K1p0BR8MH9f62Ko92w92R8JhPv3gDfKKivnR+sW8NfHU1sfumXgCV65cypDx+7Ue27LlE3ZuXM67H8xmhf9iVh/Vn+VBZzyyfPaTTDrxpCS/6a6TRNEx3UTkSOAmVScnt4hcD6Cqv4op87Jb5j0RCQBbgHJgWmzZ2HLJzjlp0iSdO3dut+r7v3fdwn37nA3AsOhmtviGJy1/wabX0Jb1ByJEfD5UfER8fhqCAQLRCMFIhPxQiLDfT1FjM3WFeTQGglTnF1IdKGJk3Q5eHXB063NOrp/PiOpdRH1CyO+noqg/JaEmKgr6U+nv3/rpssXIyAYO3PUpZbX1LBy+BxHxsTy4N/s2f8yyvH12q/MXts2ipKGRiv79qMovJD8SZntBfzYERlAr7RcijY6sY8+6zdTkFTC/4KDdngtg79Byqvz92eobxgFNS1GEJfkTE/7OWur1hYpZlNbWM3fkntT6i9jgH00/raJGSuP+XLHWUCf92CO8hrUxAamjz4WWc/CWNYg64yKi4FPFpz5UYXtZEcVNjUR8PvyRKH63TFPQT9jvI785THVhAQXhMDPKJzMyvIkzFy1g++B+NOUFCPuEiPiI+IWoL0BdMI8xO6p4ZKzzor9g0wxKqxuoLC1u151WUh+iuD7EypEDqc4vZFNeOYdvX8GodVvZOmIQKoKKk74z6vORF47SFAxQEGpmY+kAVIR9Nm6hvjDoTlUWUMWnSkldI9sGlvJc+XHtfhfl0W2csn4BUb+fuvwAAhQ0h9lYWsbC4n1bPyC1OLXybUqammkIBsmLhCluCFFU38zm8v6s7V/OiLodjF9TzV2HnMkFm17hqRGnxf0b7BFeS/9ILYvz2970pu54i5LGtm6ssN+HilDYFKKqqJD8cJgdxcUUhEP0a2zEH4miPh9RgaZgkNq8PPyqvF56ZLtznbfldQobm2nMDxIVIRBVwj4fFSUlvFNyeLuyRVrHKRWz6V/XxI7+xdTkF6BATbCQ0dU7eG3w5NbtHABOrH6P8qoawgE/TcEAIb+f2mABAY2wuGQClTKwtexe4ZUcvG0twUiEtQMHEhUfdf4CBjbX8Ga/tvG6CaEVrAhO4OxtM2n2B9hRWMLOvH4Maq5mdtGhcX+fJ1S/z4idu9gyoD9bCgeQFw2zLTiQsQ2bebdkUmu5c7bO5NgQfPVrP4z7PJ0RkXmqOmm34xkMIBcAU1X1cvf+14ApqnplTJklbpkN7v1VwBTgJuB9Vf2Xe/w+4EVVfSrOea7AyeXFmDFjDlu7NnnTM5Ff/Okm7j7Q2TjquJoPeLPfFESjqGTfMNKgaAX71a1u98+ZzBc3v8Z/hp1AVDrv4plcv4CPC8dRLWUJyxRoA41SyH7Ny8iPNrcGmGKtpU5K2pUt053scl9sZ21/m+cHH+OpzlPq5/OB+6LK10aapPOZXj0p3v/GyMgGNvpHpfg8EdTD3y2ozYTcrryjaucyt/gAmqXzFDN7h5ZT6yve7cNJMl/Z+hoPD3U+/p+/aSZPjzgRgAKtpzFmV890GqA7272JJ3PO1plMH3K8p9fyXuFVrPOPbv1dxvJpZLfXkF/DRKR9R0+R1u0WoIPazNmb3+LpEc7vrbP3lnO2vsFzQ47Z7bnjOaXyXQ5dv5Vrr76507KJJAogGevCIv5UpY7RLFEZLz/rHFS9B7gHnBZIVyoYyxeTTHFoVQ30gwnhVSwPTqBQ61u7sgCunvUAPndEUSTqfKLVKP5IMdHQLlTyiUqIvOJCwuEwYRE0L0igsZFoQYCCyjqqhpXTXJzHusEDWVsyhFOXzKOsVqiLhNCiAHk7dlI/vJz8hkZqBpSxsbyU/FCE/ww7kbJoNVcEyzjk3YeJBPxEA0FqSovot6uWUH2I1QeMo7S+iblDxzGsaQfXHXU64+7/m7PQy+8nEFaaC0oI1O9g++gxlO6ER6YcwjbfUCYvW8ExjR9RVVZCJOCHaJRPxoxg4rp67tvX+ZQ7OryRFcHxBKIRrmouZOHCh9kxdAAD125h4WH7sb5wCJ8GnGb7yNBWduU5L/RvDh/J+Jn/RAJBmvLyCBUEUZ8wu3AcxYObGLtjJ4+NcT7NlzY1gPsrb5IC9g4tZ3mwbZuafZs/5uQF87lv8rk0SBHDo5v40tuvEhV1ug7dT/T4fAQiigaL8NdU01hSQrShkaBfUF8U/D6iDWEk6EeDfkT8VPsLCQ3I47G9pgLw3fceIz/UTH5UCIZ9+NRHU30Vof4FPH3sCa2too3+UeRrI9955ymkoR7x+xGFuvIBhFH6b68n4i9g88B+PHmA89znbX6NUcs3kd9UC/5Con5o8Ofh80eRgiLyahppLMrnr4efT7HW8q1ZT5BfkE9jkyJBH7UD+lNQp8j2Hfxl6jcA2G/tBo7Ztpj6/HyC2oQ0hQgEhZ1DhtCUF6CoromPx4xgVv8jmFC1jYPe2IPQoAcJRxqRUDOR/CKay/rTUFJA0eZdzGoYxue1jpnHjyEqwgXjj+CRaucNcOjmCo7qPxe/RjhkwSaqh+fRb3sldx9xAc1SwGUf/Zfiip20vaSFQLiBiE8I55VAWIkW5TEoWEzz9q005eURDgoSBQIBfA31RAJC04D+bB9Uyn+Gncylc/7DsXvswYcfP48GBX9dCCJRwkEfwQEDqMrPo7i2jr9Mdsawvj58L8bO+idhfGheEIlGEZ+ffj5he1E+6vMTQVg1qpwRc9ZwYcE8mkN1RPGD34f6hLKQj0Z/lKp+BdT0K+TTIUOYtOAj9t1zIss+/QhQonl5hPPyKNi2g7rCZhpH7kVNcQH/GXoiIcnjt+d+h5G/uxH1QaChiWhREaGiIA1l/XhvxDiOXLOK+/d2Ok33qa5nz48fJioQzstHw1EaBpfQr3I7O4cMIxT0kxdW3h+9F2NXb+Ha627p7ltfUpkMIBuA0TH3RwGbEpTZ4HZhleJsbuXlZ1PKFzOrqaza6bNt8BXw3YX/xd/cSNXgoTy05wkAnHfeN9jnoNRs4PjXp05n/NbBfP/SWxg0aK+kZX/9tz9xeP8FHL9xB6dccR2nnHFR62P/fPIOxn3uIKZMPoEr7vwzm6r9XPLGm7xfUs4eU/fhxzffkfS5BzxwJw2hj7jmuv8DoGpXJQtmz2LSESfxoxfvZYmvjPM3zWJDWT8a3Ey8QQ1zxjlf5YxzvrTb8/32qbv57aApHFS7g6XuB8XJh5/F0cec167cJ1t3Unf/vZTuamDcsNGMem0WNUOEop1beeX4tnJjays4as0aqgf7aC4IsH+VcM20X9P0l19w7/7nEcHP9b/ouMCwe55/eSL1kUIaNuQT3NXMDdP+D1+C8ZXz1yznj2/+l0DUz6ryARy7dAk/vuF3cctGIhEe/fc5DCPK4pkj2OPAKtZuKuWum5JPcAg3NbH16b8wev0Wpt38590er6paz7Lln/IXd3Zt3tYarr01eSqM1etW88x7z1O+/TU0UMYl199JfuHuExwqK1dz9IJT+ftHl/DhW0M5YY+FHHnm/5D3+vs0UYAvHOHY58ehQL+mAZRuHsz37jqeuc/dz3slk/A1ws9u8paWozOP/Oer5C1bzZlHDuf4Y/+HzycoN//pf7Fm+UrO2PEmYzZXcPRVt3L08Wcnfe59fvY0jesKePricg47OPmCvw9feIplry7iot87/2/nxSmzdeuHLPrwZt54a28Gl2/j9BXzKDzxYH52Y/z/DYDGml28+94r5EVDfO/Sa8nPj9/qbmraBQgH3PAWzSuEQOdLq7otkwFkDjBBRMYBG3EGxb/cocx04FKcrMAXAK+rqorIdOAREfk9ziD6BGB2OisrMft+FO9yXom7fP0p91VSNqScvwXaZmHttU/8wbXuOHDY0dTUzOs0eAD85DtX85MEj339wrZZWQ9cO42Vaxdzyl3rOHL4B56ydH77G+0XT5WWDeCE05wuvSOai3h6jI+Rr2/hgF0VPPdFZwZKx8HUWNdd8G2+W1HJE7PbFrQF83fvFvjc0IH8cdpPAWf20z9f+CuF2w+lufATvr7iBf45wflEJqrcdm2cq29wWo5NHrprvJpy2H/YunUBeb7/h38w+Hw/TVh2r7F785exP+buad9j7w9DfPf3iRds+v1+vnrBcyxc+Cirml7k4bkXcd2UzvMjBfLzufPL1yV8vLR0NJMOGcqVf/obwYY8dPz4Tp9zzzF7ct2Yq4Dkg/4DBuzJgQc8y+XczaTB9zChn/M3z6PJCSChCHlu6+Ls6/dly/KPEfETeaOEs/efSb+a1HUBT97/G5Tl/y8HH/jvpOUOPf+rHAp0Zd33zUf/Cr9E2GfCq52WPfDMCzjwzAuSlhk69EBOO/Vp3l7/J2pnhPiwuPP/z4J+ZXz1UaUgsJL8qYm7bPPzne7lC0c9z+rq4Ry/T/qSVWYsgKhqWESuBF4G/MA/VHWpiNwCzFXV6cB9wEMishKn5XGx+7NLReQJ4CMgDHw/nTOwACQmF1awxllpO655Pd+7+o8AXDtzYdvjCaatdscxxyRfLNRdo4ftyS1TvsuoAVM/83P9z9e+y6XRKCMKnWmFR9Y5ExWCkeR/kuLyARw8fG/wOIszL6+IsM/J+BoJKsHmthT7tXnxX1A+91+8idT9TQYP3pvBg/dGpJghQ7x9WPj2bckXiMU66KCLWb/hQY475Rr6F1/R3Wq2EwjksXB9A5X+AI9+O/4Ad3eVl0+kX8kEDuY5mpqcxbR52gwCEV/bmMDQvQ9g6N5ODqdTKpUPFk5kwtDXU1aP8eNPYfz4U1L2fLEKo8OACMXFqc02MXViHhXNM5g67jXg2k7Lf+v+izot0+KovRZxUuErBIJf+ww1TC6TLRBU9QWcjapij90Yc7sRiDvZXlV/CaQnL3Ucvpjkc6VDRnDtm3PJ314Kn/39NyPy84v52rlvdTqV2Cufz9c68NeSeLJ9Asr4Djl0Mv/38L3siDYBnXf71Q5soLAuSrTYR6C57fkr8+KvHxk2xJktFzuDJlUmTkzUSfLZiAjFxQegugJ/IHWDznPdTAmDB3lbed0VAwZ8jppaKHFT2OxTt4Z3+pXjI8jx5+4i3Nx+U6miaJSjmoL4/b0ja8PZX3gu4XqOz6KwoIgL9n6WSCTxNOxuEyd1ybixZ6T+uV0ZDSC9icQEkEu+9UMeveqP1Je9R0t8u/GDu7llyrczVLvuSVXwaPGzxU9Qp8obn3PeqIZWectx9M2vfMvzORZO2YPAxsXsKWWEYj7d7vLHT1d98nEncPPaRo6pmYOXAJUtjjn6FubOHcZhh6buf+rrEx+lLL8Kn++slD1niwkTTiWY9yhjRjsTdcpeG8aRI+awT79i9p8ar7PIfT1JaveKSadUv14AAn6n5RwKpT6P3sR9f8+OncsYNcrbbMzuyL45qNkqZhC9uF8/Ln/gBn7wx/taj31v2l2ZqFVWuerq25h2ze1szXNeDMPXpr7v9fDywxgW2MEll07DP3gCwyJO6ofJqz6OW37PMeMpe2U1NfN7V8rwvLxCjjrqRyndRe6EUSM4Ykx6ungAxu4xuXUywSFNRZzy6V6U5if41O72OPuDqW8Z9ibhiLP+xe+fkPLnHjfuOCYdlt4PtdYC8UiiveeTUqadt2AB7+23nf6f7r4X9md19UnHwEnOWpFvn3QsT975AQXhjcwrGBq3fCAQ4N4vDmDCHkekvC69zUkn/brHziXuwHmi2WnitkD6lyZfkNvXHbD/xbz19iqOO7bz8Y9sZC0Qj1K//U7fNWIxfP6/owikeWeUUcNH88YPjqcsP8SPT9t9ZX2LYw8/kWFDRqS1LiY+ny/RIkf3FZWGja16k0AgjxNPuJmiovhZFrKdtUA88mfhivPs5QwIhoPpz1Y7bEg5C28+N+3nMd0jCRNYOoEj7DHJpMlO9q7oUV6CaaImDnV+V+H81G9gY3qXROPOLV1YEQsgvZoFEI8ChT2Tw6cvEJwA0lyY290TBkiwT31LC0TD9j/Sm1kA8aiof+LkgaY9dVOyRwvszSFnubOsiooGJyjgPB4Np3X9r0kzCyAeDR7e+SDsXuHUbxPaG0V9ztRM7SWLxEzqtewsGciLv43q4V8eSWHDYo74xqU9WS2TYjaI7tHg8pHQybq4Lz1Tgvg3w4m9Z8FaOvgjnxDiCJoCqd+D2fQOos5uOJJgGu8Bp3+RA9K+j6hJNwsgHg0dPgZW7EhaJqBBJzNXjptxciH1jas4d0LqVzyb3sJpgfhTtI2wyU7WheXRkEGp30+4r2qKHEz1imbOPfK4zgubPqmlC0ts+nufZi0Qj4rz05DsrI966tLjqQkfQ1me/c5yVksA8VsA6cssgHiUKCWD2Z3f56Msz35fuUzcLAT2uunb7K+bQrV5jxOOzMh0NYzJPLcF4kthMkiTfeyvm0I//cVPoDqtO+sa0yu0rDRPRwp0kz0sgKTSoL2cL2Nynbv5Ujo2YTLZw7qwuqBA07e3sDF9ibQmS7R8aH2ZtUC64JrpmyEQzvmFgsZ0zgkg0bAtjOrLLIB0QbDR8mEZ44nbdRUOpX5TMZM9rAvLGJNyQTYCUDxgQIZrYtLJAogxJuXOu+0S9j9kCROOOS3TVTFpZAHEGJNypcNGcvy3f5Dpapg0swBijDGmWyyAGGOM6RabhdUFm/MeoSgaBU7KdFWMMSbjLIB0wehLr2drXVOmq2GMMVnBAkgXfPswS1NijDEtbAzEGGNMt1gAMcYY0y0ZCSAiMlBEZojICvd73OWqInKpW2aFiFwac/yXIrJeRGp7rtbGGGNiZaoFMg14TVUnAK+599sRkYHAL4ApwGTgFzGB5ln3mDHGmAzJVAA5B3jQvf0gcG6cMqcDM1R1p6pWAjOAqQCq+r6qbu6RmhpjjIkrUwFkaEsAcL8PiVNmJLA+5v4G91iXiMgVIjJXROZWVFR0q7LGGGN2l7ZpvCLyKjAszkP/6/Up4hzr8vZmqnoPcA/ApEmTbHs0Y4xJkbQFEFU9JdFjIrJVRIar6mYRGQ5si1NsA3BCzP1RwBspraQxxphuk0zsWSwivwF2qOptIjINGKiqP+lQZiAwDzjUPTQfOExVd8aUqVXVki6ctwJY281qDwa2d/Nns1lfvK6+eE1g19Xb9KXr2kNVyzsezFQAGQQ8AYwB1gEXqupOEZkEfEdVL3fLfRP4mftjv1TV+93jvwa+DIwANgF/V9Wb0lznuao6KZ3nyIS+eF198ZrArqu36avXFSsjAaQ36qv/DH3xuvriNYFdV2/TV68rlq1EN8YY0y0WQLy7J9MVSJO+eF198ZrArqu36avX1cq6sIwxxnSLtUCMMcZ0iwUQY4wx3WIBpBMiMlVEPhGRle6alV5DRP4hIttEZEnMsbiZkMXxZ/c6PxSRQxM/c2aJyGgRmSkiy0RkqYhc7R7v1dcmIgUiMltEFrnXdbN7fJyIfOBe1+Mikucez3fvr3QfH5vJ+icjIn4RWSAiz7n3+8I1rRGRxSKyUETmusd69f9gV1kASUJE/MCdwBnAROASEZmY2Vp1yQO4CShjJMqEfAYwwf26Arirh+rYHWHgR6q6L3AE8H3379Lbr60JOElVDwIOBqaKyBHA7cAf3OuqBC5zy18GVKrqeOAPbrlsdTWwLOZ+X7gmgBNV9eCY6bq9/X+wa1TVvhJ8AUcCL8fcvx64PtP16uI1jAWWxNz/BBju3h4OfOLevhu4JF65bP8C/guc2peuDSjCyb4wBWc1c8A93vo/CbwMHOneDrjlJNN1j3Mto3DeTE8CnsPJc9err8mt3xpgcIdjfeZ/0MuXtUCSS0lG4CyTKBNyr7xWt4vjEOAD+sC1uV09C3Hyw80AVgG7VDXsFomte+t1uY9XAYN6tsae/BH4CRB17w+i918TOMldXxGReSJyhXus1/8PdkXakin2ESnJCNxL9LprFZES4GngGlWtFol3CU7ROMey8tpUNQIcLCJlwL+BfeMVc79n/XWJyOeBbao6T0ROaDkcp2ivuaYYR6vqJhEZAswQkY+TlO1N1+WZtUCS2wCMjrk/Cif3Vm+21c2ATIdMyL3qWkUkiBM8HlbVZ9zDfeLaAFR1F0726SOAMhFp+bAXW/fW63IfLwV2kl2OBs4WkTXAYzjdWH+kd18TAKq6yf2+DSfYT6YP/Q96YQEkuTnABHfGSB5wMTA9w3X6rKYDLfvLX4ozftBy/OvubJEjgCrN0l0fxWlq3AcsU9XfxzzUq69NRMrdlgciUgicgjPwPBO4wC3W8bparvcC4HV1O9izhaper6qjVHUszuvndVX9Cr34mgBEpFhE+rXcBk4DltDL/we7LNODMNn+BZwJLMfpi/7fTNeni3V/FNgMhHA+AV2G05/8GrDC/T7QLSs4M85WAYuBSZmuf5LrOgan+f8hsND9OrO3XxtwILDAva4lwI3u8T2B2cBK4Ekg3z1e4N5f6T6+Z6avoZPrOwF4ri9ck1v/Re7X0pb3ht7+P9jVL0tlYowxplusC8sYY0y3WAAxxhjTLRZAjDHGdIsFEGOMMd1iAcQYY0y3WAAxvZKIRNwsqEtE5NmW9RO9gYhcIyJFXS0nIi+k4jpF5CYRua6LP/OGiPTp/b1N11kAMb1VgzpZUPfHWan8/VQ8qZuBOd2uwUmW2KVyqnqmOivUjckKFkBMX/AeMYnpROTHIjLH3XehZU+NsSLysYg86B5/quXTvbuvw40i8jZwoYjsJSIvuUny3hKRfdxyF7otnkUi8qZ7zC8iv4k537fd4ye4n9qfcs/7sLsK+QfACGCmiMx0y94lInOl/R4g8cqtEZHB7u1r3bosEZFrYq5xmYjc6z7XK+6K9oTcOt4uzj4ky0XkWPd4oYg85l7T40BhzM+cJiLvich8EXlSREpEpFScfXM+55Z5VES+9Zn+qib7ZXolo33ZV3e+gFr3ux9n5fJU9/5pwD04K399OOnDj8NJa684CfAA/gFc595eA/wk5rlfAya4t6fgpNMAZwXxSPd2mfv9CuDn7u18YC4wDmfVdRVOziMfTpA7JuZ8g2PONzDmWt4ADkxQbg0wGDjMrUsxUIKzEvoQ9xrDwMFu+SeAr8b53d0Uc+1vAL9zb58JvOrevhb4h3v7QPd5J7nnfxModh/7KW0r5k91r/Ni4KVM/4/YV/q/rAVieqtCcdKe7wAG4qQ+ByeAnIaTEmQ+sA/OJj4A61X1Hff2v3BSorR4HFoz/B4FPOk+/904+zoAvAM84H6ybunqOg0nx9FCnJTyg2LON1tVN6hqFCfdytgE1/IlEZnv1nk/nM3LkjkG+Leq1qlqLfAMcKz72KequtC9PS/JOWO1JKOMLX8czu8IVf0QJ70KOMkdJwLvuNd8KbCHW24GTmC7E7jcw3lNL2fp3E1v1aCqB4tIKU4r4/vAn3FaHr9S1btjC4uzb0jHvD2x9+vc7z6cvSoO7nhCVf2OiEwBzgIWisjB7vmuUtWXO5zvBJwdBltEiPN6E5FxwHXA4apaKSIP4OSDSiZh3vo450zahdXhZzrWMV6eIwFmqOoluz0g4sNJP9+AE9Q3eDi36cWsBWJ6NVWtAn4AXCdOiveXgW+6LQlEZKQ4+zUAjBGRI93blwBvx3m+auBTEbnQ/XkRkYPc23up6geqeiPOTnmj3fN91z03IrK3m501mRqgwPfgygAAASxJREFUn3u7P07wqhKRoThbn8YrF+tN4FwRKXLPdR7wVifn7Ko3ga8AiMj+ON1YAO8DR4vIePexIhHZ233shzjZgy8B/tHyOzF9l7VATK+nqgtEZBFwsao+JCL7Au+Js8FULfBVnE/Xy4BLReRunGypifal/gpwl4j8HAji7GOxCPiNiEzA+RT+mnvsQ5xun/ninLACOLeTKt8DvCgim1X1RBFZgDOOsRqnmyxuuZjrne+2VGa7h/7u/g7GdnLerrgLuF9EWjIez3bPXSEi3wAeFZF8t+zP3d/15cBkVa1xJxn8HPhFCutksoxl4zU5wX1zfU6dab/GmBSwLixjjDHdYi0QY4wx3WItEGOMMd1iAcQYY0y3WAAxxhjTLRZAjDHGdIsFEGOMMd3y/wHT8o4wfMqoQAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "for i in range(50):\n",
    "    plt.plot(np.arange(len(soap_raw[0])), data.iloc[i])\n",
    "plt.xlabel('Representation Index')\n",
    "plt.ylabel('Value (arbitrary units)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Splitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "shuffled_names = list(entry_names)\n",
    "np.random.shuffle(shuffled_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "training_fraction = 0.8\n",
    "n_train = int(len(df) * training_fraction)\n",
    "\n",
    "training_names = shuffled_names[:n_train]\n",
    "testing_names = shuffled_names[n_train:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2012, 567)\n",
      "(2012,)\n"
     ]
    }
   ],
   "source": [
    "training_inputs = data.loc[training_names].values\n",
    "print(training_inputs.shape)\n",
    "training_outputs = df.loc[training_names]['Formation Energy']\n",
    "print(training_outputs.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(503, 567)\n",
      "(503,)\n"
     ]
    }
   ],
   "source": [
    "testing_inputs = data.loc[testing_names].values\n",
    "print(testing_inputs.shape)\n",
    "testing_outputs = df.loc[testing_names]['Formation Energy']\n",
    "print(testing_outputs.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'rest_framework'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-50fbbd75fc36>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mapi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrest\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mQueryEngine\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_KRR\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'CdTe'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/dev/materialsweb/api/rest.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0menviron\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msetdefault\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"DJANGO_SETTINGS_MODULE\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"materialsweb2.settings\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mdjango\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msetup\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpymysql\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mdba\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0msimulation\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmaterials\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mentry\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mEntry\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/django/__init__.py\u001b[0m in \u001b[0;36msetup\u001b[0;34m(set_prefix)\u001b[0m\n\u001b[1;32m     22\u001b[0m             \u001b[0;34m'/'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFORCE_SCRIPT_NAME\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFORCE_SCRIPT_NAME\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m         )\n\u001b[0;32m---> 24\u001b[0;31m     \u001b[0mapps\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpopulate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mINSTALLED_APPS\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/django/apps/registry.py\u001b[0m in \u001b[0;36mpopulate\u001b[0;34m(self, installed_apps)\u001b[0m\n\u001b[1;32m     89\u001b[0m                     \u001b[0mapp_config\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mentry\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     90\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 91\u001b[0;31m                     \u001b[0mapp_config\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAppConfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcreate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mentry\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     92\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mapp_config\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlabel\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapp_configs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     93\u001b[0m                     raise ImproperlyConfigured(\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/django/apps/config.py\u001b[0m in \u001b[0;36mcreate\u001b[0;34m(cls, entry)\u001b[0m\n\u001b[1;32m     88\u001b[0m             \u001b[0;31m# which may specify an app config class with default_app_config.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m             \u001b[0;31m# Otherwise, entry is a path to an app config class or an error.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 90\u001b[0;31m             \u001b[0mmodule\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mimport_module\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mentry\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     91\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     92\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mImportError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/importlib/__init__.py\u001b[0m in \u001b[0;36mimport_module\u001b[0;34m(name, package)\u001b[0m\n\u001b[1;32m    125\u001b[0m                 \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    126\u001b[0m             \u001b[0mlevel\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 127\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_bootstrap\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_gcd_import\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mlevel\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpackage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    128\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    129\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'rest_framework'"
     ]
    }
   ],
   "source": [
    "from api.rest import QueryEngine as q\n",
    "\n",
    "model=q.get_KRR('CdTe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_predictions = model.predict(training_inputs)\n",
    "plt.scatter(testing_predictions, training_outputs, alpha=0.2)\n",
    "plt.plot([-0.4, 0.4], [-0.4, 0.4], 'k--')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(testing_predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's do ML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SVR().fit(training_inputs, training_outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_predictions = model.predict(testing_inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sqrt(np.mean(np.subtract(testing_predictions, testing_outputs)**2)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(testing_predictions, testing_outputs, alpha=0.2)\n",
    "plt.plot([-0.4, 4], [-0.4, 4], 'k--')\n",
    "plt.axis('equal')\n",
    "plt.show()\n",
    "import pickle\n",
    "pickle.dump(model,open('Li-Ge.sav','wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalizing Representations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often, it is helpful to normalize the inputs for machine learning. Earlier, we observed that the ranges across different features in the representation can vary dramatically. Here, we normalize all 576 features independently such that they are centered at 0 and have a variance of 1. This is not always necessary, depending on the representation chosen. It does not appear to help with SOAP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "scaler = StandardScaler().fit(training_inputs)\n",
    "training_inputs = scaler.transform(training_inputs)\n",
    "testing_inputs = scaler.transform(testing_inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "for i in range(50):\n",
    "    plt.plot(np.arange(len(soap_raw[0])), training_inputs[i])\n",
    "plt.xlabel('Representation Index')\n",
    "plt.ylabel('Value (arbitrary units)')\n",
    "plt.show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SVR().fit(training_inputs, training_outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing_predictions = model.predict(testing_inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sqrt(np.mean(np.subtract(testing_predictions, testing_outputs)**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(testing_predictions, testing_outputs, alpha=0.2)\n",
    "plt.plot([-0.4, 0.4], [-0.4, 0.4], 'k--')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "pickle.dump(model,open('Cd-Te.sav','wb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random search to select best hyperparameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Almost every machine learning algorithm is associated with hyperparameters -- parameters that are not optimized by the algorithm itself, like weights. Examples include the choice of and options associated with representations, depth of the neural network, and regularization strength in various solvers. Here, we optimize Gamma and C, which are hyperparameters of support vector regression.\n",
    "\n",
    "The C parameter trades off misclassification of training examples against simplicity of the model.\n",
    "\n",
    "Gamma defines an overall scale factor for the support vector machine's notion of distance between two points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cScale = 5  # scale for exponential distribution\n",
    "gScale = 0.001  # scale for exponential distribution\n",
    "n_iter = 5  # number of iterations for random search. Realistically ~100\n",
    "    \n",
    "param_dist = {'C': scipy.stats.expon(scale=cScale), \n",
    "              'gamma': scipy.stats.expon(scale=gScale), \n",
    "              'kernel': ['rbf']}\n",
    "\n",
    "# using random search cross validation to select best hyperparameters,\n",
    "# performing 5-fold cross validation\n",
    "cv = CV(SVR(), \n",
    "        param_distributions=param_dist, \n",
    "        cv=5, \n",
    "        scoring='neg_mean_squared_error', \n",
    "        n_iter=n_iter, \n",
    "        n_jobs=n_cores, \n",
    "        refit=False) \n",
    "\n",
    "search=cv.fit(training_inputs, training_outputs) #trains the model using training set alone\n",
    "\n",
    "C = search.best_params_['C']\n",
    "gamma = search.best_params_['gamma']\n",
    "print(search.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SVR(C=C, gamma=gamma)\n",
    "model.fit(training_inputs, training_outputs)\n",
    "testing_predictions = model.predict(testing_inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sqrt(np.mean(np.subtract(testing_predictions, testing_outputs)**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(testing_predictions, testing_outputs, alpha=0.2)\n",
    "plt.plot([-0.4, 0.4], [-0.4, 0.4], 'k--')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "pickle.dump(model,open('Cd-Te.sav','wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = pickle.load(open('CdTe.sav','rb'))\n",
    "testing_predictions = clf.predict(testing_inputs)\n",
    "print(np.sqrt(np.mean(np.subtract(testing_predictions, testing_outputs)**2)))\n",
    "plt.scatter(testing_predictions, testing_outputs, alpha=0.2)\n",
    "plt.plot([-0.4, 0.4], [-0.4, 0.4], 'k--')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore', category=DeprecationWarning)\n",
    "warnings.filterwarnings('ignore', category=UserWarning)\n",
    "import nglview as nv\n",
    "import pytraj as pt\n",
    "\n",
    "p_traj = pt.load(nv.datafiles.TRR, top=nv.datafiles.PDB)\n",
    "p_view = nv.show_pytraj(p_traj)\n",
    "p_view._display_image()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even with as few as 5 random iterations, we see a huge effect on performance."
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}