machine_learning_lecture2_2020.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "ml_lecture2_2020.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "dbPgZ-eSVsIK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Install packages first\n",
        "!pip install qml\n",
        "!pip install ase"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ch_WRVLhRVDy",
        "colab_type": "text"
      },
      "source": [
        "# Lecture 2: Introduction to Machine Learning in Python, May 2020\n",
        "By Dr. Anders Christensen `anders.christensen @ unibas.ch`\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aVYtPM6QyYfJ",
        "colab_type": "text"
      },
      "source": [
        "## Part 1: Non-linear regression?\n",
        "Last time we saw *linear least squares regression*:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RMwRKix_QqDl",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "# plt.rcParams['figure.figsize'] = [8, 6]\n",
        "import numpy as np\n",
        "\n",
        "# X-values\n",
        "x = np.arange(0,20.0, 0.2)\n",
        "\n",
        "# Y-values: Y = 1.2*X + random noise \n",
        "y = 1.2 * x + np.random.normal(scale=2.0, size=len(x))\n",
        "\n",
        "plt.scatter(x,y)\n",
        "plt.plot(x, x*1.2, color=\"g\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WzeQ4TgjfLc_",
        "colab_type": "text"
      },
      "source": [
        "Approximate a target function as a weighted sum of the features:\n",
        "\\begin{equation}\n",
        "y(\\mathbf{x}) = x_1 \\alpha_1 + x_2 \\alpha_2 + \\dots + x_n \\alpha_n\n",
        "\\end{equation}\n",
        "In matrix notation:\n",
        "\n",
        "\\begin{equation}\n",
        "\\mathbf{y} = \\mathbf{X} \\mathbf{\\alpha}\n",
        "\\end{equation}\n",
        "\n",
        "Minimze the error:\n",
        "\\begin{equation}\n",
        "\\mathbf{\\hat{\\alpha}} = \\text{arg min} || \\mathbf{y}^\\text{ref} - \\mathbf{X}\\mathbf{\\alpha}||^2\n",
        "\\end{equation}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9uYJbdwGfd61",
        "colab_type": "text"
      },
      "source": [
        "## What a about linear regression for $\\sin \\left(x\\right)$?\n",
        "\n",
        "Same example as above, just with $y(x) = \\sin \\left(x\\right)$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RyL6ZS4tfeli",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "\n",
        "x = np.arange(0,6.6, 0.6)\n",
        "y = np.sin(x) #+ (np.random.random(size=len(x)) - 0.5) * 0.5\n",
        "\n",
        "print(x.shape)\n",
        "print(y.shape)\n",
        "\n",
        "xplot = np.arange(0,6.6, 0.01)\n",
        "\n",
        "plt.scatter(x,y, label=\"Training\")\n",
        "plt.plot(xplot, np.sin(xplot), color=\"g\", label=\"sin(x)\")\n",
        "\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lMFTEmJmSMq5",
        "colab_type": "text"
      },
      "source": [
        "Just like yesterday, we can use `numpy.linalg.lstsq()` to solve the optimal coefficients by minimizing the error:\n",
        "\n",
        "\\begin{equation}\n",
        "\\mathbf{\\hat{\\alpha}} = \\text{arg min} || \\mathbf{y}^\\text{ref} - \\mathbf{X}\\mathbf{\\alpha}||^2\n",
        "\\end{equation}\n",
        "Closed-form solution:\n",
        "\\begin{equation}\n",
        "\\mathbf{\\hat{\\alpha}} = \\left(\\mathbf{X}^\\top\\mathbf{X} \\right)^{-1}\\mathbf{X}^\\top\\mathbf{y}^\\text{ref}\n",
        "\\end{equation}\n",
        "However, Numpy uses a singluar-value decomposition."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7fLArAZOgYWp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_reshape = x.reshape(len(x),1)\n",
        "alpha, sing, rank, vecs = np.linalg.lstsq(x_reshape, y, rcond=None)\n",
        "\n",
        "print(alpha)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IedrklG3THXJ",
        "colab_type": "text"
      },
      "source": [
        "Making new predictions:\n",
        "*   Test set of new $x$-values\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xAmLtPPpg_hg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_test = np.arange(0.0, 6.6, 0.6)\n",
        "print(x_test)\n",
        "print(x_test.shape)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J8ZuXyFxTg9V",
        "colab_type": "text"
      },
      "source": [
        "Predictions are made with the same equation:\n",
        "\\begin{equation}\n",
        "\\mathbf{y} = \\mathbf{X} \\mathbf{\\alpha}\n",
        "\\end{equation}"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vK4vph3sTfbU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_test = x_test.reshape(len(x_test),1)\n",
        "y_test = np.dot(x_test, alpha)\n",
        "\n",
        "print(y_test)\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7E1Jg9umTfwk",
        "colab_type": "text"
      },
      "source": [
        "Making a plot of the predictions:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "C9l-otXyhb8y",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.plot(xplot, np.sin(xplot), color=\"g\", label=\"sin(x)\")\n",
        "plt.scatter(x, y, label=\"Training\")\n",
        "plt.scatter(x_test, y_test, color=\"r\", label=\"Test\")\n",
        "\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "stvSrhLWiNcR",
        "colab_type": "text"
      },
      "source": [
        "## Part 1½: Kernel ridge regression\n",
        "The exercise today is about coding kernel ridge regession!\n",
        "\n",
        "The idea is to describe the target function in terms of basis functions (\"kernel functions\") placed on the training points.\n",
        "\n",
        "Plot of basis functions for the $\\sin\\left(x\\right)$ curve:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2jMjypubgW6a",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import qml.kernels\n",
        "import qml.math\n",
        "\n",
        "sigma = 0.3\n",
        "K = qml.kernels.gaussian_kernel(x_reshape, x_reshape, sigma)\n",
        "K[np.diag_indices_from(K)] += 1e-10\n",
        "alphas = qml.math.cho_solve(K, y)\n",
        "\n",
        "xgauss = np.arange(0.0,6.0, 0.01)\n",
        "fit_curve = np.zeros(len(xgauss))\n",
        "\n",
        "for alp, xval in zip(alphas, x):\n",
        "  ygauss = np.exp(-(xgauss-xval)**2/(2*sigma**2)) * alp\n",
        "  fit_curve += ygauss\n",
        "  # plt.plot(xgauss, ygauss, color=\"C0\")\n",
        "\n",
        "plt.scatter(x,y, label=\"Training\")\n",
        "plt.plot(xgauss, np.sin(xgauss), color=\"g\", label=\"sin(x)\")\n",
        "plt.plot(xgauss, fit_curve, color=\"C3\", label=\"KRR\")\n",
        "\n",
        "# plt.ylim([-1.5,1.5])\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hQfEiCtqll6-",
        "colab_type": "text"
      },
      "source": [
        "$y\\left(x\\right)$ is now a sum of weighted basis functions:\n",
        "\\begin{equation}\n",
        "\\hat{y}\\left(\\tilde{\\mathbf{x}}\\right) = \\sum_i \\kappa\\left(\\mathbf{x}_i, \\tilde{\\mathbf{x}} \\right) \\alpha_i\n",
        "\\end{equation}\n",
        "\n",
        "As in linear regression, $\\mathbf{\\alpha}_i$ are our regression weights.\n",
        "\n",
        "In matrix form:\n",
        "\\begin{equation}\n",
        "\\mathbf{y} = \\mathbf{K}\\mathbf{\\alpha}\n",
        "\\end{equation}\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "\n",
        "For example, Gaussian kernel:\n",
        "\\begin{equation}\n",
        "\\kappa\\left(\\mathbf{x}, \\tilde{\\mathbf{x}} \\right) = \\exp \\left(-\\frac{\\|\\mathbf{x} - \\tilde{\\mathbf{x}} \\|^2}{2\\sigma^2} \\right)\n",
        "\\end{equation}\n",
        "Alternative definition: (used in the exercises today)\n",
        "\\begin{equation}\n",
        "\\kappa\\left(\\mathbf{x}, \\tilde{\\mathbf{x}} \\right) = \\exp \\left(-\\gamma \\|\\mathbf{x} - \\tilde{\\mathbf{x}} \\|^2 \\right)\n",
        "\\end{equation}\n",
        "The Gaussian kernel is always number between 0 and 1:\n",
        "*   0 if $x$ and $\\tilde{x}$ are infinitely apart\n",
        "*   1 if $x$ and $\\tilde{x}$ are the same\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IbAy8esgp5Q5",
        "colab_type": "text"
      },
      "source": [
        "Fit the best $\\alpha$ coefficients:\n",
        "\n",
        "\\begin{equation}\n",
        "\\mathbf{\\hat{\\alpha}} = \\text{arg min} || \\mathbf{y}^\\text{ref} - \\mathbf{K}^\\mathrm{train}\\mathbf{\\alpha}||^2\n",
        "\\end{equation}\n",
        "Where $\\mathbf{y}^\\mathrm{ref}$ are the training lables and $\\mathbf{K}^\\mathrm{train}$ is the pairwise kernel matrix for all the points in the training set, defined as:\n",
        "\\begin{equation}\n",
        "\\mathbf{K}^\\mathrm{train}_{ij} = \\kappa\\left(\\mathbf{x}_i, \\mathbf{x}_j \\right)\n",
        "\\end{equation}\n",
        "\n",
        "Closed-form solution for regression coefficients:\n",
        "\n",
        "\\begin{equation}\n",
        "\\alpha = \\left(\\mathbf{K}^\\mathrm{train} + \\mathbf{I}\\lambda\\right)^{-1}\\mathbf{y}^\\text{ref}\n",
        "\\end{equation}\n",
        "\n",
        "$\\lambda$ is a small number to be added to the digonal for numerical reasons (regularization and numerical stability).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VJREUUq0tHVJ",
        "colab_type": "text"
      },
      "source": [
        "First, let's define the kernel function"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zxmUjYrNjR08",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def kernel(xi, xj):\n",
        "\n",
        "  sigma = 0.3\n",
        "  k = np.exp(-np.linalg.norm(xi - xj)**2 / (2 * sigma**2))\n",
        "  \n",
        "  return k"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DdrxS1ZZtLoa",
        "colab_type": "text"
      },
      "source": [
        "Next, let's calculate the kernel matrix:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TKBB2xPwtD5k",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "K = np.zeros((len(x),len(x)))\n",
        "\n",
        "for i in range(len(x)):\n",
        "  for j in range(len(x)):\n",
        "    K[i,j] = kernel(x[i], x[j])\n",
        "\n",
        "np.set_printoptions(linewidth=666)\n",
        "print(K)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bskck9aIsK_M",
        "colab_type": "text"
      },
      "source": [
        "With the Kernel matrix above, we can now find the regression coefficients:\n",
        "\\begin{equation}\n",
        "\\alpha = \\left(\\mathbf{K}^\\mathrm{train} + \\mathbf{I}\\lambda\\right)^{-1}\\mathbf{y}^\\text{ref}\n",
        "\\end{equation}\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gD_ZO-tFlXtJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "alphas = np.matmul(np.linalg.inv(K + np.eye(11)*1e-10), y)\n",
        "print(alphas)\n",
        "print(len(alphas))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LgI608JesRnS",
        "colab_type": "text"
      },
      "source": [
        "Another way to solve the equation is via a so-called Cholesky decomposition which is more numerically stable:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E-0dEmyHsoYc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import qml.math\n",
        "\n",
        "# Add to diagonal\n",
        "K[np.diag_indices_from(K)] += 1e-10\n",
        "\n",
        "# Solve\n",
        "alphas = qml.math.cho_solve(K, y)\n",
        "\n",
        "print(alphas)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XXaO4U3TB--Q",
        "colab_type": "text"
      },
      "source": [
        "### Predictions:\n",
        "Remembering again:\n",
        "\n",
        "\\begin{equation}\n",
        "\\hat{y}\\left(\\tilde{\\mathbf{x}}\\right) = \\sum_i \\kappa\\left(\\mathbf{x}_i, \\tilde{\\mathbf{x}} \\right) \\alpha_i\n",
        "\\end{equation}\n",
        "In matrix form:\n",
        "\\begin{equation}\n",
        "\\mathbf{y} = \\mathbf{K}\\mathbf{\\alpha}\n",
        "\\end{equation}\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9wKl_636l-43",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Test points\n",
        "x_test = np.random.random(size=(20,1)) * 6 \n",
        "# x_test = np.arange(-5, 11.5, 0.5)\n",
        "\n",
        "\n",
        "# Zero kernel\n",
        "K_test = np.zeros((len(x_test),len(x)))\n",
        "\n",
        "# Calculate pair-wise Gaussian kernel function\n",
        "for i in range(len(x_test)):\n",
        "  for j in range(len(x)):\n",
        "    K_test[i,j] = kernel(x_test[i],  x[j])\n",
        "\n",
        "# Make prediction\n",
        "y_test = np.dot(K_test, alphas)\n",
        "\n",
        "\n",
        "# Plot everything\n",
        "plt.scatter(x,y, label=\"Training\")\n",
        "plt.plot(xplot, np.sin(xplot), color=\"g\", label=\"True curve\")\n",
        "plt.scatter(x_test, y_test, color=\"r\", label=\"Test\")\n",
        "plt.grid(True)\n",
        "plt.legend()\n",
        "plt.show()\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZcIVOmKJPO1N",
        "colab_type": "text"
      },
      "source": [
        "## Note on Hyperparameters!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p4wxCGxIF-LM",
        "colab_type": "text"
      },
      "source": [
        "## Everything repeated, but with Scikit-learn:\n",
        "\n",
        "Everything we've just coded is of course all available in Scikit-learn:\n",
        "\n",
        "https://scikit-learn.org/stable/modules/generated/sklearn.kernel_ridge.KernelRidge.html"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jsbn6R7knbpL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from sklearn.kernel_ridge import KernelRidge\n",
        "\n",
        "# Gamma instead of sigma!\n",
        "gam = 1 / (2.0 * sigma**2)\n",
        "\n",
        "# Make the KRR object\n",
        "krr = KernelRidge(alpha=1e-10, kernel=\"rbf\", gamma=gam)\n",
        "# krr = KernelRidge(alpha=1e-10, kernel=\"linear\")\n",
        "\n",
        "# Fit the machine\n",
        "krr.fit(x_reshape, y)\n",
        "\n",
        "# Prediction\n",
        "y_scikit = krr.predict(x_test)\n",
        "\n",
        "# Plot everything\n",
        "plt.scatter(x,y, label=\"Training\")\n",
        "plt.plot(xgauss, np.sin(xgauss), color=\"g\", label=\"True curve\")\n",
        "plt.scatter(x_test, y_scikit, color=\"r\", label=\"Scikit-learn\")\n",
        "\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q7LPV4bT0dtf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HOt3TAiiHEKF",
        "colab_type": "text"
      },
      "source": [
        "What about hyper parameters?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FDmEmTt2osrF",
        "colab_type": "text"
      },
      "source": [
        "## Part 2: Molecules in and Machine Learning?\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5n75MYKOLltw",
        "colab_type": "text"
      },
      "source": [
        "### Molecule data files"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vNcRITNkLo3D",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Download Coordinate file for ethanol\n",
        "!wget -O ethanol.xyz https://raw.githubusercontent.com/andersx/ml-intro/master/ethanol.xyz"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QRY2Cg2XLoUI",
        "colab_type": "text"
      },
      "source": [
        "Molecules are often stored in the \".xyz\" format. For example, a conformation for ethanol is stored in the above file:\n",
        "\n",
        "```\n",
        "9\n",
        "\n",
        "C         -0.69924       -0.01497        0.32426\n",
        "C         -1.90198       -0.93192        0.56363\n",
        "O          0.50350       -0.72657        0.14397\n",
        "H          0.30059       -1.64320       -0.17552\n",
        "H         -2.08482       -1.56542       -0.33061\n",
        "H         -2.80809       -0.31867        0.75417\n",
        "H         -1.71524       -1.58372        1.44281\n",
        "H         -0.89079        0.62463       -0.56368\n",
        "H         -0.58049        0.63584        1.21559\n",
        "```\n",
        "\n",
        "Such files can for example be visualized with ASE (Atomic Simulation Environment).\n",
        "https://wiki.fysik.dtu.dk/ase/"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AkWjcUR92Weg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip install ase"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-r-sbAvPMJCK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import ase.io\n",
        "import ase.visualize\n",
        "\n",
        "ethanol = ase.io.read(\"ethanol.xyz\")\n",
        "ase.visualize.view(ethanol, viewer=\"x3d\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xT-O7OtWM6Y2",
        "colab_type": "text"
      },
      "source": [
        "# Making representations for molecules?\n",
        "\n",
        "We need to calculate the pair wise kernels between molecules:\n",
        "\\begin{equation}\n",
        "\\kappa\\left(\\mathbf{x}, \\tilde{\\mathbf{x}} \\right) = \\exp \\left(-\\frac{\\|\\mathbf{x} - \\tilde{\\mathbf{x}} \\|^2}{2\\sigma^2} \\right)\n",
        "\\end{equation}\n",
        "\n",
        "Also sometimes called the *similarity kernel*: 1 if the two molecules are identical, 0 if they are completely different.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yhEY-a8iNphG",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "Problem: how do we get an input representation? Which features are good?\n",
        "\n",
        "Anything could in principle be used:\n",
        "*   Numbers of atoms\n",
        "*   Physical observables\n",
        "*   Name string (e.g. SMILES)\n",
        "*   Coordinates\n",
        "*   Nuclear charges\n",
        "*   Bonding information\n",
        "*   Etc ...\n",
        "\n",
        "Some desireable properties:\n",
        "*   Rotational, translations invariance?\n",
        "*   Permutational invariance\n",
        "*   Uniqueness (injectivity)\n",
        "\n",
        "This kind of problem in machine learning is called \"*feature engineering*\".\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rolBX0MjN_Mv",
        "colab_type": "text"
      },
      "source": [
        "There are *many* representations available for molecules. Two examples are follow, which will also be used in the exercies for today."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9-d5GObuo339",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip install qml"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qdb6EYL_hPZf",
        "colab_type": "text"
      },
      "source": [
        "### The Coulomb Matrix:\n",
        "Proposed by Rupp et al. (2012) Phys Rev Lett https://doi.org/10.1103/PhysRevLett.108.058301\n",
        "\n",
        "\\begin{equation}\n",
        "x_{ij}^\\text{CM} = \n",
        "  \\begin{cases}\n",
        "    0.5Z_i^{2.4} & \\text{for } i=j \\\\\n",
        "    \\frac{Z_i Z_j}{R_{ij}} & \\text{for } i \\neq j \n",
        "  \\end{cases}\n",
        "\\end{equation}\n",
        "\n",
        "The above takes care of translational and permutational invariance.\n",
        "\n",
        "Furthermore, the rows and columns are sorted to ensure permutational invariance.\n",
        "\n",
        "Example for ethane:\n",
        "\n",
        "\n",
        "\n",
        "```\n",
        "73.51669472 \n",
        "34.06515641 36.8581052\n",
        "19.58838249 23.5104435  36.8581052\n",
        " 8.06700497  3.03799493  2.46941255  0.5\n",
        " 3.91313164  5.40535191  2.78865961  0.35565997  0.5\n",
        " 3.87121997  5.40076757  2.76284181  0.38595636  0.55366084  0.5\n",
        " 2.9519494   2.75461693  5.40414993  0.38673507  0.39949749  0.32304269  0.5\n",
        " 2.89651793  2.75224092  5.4003281   0.41811064  0.32445499  0.39915951  0.55199418  0.5\n",
        " 2.35852226  2.76046537  5.40252013  0.28533493  0.40534695  0.39832789  0.5530945   0.55433763  0.5\n",
        "```\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4TXpPyqmSuWu",
        "colab_type": "text"
      },
      "source": [
        "Example of how to generate with QML:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Qbirf6c-SyUx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import qml\n",
        "mol = qml.Compound(xyz=\"ethanol.xyz\")\n",
        "mol.generate_coulomb_matrix(size=12)\n",
        "\n",
        "np.set_printoptions(linewidth=100)\n",
        "print(mol.representation)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vhez3GIMU1jN",
        "colab_type": "text"
      },
      "source": [
        "### Bag-of-Bonds:\n",
        "Proposed by Hansen et al. (2015) J Phys Chem Lett https://doi.org/10.1021/acs.jpclett.5b00831\n",
        "\n",
        "Same content as the Coulomb Matrix, but items are grouped differently, so only identical terms will be compared.\n",
        "\n",
        "More accurate for molecules? (Exercise of today)\n",
        "\n",
        "<img src=\"https://i.imgur.com/eiTZ5g7.png\">"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FfFN-C3nUxXR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "mol.generate_bob(asize={\"O\":2, \"C\": 3, \"H\": 8})\n",
        "np.set_printoptions(linewidth=100)\n",
        "print(mol.representation)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M4jcBu0Vb0mC",
        "colab_type": "text"
      },
      "source": [
        "## Exercises set #2: The QM7 dataset.\n",
        "Rupp et al. (2012) Phys Rev Lett https://doi.org/10.1103/PhysRevLett.108.058301\n",
        "\n",
        "The QM7 dataset contains XYZ structures for 7101 small molecules, with up to 7 atoms of type CNO, saturated with H.\n",
        "\n",
        "Some of them look like this:\n",
        "![alt text](https://i.imgur.com/vQymf39.png)\n",
        "\n",
        "Attempt to map all organic molecules with up to 7 CNO-atoms.\n",
        "\n",
        "For each molecule you will be given the atomization energy calculated using a QM method (PBE0/def2-TZVP).\n",
        "\n",
        "\n",
        "*   Instead of the raw 7101 XYZ files, you will be given the Coulomb Matrix features and Bag-of-Bonds features for each molecule, along with the atomization energies."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o8NJYnarbwGV",
        "colab_type": "text"
      },
      "source": [
        "## How to determine which machine learning method is better than the other?\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4mCN3IbPhRJt",
        "colab_type": "text"
      },
      "source": [
        "For a given training/test split, one could calculate, for example, the mean-absolute-error (MAE) or root-mean-squared-error (RMSE) of the predictions, compared to the true values.\n",
        "\n",
        "\\begin{equation}\n",
        "\\text{MAE} =\\frac{1}{N} \\sum_{i=1}^{N}|y_i - \\tilde{y}_i|\n",
        "\\end{equation}\n",
        "\n",
        "\\begin{equation}\n",
        "\\text{RMSE} = \\sqrt{\\frac{1}{N}\\sum_{i=1}^{N}\\left(y_i - \\tilde{y}_i\\right)^2}\n",
        "\\end{equation}\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WZM7Qw4I720a",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "y_true = np.array([0.00, 1.32, 2.64, 3.96])\n",
        "y_pred = np.array([0.30, 1.14, 3.18, 2.70])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oYTotr_u8cQm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "diff = y_true-y_pred\n",
        "print(diff)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q6h5NMDK9xBo",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Mean-absolute-error\n",
        "mae = np.mean(np.abs(diff))\n",
        "print(mae)\n",
        "\n",
        "# Root-mean-squared-error\n",
        "rmse = np.sqrt(np.mean(diff**2))\n",
        "print(rmse)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5XkTs7GBih44",
        "colab_type": "text"
      },
      "source": [
        "## Learning Curves:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ameti8vx6a-6",
        "colab_type": "text"
      },
      "source": [
        "What is learning? Able to make better prediction with more training data.\n",
        "\n",
        "The error decays according to a power law with the training set size:\n",
        "\\begin{equation}\n",
        "\\mathrm{Error} \\propto \\frac{a}{N^b}\n",
        "\\end{equation}\n",
        "Which is the same as:\n",
        "\\begin{equation}\n",
        "\\log\\left(\\mathrm{Error}\\right) \\propto \\log\\left(a\\right) - b \\log\\left(N\\right)\n",
        "\\end{equation}\n",
        "On a log-log scale the error is linear with the training set size!\n",
        "\n",
        "It turns out, for the same dataset, these are most always parallel (same $b$-value). \n",
        "\n",
        "<img src=\"https://i.imgur.com/4xt7TMj.png\" width=\"500\">\n",
        "\n",
        "We can compare machine learning models (for example based on different representations), by looking at the learning curve!\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cRoWAzIt5s8P",
        "colab_type": "text"
      },
      "source": [
        "## Training, Test, and Validation Splits:\n",
        "\n",
        "One common way to split data is as follows:\n",
        "\n",
        "![alt text](https://cdn-images-1.medium.com/freeze/max/1000/1*ZiYvylk60EY2XG7ck1lqJA.png)\n",
        "\n",
        "This is a method to avoid overfitting (i.e. fitting parameters to your \"test\" data)\n",
        "\n",
        "Procedure for calculating the error for a model:\n",
        "1.   Train model on **Training Set**\n",
        "2.   Minimize error for **Validation Set**, i.e. optimize hyperparameters + goto 1.\n",
        "3.   Calculate error for **Test Set**\n",
        "\n",
        "The last error on the Test Set (from 3.) is your error.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yIe6oQ6pI6dj",
        "colab_type": "text"
      },
      "source": [
        "Example how to split 100 values into three\n",
        "\n",
        "*   50 Training Set\n",
        "*   10 Validation Set\n",
        "*   40 Test Set\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kBemuPD7ILyO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "data_set = np.arange(0,100)\n",
        "print(data_set)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "duSiG-vlJqaB",
        "colab_type": "text"
      },
      "source": [
        "Remembering the Numpy *slice-notation*:\n",
        "\n",
        "The general syntax is\n",
        "\n",
        "`array[row]`\n",
        "\n",
        "or\n",
        "\n",
        "`array[row, column]`\n",
        "\n",
        "Special notation:\n",
        "\n",
        "    n   = what is in index n\n",
        "    :n  = up to n\n",
        "    n:  = n and onwards\n",
        "    n:m = from n to m \n",
        "    :   = everything\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bOo4JSZWJMl2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# First 50\n",
        "training_set = data_set[:50]\n",
        "print(training_set)\n",
        "# Next 10\n",
        "validation_set = data_set[50:60]\n",
        "print(validation_set)\n",
        "# Next 40\n",
        "test_set = data_set[60:100]\n",
        "print(test_set)"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}