{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPyOpt: Modular Bayesian Optimization \n",
    "\n",
    "### Written by Javier Gonzalez, Amazon Reseach Cambridge\n",
    "\n",
    "\n",
    "*Last updated, July 2017.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [Introduction Bayesian Optimization GPyOpt](./GPyOpt_reference_manual.ipynb) we showed how GPyOpt can be used to solve optimization problems with some basic functionalities. The object \n",
    "\n",
    "```\n",
    "GPyOpt.methods.BayesianOptimization\n",
    "```\n",
    "is used to initialize the desired functionalities, such us the acquisition function, the initial design or the model. In some cases we want to have control over those objects and we may want to replace some element in the loop without having to integrate the new elements in the base code framework. This is now possible through the modular implementation of the package using the\n",
    "\n",
    "```\n",
    "GPyOpt.methods.ModularBayesianOptimization\n",
    "```\n",
    "\n",
    "class. In this notebook we are going to show how to use the backbone of GPyOpt to run a Bayesian optimization algorithm in which we will use our own acquisition function. In particular we are going to use the Expected Improvement integrated over the jitter parameter. That is\n",
    "\n",
    "$$acqu_{IEI}(x;\\{x_n,y_n\\},\\theta) = \\int acqu_{EI}(x;\\{x_n,y_n\\},\\theta,\\psi) p(\\psi;a,b)d\\psi $$\n",
    "where $p(\\psi;a,b)$ is, in this example, the distribution [$Beta(a,b)$](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.random.beta.html).\n",
    "\n",
    "This acquisition is not available in GPyOpt, but we will implement and use in this notebook. The same can be done for other models, acquisition optimizers etc.\n",
    "\n",
    "As usual, we start loading GPy and GPyOpt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import GPyOpt\n",
    "import GPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will use the Branin function as a test case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# --- Function to optimize\n",
    "func  = GPyOpt.objective_examples.experiments2d.branin()\n",
    "func.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we are won't use the pre implemented wrapper, we need to create the classes for each element of the optimization. In total we need to create:\n",
    "\n",
    "* Class for the **objective function**,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "objective = GPyOpt.core.task.SingleObjective(func.f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Class for the **design space**,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "space = GPyOpt.Design_space(space =[{'name': 'var_1', 'type': 'continuous', 'domain': (-5,10)},\n",
    "                                    {'name': 'var_2', 'type': 'continuous', 'domain': (1,15)}])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Class for the **model type**,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = GPyOpt.models.GPModel(optimize_restarts=5,verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Class for the **acquisition optimizer**,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aquisition_optimizer = GPyOpt.optimization.AcquisitionOptimizer(space)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Class for the **initial design**,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "initial_design = GPyOpt.experiment_design.initial_design('random', space, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Class for the **acquisition function**. Because we want to use our own acquisition, we need to implement a class to handle it. We will use the currently available Expected Improvement to create an integrated version over the jitter parameter. Samples will be generated using a beta distribution with parameters a and b as it is done using the default [numpy beta function](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.random.beta.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from GPyOpt.acquisitions.base import AcquisitionBase\n",
    "from GPyOpt.acquisitions.EI import AcquisitionEI\n",
    "from numpy.random import beta\n",
    "\n",
    "class jitter_integrated_EI(AcquisitionBase):\n",
    "    \n",
    "    analytical_gradient_prediction = True\n",
    "    \n",
    "    def __init__(self, model, space, optimizer=None, cost_withGradients=None, par_a=1, par_b=1, num_samples= 10):\n",
    "        super(jitter_integrated_EI, self).__init__(model, space, optimizer)\n",
    "        \n",
    "        self.par_a = par_a\n",
    "        self.par_b = par_b\n",
    "        self.num_samples = num_samples\n",
    "        self.samples = beta(self.par_a,self.par_b,self.num_samples)\n",
    "        self.EI = AcquisitionEI(model, space, optimizer, cost_withGradients)\n",
    "    \n",
    "    def acquisition_function(self,x):\n",
    "        acqu_x = np.zeros((x.shape[0],1))       \n",
    "        for k in range(self.num_samples):\n",
    "            self.EI.jitter = self.samples[k]\n",
    "            acqu_x +=self.EI.acquisition_function(x)           \n",
    "        return acqu_x/self.num_samples\n",
    "    \n",
    "    def acquisition_function_withGradients(self,x):\n",
    "        acqu_x      = np.zeros((x.shape[0],1))       \n",
    "        acqu_x_grad = np.zeros(x.shape)\n",
    "        \n",
    "        for k in range(self.num_samples):\n",
    "            self.EI.jitter = self.samples[k]       \n",
    "            acqu_x_sample, acqu_x_grad_sample =self.EI.acquisition_function_withGradients(x) \n",
    "            acqu_x += acqu_x_sample\n",
    "            acqu_x_grad += acqu_x_grad_sample           \n",
    "        return acqu_x/self.num_samples, acqu_x_grad/self.num_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we initialize the class for this acquisition and we plot the histogram of the used samples to integrate the acquisition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "acquisition = jitter_integrated_EI(model, space, optimizer=aquisition_optimizer, par_a=1, par_b=10, num_samples=200)\n",
    "xx = plt.hist(acquisition.samples,bins=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Finally we create the class for the **type of evaluator**,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# --- CHOOSE a collection method\n",
    "evaluator = GPyOpt.core.evaluators.Sequential(acquisition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With all the classes on place,including the one we have created for this example, we can now create the **Bayesian optimization object**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bo = GPyOpt.methods.ModularBayesianOptimization(model, space, objective, acquisition, evaluator, initial_design)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we run the optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "max_iter  = 10                                            \n",
    "bo.run_optimization(max_iter = max_iter) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the acquisition and the diagnostic plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bo.plot_acquisition()\n",
    "bo.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
}