{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPyOpt: Integrating the model hyperparameters\n",
    "\n",
    "### Written by Javier Gonzalez and Zhenwen Dai, University of Sheffield.\n",
    "\n",
    "## Reference Manual index\n",
    "\n",
    "*Last updated Friday, 11 March 2016.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrating the model hyperparameters of a GP is easy with GPyOpt. To see how to use this feature we start loading GPyOpt. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline  \n",
    "import GPyOpt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We minimize the use the Six-hump camel function \n",
    "\n",
    "$$f(x_1,x_2) = \\left(4-2.1x_1^2 = \\frac{x_1^4}{3} \\right)x_1^2 + x_1x_2 + (-4 +4x_2^2)x_2^2,$$\n",
    "\n",
    "in $[-3,3]\\times [-2,2]$. This functions has two global minimum, at $(0.0898,-0.7126)$ and $(-0.0898,0.7126)$. As in the previous case we create the function, which is already in GPyOpt. In this case we generate observations of the function perturbed with white noise of $sd=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# create the object function\n",
    "f_true = GPyOpt.objective_examples.experiments2d.sixhumpcamel()\n",
    "f_sim = GPyOpt.objective_examples.experiments2d.sixhumpcamel(sd = 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in previous cases we create the bounds and plot the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bounds =[{'name': 'var_1', 'type': 'continuous', 'domain': f_true.bounds[0]},\n",
    "         {'name': 'var_2', 'type': 'continuous', 'domain': f_true.bounds[1]}]\n",
    "f_true.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we are going to use **Hibrid Monte Carlo Sampling** over the parameters of the GP. If ``model_type = 'GP_MCMC'`` is selected, samples from the hyperparameters are generated when the model is updated (instead of doing the default MLE update). To  use the samples to integrate Expected improvement we need to select the option ``acquisition_type='EI_MCMC'``. Similar syntaxis works for the rest of available acquisitions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Four parameters can be controlled in the sampling generation (apart from the hyperparameter priors that can be directly changed in the model):\n",
    "\n",
    "* ``n_samples``: number of generated samples.\n",
    "* ``n_burning ``: number of initial samples not used.\n",
    "* ``subsample_interval``: interval of subsampling from HMC samples.\n",
    "* ``step_size``: size of the gradient steps in the HMC sampler.\n",
    "* ``leapfrog_steps``: number of gradient steps before each Metropolis Hasting step. \n",
    "\n",
    "All this attributes are handled as extra arguments to the ``BayesianOptimization`` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This make take some time, in the initialization samples form the hyperparametes are generated\n",
    "myBopt_mcmc = GPyOpt.methods.BayesianOptimization(f_sim.f,\n",
    "                                               domain=bounds,\n",
    "                                               model_type= 'GP_MCMC',\n",
    "                                               acquisition_type='EI_MCMC',                                                \n",
    "                                               normalize_Y = True,\n",
    "                                               n_samples=5) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This make take some time,\n",
    "max_iter = 10\n",
    "myBopt_mcmc.run_optimization(max_iter)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot again the samples generated in the last iteration to see how the optimization affected to out beliefs about the model hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "labels = ['rbf variance', 'rbf lengthscale','noise variance']\n",
    "\n",
    "plt.figure(figsize=(14,5))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(myBopt_mcmc.model.hmc_samples,label = labels)\n",
    "plt.title('Hyperparameters samples',size=25)\n",
    "plt.xlabel('Sample index',size=15)\n",
    "plt.ylabel('Value',size=15)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "from scipy import stats\n",
    "xmin = myBopt_mcmc.model.hmc_samples.min()\n",
    "xmax = myBopt_mcmc.model.hmc_samples.max()\n",
    "xs = np.linspace(xmin,xmax,100)\n",
    "for i in range(myBopt_mcmc.model.hmc_samples.shape[1]):\n",
    "    kernel = stats.gaussian_kde(myBopt_mcmc.model.hmc_samples[:,i])\n",
    "    plot(xs,kernel(xs),label=labels[i])\n",
    "_ = legend()\n",
    "plt.title('Hyperparameters densities',size=25)\n",
    "plt.xlabel('Value',size=15)\n",
    "plt.ylabel('Frequency',size=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, of course, you can check how the optimization evolved when you integrate out the acquisition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myBopt_mcmc.plot_convergence()"
   ]
  }
 ],
 "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": 0
}