QML_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QML Tutorial\n",
    "\n",
    "by Anders S. Christensen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial is a general introduction to kernel-ridge regression with Quantum Machine Learning (QML).\n",
    "\n",
    "QML is freely available under the terms of the MIT open-source license. However, we ask that you properly cite the relevant, underlying work in your published work (provided below).\n",
    "\n",
    "##### QML Citation: \n",
    "AS Christensen, FA Faber, B Huang, LA Bratholm, A Tkatchenko, KR Müller, OA von Lilienfeld (2017) \"QML: A Python Toolkit for Quantum Machine Learning\" https://github.com/qmlcode/qml\n",
    "\n",
    "##### Theory Citation:\n",
    "[Rupp et al, Phys Rev Letters, 2012](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.108.058301)\n",
    "\n",
    "[Montavon et al, New J Phys, 2013](http://iopscience.iop.org/article/10.1088/1367-2630/15/9/095003/meta)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory \n",
    "Regression model of some property, \n",
    "\n",
    "\\begin{equation*}\n",
    "y\n",
    "\\end{equation*}\n",
    "\n",
    "for some system, \n",
    "\n",
    "\\begin{equation*}\n",
    "\\widetilde{\\mathbf{X}}\n",
    "\\end{equation*}\n",
    "\n",
    "this could correspond to e.g. the atomization energy of a molecule:\n",
    "\n",
    "\\begin{equation*}\n",
    "y\\left(\\widetilde{\\mathbf{X}} \\right) = \\sum_i \\alpha_i \\  K\\left( \\widetilde{\\mathbf{X}}, \\mathbf{X}_i\\right)\n",
    "\\end{equation*}\n",
    "\n",
    "E.g. Using Gaussian kernel function with Frobenius norm:\n",
    "\n",
    "\\begin{equation*}\n",
    "K_{ij} = K\\left( \\mathbf{X}_i, \\mathbf{X}_j\\right) = \\exp\\left( -\\frac{\\| \\mathbf{X}_i - \\mathbf{X}_j\\|_2^2}{2\\sigma^2}\\right)\n",
    "\\end{equation*}\n",
    "\n",
    "Regression coefficients are obtained through kernel matrix inversion and multiplication with reference labels\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\alpha} = (\\mathbf{K} + \\lambda \\mathbf{I})^{-1} \\mathbf{y}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load QML into Jupyter Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Printing version info for help reporting bugs\n",
      "Python version: 3.6.5 (default, Apr 26 2018, 00:14:31) \n",
      "[GCC 4.8.5 20150623 (Red Hat 4.8.5-11)]\n",
      "QML version: 0.4.0\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import sys\n",
    "print(\"Printing version info for help reporting bugs\")\n",
    "print(\"Python version:\", sys.version)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "try:\n",
    "    import qml\n",
    "    from qml.kernels import gaussian_kernel\n",
    "    from qml.kernels import laplacian_kernel\n",
    "    from qml.math import cho_solve\n",
    "    print(\"QML version:\",qml.__version__)\n",
    "except ImportError:\n",
    "    print(\"Failed to find QML\")\n",
    "    print(\"Please follow instructions here: http://www.qmlcode.org/installation.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial exercises\n",
    "\n",
    "If you got this notebook without the repo, please clone the following GIT repository to access the necessary scripts and QM7 dataset (atomization energies and relaxed geometries at PBE0/def2-TZVP level of theory) for ~7k GDB1-7 molecules. \n",
    "[Ruddigkeit et al, J Chem Inf Model, 2012](https://pubs.acs.org/doi/abs/10.1021/ci300415d)\n",
    "[Rupp et al, Phys Rev Letters, 2012](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.108.058301)\n",
    "\n",
    "git clone https://github.com/qmlcode/tutorial.git\n",
    "\n",
    "Additionally, the repository contains Python3 scripts with the solutions to each exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1: Representations\n",
    "\n",
    "In this exercise we use QML to generate the Coulomb matrix and Bag of bonds (BoB) representations. [Montavon et al, New J Phys, 2013](http://iopscience.iop.org/article/10.1088/1367-2630/15/9/095003/meta)\n",
    "In QML data can be parsed via the <font color=red>Compound</font> class, which stores data and generates representations in Numpy's ndarray format.\n",
    "If you run the code below, you will read in the file <font color=red>qm7/0001.xyz</font> (a methane molecule) and generate a coulomb matrix representation (sorted by row-norm) and a BoB representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[36.8581052   5.50857022  0.5         5.50857007  0.56221501  0.5\n",
      "  5.5085695   0.56221605  0.56221611  0.5         5.50856526  0.56221669\n",
      "  0.56221777  0.56221405  0.5       ]\n",
      "[36.8581052   0.          0.          5.50857022  5.50857007  5.5085695\n",
      "  5.50856526  0.          0.          0.          0.          0.\n",
      "  0.          0.5         0.5         0.5         0.5         0.\n",
      "  0.56221777  0.56221669  0.56221611  0.56221605  0.56221501  0.56221405\n",
      "  0.          0.          0.          0.        ]\n"
     ]
    }
   ],
   "source": [
    "# Create the compound object mol from the file qm7/0001.xyz which happens to be methane\n",
    "mol = qml.Compound(xyz=\"qm7/0001.xyz\")\n",
    "\n",
    "# Generate and print a coulomb matrix for compound with 5 atoms \n",
    "mol.generate_coulomb_matrix(size=5, sorting=\"row-norm\")\n",
    "print(mol.representation)\n",
    "\n",
    "# Generate and print BoB bags for compound containing C and H\n",
    "mol.generate_bob(asize={\"C\":2, \"H\":5})\n",
    "print(mol.representation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The representations are simply stored as 1D-vectors.\n",
    "Note the keyword <font color=red>size</font> which is the largest number of atoms in a molecule occurring in test or training set. \n",
    "Additionally, the coulomb matrix can take a sorting scheme as keyword, and the BoB representations requires the specifications of how many atoms of a certain type to make room for in the representations.\n",
    "\n",
    "Lastly, you can print the following properties which is read from the XYZ file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.041682 -0.0562   -0.071481]\n",
      " [ 2.130894 -0.056202 -0.071496]\n",
      " [ 0.678598  0.174941 -1.072044]\n",
      " [ 0.678613  0.694746  0.62898 ]\n",
      " [ 0.678614 -1.038285  0.228641]]\n",
      "['C', 'H', 'H', 'H', 'H']\n",
      "[6 1 1 1 1]\n",
      "qm7/0001.xyz\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "# Print other properties stored in the object\n",
    "print(mol.coordinates)\n",
    "print(mol.atomtypes)\n",
    "print(mol.nuclear_charges)\n",
    "print(mol.name)\n",
    "print(mol.unit_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2: Kernels\n",
    "\n",
    "In this exercise we generate a Gaussian kernel matrix, \n",
    "\n",
    "\\begin{equation*}\n",
    "\\mathbf{K}\n",
    "\\end{equation*}\n",
    "\n",
    "using the representations\n",
    "\n",
    "\\begin{equation*}\n",
    "\\mathbf{X}\n",
    "\\end{equation*}\n",
    "\n",
    "which are generated similarly to the example in the previous exercise:\n",
    "\n",
    "\\begin{equation*}\n",
    "K_{ij} = \\exp\\left( -\\frac{\\| \\mathbf{X}_i - \\mathbf{X}_j\\|_2^2}{2\\sigma^2}\\right)\n",
    "\\end{equation*}\n",
    "\n",
    "QML supplies functions to generate the most basic kernels (E.g. Gaussian, Laplacian). In the exercise below, we calculate a Gaussian kernel for the QM7 dataset.\n",
    "In order to save time you can import the entire QM7 dataset as  <font color=red>Compound</font> objects from the file  <font color=red>tutorial_data.py</font> found in the tutorial GitHub repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[73.51669472 19.91549544 73.51669472 ...  0.          0.\n",
      "   0.        ]\n",
      " [73.51669472 24.77884434 53.3587074  ...  0.          0.\n",
      "   0.        ]\n",
      " [73.51669472 28.39419169 73.51669472 ...  0.          0.\n",
      "   0.        ]\n",
      " ...\n",
      " [73.51669472 39.57492699 36.8581052  ...  0.          0.\n",
      "   0.        ]\n",
      " [73.51669472 27.22906335 73.51669472 ...  0.          0.\n",
      "   0.        ]\n",
      " [73.51669472 39.78475661 36.8581052  ...  0.          0.\n",
      "   0.        ]]\n",
      "[[1.         0.99795995 0.99864847 ... 0.99905645 0.99812809 0.99747646]\n",
      " [0.99795995 1.         0.99795162 ... 0.99828487 0.99972878 0.9976251 ]\n",
      " [0.99864847 0.99795162 1.         ... 0.99885073 0.9981792  0.99752914]\n",
      " ...\n",
      " [0.99905645 0.99828487 0.99885073 ... 1.         0.99857214 0.9987595 ]\n",
      " [0.99812809 0.99972878 0.9981792  ... 0.99857214 1.         0.99788663]\n",
      " [0.99747646 0.9976251  0.99752914 ... 0.9987595  0.99788663 1.        ]]\n"
     ]
    }
   ],
   "source": [
    "# Import QM7, already parsed to QML\n",
    "from tutorial_data import compounds\n",
    "\n",
    "# For every compound generate a coulomb matrix or BoB\n",
    "for mol in compounds:\n",
    "\n",
    "    mol.generate_coulomb_matrix(size=23, sorting=\"row-norm\")\n",
    "    # mol.generate_bob(size=23, asize={\"O\":3, \"C\":7, \"N\":3, \"H\":16, \"S\":1})\n",
    "\n",
    "# Make a big 2D array with all the representations\n",
    "X = np.array([mol.representation for mol in compounds])\n",
    "\n",
    "# Print all representations\n",
    "print(X)\n",
    "\n",
    "# Run on only a subset of the first 1000 (for speed)\n",
    "X = X[:1000]\n",
    "\n",
    "# Define the kernel width\n",
    "sigma = 1000.0\n",
    "\n",
    "# K is also a Numpy array\n",
    "K = gaussian_kernel(X, X, sigma)\n",
    "\n",
    "# Print the kernel\n",
    "print(K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3: Regression\n",
    "\n",
    "With the kernel matrix and representations sorted out in the previous two exercise, we can now solve the \n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\alpha}\n",
    "\\end{equation*} \n",
    "\n",
    "regression coefficients:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\alpha} = (\\mathbf{K} + \\lambda \\mathbf{I})^{-1} \\mathbf{y}\n",
    "\\end{equation*}\n",
    "\n",
    "One of the most efficient ways of solving this equation is using a Cholesky-decomposition.\n",
    "QML includes a function named <font color=red>cho_solve()</font> to do this via the math module <font color=red>qml.math</font>.\n",
    "In this step it is convenient to only use a subset of the full dataset as training data (see below).\n",
    "The following builds on the code from the previous step.\n",
    "To save time, you can import the PBE0/def2-TZVP atomization energies for the QM7 dataset from the file <font color=red>tutorial_data.py</font>.\n",
    "This has been sorted to match the ordering of the representations generated in the previous exercise.\n",
    "Extend your code from the previous step with the code below:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.         0.99987237 0.99991548 ... 0.999941   0.9998829  0.99984209]\n",
      " [0.99987237 1.         0.99987185 ... 0.99989272 0.99998305 0.9998514 ]\n",
      " [0.99991548 0.99987185 1.         ... 0.99992813 0.9998861  0.99984539]\n",
      " ...\n",
      " [0.999941   0.99989272 0.99992813 ... 1.         0.9999107  0.99992242]\n",
      " [0.9998829  0.99998305 0.9998861  ... 0.9999107  1.         0.99986778]\n",
      " [0.99984209 0.9998514  0.99984539 ... 0.99992242 0.99986778 1.        ]]\n",
      "1000\n",
      "1000\n",
      "[ 2.27792067e+09 -1.02252571e+09 -5.08751603e+08  1.22295207e+09\n",
      " -6.57718637e+08  2.57277760e+09 -1.57445226e+08 -5.20879350e+08\n",
      "  1.25864418e+09  1.59068718e+07 -1.61319348e+09 -5.47179964e+08\n",
      " -8.83818584e+08 -7.08245848e+08  4.76566462e+08 -1.19741800e+09\n",
      " -1.03916151e+09  1.22268511e+08 -1.43150882e+09 -1.44851064e+09\n",
      " -8.65609117e+08  2.19972771e+09 -9.37737176e+08  5.35442253e+08\n",
      " -3.31995869e+06 -1.44634843e+09 -2.03426603e+09 -2.18121529e+09\n",
      "  1.80982079e+09 -1.52603827e+09  1.38844883e+09 -1.69434225e+09\n",
      " -6.85272902e+08  3.00409687e+08  8.34188323e+08 -9.51008774e+08\n",
      "  1.23761690e+09 -5.47571330e+08  4.89914734e+08 -2.33748001e+09\n",
      "  1.13305932e+09 -1.83517743e+08 -1.29375497e+08  6.55694611e+08\n",
      "  8.19053473e+08  3.80641967e+08 -2.07306207e+09  7.92872484e+08\n",
      "  2.11983106e+09  6.17835756e+07  2.08786005e+09 -4.96667410e+08\n",
      "  2.17877772e+08  1.09601270e+09  1.55395466e+09 -2.30753575e+08\n",
      "  1.96164846e+08  1.01624124e+09 -1.54563315e+08 -1.96483182e+09\n",
      " -7.42128581e+06 -2.90216879e+08  3.07496741e+08 -2.50825726e+09\n",
      "  6.91418461e+08  7.93575619e+08 -9.46967098e+08  9.12071336e+08\n",
      " -2.03777935e+08  1.02308117e+09  2.71768060e+07  8.61292766e+08\n",
      " -1.64511096e+09 -2.24421351e+09 -3.44101744e+07 -4.99974952e+08\n",
      "  5.44165412e+08  2.33133835e+08 -4.57490727e+08 -1.14693378e+09\n",
      " -7.93426598e+08 -6.38744659e+08  5.80791660e+08  2.09444083e+09\n",
      " -6.77986513e+08 -3.66926238e+08 -5.06519466e+08  5.70024716e+08\n",
      "  1.14588288e+09 -3.80431613e+08  8.37248488e+08  1.55751945e+09\n",
      " -1.63724334e+09 -1.16648077e+08  6.14020377e+08  1.26539280e+09\n",
      " -7.68801300e+08 -8.66526417e+08 -3.28666406e+07  1.09925263e+09\n",
      " -8.93347374e+08 -4.48077253e+08  3.85801890e+08 -9.38287541e+08\n",
      "  1.34857736e+09 -5.92431829e+07 -7.63394858e+08 -2.40781532e+09\n",
      " -6.76005251e+07 -5.72835729e+08  1.81267199e+09 -9.91418714e+08\n",
      "  1.16105053e+09 -1.17429281e+09 -1.49205490e+09  1.47097380e+09\n",
      " -2.68122581e+07 -7.08110101e+08  3.82067735e+08  6.26298441e+08\n",
      "  2.19213088e+08 -1.95650988e+09 -6.82238442e+08  5.38923626e+08\n",
      " -3.07624339e+07 -7.61673671e+08 -1.24363932e+09  1.86214791e+09\n",
      " -2.29564433e+09  1.51968830e+09 -5.13922808e+08  2.44485036e+08\n",
      " -1.66720008e+09  1.98909380e+08 -3.74166788e+08 -1.25425081e+09\n",
      " -2.24255250e+06 -1.17051543e+09 -8.53505199e+08  1.11994302e+09\n",
      " -1.11698623e+09  1.84196985e+09  2.56026515e+07  1.00353683e+09\n",
      " -4.27941383e+08  1.84628697e+09  1.86876638e+08  1.03439693e+09\n",
      " -2.59283720e+09 -1.20135838e+08  3.36970879e+07 -6.21579483e+08\n",
      " -4.19706426e+07 -2.23077559e+09 -3.33421083e+08  7.47176831e+08\n",
      " -1.57288636e+08 -2.42635125e+08  3.38703180e+08 -9.98253395e+08\n",
      " -1.33123835e+08 -1.80791278e+08  3.09868111e+08  1.21365146e+09\n",
      "  3.58589336e+08  9.35308294e+08  1.56713948e+09  4.69299621e+08\n",
      "  1.30212988e+09 -6.93268917e+08 -8.04694152e+08  3.49086615e+09\n",
      " -1.74479765e+09  5.67231057e+08  1.40291740e+09 -1.64021691e+09\n",
      " -5.45334377e+07 -1.15133793e+08  1.43744589e+09  6.22912780e+08\n",
      " -2.07149208e+09 -1.53992753e+08  1.61280684e+08  1.06600543e+09\n",
      "  1.14173545e+08  9.12909467e+07  7.01456972e+08 -1.26733377e+09\n",
      " -1.25188944e+09 -1.21440876e+09 -9.84911439e+07 -5.60881803e+08\n",
      " -9.47765781e+08  2.64543565e+09  1.82260864e+09 -1.16709132e+09\n",
      "  2.89669825e+08  6.97515559e+08  1.64028086e+09 -6.02827968e+07\n",
      "  6.14236631e+08 -1.26213091e+09  1.43326480e+09 -6.80766480e+08\n",
      "  4.34103025e+08  1.17918880e+09 -4.65186088e+08 -2.98595095e+08\n",
      " -9.24794781e+08 -1.40636192e+09  1.35402828e+09  3.22378732e+08\n",
      " -5.24314055e+08  2.10066483e+09  3.12916775e+08  4.46327718e+08\n",
      " -1.05423736e+09 -1.03250734e+09  1.29635898e+09 -3.74698117e+08\n",
      " -2.09776568e+08  5.60691253e+08 -3.22841646e+08  3.17757865e+08\n",
      " -1.59630359e+09  6.25558997e+08 -1.88936971e+09  5.14719810e+07\n",
      " -2.44286494e+09  1.48812254e+09 -9.45248895e+08  1.21719157e+09\n",
      "  5.78830882e+08  2.92222346e+08 -2.70790980e+08 -7.42832084e+08\n",
      "  5.42494716e+09 -6.84552218e+08 -2.36992712e+08 -2.40277211e+09\n",
      "  1.47943864e+09 -1.67651780e+08 -1.20908518e+09 -2.15821083e+08\n",
      "  1.31928947e+08 -2.22954026e+08 -1.48721501e+09 -1.29377448e+08\n",
      "  1.78585728e+09 -3.51954404e+09 -6.03992744e+08  8.34605917e+08\n",
      "  8.13809715e+07 -1.67933932e+09  5.69296883e+09  1.81603912e+09\n",
      "  9.34539070e+08 -1.49732553e+08 -9.98814388e+08  3.15486615e+09\n",
      " -1.22009270e+08 -3.86789668e+08  7.38856434e+08 -1.29908241e+09\n",
      "  1.21035328e+09  2.96208800e+08 -1.17756725e+08 -1.34953496e+09\n",
      "  1.74877741e+07 -2.87548381e+09  9.43639879e+08  5.77943042e+08\n",
      "  1.16910416e+09  1.85346426e+09 -3.64812635e+08  6.97711118e+08\n",
      " -7.05983448e+08  5.11444469e+08  7.47318485e+08  8.89263063e+08\n",
      "  2.17356778e+08  1.97479609e+08 -1.55197015e+09 -1.04893125e+09\n",
      " -3.03307132e+08 -3.38385390e+08 -8.36061689e+08 -3.13249371e+08\n",
      " -2.45166603e+09 -1.56447633e+09  1.07626445e+09 -4.08970203e+07\n",
      "  6.46185012e+08 -1.68097317e+09  3.09021041e+09 -1.99751080e+08\n",
      "  5.47378928e+08 -1.56363857e+09 -1.19799112e+09  4.33260617e+08\n",
      " -1.68125109e+08 -2.28169924e+08 -1.59502352e+09  2.40230367e+08\n",
      " -2.22566062e+09  1.55678826e+09 -5.94832435e+08 -2.47376975e+07\n",
      " -5.78205404e+08  3.05277298e+09  1.05653036e+09 -3.78923525e+08\n",
      " -9.57925106e+08 -2.15868502e+09 -8.93848144e+08 -7.68804843e+08\n",
      "  2.21547640e+09  1.72425028e+09  7.13503986e+08 -3.79701118e+08\n",
      "  1.35115087e+09 -3.52306153e+08  3.05004874e+08  5.49969826e+08\n",
      " -2.75511244e+09 -1.73644339e+09  2.47032914e+08 -1.32295305e+09\n",
      " -3.68632924e+08 -1.60301629e+08  1.36964829e+09 -1.28713855e+09\n",
      "  1.12355745e+09  1.14968679e+09  1.49962441e+09 -1.31513482e+09\n",
      "  9.82359113e+08  3.38980451e+08 -1.16739778e+09 -3.57117523e+08\n",
      " -1.83299027e+08 -5.32032962e+08 -1.16304719e+09  2.03361727e+08\n",
      "  1.28872088e+09  1.04218304e+09  1.45771184e+08  1.66756001e+09\n",
      " -1.54490175e+09  3.55274605e+08  3.51988653e+08 -5.89134097e+08\n",
      " -1.38647185e+09  1.47940768e+09  3.16580376e+08  5.05651414e+08\n",
      "  1.10538071e+08  1.44218368e+09  8.57624226e+08 -8.39863468e+08\n",
      "  3.56545409e+07 -2.44122310e+08  1.50501771e+09 -8.03658769e+08\n",
      " -8.97109353e+08  2.56872565e+09  4.82383722e+08 -4.68892052e+08\n",
      "  8.47556250e+08  1.35131484e+09 -9.87133362e+07  5.80697294e+08\n",
      " -1.63006358e+08  1.59513307e+09  2.98368267e+08 -9.38288717e+08\n",
      "  1.69124916e+09  8.61583570e+08 -5.86509461e+08  1.62107049e+08\n",
      " -1.66517849e+09  5.35511212e+08  4.80544105e+08  1.05062510e+09\n",
      "  4.87996578e+08 -7.57424374e+08  1.02570852e+09  3.38718107e+08\n",
      " -1.11645757e+09  2.85827549e+08  2.35136892e+08 -4.18046702e+08\n",
      "  1.54407859e+09 -1.94579470e+09  1.37656825e+09  4.40966186e+08\n",
      "  2.75942541e+08 -1.04215659e+09  7.65858502e+08 -2.43985773e+09\n",
      "  8.04421943e+07 -1.07714650e+09  7.71168569e+08 -5.89371469e+08\n",
      "  6.34567787e+07 -9.14946024e+07 -9.81760778e+08 -1.21103065e+09\n",
      " -2.69264839e+08 -1.94898499e+09 -4.82008370e+08 -7.95904655e+08\n",
      "  5.35001762e+08  1.92578914e+09 -2.28660800e+09  4.20822095e+08\n",
      " -2.00682788e+09 -5.16706370e+08 -2.73295985e+08  1.69943041e+08\n",
      "  1.33097324e+09 -5.14595409e+08  9.02965174e+07 -2.88769594e+09\n",
      "  4.63750676e+09  3.35176113e+08 -9.06594287e+08  3.26480482e+08\n",
      " -4.06273179e+08 -8.11607535e+08 -1.93964056e+09  6.35689986e+08\n",
      " -6.90730790e+08 -3.73807760e+08  1.41527485e+08 -1.84356073e+08\n",
      "  9.39134879e+08 -1.30243557e+08 -1.90106390e+09  8.06162450e+07\n",
      "  1.18209747e+09 -1.54116248e+09  5.18747727e+08  1.36497764e+09\n",
      " -8.31086804e+08 -5.77001313e+08  2.86506497e+08  1.20221072e+09\n",
      "  1.21730326e+09  3.56668963e+08 -6.09910124e+08 -1.03483196e+09\n",
      " -9.33531472e+08  8.36911854e+07  8.75969944e+08 -2.49468152e+07\n",
      "  4.11383600e+08  4.69463446e+08 -9.53774858e+08  1.33891111e+09\n",
      "  3.48160253e+07 -1.31525646e+09  1.04046235e+09 -7.12162069e+08\n",
      " -6.28045627e+08 -8.30541699e+08  2.15112783e+09  4.41490824e+08\n",
      " -1.71354005e+08  2.48363572e+09 -4.85440505e+08  3.74859589e+08\n",
      "  1.03183206e+08  1.62632277e+09 -6.68946412e+08  5.58019433e+08\n",
      " -1.05979737e+09  8.17092192e+08  1.03218821e+09 -2.28336157e+09\n",
      " -2.77353415e+09  2.11454414e+08  1.65186774e+09  7.19845665e+08\n",
      "  1.47646382e+09 -3.48070097e+08 -2.71595411e+08  8.38768677e+07\n",
      " -4.12458538e+08  1.55900144e+09  1.61607274e+09  1.55868418e+09\n",
      "  5.53470106e+08  9.17003441e+08 -1.22043886e+09  3.00643551e+09\n",
      "  1.04364205e+09  1.16834008e+09  1.98603533e+09  4.32248240e+08\n",
      " -5.43649428e+08  4.11483176e+08 -1.28233221e+08 -4.69054727e+08\n",
      " -1.32598195e+08 -1.00369848e+09 -1.99768490e+09  9.22047578e+08\n",
      " -1.75229042e+08 -5.67113168e+08 -1.66341275e+06  6.53740622e+08\n",
      "  5.95546344e+07  1.18911087e+09 -1.36938574e+09  1.63800605e+09\n",
      "  2.55291375e+08 -1.76868050e+09 -2.00446992e+09 -2.39044450e+08\n",
      "  3.58554078e+08  2.06195665e+08 -2.03284825e+08 -1.93574346e+09\n",
      " -3.07865252e+07  9.82165492e+07  2.61220316e+09  1.90951958e+09\n",
      " -1.76424179e+09 -6.79294908e+08  2.62944125e+09  2.85746900e+08\n",
      " -1.72472663e+09 -3.36301206e+08 -3.03105761e+08  9.43398427e+08\n",
      "  1.63087831e+09  9.79010139e+07 -7.75502763e+08  6.20691762e+08\n",
      "  1.50545959e+09  4.71401816e+08  3.11941233e+09  3.18521847e+08\n",
      " -9.13203309e+08  5.97881937e+08  1.01902895e+08 -4.85121921e+08\n",
      "  8.67488371e+08  1.86587225e+08 -7.40435813e+08 -4.55674240e+08\n",
      "  2.12350140e+08  8.25719005e+08  4.25671778e+08 -8.73402909e+08\n",
      " -1.09679171e+09  1.10045358e+09  4.84848153e+08  7.86266329e+08\n",
      "  9.39736019e+08 -8.42425800e+08  8.32203802e+08 -7.29910993e+08\n",
      " -1.76571371e+09  9.41047172e+08 -1.63028365e+09 -7.96617704e+08\n",
      " -2.64254820e+09 -1.76511789e+09  1.51677397e+09  9.97665301e+08\n",
      "  3.02657779e+08  1.05506745e+09 -5.72830012e+08 -2.15273098e+09\n",
      "  1.83249977e+09  1.68388221e+09 -1.29479583e+09  1.37513229e+09\n",
      "  1.27281193e+09  9.94991617e+08  6.16441055e+08 -1.68104921e+09\n",
      " -1.07612450e+09  1.53389257e+09 -1.28573948e+08  1.74674645e+07\n",
      "  4.98754150e+08 -1.85073867e+09 -8.39739077e+08  1.17438591e+09\n",
      "  1.29322607e+09  2.52911669e+08  6.83826681e+07 -1.28319127e+09\n",
      " -1.15472688e+09 -1.85297759e+09  1.01514001e+09  4.02189697e+06\n",
      "  2.04493010e+08  6.40624800e+07 -1.33744737e+09 -1.30005787e+09\n",
      " -5.94421627e+08 -6.19308272e+08 -2.57047273e+09 -5.61799930e+08\n",
      " -2.96854103e+08 -1.85448494e+09  7.85479051e+08 -2.58150635e+09\n",
      "  8.39906918e+08  6.04772955e+08 -2.04372101e+09  1.78394397e+09\n",
      "  8.63686321e+07  2.15392977e+08  3.31967168e+08  7.08584980e+08\n",
      " -2.22296669e+09 -1.85858879e+09 -4.93889048e+08 -6.82069262e+08\n",
      "  8.84794859e+08 -1.20660252e+09  6.23236585e+08  1.57438154e+09\n",
      "  1.68319287e+09 -9.40582230e+07  1.26042933e+08 -1.92698497e+09\n",
      " -7.23984313e+08  1.52235216e+09 -8.92726202e+08 -9.09092083e+07\n",
      "  3.53004810e+08  4.81475443e+08 -2.00522127e+09 -9.03514220e+08\n",
      "  1.42554105e+09  9.60524649e+08 -3.73811121e+08  7.85333555e+08\n",
      " -4.25656358e+08 -1.60964132e+09 -1.67735796e+08 -1.78609816e+09\n",
      "  1.45463946e+09  6.51568696e+08  3.46508202e+08 -6.97319507e+08\n",
      " -3.72550033e+08  1.18818204e+09 -2.41349819e+08 -2.39411990e+08\n",
      " -2.75649160e+08 -1.12639203e+09  1.94845352e+08  9.18166109e+08\n",
      " -1.43635613e+09 -1.42360602e+08 -5.04869201e+08  3.84733667e+08\n",
      " -6.71441775e+08  8.35452322e+08 -1.87208189e+08  6.27073382e+08\n",
      "  5.47492899e+08 -7.98319978e+08 -1.11789905e+09 -6.41160847e+08\n",
      " -9.15614948e+07  4.49000774e+08  1.20848941e+08 -1.66933985e+07\n",
      "  1.65343205e+09  1.35392106e+09  1.05128807e+09  5.24076452e+08\n",
      " -2.94207883e+08  2.43281013e+08 -1.25853274e+09  1.12416683e+09\n",
      "  8.45030194e+08 -8.51189499e+08 -1.45836860e+09  1.85035636e+09\n",
      " -1.43331747e+09  1.79649557e+09 -1.06112840e+09 -9.01813159e+08\n",
      " -1.46375110e+09 -1.19434727e+09 -1.96839530e+09 -4.00879843e+08\n",
      " -1.34423141e+09 -1.62740838e+09 -1.27212443e+09 -6.30634662e+08\n",
      "  4.88995776e+07 -3.36723894e+08  2.71824600e+09  6.76316299e+07\n",
      " -9.31554640e+08  1.38278531e+09  9.87857342e+08 -2.41082209e+09\n",
      "  6.50548151e+08 -1.88094282e+09  3.52127630e+09 -7.29116141e+08\n",
      "  6.44538466e+08  6.40392740e+08  5.99280789e+07  2.66464883e+09\n",
      "  5.20842462e+08  5.75005860e+08  1.82407014e+09 -2.62475906e+09\n",
      " -2.35050920e+09 -5.87404574e+07  2.49646046e+08 -1.40916115e+09\n",
      " -2.63808384e+08  1.83233764e+09 -1.30630049e+09  9.36309699e+07\n",
      "  5.75761283e+08  1.70974950e+09  4.27104966e+08 -3.29262335e+08\n",
      " -8.68018481e+08  3.61147926e+08 -4.94772368e+08 -6.08676980e+08\n",
      " -2.19695716e+09 -4.77804362e+08 -1.01663045e+09  5.11904915e+07\n",
      " -4.36715579e+08 -3.94038963e+08  4.92320120e+08  8.38537713e+08\n",
      " -1.44915806e+09 -1.09169128e+09  1.43279747e+09  2.20967492e+09\n",
      "  4.12193824e+08 -2.50295196e+09 -1.98896568e+09 -7.60883833e+08\n",
      "  9.48267776e+08  2.03885637e+09  2.55227748e+08 -1.26585947e+09\n",
      " -1.13515856e+08  9.70492913e+08 -1.19825116e+09 -1.06102230e+09\n",
      " -6.34430728e+07  5.03776089e+08 -1.45575795e+09 -1.32104648e+09\n",
      "  1.45490565e+09  6.51620806e+08 -9.98253613e+08 -2.44827794e+09\n",
      " -2.26520947e+08 -1.62419825e+09  3.32838887e+09  2.02396748e+09\n",
      "  1.09546110e+09 -1.41077234e+08  1.64691995e+09  1.51250439e+09\n",
      " -1.02026279e+09 -4.36454713e+08  2.58874692e+09  1.14273696e+09\n",
      " -1.77722554e+08 -2.25539650e+08 -2.56915702e+08 -6.85091885e+08\n",
      " -1.49431091e+09 -2.03338625e+09  1.48200271e+09  2.04184210e+09\n",
      " -6.59664237e+08 -3.47500966e+08 -1.34063205e+09 -2.04852707e+09\n",
      "  1.93327350e+09  7.00589138e+08  8.59039164e+08 -1.38436501e+09\n",
      " -3.25863446e+08 -6.21591626e+08 -1.53856500e+09  2.16289296e+08\n",
      " -6.74567248e+08  2.12068286e+09  7.15569759e+08 -7.62724832e+08\n",
      "  2.60521612e+07  4.43854613e+08  8.26880379e+08 -1.84851175e+09\n",
      "  2.80066452e+08 -1.24895354e+09  4.50286138e+08  7.40946708e+08\n",
      " -5.28834802e+08 -1.60322917e+09 -4.38304327e+08 -1.91738921e+09\n",
      " -4.49391072e+08  9.42793979e+07  2.69539268e+07  1.42559323e+09\n",
      "  3.32956352e+07  7.14810699e+08 -4.98639064e+08  2.92955049e+08\n",
      " -1.48326793e+09  4.11799772e+08 -5.61423762e+08  4.07101971e+08\n",
      " -1.36929345e+09 -1.21552628e+08  1.55600218e+09  1.14792042e+09\n",
      "  5.89276368e+08  4.73463946e+08  1.50158348e+09 -2.85682724e+08\n",
      " -9.38559166e+08  6.86133474e+08  4.09221014e+08 -4.43496956e+08\n",
      " -9.03646881e+08 -2.21589013e+09 -1.16625629e+08 -8.98293608e+08\n",
      " -1.99270237e+08  6.15077136e+08  3.87719344e+08 -7.04535685e+07\n",
      "  7.60160523e+07 -1.21107380e+09 -1.01034972e+09  5.65784178e+08\n",
      "  6.78096624e+08 -7.27185006e+07 -1.65826517e+09  1.86585723e+08\n",
      "  1.65204350e+09  2.30828150e+09 -9.44249886e+08  5.48146300e+05\n",
      "  1.64125396e+09 -2.55963462e+09  2.76314536e+07 -6.97147057e+08\n",
      "  2.26399760e+08  1.88353868e+09  6.24741117e+08  1.96796503e+09\n",
      "  1.74803007e+08  1.83045758e+09 -1.23742609e+09  1.23832133e+09\n",
      " -4.69262972e+07  7.69278536e+08 -6.66711871e+08 -1.42224616e+09\n",
      "  2.06699559e+09  8.17888076e+08  1.30251629e+09 -2.52578500e+08\n",
      " -5.41166900e+08 -3.26019442e+09  1.00212526e+09 -6.78189488e+08\n",
      " -1.76213895e+09  1.56236099e+09 -7.10452450e+08  1.15549238e+09\n",
      "  1.16657774e+08  3.26745881e+08 -1.31287242e+09  2.31410169e+09\n",
      " -5.82907201e+08  9.37711426e+08  1.16840455e+09  9.54019888e+07\n",
      "  5.85145487e+08  3.61969784e+08  1.42393215e+09  1.41875793e+08\n",
      " -5.02740415e+08  1.98398044e+09 -1.05472319e+09  5.11302324e+08\n",
      "  3.02772097e+09  1.20262740e+08  7.24051526e+08  1.34459510e+09\n",
      "  6.78746677e+08 -3.27359631e+08  7.23629981e+08 -4.99837380e+08\n",
      "  7.95329161e+08  1.22113886e+09  1.68188473e+08 -1.24126676e+09\n",
      "  7.21052741e+08 -9.37759781e+08  1.14709700e+09  1.59732526e+09\n",
      "  1.11240050e+08 -2.57152912e+08  2.36599012e+08 -6.87713852e+08\n",
      "  2.18088601e+08  1.11664237e+08 -4.51805760e+08  2.39235386e+08\n",
      "  1.44893845e+09  2.18869442e+09 -5.32507398e+08 -5.55847765e+08\n",
      "  6.95955283e+08 -1.36467628e+09 -9.00079540e+08 -1.50108772e+09\n",
      " -2.55307572e+09  1.10616780e+09  2.30877362e+09  3.44568276e+08\n",
      "  1.28221768e+09  1.43845202e+09 -1.79922342e+08  2.64412710e+09\n",
      " -2.56092979e+09  1.85322722e+09 -1.24305865e+09 -6.17222599e+08\n",
      " -7.36710331e+08  1.11144873e+09 -2.82782324e+09  1.36016404e+08\n",
      " -1.03646387e+09  2.15612421e+08  8.25772312e+08 -9.08423387e+08\n",
      "  7.63260332e+08  1.01668425e+09 -4.59020186e+08 -2.87794062e+09\n",
      "  2.78019854e+08 -1.17864588e+09 -6.39781027e+08 -9.37100038e+08\n",
      "  3.23851708e+08  1.17495814e+09  1.05914797e+09 -1.67341299e+09\n",
      "  3.61312043e+08 -3.29596630e+08  1.87922493e+09  2.24231619e+09\n",
      "  4.18780570e+08  1.23806393e+09  1.32664834e+08 -1.00127818e+09\n",
      " -6.82360606e+08 -4.36574366e+08 -8.21068709e+08 -6.54151272e+08\n",
      " -1.09325830e+09  5.50385373e+08  6.56047314e+08 -6.36702494e+08\n",
      "  4.07663683e+08 -5.26990119e+08  3.41793873e+08  4.05668893e+08\n",
      "  1.39689452e+08 -1.04849534e+09 -6.19158796e+08 -3.75627834e+09\n",
      "  1.77369721e+09 -4.56713720e+08  9.09141349e+08 -2.37914918e+09\n",
      "  1.76619313e+09 -5.43908179e+08  3.97552199e+08  1.28200817e+09]\n"
     ]
    }
   ],
   "source": [
    "from tutorial_data import energy_pbe0\n",
    "\n",
    "# Assign 1000 first molecules to the training set\n",
    "X_training = X[:1000]\n",
    "Y_training = energy_pbe0[:1000]\n",
    "\n",
    "sigma = 4000.0\n",
    "K = gaussian_kernel(X_training, X_training, sigma)\n",
    "print(K)\n",
    "\n",
    "# Add a small lambda to the diagonal of the kernel matrix\n",
    "K[np.diag_indices_from(K)] += 1e-8\n",
    "\n",
    "print(len(K))\n",
    "print(len(Y_training))\n",
    "\n",
    "\n",
    "# Use the built-in Cholesky-decomposition to solve\n",
    "alpha = cho_solve(K, Y_training)\n",
    "\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4: Prediction\n",
    "\n",
    "With the\n",
    "\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\alpha}\n",
    "\\end{equation*}\n",
    "\n",
    "regression coefficients from the previous step, we have (successfully) trained the machine, and we are now ready to do predictions for other compounds.\n",
    "This is done using the following equation:\n",
    "    \n",
    "\\begin{equation*}\n",
    "y\\left(\\widetilde{\\mathbf{X}} \\right) = \\sum_i \\alpha_i \\  K\\left( \\widetilde{\\mathbf{X}}, \\mathbf{X}_i\\right)\n",
    "\\end{equation*}\n",
    "\n",
    "In this step we further divide the dataset into a training and a test set. Try using the last 1000 entries as test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "254.2353172018013\n"
     ]
    }
   ],
   "source": [
    "# Assign 1000 last molecules to the test set\n",
    "X_test = X[-1000:]\n",
    "Y_test = energy_pbe0[-1000:]\n",
    "\n",
    "# calculate a kernel matrix between test and training data, using the same sigma\n",
    "Ks = gaussian_kernel(X_test, X_training, sigma)\n",
    "\n",
    "# Make the predictions\n",
    "Y_predicted = np.dot(Ks, alpha)\n",
    "\n",
    "# Calculate mean-absolute-error (MAE):\n",
    "print(np.mean(np.abs(Y_predicted - Y_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5: Learning curves\n",
    "\n",
    "Repeat the prediction from Exercise 2.4 with training set sizes of 1000, 2000, and 4000 molecules.\n",
    "\n",
    "Note the MAE for every training size.\n",
    "\n",
    "Plot a learning curve of the MAE versus the training set size.\n",
    "\n",
    "Generate a learning curve for the Gaussian and Laplacian kernels, as well using the coulomb matrix and bag-of-bonds representations.\n",
    "\n",
    "Which combination gives the best learning curve? Note you will have to adjust the kernel width (sigma) underway."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6: Delta learning\n",
    "\n",
    "A powerful technique in machine learning is the delta learning approach. Instead of predicting the PBE0/def2-TZVP atomization energies, we shall try to predict the difference between DFTB3 (a semi-empirical quantum method) and PBE0 atomization energies.\n",
    "Instead of importing the <font color=red>energy_pbe0</font> data, you can import the <font color=red>energy_delta</font> and use this instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from tutorial_data import energy_delta\n",
    "\n",
    "Y_training = energy_delta[:1000]\n",
    "Y_test = energy_delta[-1000:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally re-draw one of the learning curves from the previous exercise, and note how the prediction improves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}