NeuralNetwork - TotalEnergy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Total Energy Prediction - Neural Network\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this notebook we will machine-learn the relationship between a molecular descriptor and total energy using neural networks.\n",
    "\n",
    "The energy of ~134k molecules was calculated at the CCSD level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- INITIAL DEFINITIONS ---\n",
    "from sklearn.neural_network import MLPRegressor\n",
    "import numpy, math, random\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.sparse import load_npz\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pick a descriptor. Allowed types are:\n",
    "\n",
    "1. cnt: atom counts\n",
    "2. bob: bag of bonds\n",
    "4. soap: smooth overlap of atomic positions; choose from:\n",
    "    1. soap.sum - all atoms summed together\n",
    "    4. soap.mean - mean of all atom SOAP\n",
    "    4. soap.centre - computed at the central point\n",
    "5. mbtr: many-body tensor representation\n",
    "6. cm: Coulomb matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TYPE is the descriptor type\n",
    "TYPE = \"cm\"\n",
    "\n",
    "#show descriptor details\n",
    "print(\"\\nDescriptor details\")\n",
    "desc = open(\"./data/descriptor.\"+TYPE.split('.')[0]+\".txt\",\"r\").readlines()\n",
    "for l in desc: print(l.strip())\n",
    "print(\" \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and load the databases with the descriptors (input) and the correct charge densities (output). Databases are quite big, so we can decide how many samples to use for training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load input/output data\n",
    "trainIn = load_npz(\"./data/energy.input.\"+TYPE+\".npz\").toarray()\n",
    "trainOut = numpy.load(\"./data/energy.output.npy\")\n",
    "trainIn = trainIn.astype(dtype=numpy.float64, casting='safe')\n",
    "\n",
    "# decide how many samples to take from the database\n",
    "samples  = min(trainIn.shape[0], 9000)\n",
    "vsamples = min(trainIn.shape[0]-samples,1000)\n",
    "print(\"training samples:   \"+str(samples))\n",
    "print(\"validation samples: \"+str(vsamples))\n",
    "print(\"number of features: {}\".format(trainIn.shape[1]))\n",
    "\n",
    "# split between training and validation\n",
    "validIn = trainIn[samples:samples+vsamples]\n",
    "validOut = trainOut[samples:samples+vsamples]\n",
    "\n",
    "trainIn  = trainIn[0:samples]\n",
    "trainOut = trainOut[0:samples]\n",
    "\n",
    "# shift and scale the inputs\n",
    "train_mean = numpy.mean(trainIn, axis=0)\n",
    "train_std = numpy.std(trainIn, axis=0)\n",
    "train_std[train_std==0] = 1\n",
    "for a in range(trainIn.shape[1]):\n",
    "    trainIn[:,a] -= train_mean[a]\n",
    "\n",
    "for a in range(trainIn.shape[1]):\n",
    "    trainIn[:,a] /= train_std[a]\n",
    "# also for validation set\n",
    "for a in range(validIn.shape[1]):\n",
    "    validIn[:,a] -= train_mean[a]\n",
    "for a in range(validIn.shape[1]):\n",
    "    validIn[:,a] /= train_std[a]\n",
    "    \n",
    "\n",
    "# show the first few descriptors\n",
    "print(\"\\nDescriptors for the first 5 molecules:\")\n",
    "print(trainIn[0:5])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we setup a multilayer perceptron of suitable size. Out package of choice is scikit-learn, but more efficient ones are available.<br>\n",
    "Check the scikit-learn <a href=\"http://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPRegressor.html\">documentation</a> for a list of parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup the neural network\n",
    "nn = MLPRegressor(hidden_layer_sizes=(1000,200,50,50),  activation='tanh', solver='lbfgs', alpha=0.01, \n",
    "                  learning_rate='adaptive')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training\n",
    "Now comes the tough part! The idea of training is to evaluate the ANN with the training inputs and measure its error (since we know the correct outputs). It is then possible to compute the derivative (gradient) of the error w.r.t. each parameter (connections and biases). By shifting the parameters in the opposite direction of the gradient, we obtain a better set of parameters, that should give smaller error.\n",
    "This procedure can be repeated until the error is minimised.\n",
    "\n",
    "\n",
    "It may take a while..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use this to change some parameters during training if the NN gets stuck in a bad spot\n",
    "nn.set_params(solver='lbfgs')\n",
    "\n",
    "nn.fit(trainIn, trainOut);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the ANN quality with a regression plot, showing the mismatch between the exact and NN predicted outputs for the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate the training and validation error\n",
    "trainMLOut = nn.predict(trainIn)\n",
    "validMLOut = nn.predict(validIn)\n",
    "\n",
    "print (\"Mean Abs Error (training)  : \", (numpy.abs(trainMLOut-trainOut)).mean())\n",
    "print (\"Mean Abs Error (validation): \", (numpy.abs(validMLOut-validOut)).mean())\n",
    "\n",
    "plt.plot(validOut,validMLOut,'o')\n",
    "plt.plot([min(validOut),max(validOut)],[min(validOut),max(validOut)]) # perfect fit line\n",
    "plt.xlabel('correct output')\n",
    "plt.ylabel('NN output')\n",
    "plt.show()\n",
    "\n",
    "# error histogram\n",
    "plt.hist(validMLOut-validOut,50)\n",
    "plt.xlabel(\"Error\")\n",
    "plt.ylabel(\"Occurrences\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Compare descriptors\n",
    "Keeping the size of the NN constant, test the accuracy of different descriptors with the same NN size, and find the best one for these systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DIY code here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Combine with Principal Component Analysis - Advanced\n",
    "Reduce the descriptor size with PCA (check the PCA.ipynb notebook) and train again. Can you get similar accuracy with much smaller networks?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}