GeneticAlgorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Genetic Algorithm\n",
    "\n",
    "GA is a simple optimisation algorithm that mimics evolution. It starts with a random population of elements, and evaluates how <i>fit for survival</i> they are. Elements are then randomly selected in pairs to mate and produce a new generation of offspring. The selection ensures that elements with higher fitness are more likely to be chosen for mating, thus increasing the chances that their offspring retains the characteristics that made their parents fit in the first place. Random mutation can also appear.\n",
    "\n",
    "If an optimisation problem is too complex to be solved via gradient methods, GA can still give an answer.\n",
    "\n",
    "In this tutorial we are going to use GA to find a set of parameters that fit a signal.<br>\n",
    "More details about the GA are given in the following sections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IMPORTS\n",
    "import os\n",
    "import math\n",
    "import random\n",
    "import numpy\n",
    "from GA import GAEngine\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of the problem\n",
    "We are given an compound oscillatory signal $W(t)$ at given times and we would like to find its components. We are too primitive to use FFT or non-linear fitting. We are dumb and strong, so we brute-force it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Definition of element\n",
    "\n",
    "Each $i$-th element of the genetic population represents a possible candidate signal $W_i(t) = \\sum_{k=1}^n A_k \\cos\\left( 2\\pi k t\\right)$. Thus, the element can be represented by the list of amplitudes for each wave component: $\\mathbf{e} = \\{ A_1, A_2, ... A_n\\}$.\n",
    "\n",
    "The training set is a collection of $\\left( t,W(t) \\right)$ points corresponding to the signal we are trying to fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set a decent population size to work with\n",
    "populationSize = 100\n",
    "\n",
    "# set the amount of plane waves\n",
    "nWaves = 5\n",
    "\n",
    "# limit amplitudes to [-2, 2]\n",
    "scale = 2.0\n",
    "\n",
    "# create a GA engine\n",
    "ga = GAEngine(popSize=populationSize, dnaSize=nWaves, scale=scale)\n",
    "\n",
    "ga.tau = 0.1\n",
    "\n",
    "# print the first 3 elements to see what we are dealing with!\n",
    "print(ga.population[0:3])\n",
    "\n",
    "# load a training set\n",
    "trainT = numpy.load(\"./data/ga-train-T.npy\") # time values\n",
    "trainY = numpy.load(\"./data/ga-train-Y.npy\") # corresponding signal values\n",
    "\n",
    "# plot the training set\n",
    "plt.plot(trainT, trainY, 'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##### Mixing elements\n",
    "\n",
    "After the population is evaluated, all elements are sorted by descending fitness (best fit first). Pairs of elements are selected at random by their index, using an exponential distribution, and mixed to produce an offspring population.\n",
    "The simplest mixing function would be the one that assign each DNA component of the offspring from either parent at \n",
    "random. This function is already defined in the GAEngine.\n",
    "\n",
    "We can control the PDF for the selection using one parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ga.tau = 1\n",
    "selected = ga.TrySelection(10000)\n",
    "plt.hist(selected, populationSize, density=True)\n",
    "plt.xlabel('element index')\n",
    "plt.ylabel('Probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Large values of $\\tau$ will cause elements with low fitness to be almost never selected. This will make the GA \"converge\" faster, halting evolution.\n",
    "\n",
    "Small values of $\\tau$ will eventually give equal chance to all elements regardless of their fitness. This will slow down the evolution.\n",
    "\n",
    "You will have to find a good compromise!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##### Mutating elements\n",
    "\n",
    "Occasionally, mutation can appear when mixing two elements. There are several ways to implement mutations, but here we stick to simplicity. A mutation adds a random value to one random gene of the element. <br>\n",
    "Mutations should be rare or they will mess up the evolution!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the mutation rate\n",
    "ga.mutationRate = 0.01 # this means 1% of the elements will get it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##### Evaluating fitness\n",
    "\n",
    "This function takes in the element descriptor, performs the calculations necessary to estimate how well it performs, and returns its fitness.<br>\n",
    "It is actually easy to calculate its badness (mean square error) and flip the sign. This way, the perfect element will have 0 fitness, while any mismatch with the training will give a negative value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EvaluateFitness(element):\n",
    "\n",
    "    # prepare a 0 array\n",
    "    y = numpy.zeros(len(trainT))\n",
    "    freq = 1 # starting frequency\n",
    "    \n",
    "    for amp in element: # loop through the amplitudes\n",
    "        y += amp * numpy.cos(2*numpy.pi*freq*trainT) # add the wave\n",
    "        freq += 1 # increment the frequency for the next wave\n",
    "    \n",
    "    # compute square error\n",
    "    y = (y-trainY)*(y-trainY)\n",
    "    \n",
    "    return -numpy.mean(y)\n",
    "\n",
    "# tell the engine to use this evaluation function\n",
    "ga.Evaluate = EvaluateFitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Evolving...\n",
    "\n",
    "Let the population evolve!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of generations to compute\n",
    "nGens = 30\n",
    "\n",
    "stats = numpy.zeros((nGens,3))\n",
    "\n",
    "# start the evolution loop!\n",
    "for g in range(nGens):\n",
    "    stats[g] = ga.Evolve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evolution plot\n",
    "plt.plot(stats[:,0],'o-', label='best')\n",
    "plt.plot(stats[:,1],'o-', label='average')\n",
    "plt.plot(stats[:,2],'o-', label='worst')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# best element info\n",
    "print(\"best model amplitudes:\", ga.best)\n",
    "print(\"best model fitness:\", ga.bestFit)\n",
    "t = numpy.arange(0,1,0.01)\n",
    "y = t * 0\n",
    "freq = 1 # starting frequency\n",
    "\n",
    "for amp in ga.best: # loop through the amplitudes\n",
    "    y += amp * numpy.cos(2*numpy.pi*freq*t) # add the wave\n",
    "    freq += 1 # increment the frequency for the next wave\n",
    "        \n",
    "plt.plot(trainT, trainY, 'o', label='training points')\n",
    "plt.plot(t, y, '-', label='best model')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Final remarks\n",
    "Do not be fooled by the speed of this example... in real applications, training an AI with GA can take up to several weeks on supercomputers.<br>\n",
    "You can see one in action <a href=\"http://www.nanolayers.com/nano_copter.php\">here!</a> Try to beat it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Beast Exercise\n",
    "\n",
    "Train neural networks to output the energy of a molecule using ACSFs as input. \n",
    "There should be one network for each atom type and ACSF contribution. For example, all H atoms will have 5 NNs that process the ACSFs calculated for that atom with respect to C, N, O, or F atoms.\n",
    "Additionally there will be NNs for the 3-body parts: H-HH, H-HC, H-HN, ... total 15 NNs\n",
    "In total, each atom type should have 5 + 15 = 20 NNs, giving a total of 100 NNs and lots of parameters to train!\n",
    "\n",
    "Some serious python hacking is required!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load structures and energies from the database\n",
    "# the file is data/structures.xyz and it contains lots of molecules in the following xyz format:\n",
    "# \n",
    "# natoms_mol1\n",
    "# energy_mol1\n",
    "# type1 x y z\n",
    "# type2 x y z\n",
    "# ...\n",
    "# natoms_mol2\n",
    "# energy_mol2\n",
    "# type1 x y z\n",
    "# ...\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use DScribe to compute ACSFs\n",
    "# ...\n",
    "\n",
    "# the ACSFs should be divided into their logical contributions (x-H, x-C, ... x-HH, x-HC, ...)\n",
    "# the length of each contribution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPRegressor\n",
    "\n",
    "# setup model networks for all possible atoms\n",
    "# there are only H, C, N, O, F atoms in the database\n",
    "#\n",
    "# H -> H (2-body)\n",
    "# H -> C (2-body)\n",
    "# ...\n",
    "# H -> H-H (3-body)\n",
    "# H -> H-C (3-body) (no need for H -> C-H, it is symmetric!)\n",
    "# ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sklearn network weights and neuron biases are stored in\n",
    "# model.coefs_ and model.intercepts_\n",
    "# these are assigned only after calling model.fit(x, y), does not matter if the training succeeded\n",
    "# as long as x,y have the correct shape\n",
    "#\n",
    "# the GA DNA is the concatenation of all these numbers\n",
    "# \n",
    "# calculate the total length of the DNA\n",
    "# ...\n",
    "\n",
    "# write a function that converts the DNA (list of numbers) into a set of model.coefs_ and \n",
    "# model.intercepts_ for a sklearn model\n",
    "# ...\n",
    "\n",
    "# write an evaluation function for the fitness of a GA element\n",
    "# ...\n",
    "\n",
    "# evolve the population\n",
    "# ..."
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}