demo_constrained_bo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Optimization Demo\n",
    "*Joachim van der Herten*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Bayesian Optimization, the typical goal is optimizing an *expensive*-to-evaluate function $f: \\mathcal{X} \\rightarrow \\mathbb{R}^p$: \n",
    "$$ \\underset{\\boldsymbol{x}}{\\operatorname{argmin}} f(\\boldsymbol{x})$$\n",
    "over the optimization domain $\\mathcal{X} \\subset \\mathbb{R}^d$ by evaluating $f$ on a sequence of decisions $\\boldsymbol{x}_0, \\boldsymbol{x}_1, ... \\boldsymbol{x}_{n-1}$. Each decision $\\boldsymbol{x}_i \\in \\mathcal{X}$ itself is the result of an optimization problem\n",
    "$$\\boldsymbol{x}_i = \\underset{\\boldsymbol{x}}{\\operatorname{argmax}} \\alpha_i(\\boldsymbol{x})$$\n",
    "The function $\\alpha$ acts as a surrogate for the expensive function, and is cheap-to-evaluate. It is typically a mapping of the predictive distribution of a model built on all decisions and their corresponding (noisy) evaluations. The mapping introduces an order in $\\mathcal{X}$ to obtain a certain goal. The typical goal within the context of Bayesian Optimization is typically to obtain *optimality* or *feasibility* with a small number of evaluations of $f$. As we can have several functions $f$ representing objectives and constraints, Bayesian Optimization may invoke several models and mappings $\\alpha$. The mappings are typically referred to as *acquisition functions*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import rc\n",
    "#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "rc('text', usetex=True)\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import GPflow\n",
    "import OptDemo\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates the optimization of an analytical function using the well known Expected Improvement (EI) function. In addition, a black-box constraint function is added to the optimization. The feasible regions are learnt jointly with the optimal regions by considering a second acquisition function known as the Probability of Feasibility (PoF), following the approach of Gardner et al. (2014)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f9ef7b8a810>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAAGrCAYAAAB0e2vOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3U+MnVWe3vHn9DCeJoFQtuNpqek4obywFBGccVc2rFpK\nNepZkEXikVcsu7IzQqC2smFBFpEJCOHsPEukSFZ6VkRq0tgSK1YFLTO9sSIXCd29SBi7ioaEtGfw\nyaLeF19f3z/vv/Pvd74fqURd31v3HnSr7n3uc855X+e9FwAAAPLxndQDAAAAwIMIaAAAAJkhoAEA\nAGSGgAYAAJAZAhoAAEBmCGgAAACZIaABAABkhoAGAACQGQIaAABAZghoAAAAmXkk9QDGOnrsO/77\nP+j3v3H7m8eCjOX3d7876ue/+cMfTTKOP/rDJHfT/3H/33SnDfvmu26y+5rSN38S/jH+2Z9+L/yD\nAEjir//3/xp9H2Ne46d4nf7O1387+j5yce/RPx78s33fp9r3j7u/+e3feO9PrLt98QHt+z94RP/5\nv3Z/Q3tn/9lgY/nlZ6dH38dXnz4x+j4ev5W+GD16827vn9k/fSTASKb15al7wR/jjqS9Cy8HfxwA\n6WxefnPUzw99nR/y2jzv0V//bvR95OTrp58c/LN937e+PHVPn774yv/sctviAxoe9uWpe8lDWglh\nK2ftizdBDbBp78LLo0PaEPunj4wOaV8//aSpkPbor383OKQdvXm31/tdn/fm9FVLJO/sP5t9ezal\nGE0PwkvxAg4gjr0LLw/+EDbmNZ4P0A/LMXCaD2ihgxnqkqKZJKQBtqUIaWONmRbM1dCQNsW08SKm\nA1qJweyxp76Y7L5o0ezYvPwmQQ0wbEybNgQt2mI5hTRTAa1ty2jN7iOk2UJIA2zrG9Jo0aaXS0gr\nfpPA7W8eI4ytkcOmAUyHDQSAbe3fdugPZFNsGEA4vGtnaMppzhZNmj20aYBtXT+EpdwwQIv2oCkD\nLwENKBhr0wDbuq5N40P49FKHNAJapmjR0AchDbAt5JIGWrTlUoY0AlplCGl20aYBtq1r03h9DyPV\nMdIIaBkL0aJJ/BGPlfuGC0IaYFuINo0WbbUhIW1si5b3Ow2CIaTZRpsG2LYspHHYjbyMCWkEtMyF\natGkwz9kgppthDTArmVTnkNf1zl47Wqx16MR0Cby3Mmbwe47ZEiTaNOso00DbJtyypOpztVirkcj\noEESIa2v3NehLUJIA+yab9N4TQ8n1nq0ZO8yzrlL665zzu3EG1HeQrdoEn/QNaBNA2ybIqTRoq0X\nI6QlCWhN8Dq34iY7zrlbkvYiDakIsUIaQQ0AytW2aWOaftajhdEnpCUJaN77K1odvn7qvT/lvb8W\na0xTCLkOrRUjpEm0aV2UOM0JoB433n4p2WPToo2X6zvMMefctnPuZ6kHkqOYIY2gBgDluvH2S4OD\nGlOd64UMaVkGNO/9laY9O+6c2049nhzFCmkSbRoAlC5lm2ZdqJCWXUBzzu0459r1abclbS65za5z\nbvf/7v8h7gDXiDHN2Yod0ghqDyt1mpONAkB9hrRptGjpZPPu4pzbaL7dk9SuPTsuaXf+tk3DtuW9\n3/p7R/8k1hA7sxrSJNo0S9jRCaALNgysF6JFS7WL85ykrbnDaFyXpGZqc7u5zW3v/ccpxliSFCGN\noHZfqS1ai5AGIKRaWrSpQ1qqXZw/994fbXZztv/2w7nrf+69fz3F+KYQs0WT4oc0iTbNEto0AKsw\n1Rlf2R/9M1dLSCOold+itQhpAJZhqnO9KVs0G+8qGashpEm0aZKtkEZQA+xKtaOzlhZtqpBm4x0l\nczWFNIKaHYQ0wK6hx0djqrObKUIaAS2SWkKaVHebZqVFa9GmAbalCGnoxta7SeZShDTatPishTSJ\nNg2wbMzZBoagRevG3jtJ5mKHNIk2LQVCGoDS9AlpTHWGZ+9dpAA1hrQag5rVkEZQA+yK3aZZN6ZF\ns/cOUojaQppUb5tmESENsK1LSKNF62ZoSCs+oP3+7ndTD2Gw507erGrzgFRfm2axRWvRpgG2dWnT\nCGnh2H33KEhNmwdatYU060ENgF1MeY43pEUz8a7xy89Opx7CaEx52mc9pBHUALtWtWm0aGHYfcco\nUK1TnjWxHNIk2jTAOkLacH1bNDPvFhZatFaNIa2moFZDSCOoAXaF2ulJSHuQ7XeKgsVu01KHNKmu\nNs36ujSJNg2wbj6kcYaBadl+hzCAkGZbDSGNoAbYNd+mMdU5HdvvDkYQ0myzHtIk2jTAuimnPAlp\nh0y9M1hahzYv5pRnDiGtNrVMeRLUALvakDbFVCchzVhAq0EtIa22Fq1lPaRJtGmAZZwqajr23w0M\nIqTZRpsGoHQfvHdx9H3U3qLZfhcwjJBmn/WQJtGmAZZ98N7F0UGt5pBm/x3AMEKafbW0aQDsIqQN\nU/wr/z994vv65PnXUg8jmVibB1KHtNpZD2pMeQK2TdGm1cbMK37NIU2K06alDGk1t2izLIc0iTYN\nsG5oSKuxRXPe+9RjGGVra8vv7u4+8G+v3DifaDTpxTjUyFefPhH8MRaxHk76shpa2+eZnWCAbT/6\nyaXeP9P3fJY5eu93/+kj7/3WutuZfMd748zV1ENIJvZ5PGOyGkiGsj7teebFt1IPAUBAQ9q0mpo0\nkw3arFrbtNBNGi1afqwE2EXPMW0aYFvfNq3kJq3qBm1WrW1a6CYt1Xo0KyEkBAuN2rLx06YBtrGB\n4GFlv5p3REgLg52debIQ1ADUp89OzxqmOnkVN87imjRatG5KC2mljRdAGIS0Q7wiViBkSKNFy1sp\nbVoJYwQQT9c2zXJI41URoxHS8ldKUFuFdWhAfWoOaWW/YvdQ6zq0lsWpTvTXBrWcwlqfsRDSgPrU\nGtLyeZWOgJBmZ6qTdWjj5RbUuiKkAVjEWkgr79V5pDfOXK06qNGkYV7KoDb0cc+8+BZBDYBpj6Qe\nQCpvnLk6yUFsXzj6YefbvrP/7OjHy9ljT32R7AC2GK8NS7HaySlC4ZkX3+IgtgC+9fXTTxZ9ENtZ\nyc4k4Jy75L1fOLHsnDsn6UDSpvf+yqr7WXcmgS6GBLU+wWyR1GEt1JkGYga0EqfnShMirIV63ghq\ngG19zjaQc0jL+kwCzrkdSeeWXHdWkrz312Yvh9R3ynNsOGvvY4r7GYqpTnQx9aaCkKGaKU/Atj5n\nG7CwHi1lg/a+9/7HC/79kqT3vffXnHPbks56719fdj9TNGizVrVpIQNVikaNFg1D9W3WYj9PtGmA\nbV3btBybtKwbtDU2JN2ZuXw85oMva9NCt10p2jRaNAw126x1+YqNNg2wrYazDeQY0JKb3+kZKzyl\nnPKcEgeuRQ7Y6QnYZv1sAzkGtANJx5rvNyTdTjWQFIfjiB3SaNFgHSENsM1qSMsmoDnnNppvr0ra\nbL7flHRtwW13nHO7zrndzz//POi4aghpgHW0aYBtXdu0kqTaxXlO0lazm7N1XZK89x83t9mWdNBe\nnuW9v+K93/Leb504cSL4eM+c/I3OnPxN8MeZFTOkhWjRmOZEjghpgG2rQlppLVqSgOa9/7n3/ujs\nMc689z+c+f6K9/7aumOgxWY5pAG1oE0DbLMS0rKZ4ixF7DYtVkhjLRpqQ0gD7Fo15VlKSCOgDWQx\npAG1oU0DbCs5pBHQRkixNq0krENDKQhpgF3L2rTcQxoBbQIxQlqMFo1pTtSMkAbYVlpII6BNhDYN\nKB9TnoBti9q0XEMaAW1iIUMaa9GAOAhpgG0lhDQCWgAlt2lMcwKHaNMA2+bbtNxCGgEtoBAhrbQW\njY0CKB0hDbAt15BGQAus5DYNwCHaNMC22TYtl5BGQIuEkAaUj5AG2PbBexf16K9/l3oYkghoUU3V\nppU2zQlYQpsG2PaL317WL357Och992nnCGgJ5N6msVEAWI+QBtg2dUjrO3VKQEsk95AGYD3aNMC2\nkG3aOgS0hNhAANhASANsGxvShmw8IKBlgJAGlI82DbBtaEgbuiuUgJYJ2jTABkIaYFfMKU8CWmYs\nhjQOVova0KYBtnUNaWOOqUZAyxBtGmADIQ2wK3SbRkADgIAIaYBtoUIaAQ0AAmPKE7AtRJtGQAOA\nSAhpgG2zIW3sOT0JaAAAABOZqk0joGWMzQIAAJTpg/cujvp5AloB5kPaO/vPJhoJAADo6oP3Lg4O\nagS0QtCkATawDg2oz5CQRkArCFOegA2ENKA+fds0AlqB3jhzNfUQAIzEoTcArEJAK9QbZ65WH9Qe\nv8WvL8pHSAOwCO9whas9pAEW0KYBmEdAM4A2DbCBkAagRUAzhJAGlI+QBkAioJlDmwaUjylPAAQ0\nowhpQPkIaUC9CGiGDW3TfvnZ6QCjATAEbRpQJwJaBWjTgPIR0oC6ENAqkXJt2lefPpHkcQFraNOA\nehDQKkObBpSPkAbYlySgOefOOee2nXM7S66/1Px34fUYh52eQPlo0wDbogc059xZSfLeX5u9PGfH\nOXdL0l7MsdWGkAaUj5AG2JSiQTsv6aD5fk/S9oLb/NR7f6oNcQhnPqSxgxMoD20aYE+KgLYh6c7M\n5eMLbnOsmQL9WaQxVY0pT8AGQhpgR5abBLz3V5r27Lhz7qGGzTm345zbdc7tfv755wlGaBMhDSgf\nbRpgQ4qAdiDpWPP9hqTbs1c24etcc/G2pM35O2gC3Jb3fuvEiRNBB1ubT55/TZ88/1rqYQAYiZAG\nlC1FQLuq+6FrU1K7WWCj+be99t90OP25G3V0kCRCGmAAbRpQrugBzXv/sSQ1U5cH7WVJ15vrr0na\nblq02zPXIzLaNMAGQhpQHue9Tz2GUba2tvzuLiVbaM+8++rgnw11JoHHb2W5hBLI1o23X0o9BKB6\nzrmPvPdb627HOxw6oU0DyseUJ1AOAhp6IaQB5SOkAfkjoKE32jSgfLRpQN4IaBiMkNbf0Zt3F34B\nqRDSgDwR0DAKbVp3q4IYIQ0p0aYB+SGgYRKEtNW6BDBCGlIjpAH5IKBhMrRpi/UJXoQ0pEabBuSB\ngIbJzYa0UMdAK8WQwEVIQw4IaUBaBDQEQZs2DiENOaBNA9IhoCGovQsvpx5CMmNDFiENAOpFQAMy\nRkgDgDoR0AAAADJDQAMyR4sGAPUhoKFYX566l3oIQBXYKADER0BDcHsXXq56s8AUaNGQGjs6gbgI\naIiGkAaUj5AGxEFAQ1S0acPRoiEXtGlAeAQ0JEFIA8pHSAPCIaAhGdq0/mjRkBvaNCAMAhqSsxrS\n9k8fST0EIBpCGjAtAhqyMLRN41AbQD5o04DpENCQFatt2pSY5kTuCGnAeAQ0ZMdSSGOaE7WiTQPG\nIaAhS2wgAGwgpAHDENCQNUIaUD7aNKA/Ahqyt65Ny32jANOcwCFCGtAdAQ3FoE27j40CKBUhDeiG\ngIailLo2jRYNuI8pT2A9AhqKVGJIA/AgQhqwnPPepx7DKFtbW353dzf1MJDQ5uU3JUmP38r/88aU\nU5O0crDkxtsvpR4CEIVz7iPv/da62+X/jgassXfh5SLCGYDlaNOAB9GgwZQSXuSnatFo0GAVbRos\no0FDlXhhB8pXwgctIDQaNJiV84v8FC0aDRpqwIcuWEODhurl/MJOuAK6yfmDFhASDRqqkOuL/NAm\njYCHGuX8oQvoKusGzTl3zjm37ZzbGXI90FeuL+wELaC7XD9oASFED2jOubOS5L2/Nnu56/XAUIQ0\noHychQC1SNGgnZd00Hy/J2m75/XAYDfefinLoNYnpBHoANo02JcioG1IujNz+XjP64HRSg9pAGjT\nYNsjqQcwRLM2bUeSTp48mXg0KFUb0nJ6gW9D2vzmAcIbYli1aSXn38EzL76V5YcuYIwUDdqBpGPN\n9xuSbve8Xt77K977Le/91okTJ4INFHXI8YV9//SRB76A0NbtKJ7yPLIh5PRBC5hCioB2VdJm8/2m\npHYzwMaq64GQcl2bBsTQNXwdvXk366DGlCcsiR7QvPcfS5JzblvSQXtZ0vU11wPBEdKA9XIOaRJt\nGmzgQLXAArzAoxZjwlYJ0+986EJusj5QLQCgfLk3aRIftlAuAhoAwDTWpqFEBDQAwGAltGhAiQho\nwAKsW0ENpgpXhDRgegQ0YAlCGtAdIQ2YFgENWIHjowHdEdKA6RDQgA4IaUD52CyAkhDQgI5o04D1\nSmjRNi+/mXoIwFoENKAnQhpQvs3LbxLUkDUCGjAAbRqwXAktWouQhlwR0IARCGlAeR6/9eBbHyEN\nOSKgASPRpgEPK6lFk5jyRH4IaMBECGnAg0oLaRJtGvJBQAMmRJsGlI82DTkgoAEBENJQgv3TR4I/\nRq4t2vw6tEUIaUiJgAYEQpsGlI82DakQ0IDACGlA+QhpiI2ABkRAm4aa5TrN2RdtGmIioAEREdKA\nfHRZh7YIIQ0xENCAyAhpQPlo0xAaAQ1IgClP5CLGTk7JzjTnPEIaQiGgAQkR0oDy0aYhBAIakBht\nGmADIQ1TIqABmSCkAfEN3SiwDG0aptLpN9M595Rz7hXn3D9pLv/rkIMCakWbhhRYhzY9QhrG6vrR\n4d9Iui7pL5xz/1zSj8MNCQAhDSgfIQ1jPNLxdr+SdMt7/yvn3L+UtBlwTAB0P6SdefGtxCMBMNQz\n774qSfrk+dcSjwSl6dqg7UnaliTv/XVJrwcbEYAH0KYB5WuDGtCV896vvoFz/8B7//tI4+lta2vL\n7+7uph4GEAVtGkKKtUYs1pq3rr48dS/YfT/21BcP/RttWt2ccx9577fW3a5Lg/bvmnVncs79Wfs9\ngPho04DpTb2Tcx3aNHTR5bdyV9Jm06T9StKxwGMCsAI7PYHyPfPuqwQ1rNQloG1KOi7pdefcf5N0\nNuyQAHRBSMPUcpt6rAEhDct02cW5573/K0l/KXEMNCAn7PQEysdOTyyytkHz3v/VzAFq/0zSqcBj\nAtATbRqmQouWDm0aZnVaGem9/x/Nf3/lvf+PQUcEYBBCGlA+1qahxbk4AUPYQIAp0KKlR0hDkoDm\nnDvnnNt2zu0suf5S89+F1wNYjZAGlI+QVrfoAc05d1aSvPfXZi/P2XHO3dLhGQwADEBIwxi0aEBa\nKRq085IOmu+/PYXUnJ9670+1IQ4AAKAmKQLahqQ7M5ePL7jNsWYK9GeRxgQAmFNTixb7bALAOln+\nRnrvrzTt2XHn3EMNm3Nuxzm365zb/fzzzxOMEACA8FiHVq8gAa0JUPNf55qrD3T/dFEbkm4v+Nn2\ntrd1eCaDBzQBbst7v3XixIkQ/wuACaxDw1g1tWi5IqTVqcuZBHrz3l9ZcfVVSe1Z3DcltZsFNrz3\nBzpcl7bbXH+8vR7AMJxtADkh8A3D2QbqE32K03v/sSQ1U5cH7WVJ15vrr0nablq02zPXAxiBNg1D\nEaryQZtWD+e9Tz2GUba2tvzu7u76GwL4Fm0ahjh68+7o+8g57H156t7k9/nYU19Mfp/PnbwpSXrj\nzNXJ7xvhOec+8t5vrbtdlpsEAIRFm4YhxoarnMNZiV65cT71EBAQAQ2oFCENQxCy8vLKjfMENaMI\naEDFOHcnhhgS0gh202inN+cR0uwhoAEgpKE3Ald+aNNsIaABkESbhv66hjTCXFyENBsIaAAeQEhD\nH/unjywNYKuuQ1i0aeUjoAF4CG0a+mrD2OwX0iOklYuABmApQhpQPkJamQhoAFaiTQPKx5RneQho\nADohpAHlI6SVg4AGoDPaNKB8tGllIKAB6I2QBsS37CC1QxHS8kZAAzAIbRpQPtq0fBHQAIxCSAPK\nR0jLDwENwGi0aUD5aNPyQkADMBlCGlC+G5/9o9RDgAhoACZGSAPK9cLRDyUdhjSCWloENACTY8oT\nsIGQlg4BDUAwhDSgHG17No82LQ0CGoCgaNMAGwhpcRHQAERBSAPKR5sWDwENQDS0aYANhLTwCGgA\noiOkIUeP3+ItsQ9CWlj8NgJIgjYNNfjq0ydSDyEopjzDIaABSIqQBpSPkDY9AhqA5GjTgLSWHWKj\nD9q0aRHQAGSDkAYs98vPTqceQieEtGkQ0AAAqNw7+89Oen+0aeMR0AAAQBCEtOEIaAAAIBhC2jAE\nNAAAgMwQ0ABkhR2dAEBAA5ApQhqAmhUf0P77J79JPQQAgdCmATawDq2/4gOaJP35Dy6kHgKAgAhp\nKNmUp3sKeSy0qQ+1MY+Q1k+ygOacu7TiunPOuW3n3E7X+/vzH1wgqAGG0aYB5eP4aN0lCWhN8Dq3\n5LqzkuS9vzZ7uStCGmAbIQ0oHyFtvSQBzXt/RdLekqvPSzpovt+TtN33/mnTANto04AwQk9zzqJN\nWy3HNWgbku7MXD4+9I4IaYBthDTUqJRzcnZFSFssx4A2Kdo0wDbaNKB8tGkPCxLQnHM7C74Wrjlb\n4EDSseb7DUm3l9z/rnNu9+69rzvdKSENsI2QhpxNuZMztJjTnPMIafc9EuJOmzVmvTjnNrz3B5Ku\nStpq/nlT0rUl939Fkp448j3f9THakPaL317uOzwABWhD2pkX30o8EgBDtSHtzMm6j3OaahfnOUlb\nc4fRuC5J3vuPm9tsSzpoL0/l66ef1I9+svQIHwAMoE2DdaHXoaVs0Vq1t2nO+84FVJaeOPI9/+yf\nnu/1M18//eS333/w3sWphwQgEzRp6OvLU/eC3fdjT30x6f09d/LmpPc374WjHwa9/z4stWnOuY+8\n91vrbmd+k8Aij/76d99+T5sG2MUGAuSkpHVoUh4tWqvGNq3KgCY9HNIIaoBdhDR09fitct4WYxxu\nI7eQ9sqNfjNmJSvnNzGA2ZAm0aYBltGmAcPkEtLacdQS0qoOaIsQ0gDbCGmwxNpBa7t65cZ580Gt\nyk0C82Y3DcxiAwFgG5sIsEjIjQJSeZsFWik3Daxq8d44czXiSMZjk0AP81OdLdo0wDbaNKC7XKY6\n51lt0mjQZixr0iTaNMA62jS0SmvQpHgtmhS/SesTDEto02jQBljWpEm0aYB1tGmIpbTDbcyL2aT1\nfSxLbRoN2gKrmjSJNg2wjjatbqEbNKn8Fq0Vsk0bGwRzbdNo0AKiTQNso02rW0nHQpuVYkdnruvS\npPLbNBq0Jda1aC3aNMA22rQ60aL1N2WbNnXwy6lNo0EbadV6tFm0aYBttGkoScrjor2z/+zoYDXF\nfSxSYptGg7ZG1yZNok0DrKNNq0eMBk2y16LN69KqxZ4mTd2mdW3QCGgdENIAtAhpdSg5oEl5hbQc\npQxpTHEmwonXAduY8qxDrI0CoQ65UespoLp65cZ5PfPuq6mHsRIBrYOu69FmEdIAuzjxOmBDziGN\ngNbRkJAGwDZCGqZAi5bWM+++mmVQI6D1QEgDMI82DTkjpHWXW0gjoAHABAhp9sQ8YG3I0z8R0rrL\nqU0joPVEiwZgGdo05IqQ1k8OIY2ANkCfg9iyWQCoDyENQ4Q+iTohrSwEtAgIaQBQptjn5SSkoUVA\nG6jvVCdtGgAA6IqANgLHRwMATI0WDRKneppEn1NBzeK0UIBt604Ltep0QrGn1rBarFM/zQp1GqhZ\nNZ8SqmtQ/eT51yZ9XE71VADaNMC2+c0CX56698DXKvO3TREQkFboJk2iTesi1Y5OGrSJDG3RWrRp\ngF1nXnxr8oBFwxaf1RZNqrNJGxJOp2jTujZoBLQJEdIArLJ5+c0g90tYi4eQZkuKkMYUZ4F+9JNL\na9esACjX3oWXg9wv06C2xZjqlJju7CrW2QYIaBOa6iwDhDTArr0LLwcLapIIaoGlaisJafkJHdII\naJk68+JbBDXAsJAhTaJVs4iQlp+QbRpr0AIYuxZt//SRBy5z2hjAtlBr0+axVm06KYNvrDVpUh3r\n0qYKpF3XprEGLaGxU51Hb9594DJNGmBb6DatRaM2nZRhN1aTJtGm9TF1k0aDFsjULVqLNg2wLVab\nJtGoTYEmrXxTh9CvPn1i5YcuGrTEpm7RWrRpgG2x2jSJRm0KNTVpVtu0qcPnY099MckHLRq0gMa2\naNLyJk2iTQOsi9mmtWjV+ksdcmM2aZLNNi1E+GwD9PyHruwbNOfc0vMctdc553bijWh6Ux12Yxna\nNMC2mG1aK3XYKFHqUBuzSZNsrksLETrb4Dz0g1aSBq0JXhe996eWXL8v6Y6kf+u9v7bqvnJu0KTw\nLVqLNg2wjTYtf6nDLU3aOCFbtNbehZfzbtC891ck7a24yU+996fWhbMShG7RWrRpgG2p2rTUoaMk\nj9/6TjVr0iR769JCtmitPh+0cv14dMw5t+2c+1nqgeRg2YaBeRzcFrAt9FkIliGo9VNTSJNsTXnG\nCGldZRnQvPdXmvbsuHNuO/V4xpqiResa0iTaNMC6FCFNIqj1kTqksS5tuFxCWpDfIOfczoKvcz1+\ntr3tbUmbS26z65zbvXvv6ymHHkysqc4WbRpgW6qQJqVfa1WK1Gv4CGnD5RDSkh1mwzn3vvf+xzOX\nN7z3B01jttt8f0nSVe/9x8vuJ/dNArNibRiYxwYCwLYUGwhaqUNICVIH2tibByQ7GwhChM6//lf/\nPt9NAk1DtjV3GI3rktRMbW43t7m9KpyVJnaL1qJNA2yjTctbbZsHJDttWsqgyYFqI0vVorVo0wDb\naNPyVsupoVo0aQ/LukHDOH02DMyjTQNso03LG01amZ47eTN62CSgRZZqmnMeIQ2wK9XhOCR2enaR\ncsqTHZ7jxAxpBLRCjWnRpMMX0ZRTIQDCo03LW01tGiGtPwJaAjm0aO0Lw+blNwlqgGGp2zSslrpN\nwzAxpjwJaAUb26LNIqQBthHS8lZDSLPUorVChjQCWsXmXxAIaYBtKU8VhfVStWmEtHFCtWkEtESm\nmuYc26KxLR6oDyEtb9ZDmlVThzTenQGgQinaNEJadynatFghzWKL1pqyTSOgGUCLBmAoQlrerIY0\n66YIabwzJzTlbs6pQhrr0ID6ENLyFrtNixHSLLdorbEhjYCGh3DoDQChEdL6sxbSajBmyrP4gHbv\n0T9OPYRsTD3VSUgDEBIhrT+WpJRpSEgz8UxPcQLyVHI4aO0qtGkAQiKk9RdryjN0i1bDNOesvm2a\niYAmlR1aTyUtAAAK/klEQVTSphRqwwAhDbAt5WmhMIyFkIblzAQ06TCkEdTGI6QBdeJsA+UhpNll\nKqC1ag9pU54Cah5TnoBtKc/diWFSns8T4Zh9RksKaTmuQ1v3x05IA2zj0BvlCRnSQrVota1D68Ns\nQJPKCmklok0DbKNNK0+JIQ2LmQ5oUr0hbYppzq5/6IQ0wLZYIY0WbRpMd9pQxbNYa0ibQp+QRlAD\n7KJNK0uokEaLFs8jqQcQy9dPP5nlWq+Qjt68q/3TR6I+5ublN3kRBwzbu/By0A9jX566RwM0kcdv\nfYdWsoMXjn648vp39p+NNJIHVRPQpPtNWm1Bbay+f+TtizdBDbCp/dumNc9fiJD21adP6LGnvpj0\nPmNbF8pW3TZWYKvyY0pNU54hD7mxDi/egG2hPoTR+kyLRvLQC0c//PYrh/tZp9pnraaQNoWhf+Cs\nTQNsoykvw9QhrbS1aKHCVMigVm1Ak+oJaVO1aGP+wAlpgF1sIChDrU1a6KarfYypH6fOZ2sGp4eK\nhzYNsG3KkMY0J6YQI5zNP95Uj1l9QGulDGklBcQpPoER0gC7aNPyNmWLlvs0Z+xwNvVjE9BmWG7T\nUm4WWIQ2DbCNkJavGqY6U4azqcZg/1kaIGZIKzEQTvnHTUgD7KJNQwo5hLPWmLEQ0Jaw2Kbl1qK1\nCGmAbUNDGuvQwrHaouUUzlpDx1TVgWqHCHlw25ID4NQHP+TgtoBtNRzctu+H4NhnegnBwkFrY3jh\n6Ie9D3BrM0IHMHWjVnI4a4X4BPbMu69Ofp8A8mHtQ9jRm3e//Yr5s1PIoUV77uTNye4rx/ZsVt/x\npX92CjNFUEsZznKd5pz1zLuvEtQAw6ysTZvy9TRVWMshpNWkT0gr/pn55rsuyeO2Qa1v2LLQnM0K\neXRqQhpgW6khLXSQStmqlSr39mwIE2vQ9k8fSfrLPB+62vVq1sJYCm1I++T51xKPBEAI69amfXnq\nXlYtT8z3mqM375pYp4Zh8vmtNyT3HaBTv8DEOMcbbRpgWwltWooiIEabNvY1PPUBay22Z5KhgMan\nDFuWhTSCGmBXzmvTUk85pn58xJckoDnndpqvS0uuP+ec23bO7fS5X0JaOiGmIJZ9KiOkAbblFtJy\nCUchx5HTNHIfVtszKUFAc85tS7rmvb8iabO5PHv9WUny3l+bvYxp5fKCMxRtGmBbLiEtt9fK3MaD\ncFJE5k1JbSjbay7POi/pYOb6bfVAi5ZOzBatRUgD7Mp5yjMlQlodogc07/2Vpj2TpLOSdudusiHp\nzszl4/P30UyP7jrndv/u6/8TaKTIRZeQRlAD7Lrx9ktJHjfnIBRibKVNc1qe3pQSbhJopi4/9t5/\n3Pdnm5C35b3feuTRv//Q9bRo6aT8AyekAXbdePulqEEt53DWKmGMGC7Iu+nMJoDZr3NzN9v23l9c\n8OMHko41329Iuj1kDIS09Ur64+66jZs2DbAtVZuWq5Jex9FPkAPVzkxhLuSc2/Hev958v+29v+ac\n2/DeH0i6KmmruemmpGshxohwpj6ReouT8gKQ7oe0My++FeT+CT1xjDkPp/XpTSndLs5Lzrlbzrn9\nmauuS1I75dnc7mDIFGiLFg0A7KJNO0SgtCn6qZ6aw2ccXfDvP5z5fmUDh/zRogGIYeo2rdSwM9Vp\noUK9dqO/srZsDECLtlqJL0apTysCID+0abDGfEBDOuzoBBDT2J2eJX5gnVX6+PGgKgIaLZo9XVo0\ndnQCgD01bBCQKgloSCdki9bn0BsAAJSkmoBGi1Y32jQAq1iZHkz1/8Ha4OlVE9CQTg4tWouQBgCr\nlXbKJ6uqehZo0WwaEtIIaoBdte/otNIG1q6qgIbFYvwx5/iJjJAG2NVnRyeBBjnK710TGGDo+gfa\nNMC22ts0lKu6gMY0ZzqhW7Qxi1QJaYBdY4+PBqRQXUADliGkAbYtCmlWpzet/n/VhICGqHJu0SSm\nPAHraNNQiioDGtOctk1xPB5CGmDbjbdfomVC1qoMaEgrxo7OqUIaQQ2w64P3LuqD9y6mHgawEAEN\nWIOQBthmNaTREJaNgIYkSmnRWrRpgG20acgNAQ2m1+RNfX44QhpgGyENuag2oFkOJaXI8ewCXdCm\nAbbRpiEHZb5DAj1M3aK1CGmAbYQ0pERAQ1Kltmgt2jTANto0pFL2uyPQUagWrUVIA2yrLaSV/uHZ\nAp6BytW0Fi9GSCOoAXbVFtJy9c7+s6mHEAUBDcnF/KQWOqRJtGmAZUx5IhYCGhAAbRpgGyENoVUb\n0DjCcl7Tm9ZatBYhDbCLNg0hVRvQgFgIaYBtOYc0q2VEDevQCGiVyqk9a1lt0SSmPAHraNP6++Vn\np1MPIWsENFQrdkiTaNMA6whp8Vhv0QhoyEoNx96hTQNso03DFOy/G+IhOU5vppKiRWsR0gDbCGnh\nldai9RlvlQHN6qLJLghnD0sd0ghqgF20aeGVFtK6qjKgIW81THPOI6QBthHSFqtpo0DfIFnfO2HF\naM+WS9miAagDbVo4Fls0AhqylKJFI6QBiCFmSKvpg3nOIW3I2KoLaLWuP6vpjxQAcmexScvhQ27O\nIa2v6gIasErqDQMA6sGU56Gp16HlFtKGjidJQHPO7TRfl5Zcf6m9XdyR2VRqe5ZqswAhDQDKlltI\nGyL6O6BzblvSNe/9FUmbzeV5O865W5L24o4OSI+QBqA2IXZzvrP/bPKgNubxU1QUm5LaULbXXJ73\nU+/9Ke/9tXjDsqnU9qxVY4smcXw0AJhKqpA29nGd936ioQx4cOfel3TRe//x3L/v6DC8nfXev77g\n53YktdOfpyXdDD3WOf9Q0t9EfkyExXNqE8+rPTynNtX0vP5j7/2JdTdKFtCcc2clnffeL10h2axF\nez+3Js05t+u930o9DkyH59Qmnld7eE5t4nl92CMh7nTJ4v473vufz1zeXhTOmp9tb3tbi6dAAQAA\nzAoS0JoNAEs553baqUvn3Lb3/ppzbsN7f6DDqc3d5qbHJWXVngEAAISWahfnJefcLefc/sxV1yWp\nmc7cds6dk3R7fn1aJlYGUBSJ59Qmnld7eE5t4nmdk3STAAAAfTjnLq1auwxYwZkEgAWWHUQZZXDO\nnXPObXOwa1ua5/Nc6nFgOusOXF8zAtpI/HLZw5tA2Zod4u1yiW8vo3zN+mYOYG5ExwPXV4uANgK/\nXDbxJlC885IOmu/3dP/A2ADy0uXA9dUKsouzIpvNV/uGzi8XkN6GpDszl4+nGgiA5eaO+HBW0tVU\nY8kRAW0EfrkAABinWYbwcaZHbUiGgDYBfrnK0vFAyijXgaRjzfcbOjzgNYB8LTxwfe0IaGuMOSsC\n8rTuQMoo3lVJ7SljNsXBrs1ojo+51RzsnL9jAxYduD71mHLBcdBGmn2h4JfLhuZN4C8lXeRNoEzN\nB6s9SZs8h0Cemo11/0WHa0aPSfoL3kPvI6CNwC8XAAAIgYAGAACQGY6DBgAAkBkCGgAAQGYIaAAA\nAJkhoAEAAGSGgAYAAJAZDlQLAPr2+HebOjx+2r+Q9B+89werfwoAwqBBA1A959xmc3aQNpBdbcOZ\nc26zOeYhAERDQANQPe/9XvPtDyVdmzuv7jkdntMTAKIhoAGonnPubPPtpvf+YOayxLk8ASTAGjQA\nkLadc5uS3m+mM++kHhCAunGqJwBYwTn3M0nHxaYBABER0AAAADLDGjQAAIDMENAAAAAyQ0ADAADI\nDAENAAAgMwQ0AACAzBDQAAAAMkNAAwAAyAwBDQAAIDP/H4b/oaBVbFBDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9ef7b8a250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Objective & constraint\n",
    "def townsend(X):\n",
    "    return -(np.cos((X[:,0]-0.1)*X[:,1])**2 + X[:,0] * np.sin(3*X[:,0]+X[:,1]))[:,None]\n",
    "\n",
    "def constraint(X):\n",
    "    return -(-np.cos(1.5*X[:,0]+np.pi)*np.cos(1.5*X[:,1])+np.sin(1.5*X[:,0]+np.pi)*np.sin(1.5*X[:,1]))[:,None]\n",
    "\n",
    "# Setup input domain\n",
    "domain = OptDemo.domain.ContinuousParameter('x1', -2.25, 2.5) + \\\n",
    "        OptDemo.domain.ContinuousParameter('x2', -2.5, 1.75)\n",
    "\n",
    "# Plot\n",
    "def plotfx():  \n",
    "    X = OptDemo.design.FactorialDesign(101, domain).generate()\n",
    "    Zo = townsend(X)\n",
    "    Zc = constraint(X)\n",
    "    mask = Zc>=0\n",
    "    Zc[mask] = np.nan\n",
    "    Zc[np.logical_not(mask)] = 1\n",
    "    Z = Zo * Zc\n",
    "    shape = (101, 101)\n",
    "\n",
    "    f, axes = plt.subplots(1, 1, figsize=(10, 7))\n",
    "    axes.contourf(X[:,0].reshape(shape), X[:,1].reshape(shape), Z.reshape(shape))\n",
    "    axes.set_xlabel('$x_1$')\n",
    "    axes.set_ylabel('$x_2$')\n",
    "    return axes\n",
    "\n",
    "plotfx()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We proceed by assigning the objective and constraint function a GP prior. Both functions are evaluated on a space-filling set of points (here, a grid). Two GPR models are created, the EI and PoF acquisition functions are created and multiplied into a joint criterion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial evaluations\n",
    "design = OptDemo.design.FactorialDesign(4, domain)\n",
    "X = design.generate()\n",
    "Yo = townsend(X)\n",
    "Yc = constraint(X)\n",
    "\n",
    "# Models\n",
    "objective_model = GPflow.gpr.GPR(X, Yo, GPflow.kernels.Matern52(2, ARD=True))\n",
    "constraint_model = GPflow.gpr.GPR(X, Yc, GPflow.kernels.Matern52(2, ARD=True))\n",
    "\n",
    "# Setup\n",
    "ei = OptDemo.acquisition.ExpectedImprovement(objective_model)\n",
    "pof = OptDemo.acquisition.ProbabilityOfFeasibility(constraint_model)\n",
    "joint = ei * pof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now inspect our belief about the optimization problem by plotting the models, and the EI, PoF and joint mappings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAALICAYAAABW0sdqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3U+MHPed9/fPd0Txn0lpSMoPFZHWmkMgkgwdssNxDnuU\nR5vAx93R8vQcNcotAbIQ4UPOBnUNEEDj4wIBGHP38CAwniyHOuY5PBS9QAxLwoIjryPpkWKRHJs0\naVPUfHOYKrKmprq7qrqqfr+qer8Agt1d/37dU/3rT3/7V1Xm7gIAAABQzkLoBgAAAAB9QoAGAAAA\nKiBAAwAAABUQoAEAAIAKCNAAAABABQRoRMPM1szsw+Tfz81sMTNt0cxuTVjurpktzbHdJTO7W3f5\nNpjZrezzz02b+FoAAID2EaARBTNbk/SOpB+5+wVJ70v6cFKIzDnn7lsVt7ee3k6WPVdleQAYCjNb\nTgoXt5J/7wZsCwUN9AIBGrG4LOktd9+WJHff1G6I/smsBdNlKnqngXUAQK+Z2bKkn2u3/z3v7ucl\nbTa07vXZc+1VtqBRZ91AkwjQCC7pwLcLQuympNXM/ZNm9n5SDXg/s/yH6RAOM3s3mX4trRiY2Xqm\nsnI5WTatuFzOVhCSoSOrmXXfmrTeguexmKzz/cy23k0eu5aZbz0zVCW7rcvJcj+XdDLz+MxtA0BN\nP5P0TvZXPHe/2dC635k9y34lCxq11g00hQCNGCxJulPw+FYyLetSUiFZyYZP6UkQ/2Ey/R1Jl5PH\n3slUVt5393ckbbn7BXe/lFv/FUlvZdZ3s2i9U57LsqTLybzvaveLwYV0fcm63kq2fSFp42LyXFaT\ndr417TlN2TYAVLWc/OK3T9GX/aJCQfL4UvIl/8Pk/6JCRXr758kyP8/Mv5hZ/60Z29qz7oJ2z1XM\nSKZR0MBUB0I3AFBxUFbyWHZs851MZeJ9SW9q70+NFyUtZTvI5LEr6Z1ZY6Xd/WqmQ06XLVrvJNuZ\nbWxl2pc+xx8mbU9dkfR3ks7nHk+/UFTZNgCUlvxyV9gnZr/sJ/c/NLMfJZPTae8k45UvafcL/vvu\nfjWzjtXM8ovJclfSwkVaLLDdY2B+kqwnb9+2kttP1j1BdjnXbiHlQhJ8l5N59j0/d9/OFjSSaXcz\nr8kP3f188tpdntBmjAABGsG5+00zO2lmS7mAe1HTx+LdLnjsp7kO/LKkquOb06rzqrtfMrMf5tc7\nRb6Snr+fr1icKtG+/HOi6gFgbu6+ZZPPYHRRxV/2/w/tLRTcSfqka5J+nvSX708oVmy7+3vpnSSo\nvqP9xZL8Mnu2VXKIR91ixoZ2izMUNDAVQzgQi7e12/mmP+OtSlqT9NPMPCczt9/S/nB9RZlxcUkI\nLnpMkraTn/mKwmi6zJ3c/fw66vh5uq5k22vJ87imp0NHFvW0It/ktgEg72Z++EKi6Mt+at+Qu2QY\nyAXtFjY+nLCtJ8slfdkl7fb9b09pX9HwvjLqFDPK+Km7v5n+q9c0DAEBGlFIKqw/lXQ9Gf92SdKF\nfKUhGdP2oaSb+QNdkvvpmLoPJa0kj6Xj4G7p6UGJm9rt5PeNn0vasq6kAlG03jme56aka0lbrmv3\n58it5PGt5PEnVfMmtw0ABd7Wbh/5pBKdBOpJX/YLpb8gJhXmG8n6phUqVrTbj2+rXr82bd1lTHt+\nFDQwk7l76DYAc0mC5VtVzwUNAHgSBH+mp1XZ9939Pds9H/Q72v1C/9PkGJFFSR9mxgff0m7l+e+S\neRclbSZjjy/raTC9lFtuUbtFjG0lwzfc/a3s+idtKxmn/GTdyYHh2eczbbn3JV1Lnsu+55dZx/va\nLbhsJs/tXLL8up6G6Pe1O6TlybYwHgRo9F62cwzdFgAAMHwM4UBvJacZuqvdCgThGQAAdIIKNAAA\nAFABFWgAAACgAgI0AAAAUEHvL6Ry4uSCv3R2vqdx+9tjDbWmmj88Ohxku3nf/vmZ0E1o1TOHvg3d\nhEY8d/BPnW8zln20rL7uy4/+38++dvfvhm5HHc8e/o4f+s7JfY9/+2zx/DsHW25QC+zZndBNQIQO\nP/vNxGnPHajWX//h8eS+9k/fTHgzBeTfxFd/XXhU/PgzuT/Twp/3v5/v/fGLyn1w7wP0S2cP6H//\nP0/XXv4f7v5Vg62p5p9/+0qwbUvS/U+fD7r9rh079/vQTWjEX7/8Savrz+6X/67VLTWn7/vyp//j\n3/9b6DbUdeg7J/X6f/8/7Xns/pnJH64PzvYvjB546UHoJiBCr7341cRpb7zwcaV1ffD1qxOnffRl\n/YzTlsdfHA3dhH2Oflbc7xz7fH+fc/zTh3vub/6n/6VyHxzfV4gOEZ7HZSjP+Z9/+8qTf22st2+G\n8ncdiqGFZ2CSacF2WiCuMm+M4TlWk/qXoj7p3rkjunfuyFzb630Fui7C8zjd//T5wVSipf37UpXq\ndOj9cF5j3o9j8u2z00NzivCMsUmDcVE1ukzAJjxX9+DsTmEl+v6ZhcJK9JMQ/Z+qb2uUAZrwPG7p\nazCkIJ0KvX91pc/78fFb4/vhj/CMofroy9NTh3JI1arR2fWinmkhWioe0lHHqAJ0yOAcWp8DR1uG\nVo0egz7ux2MMzCmCM8YgDbuzgnTV9aG+SSFaKveLWRnBevbkOvZTpyXXnG9EDOE5VHWwj6GjK7w2\n/dGnv9XxWwtP/sWqzT74wdkdwjNG56MvT88VfuddHnu13Q8FqUAnnfKapEsTZlk3szVJ78y7rRiC\ns0R4jhmV6Pj1YT+OOSznNdkH7xyk0gxk5UPwpMo0Ybkbaf80qSJdV5AA7e4bZvbWlFnedver82wj\nluAsEZ77YMjjovss9n24T6E5q4s+GMAugnIc8l/05w3UsY6BPmlmq5KW3f29sgvFFJpThOd+oRod\nj1j34b6G5opq9cEA0Bfz/nIWZYB29w1JMrM3zWzV3Tez05OfH9cl6fh/dTTK4CyFCc+xho4+IUSH\nFeM+PJLQ/ESVPviZEycCtBAAworuU8HM0rF3knRb0lJ+HnffcPcVd185euJQtw0sifDcb7yWYcT2\nusd+IGAbqvbBzxz7TrcNBIAIRFOBNrNFd9+WtCXpRvLwKUmbk5eKE+F5GBgX3Z3Y9t+xhWZpWH0w\nALQtyKdEUt1YyZ0i6bokJT8Vribz3Hb3myHaWBfheXh4fdsV0+s7lorzkPtgAOhCqLNwXJV0NffY\nhdz03iE8DxfjotsRy/47htCcNdQ+uC2PvziqAy89CN0MABGJZghHn3GmjXFgSEdzYtl3xxacAQDN\nIEDPiarz+FCNnk8M+y/BGQAwDwL0HAjP40U1uroY9l2CMwCgCQTomgjPkKhGlxHDfktwBgA0iQBd\nEeOdkUeIniyG/ZbwDABoGgG6AsIzJmFIx14x7LMEZwBAWwjQJRCcUdbYg3QM+yzBGW3gVHYAsgjQ\nU4QKzlIcQaSOKuHl3vmdFlsS1tiGdcSyv8Yenk988ih0EwAADSBAFwgZnKV4wsgs84aVouWHFKrH\nUI2OZV+NNTgTmIeFKjSAFAE6I3RwluIJJJO0HVSy6x9KmB5ikI5pP40tPBOah+3xF0cliSANjBwB\nWgTnWUIFlKGF6SEE6Zj205iCM6F5fNIgPQkBGxi2UQfoGIKzFFcoyYopoKRtGVKQlvoRpmPcP2PY\nNwnNmGZWwMYuvmj004GXHox+Hw8WoM3ssrtfmjBtTdK2pCV332hyu7GE5lRs4SSGYDINVenuxLZv\nSuH3zyGF5lB9MJCVDWGEafRJkABtZuuS1iTt67zNbFmS3H3TzNbNbNndb86zvdhCsxRfOAkdTOo4\nfmthECFaiqMqHds+mRdyHx1ScJa674OBMhhf3i9jr0IHCdDuvmFmb02YfFHSteT2lqRVSZU67xgD\nc1ZMQaWPwTlrSEM7Uvn9o61AHdN+OA3BuXlN9sH27E4vA8+YP/hjx9lO+mPMITrGMdCLku5k7p+a\nNvMfHh2OPjCnYgssfQ/PWUMM0qlZ+82kgB3b/lZHqH10qMG5pEp9cF9lA9pYA0DMqEb3x1jfSzEG\n6EGKKcwMKTjnDTlITxLTvtWkEPvpyIPzaI25ihY7qtH9Mqa/VYwBelvSyeT2oqTb+RmS8XvrkvTs\nd+MODzGFmyEH57whjY8eG6rOwVXqgw+8EE8fNw9CdLwI0YhRNInKzBaTm1ckLSW3lyRt5ud19w13\nX3H3lQPPx9nh3f/0ecJzYMdvLYzyefdZqKoz4bl+H/zMc9/pqokAEI0g6SI5RdJKUsVIXZek9Ghv\nM1uVtN23o79jDM5jD5Fjf/59EGI/HXNwHnIfjGHi1wHEJtRZOK5Kupp77ELmdi/POxpbcMZTYxwb\n3ReMde7eUPtgAOhKjGOgeyem4CwRnqdhbHQ8CM4AgL4iQNcUW2iWCM5lUY0Oj/AMAOgzAnRFMQZn\nifBcB9Xo7nGGDQDAEBCgSyI4DxMhujtUnfc68qvPQzcBOZzKLm6czg4xIUBPEWtoThGem8GQjnZR\ndd5FYAbmQ3hGTAjQObGHZong3Baq0c0bc3gmMAPAcBGg1Y/QLBGcu0CIbkbIfTV0eCY4A8DwjTJA\n9yUwZxGeu8OQjvoIzgCAMRhFgO5jYE4RnMOhGl3NWMMzwRkAxmdwAbrPYTmL4BwHqtGzhd5XQ4Vn\ngjMAjFfvA/S3f35mMKFZCh9GUIxq9H4x7KshwjPBGQDQ+wA9JDEEEkxGNXpXLPsp4RkAEAoBOgKx\nBBKUM9ZqdEz7KeEZABASATqgmAIJqhlTNTq2/bTr8ExwBgDkBflkNLM1M1s1s/UJ0y8n/xdO77Pj\ntxae/EP/DflvGeNzIzw3Y8x9MAA0ofNPRzNbliR338zez1k3s1uStrpsW5tiDCNozlD+vjF/wSM8\nN6PpPvjws9/otRe/2vcPAIYsxBCOi5KuJbe3JK1Kupmb5213v9ppq1oQYwhBu/o6tCP2fZXw3KhO\n+uBsiP7oy9PzrAolHP1s+nv4wdl+9UlA7EIE6EVJdzL3TxXMc9LMViUtu/t7+YnJz4rrkvTMiROt\nNHIesYcRtC/2IN2nfZTw3LhG++Ajp4/N3GAapgnSzZoVmifNS5gG5hflQYTuviFJZvamma2mPzXm\npm9I0qGXv+cBmrhHn8IIupXdN0KH6T7up4TnMKr0wS+//py/8cLHT6Z98PWrE9f72otfEaLnVCU0\nz1oHQRrY7Zf+tcZyIQL0tqSTye1FSbezE5PKxp3k58Pbkpa6bV45fQwjCCu/z7QdqPu+jxKeW9Nq\nH5yG6UlBmhBdXxPhuWh9BGmMSVPHaIQI0FckrSS3lySlB7Isuvu2dsfk3Uimn0qnh9b3MIL4TNun\nyoTrIe+ToS7PPRKd9MHTgjQhurpp4fnY57P7i/tnJi9/9LMFQjQGr+mDmzsP0O5+08xWkvF12+6e\nHrxyXdIFd99MTrEkSbcz0zs15HCC+I15/+MiKe1qug9+fuFP+vGxXz+5/4v7P9gz/Y0XPiZEz2lS\neC4TnPPzTgrShOj2FAU39v3utHVWoCBjoNPxdbnHLmRud3YGjjEHFSA2hOdutNkH//jYr0uHaMzW\nRHjOLzetGo1mTQpvHFjbf1EeRNgGgjIQN8LzcBSFaLTv+KcP9z1279yRfY9NCtFUoZtVpvLJLzHt\navOc9IML0ATl7p345JHuvnIwdDPQY4x5Bsorqj4Xhef08aIQDWA+vQ/Qz/yZ0Ny1orCTfYwwjT6g\n+oy+un9mYV+IvnfuyMQQDYzVR1+eHtYYaPRT2SphOh9BGmVQfR4ehm+EkVaaZ1WdJ42BZvhGs8qE\nN4Zv9BelW8x04pNHtUIOwQizhNpHqD63pyg8cwBhfZNC7bQDAeuEZ7RjUkD+6MvThOeeowKNiQjA\naBP71zD8fufw1IrzpPBMeCjvwdmdwrNxpGG47Bk5poVnqs/tYV8Pq61hHARoFCLcoE3sX+NA5bk5\nk0K0tD8YVz1VHeEZQ5d+iWkySBOgsU9T4YYx0CgSOjwzfKNds0Iz1bj6poXorLLhmeCMscn3P/ME\nagI09ggdbgD0yx8eHy5daSY8zy8NvWWCdJn1AGM2T59EgMYTTYZnqs8owhe0cSI4Ny8fgGcFagIz\n0CwCNCQxbAPtIzyPC6G5WwRkoFsEaBCe0bpYwjPjn5v3p2+eJSwDGJ0gAdrM1iRtS1py942q09Gc\nJoINwRnTxBKe8RR9MADMp/MzqpvZsiS5+2b2ftnpaM48webuKwef/APQH/TBADC/EJckuqjdyoYk\nbUlarTgdDagSnrNhmdCMKqg+R4k+GADmFGIIx6KkO5n7pypORwPuvnJwT7ghFKNphOdo0QcDwJx6\neRChma1LWpekZ4+fqLRsDB/qsYTVWNoBdOXh62c4kLAB2T74wAvPB24NAHQvRIDelnQyub0o6XbF\n6UoOatmQpOPPn/UYQnEVIdpLWEaX+vaeHJlG++DD5894O80EgHiFGAN9RdJScntJUnqgyuK06QD6\ngfAcPfpgAJhT5wHa3W9KkpmtStpO70u6PmM6AGBO9MEAML8gY6CLzivq7hemTUd9DN9AV/pQfWYc\nNH0w+unxF0d14KUHoZsBSAozhAMdIjyjK30IzwAANIEADWB0Hr5+JnQTAAA9RoAeMKrP6Eofq8+E\naABAXQRoAKNFiAYA1EGAHiiqz+hKH6vPWYRoAEBVBOgBIjwD1RCiAQBVEKABQIRooA8ef3E0dBMA\nSSUDtJmdM7O/N7PvJ/f/ps1GAeiHvg/fyHv4+pkogzR9MPAUIRoxKFuB/lvtXqXqLTP7byS92V6T\nMA+GbwDzizBE0wd3hAt19MPjL44SpBFU2SsR/lLSLXf/pZn9SNJSi20CgODSEB3JVQvpgzt04KUH\nhLOeaOvv1LcvUk2/Dn17/iGUDdBbklYl/ZO7XzezFpuEuqg+A83LVqMDhmn64I7VCRCE7uEY+9+y\n6PkTqveaGaDN7Dl3/1TSp+lj7n691VYBiN7Qxj+XkR/a0UWgpg/uDwJGc8YeYGOU/k3Yz3eVqUD/\nxMyuuPu/mNlfSnJ3/5e2G4ZqqD4D3SsaK91CqKYPxuikIY0gHZ/HXxwlRKtcgL4hacnMtpLxd2/M\nu1EzW5O0LWnJ3TcKpl9290tmtl40HQBiVesAxOmZO+o+2L9ZmCvk8EEM9A8hulyAXtJuR/uemZ2T\ndE3SB3U3aGbLkuTum2a2bmbL7n4zN9t60sG/U3c7Y0L1GRi0QffBXVYYx/6B30cc0BmvsYfoMgF6\ny93/UdLPpEbOP3pRux8A0tMDY/Kd99vufrXMyhYeflP6J9MIT00F9NbdVw6Ochx0AFH3wX3CGE4A\nTZl5Hmh3/8fMyfv/UtL5Obe5KOlO5v6pgnlOmtmqmb1btIKkanLDzG482nlYesNHfvX5nn9DQPUZ\nGLbY++Bv7/9xzuZ0j4omgHmVOo2du/8m+f+X2j0faavSMXdm9qaZrbr7ZsH0DUl6/uBpr7udbIim\nOg0gVjH3wYde/l7tPhgA+qrseaArMbP1gofvJD8Jbks6mTy2KOl2wbLpvLfV0QUD0jDdpyBN9Rmh\nMYwjTn3sg7s29vGbAObTSoCecdT2FUkrye0lSZuSZGaL7r6t3TF5N5Lpp9LpXTnyq897FaIBIK/P\nfTAA9MHMMdBNS4/2NrNVSduZo7+vJ9M3Ja0mR4DfLjg6vHVDGiMNtI1fQvqlD30wAMSulQr0LEXV\nEXe/kLkdxdHfMVejCS2ICUM5+qUvfTAAxKrzCnTfUIkGyuFLHQBgLAjQJcQWogkqiBX7JgBgDAjQ\nJcUWooFY3X3lIEEaADBoBOieIZigLwjSAIChCnIQYV/FfFDhUNw7v1N5meO3+B4Ys2yI5kBDAMAQ\nEKDRuTohuc76CNbxIUwDAIaAAF1RyCp0H38Obzosz7NtAnVc8vszgRoA0BcEaDQuZGieJtsuwnR8\nJn1BJFgDAGJDgO6J2KvPsYbmSdL2EqTjN2vf73vAvvvKQen/Ct0KAEAVBGjMpW/BOY8g3X/zfrls\nMoDH/kUXANAMAnQNXY+DjvFDue/BOY8gPV4xvr+AKl578auZ83z05ekOWgKMBwEalQwtOOcRpAH0\nQZnQPG1+AnW7jn5W7TPkwdlhf7YOUbAAbWaX3f3ShGlrkrYlLbn7RrctQ5GhB+e8e+d3CNEYNPrg\n/qkamsusiyDdnKqhedKyhOluvfbiV/rXGssFSQhmti5pbcK0ZUly983s/bGK4eflsYXn1L3zO6N9\n7hg2+uD+aTI859fb1rrH4uhnC3OF57bXh6fS/T37r64gFWh33zCztyZMvijpWnJ7S9KqpJudNAx7\nEB53UY3G0NAH98usD/k3Xvi41Ho++PrVqdugGl3drKB77PPpn6P3z0xe/uhnC1SjG9LGl8QYU8Gi\npDuZ+6dCNWSSrg4gDFl9JjzvxeuBEYm+Dx6TaR/8b7zwcenwXGZ+KtHVTAvPxz7fmRmey8xHJXp+\nbe3XvfzLmNm6md0wsxuPdh6Gbs7gEBaLMaQD2JXtg7+9/8fQzRmsWeG5rnmWxa5JwXZWID7+aXFm\nmbYcIbq+Nr8UtjKEIxlfl3fH3a+WWHxb0snk9qKk2/kZkoNaNiTp+YOnvW47Yxaq+kxAnI0hHYhd\nl33woZe/N8g+OGZFAfjHx349cf5f3P9B4TqKhnQwlKO+SQE4H5rT+/fOHWm9TWPW9i8qrQToOkdt\nm9miu29LuiJpJXl4SdJmk21DMYJzNYRoxKzLPnjhEWM1Q5sWntPpRSEa9ZStCE+qNuenZ4P0sc93\nCsdF8x5rR/pl9D/UWDbUWTjWJK3kqiTXJcndbybzrEraTu+PSdfVZ8JzPbxu6Ks2+mB+Zu4Gwy+A\nOIQ6C8dVSVdzj13I3I72vKNdXoGwC4TA+VCJRh+10QdTHQOqYxhHf/HJH5kuq8+E52bwOmLMdg4S\nnrtUNG75F/d/MHWIRpXhG4x/nq3s/n7v3JGpAblo2qTT2vEea8cHX7869fSO03Ap75Ei9DWLSjSA\nrnzw9auFQzmqBOW6oQGT3T+zUHggYdkq87RzQqO6j7483eqBhPy1Kmh7+EYMVx1EfXwpAdCkadXg\neQLwpGWpPpc3qSJ8/8xCrSA8bRmqz/W1uU9TgR4hgl57qESjDN6DKGtaFS0Nwk1ciRDVPTi7M/Hg\n2WwgnnR6uzJBm/A8v7Yq0QTokoZSfeaDu32EaPA+Q5NmBYB5gzGV5/qmhehUnYo0wblZ6T7eZJAm\nQI8IH+pAc3g/oUttBIDselFfGnabOJUjwbldTe7vBOgShlB95sO+W1Shh4P3DmLSVJAmODcvG36r\nhmmCc/8QoEeAABAGIbp/eK+gL/IBeFagJjB3i0A8fAToGYZQfUY4hOh4EZYxJARkoFsE6IEjJIRH\niI4D7wUAQFMI0FP0vfpMYIgHIbp77P8AgLYQoCcgPKNphOj2sd8DALoQ7NPczC7PmmZm69216Km2\nwzPGi4DXvHvnd578Q3kx98EAELsgATrplNemzLJuZrckbXXUpCe6CM9Un4H5EJrnE3MfDAB9ECRA\nu/uGpnfMb7v7eXff7KpN0jAqzwSK+PE3qofQ3JxY+2AA6ItYB2SeNLNVM3u3qw12FZ7brD4TLPqD\nv1V5hOYgOu+DAaBPogzQ7r6RVD5OmdlqfrqZrZvZDTO78Wjn4dzbG0J4Rv8QCiej2hxWlT742/t/\nDNBCAAirlbNwTDjw5I67Xy25bDrvbUlL+XmSnx83JOn5g6d9nrYOYdiG1K8wduzc76dOv//p8x21\nJDzOzLFXn/bjmHXZBx96+Xtz9cEA0EetBOikc63EzBbdfVu74/JuJA+fktTaGLwuw/NYh27MCstl\nlhl6oCZEx70P91Ff+mAA6Ksg54E2szVJK2a2nunor0u64O6bZrZmZpJ0291vNr39rqvOYxy6USc4\nl1nXUMP0WEP0mIJzk++JeYXugwGg74IE6OSnwau5xy7kprdiaOE5pgDSRUBItzHEID2mEB3TftuG\nmMJykZB9MAAMwWiuRDiUsc5ZMYWQrgPDUIP00EN0TPtsk2IPzACAZg0+QIcMzmMY9xw6OAwxSA8x\nRMeyvzYl9H4PADF4/MVRHXjpQehmBDHYAB264jz0cc+xBYhj535PiI7QkIJzbPs8ACCcQQXo0KE5\nNfRxz7EGiaFVo/scokPvo02JdV8HAITV+wC9c+TZaIKzRHiOwZCq0X0L0aH3zyb0YR9HM8b883Of\nHHjpgR5/cTR0M4A9eh+g0Z0+BYuhhWhJUQdpgjMAjNNYv4gSoBs01OpzX4PFkEK0FGc1muCMvhvr\nh3/fUIWOW/ZvM5b3EwG6IYTnOBGi22tH3/V930ZzCNH9kP8bEajjNJa/CwG6AYTnuA0xREthhnQQ\nnDFUhOj+4e/VnbGE4ioI0HMa6unqhhYyhhaipe6C9BBCszS8fRrNI0QDxdL3BUH6KQL0HLoIzyHC\ny1CDxhBDtNRekCY4Y4zSgECQBvYjSD9FgK6hq6oz4bl5Qw3R0t79pU6YHkpgzhr6/oz2EKQBTEOA\nrojw3H9DDtGpIYbhKsayL6N9BGkARQjQFQx1vPMYjSFEjxXhGW3o+jRdIX8i58sCMFuQAG1m68nN\n8+5+qWD6mqRtSUvuvtFp4yboMjxTfQaqYx8ur499cEyGPv6z6edHIMcQdX4eLDNblbSZdMpLyf3s\n9GVJcvfN7P2QCM/DNNbnPUR9+1v+9cuf7PnXpT72wei3x18cffIPGIoQFeil5N+GpK3kdtZFSdeS\n21uSViXd7Kx1GV0P2SA8d4+hHP3Xh32465A8Q2/6YAwPY8oxFJ0H6NzPgcuSruRmWZR0J3P/VH4d\nyc+P65J06PBi002URHgGYhf7/htZaH6i6T74mRMnmm4iRoBzbqPvgh1EmPwseNPdK1c2kg+ADUk6\n/vxZb7JdYzlQMPbw0SWq0P0T6/4ba2gu0lQffOjl7zXaBwNAH7QSoDMHqGTdcfermfurRQevaPfA\nlZPJ7UVJt5tuX5GQwXnspxyLASG6P2IMz7EF5z72wQDQJ60E6FlHbZvZuru/l9xedfdNM1t0923t\n/py4ksxD/gJ+AAAgAElEQVS6JGmzjTZmjS08xxhAgDJi23djC86pvvXBGCeGcaDPQp2F47KZ3TKz\nu5lJ1yUp/TkxmW+7zs+LZd195SDhGU/w2sQtpr9PiLNnNCWmPhgA+irEQYSbkvYddeLuFzK3Wz3v\naAzjnAnPcWIoR5xi2Xf7GpqzYuiDAaDvRnMlwhhCc4oxz0B5hGcAiMuBlx6M/rzegw/QMQVnKVx4\njiWE9AFV6HjEsN8SnAEAeYMM0LGF5hThGSgvhv2W8AwAKDKIAB1rYM4iPPcLVeiwYthvCc8AgEk6\nPwtH0749bKGbMBNjnoHyCM8AgNgNogIds5DhOYYg0mdUobsXep8lOAMAyiBAt4jw3H+E6O6E3mcJ\nzwAwHkc/m28QBgG6JQzbAMojPAPjNNSrEb724leV5v/oy9MttQSpeQNzHgG6BaHDc+gwMjRUoYHJ\nFh49/WB6cJbCAcaramguWpYg3aymQ3MWAbphhGegmtD7LNXn5mQ/rAjTGJN5wnPRegjS82szPEsE\n6EaFDs9oD1XodhCeh4uqNMagqeBctF5CdH3TwvOxz5vpkwjQDYkhPIcOI0AVofdXwnM3CNIYqzde\n+LjUfB98/Wrh44ToeiaF56aCcyrIeaDNbD35d3nC9MvpfN22rB7C8zjwGjcn9Gs59vDcZB/8zDfl\nPpja/jkV6Nq06nPZ8Dxr3rYq3ENVNTwf//Shjn/6sNa2Ou/RzGxV0qa7b0haSu7nrZvZLUlb3bau\nmnvnd6IIzwBQVlt98LHPdxqv8ACxaio8z7MM9qoSnucJzqkQJYElSWmHvZXcz3vb3c+7+2Z3zaom\npuAcupoHVBF6fx179Vkt98GEaAAxmzc4p8zdG1lRrY2bXZN0yd1v5h5f127Hvuzu7xUsty4p/Wnx\nFUnzfCK+IOnrOZbvu7E/f4nXYOzPXwr/GvyFu3+3643SB0dh7M9f4jUY+/OXwr8GlfvgYAHazJYl\nXXT3S1PmuSzpWpuVaDO74e4rba0/dmN//hKvwdifvzTO14A+OA5jf/4Sr8HYn7/Uz9eglbNwTDjw\n5I67X83cXy3quJNl03lvq/jnRQDABPTBANCuVgJ0cnDKRGa2nv4saGar7r5pZovuvq3dnw1vJLOe\nkhTtOGgAiBF9MAC0K9RZOC6b2S0zu5uZdF2Skp8KV81sTdLt/Ni8Fkz9oBmBsT9/iddg7M9fGtFr\nQB8cnbE/f4nXYOzPX+rhaxD0IEIAAACgbzizPTByky6mAQBoH31wPxGgsceY3shmtmZmq3254mUb\nkue+Frodocy6Ih/QtTHti/TB9MF97oMJ0Ik+/xGbMqY3cnIKr3S855P7Y5McbBb1FT/bUvKKfOgI\nfTB9cNgWhUEf3N8+mACt/v8RmzKyN/JFSdvJ7S09vTIbxqPMFfnQAfrgXfTBGJle98GtnMauh5aS\nf2nn1as/ImpZlHQnc/9UqIYgjNyp3pYlXQnVFtAHjxB98Mj1vQ8mQKv/f0QA9SU/Hd/s4HRtmIA+\nGBivvvbBBOiMvv4Ryyp5dbKx2JZ0Mrm9qN0rrmGcCq/Ih+7RB48KfTBSveyDRxOg57m07VDMujrZ\nyFyRtJLcXtJIr7aWXCxjJbky3ej2j6Ir8oVu01DRB9MH59AHiz64z30wF1JJZHfevv0Rm5K8kX8m\n6dIY3sjJB/qWpKUxPF/slRyo9nPtjsM8KemtMb7vY0EfTB8cuj3oVt/7YAK0+v9HBIA+ow8G0DcE\naAAAAKACzgMNAAAAVECABgAAACogQAMAAAAVEKABAACACgjQAAAAQAWjuZAKUFZyLtYl7Z6f9IeS\nfuru22FbBQDjQB+MPqACDWSY2VJyZbS0s76SdtxmtpScrxYA0AL6YPQFARrIcPet5OYFSZvufjMz\neU3SYvetAoBxoA9GXxCggQwzW05uLrn7dua+JHFlNABoEX0w+oIx0MBeq2a2JOla8lPhndANAoAR\noQ9GL3Apb6AkM3tX0ilxQAsAdI4+GDEhQAMAAAAVMAYaAAAAqIAADQAAAFRAgAYAAAAqIEADAAAA\nFRCgAQAAgAoI0AAAAEAFBGgAAACgAgI0AAAAUAEBGgAAAKiAAA0AAABUQIAGAAAAKjgQugFATMzs\nrqQ7uYffcffNZNo5d98O0DQAGJSkT93Sbp97UtKmu18qsUy+j75Av4yumbuHbgMQjWkhmQANAM3J\n96lmdlmSpoVo+mHEgiEcAAAgBlckLYduBFAGARoAAMTgJ5KuSZKZrZvZh8m/1cDtAvZhDDSw33Uz\ne3LH3S8EbAsADFm2v73i7u+Z2bKkt9K+NwnRP8oM28guM3PcNNAGAjSw348YXwcAnSjqby9Kej9z\n/4qkv5O0MWUZoFMM4QAAADFZzN0/FaQVwBQEaAAAEJOfS3pHksxsUdKapM2gLQJyGMIB7Pdpdgy0\npEvuvjFpZgBAc5Lz7i+b2S1J29rtg7dCtwvI4jzQAAAAQAUM4QAAAAAqIEADAAAAFRCgAQAAgAoI\n0AAAAEAFvT8Lx8GFI37kwHOhmwEAtf3hm//va3f/buh21HHsxEE/eebwnsf+8Pjwvvn+9M2zU9fj\n31DPQdwWHs2e55lvJiz7550nt3cOTd7Xv53+Ntld/uDseYbMnt2ZOv3ws3v/CM8d+NO+eZ5f2PvY\nr/6fbyr3wb0P0EcOPKe/+ncXQzcDAGr7j5//r/8Wug11nTxzWP/z1f/2yf0Pvn61cL6Pvjw9dT2P\nvzjaaLuAph39bPaXvGOfTw93s9w/M3sbD87Ot42+O/DSg6nTX3vxq32PvfHCx3vu//jYr/fc/6+/\n918q98F85QcAAJihTHAtE4DnWXbs4bkpv7j/g7nX0fsKNAAgnD88Pjyx6pyaVX0GhuT+mYXKleh5\ngvfYPP7i6NQq9Edfnt5Xhf7g61f3VaH3huj/Urkd/MUAAK0pE54ZvoG+KFsBvn9moVQoLjtflW2j\nuN/54OtXZ37Zr4IKNACgcWWrzoRn9M2DszulxkNLzVWWCc97zapCS8WVaGn/cRr5ynRZwSrQZnZ5\n1jQzW++uRQAwHk33wR99eXrPvzIIz+irLgMt4blYmf6jTH9UtyodJEAnnfLalFnWzeyWpK2OmgQA\no9FkH/ynb56tNcaZ8Iy+e3B2p9Vw2/b6h6BsP1Lli31ZQYZwuPuGmb01ZZa33f1qZw0CgBEJ3QcT\nnjEkacgtO6yj7PpQTpnhHKkmQ3SsY6BPmtmqpGV3fy90YwBgZFrpgwnOGLJs8K0apgnN80n7lrJB\nuglRBmh335AkM3vTzFbdfTN0mwBgLJrsgwnNGCMCcRhdBunoAnQyNu9O8vPhbUlLE+ZZl6TDzxzv\ntoEAMGBV++ADLzz/5HHCMoAYdBGkownQZrbo7tvaPWjlRvLwKUn7Kh9JdWRDkp4/eNo7ayQADFTd\nPvjQy99zgjOAGOX7piYDdZAAbWZrklbMbD39qVDSdUkX3H3TzNbMTJJuu/vNEG0EgKGiDwYwRk1+\n2Q91Fo6rkq7mHruQmw4AaAF9MADMh0t5AwAAABUQoAEAAIAKCNAAAABABQRoAAAAoAICNAAAAFAB\nARoAAACogAANAAAAVECABgAAACogQAMAAAAVEKABAACACgjQAAAAQAUEaAAAAKACAjQAAABQQbAA\nbWaXp0xbM7NVM1vvsk0AMBb0wQBQX5AAnXTKaxOmLUuSu29m7wMAmkEfDADzCRKg3X1D0taEyRcl\nbSe3tyStdtIoABgJ+mAAmE+MY6AXJd3J3D8VqiEAMEL0wQAwQ4wBeiYzWzezG2Z249HOw9DNAYBR\nyfbB397/Y+jmAEDnYgzQ25JOJrcXJd3Oz+DuG+6+4u4rBxeOdNo4ABi4Sn3wM8e+02njACAG0QRo\nM1tMbl6RtJTcXpK0GaZFADAe9MEAUF6os3CsSVrJnSLpuiS5+81knlVJ2+l9AEAz6IMBYD4HQmzU\n3a9Kupp77ELm9kbnjQKAkaAPBoD5RDOEAwAAAOiDIBVoAIjJw9fPlJ73yK8+b7ElAIA+IEADGJ0q\ngbnMsoRqABgXAjSAUZgnNFdZN2EaAIaPAA1g0NoMztO2R5AGgOEiQAMYpK6D87TtE6YBYFg4CweA\nwQkdnvMevn4mujYBAOojQAMYjNiDasxtAwCUxxAOAIPQRDi9+8rB0vOe+ORRrW0wRhoA+o8ADaD3\n6obnKoF51rJVA/WeNpOlAaBXCNAAeq1OeJ4nOM9aZ93KNACgPwjQAHqranhuIzhP2gZBGgCGiwAN\noJeqhOcugnPRNscQohceSUc/2z0e/cHZncCtGY8DLz2ovezjL4422BJgnIIEaDNbk7QtacndNwqm\nX3b3S2a2XjQdwLi1HZ7vnZ8eBI/fKncCo1hDdFt9cBqkswjVzZknNE9aD2EaqKfz09iZ2bIkuftm\n9n7OupndkrTVZdsAxK9seL77ysFK4fne+Z0n/6rMO2v+ENXvabrug49+tvDkH+o58NKDxsLzpHW3\ntX5gqEL0aBe1W/mQdjvn1YJ53nb382kHDwBStfBcRpXQXGY9PRGsDyZIV9N1sCVEY4zq9kkhhnAs\nSrqTuX+qYJ6TZrYqadnd3+umWQCGoEx4bivs3ju/Uzi8I7KhHMH74KOfLTC0Y4a6Yfa1F7/ac/+j\nL0/X2i5DOzBETX6Bj/IgwnTMnZm9aWar+SqIma1LWpekw88cD9BCAF0rU30OGZ7z688H6chC9FRV\n+uCDR0/o2OdPX9P7Z8p9QHHg4WRlw3M+LJeZp2ygPvDSA0J0Qzg2IKy2fvUKEaC3JZ1Mbi9Kup2d\nmHTMd9z9ajJtKb+CpHPfkKTnD572VlsLILgmwnOV4Hzs3O9LzXf/0+enbq/swYYda7QPPnbqe3v6\n4Kphmmr0XmXCc5ngPGvZMkGaEF3frNCWn857oHltDxcLEaCvSFpJbi9JSg9kWXT3be2OybuRTD+V\nTgeASZoIz2VDc9Ey04L0LAEu6d1ZH3zs853SFWk0E57feOHjPfc/+PrVqeuZFaQJ0dXUDW38ItOs\nLo616DxAu/tNM1tJxtdtu/vNZNJ1SRfcfdPM1sxMkm5npgMYoVnV53nDc53gPGkd+SCdr0LHMIyj\n6T544c87Ov7pQ0nSvXNH9k1PK9LTgjRV6HImhed8aJ40rShMv/biV5XHSaNYE6GN98L8ujpQOcgY\n6KLzirr7hcztq922CECM6lymO6uL8Jxf36wQnRWg+iypvT54WpCeZezBYVb1uU54Lpp3UkUa82ky\ntI39vTCPLs/yw29rAHprWvV5Wng+du73jYfnaevu0SnuGpEG6azs2Gg0o0p4nrbMrGEhnN5uujZC\nG6d7bN+xz3ee/KsjyrNwAMA8Qzdmhecy/vrlT0rN98+/fWXidorGRmeHcTx8/UywKnRT7M/f6MCt\nz/X4/N6/1/FPH9aqRGM+Pz72a0nSL+7/IHBLgHg18YWeAA1gUOYNz2WDc37+SUE62650KEcMY6Ex\nTARnoBv8RgAgOvMeOFhkVnj+65c/qRye88vP2mY23KfPYd5x3qH5oWf3VZ+leuOgAaALTZwdiAAN\noFfqDN0oE56bUBTCy4TooSE8d6POAYFFy8w6CwensZuujQP+OIiwfffPLDz5VwcBGkBU6lZk6x6o\n11R4nrbOsYToe+eOFIbnMh9SYw8Ms0LqpJBbJURzBo72NLn/jv29MI8uXzvGQAPojaaHbrQRnqe1\nI3tQYcRXKqxk59DC1IpzmeoOgaGcj748XXjGjGwwnnXe5/z6ZqH6XN6Dsztznz2D98L8mvg7lEGA\nBhCNpqvPTYTnf3/i/5447R/u/tXU9ecPLCwK0UPFFQire/zF0ZmnjJsUolNlq8xcPKUdaQCuGuAI\nzs2q+3eoggANoBeaHOrQRHjOTp8UpMuE6CGpE5oJDnuVDdGpWedwnrZsmbagnvx+nQ9y7PfdaLMa\nTYAGEIUuq8+zzArORfOPNUR/+yzBuWllQnSqrUoy4blZ7O/hZF/7JsM0ARpA9EJUn9s2pBBdFWFi\ntiohuuntAkPVZN/DIDUAvRVD9bnMcpNCe5uXFI/Rg7M7hOcKHn9xtNNAS3gGyqMCDSC4ri4mEkv1\neUwIzPNLg20bFWlCM1BPkABtZmuStiUtuftG1ekAxmNI50mORVt9MGG5XdmwO0+YJjQD8+t8CIeZ\nLUuSu29m75edDmBYur5wytg13QfvHGRoRgjp8I46/wDML8QY6IvarWxI0pak1YrTAWCiUGOKp50T\nOjL0wQAwpxABelHSncz9UxWny8zWzeyGmd14tPOwhSYCiEFfhm/0KDxLDffB397/YwtNBIC49fIs\nHO6+4e4r7r5ycGHyJWQBoKp/uPtXlQJxmXnz54Huu2wf/Myx74RuDgB0LsRBhNuSTia3FyXdrjgd\nwEB0dfaNOtJgXHR6uioBO8LwTB8MAHMKEaCvSFpJbi9JSg9UWXT37UnTAWBe//zbVyqfym6e4Rll\nwnOAi6nQBwPAnDoP0O5+08xWzGxV0ra730wmXZd0Ycr0XptVaTvyq887agmALswKz6GuQjjWPhgA\nmhTkPNBF5xV19wvTpvdN1Z+ms/MTpoH67n/6/NQzcdSpQlcRa3DOGkMfDABt4kqEDWpqPGe6HoI0\nhizk+Oc2QnSkwzUAAC0gQM+pzRDw8PUzhGigJWngnSdIlz1AkOAMAMNCgK6pq+oZIRqobtYwjqwq\nQbrqGTUIzgAwTAToCkL95EyIBqqrEqKlZk83VyY4H7/Vy9PwAwBEgC4l5nPVAmN1/NaC7p3fmTpP\n1RA9D0IzAIwHAXqKmIIzVWignjZDNKEZAMapVM9uZufM7O/N7PvJ/b9ps1GhPXz9TFThGUCxsuH0\n/qfPNzYeOV3XrPUdv7XQWHgeWx8MALErW4H+W+2eZP8tM7sm6U1J/9RaqwIhNAPDlg29ZarSdUJ3\n2dB84pNHVVY7ij4YAPqibID+paRb7v5LM/uRdi/vOhgEZyBOJz55pLuvHJw6T5mx0EWaPENGmdBc\nMTDnDboPBoC+KRugtyStSvond79uZi02qTsEZ2AY6oboebdZxpzBOTXIPhgA+mpmgDaz59z9U0mf\npo+5+/VWW9WyJoLzrKpYVkMfoACmSANtm0G649AsaZh9MAD0XZkK9E/M7Iq7/4uZ/aUkd/d/abth\nbZgnOFcJzJOWJUgD1ZUZxpHVdJCuciBg1fd4yTPrDKYPBoChKBOgb0haMrOtZPzdG203qml1g/M8\noXnS+gjRwK4jv/q89HuzaoiW9gffWYG67hkzqryna56Ksvd9MAAMTZkAvSRpW9J7ZnZO0jVJH8yz\nUTNbS9a55O4bBdMvu/slM1svml5WneDcdGguWj8hGuhe0+dj7iA4p3rbBwPAUJX5RNly95+5+//g\n7v+ddg9mqc3MliXJ3Tez93PWzexW3W3VOY/z3VcOth6es9sCUE0MXzxPfPLoyb9Zjvzq8yf/5tS7\nPhgAhm5mBdrd/9HMvu/uv0nG352fc5sXtVtBkZ4eWX4zN8/b7n61ykpjGaYBoLwqwzikpyG66/dt\nh9XmffrSBwPAmJQ6jZ27/yb5/5faPR/pPBYl3cncP1Uwz0kzW5W07O7v5Sea2bqkdUk6dHgxyqEa\nZbYfQ0UN6KM6Y6Krrr+qpoNzVsx98DMnTszZHADon7Lnge5UOubOzN40s9X0p8bc9A1JOv78Wa+y\n7tDBGcBeVavQqaar0XW/0LYZnEOp0gcfevl7lfpgABiCVgJ0Up3Iu5P8JLgt6WTy2KKk2wXLpvPe\nVkNX3Jr3QzZ/BH/TByQBqCcffMu815v49Sfm4BxjHwwAQ9JKgJ5x1PYVSSvJ7SVJ6YEsi+6+rd0x\neTeS6afS6fOoGp7LnD+260Ad84c1MK+6VegibQ6N6sv7MLY+GACGpvMyqrvflKRkfN12el/S9WT6\npqTV5DRLtzPTK6t6Zo1753dqX3yh68sIA0MTczht6GwaUeiyDwaAoQoyBrqoOuLuFzK35z76u2pw\nbsK98zsM7QDm0GQluglDCc15XfTBADBkg0x7ZcPzPBXnaesEUF8MoXVIFWcAQPOiPAvHPKqE577g\ngxxjE6ISzfsMAFDWoAJ0mfDcp+AMjFlXIZrgDACoalABOhaMhQaakYbbpoM0oRkAMI/BBOihVp/5\noAf2vg/qhmneSwCApgwmQM/Sx/AMYD+CMAAgtNEEaAAAgC4ceOlB6Xkff3G0xZagLQToQMpcLY1K\nGwAA8asSmGctS6Duh8EE6BOfPKp8yW4AAPro6Gf1D1R/cJYhjU2ZJziXWSdhOl6DCdAx4QwcAIA2\nzBOci9ZBmK6njeA8bTsE6fgMKkD3pQrN8A0AQBVNBOdp6yVIl9NVcJ60XYL0fJp8HwUL0GZ22d0v\nTZi2Jmlb0pK7b3TbMgAYvqb6YHuW4NWmtoLzpO0QpCerEp5fe/Gr0vN+9OXpSm0gRFfT1nsoyFgD\nM1uXtDZh2rIkuftm9n5ZZaq7ADBmTffBB156sO8f5nP0s4VaH/zHPt/Z86/OdrFf2X36tRe/qhSe\ns8uUXY73Vzl130NlBXmnJBWNrQmTL2q38qFkntVOGtWQWeOfGb4BILQu+mCCdH1VPvRnBeY6gZoQ\nvVeZ/bhOcJ5nPby3putiH45xDPSipDuZ+6emzbzw8Jt9j/VlLDQARKhSHzwLPzm3o051Obvc/TPT\nA8bRzxYYzqHy4bnIGy98XGobH3z9auE6qwztwFNdfQHs5ddMM1s3sxtmduPRzsNoKracfQPAGGT7\n4G//8MeZ81MtK6/Mh/+08Hz804d7/tVZR5W2jF1ReH7jhY9Lh+dp88+qRPO+2q/qPlt3qJPUUgU6\nGV+Xd8fdr5ZYfFvSyeT2oqTb+RmSnx83JOn5g6dd2h328PD1MxNXevzWQvDLeTN8A0AXuuyDD58/\n42XaRCV6trrheVpQTqfdO3ekcF2zKtGoZlJw/vGxX++5/4v7P+iiOaNSNjzXDcx5rQToOmfOMLNF\nd9+WdEXSSvLwkqTNOm3oehgH1WcAsYihD0bzqobn/HxFIRqTzarwlhmrnA/O2cfzIfqNFz7eN5yD\noRzNaio8S+HOwrEmaSVXJbkuSe5+M5lnVdJ2er/vqD4DiMUY+2AUh+1ZgYJhHO2ZFK7RD0EOIkx+\nRryae+xC5navzv1M9RlAnwytDwaAro0q+YUKulSfAQBtqDIso2jeWWOgx3wmjibG7E8b61xmHDTD\nN5rV5Jj/GE9j1ytUnwEATXhwdmfqkIn7ZxYKh1zcO3dk6ljoSSGbAwjn89GXp/eNg/7g61f3HUhY\n9oDBotPZoXmT3kdVEaDnUCY8c2VEAGimmofpIbppY64+z6MoRJdZJo/qczWzvoBmNfHlcRBfPx++\nfmbqKeyymqoYN1l5ZvgGgCEjPJdXJrTeP7MwVwCYd/kxmbXvTgq5H3z9aqmK8qT5yoRn3lf7dfml\nr/cV6J0jz3a+zbLhmeozgDHjA76espW0NASX/Tm6Smim+lxeGnaLTmtXZ1gG4Xk+VSrR8+h9gK5j\nnouqNB2eqT4DGBo+3OfX9c/R+W3jqcdfHC111b9pQboMhmw0J92H2wzSgw3Qsy6iUjVEc7AgAOzn\n3ywQmFvSVSUtuz0UKxuipb1BeFaYrhqaea9V02aQHmyALiMNxdOCdJ3gTPUZANCELippBOdy0vBa\nNkhLzVaVCc/1Zffxpt5Low7QqSary4x7BgA0rY0gTXCup0o1uqntoTlN7fcE6AZVCc9UnwEAVeU/\n/KsGakJzM+pUo+tuA3EiQDeEyjMAoGsE4rCyIbeJME1o7o/BBugTnzyaeSBhk9uqguozAADDMin8\nFgVrgnL/DTZAd4XwDAAAJiEsD1Owc7OZ2eVZ08xsfZ5ttD2sgmEbAPqqiz4YAIYqSIBOOuW1KbOs\nm9ktSVvzbquNkHvik0e11kv1GUAMuuyDAWCIggRod9/Q9I75bXc/7+6bTWyvbuBtcj2EZwCx6LoP\nBoChifXyeifNbNXM3m1ypXUD8LwBnPAMoGda6YMBYCiiPIgwqY7IzN40s9V8FST5+XFdkg4dXqy8\n/i7HLhOeAfRNlT74mRMnArQQAMJqJUBPOPDkjrtfLblsOu9tSUv5eZLOfUOSjj9/1udsLgAMSpd9\n8KGXv0cfDGB0WgnQafWiCjNbdPdt7Y7Lu5E8fEpSL8fgUXkGEAp9MAC0K9RZONYkreSqJNclKfmp\ncDWZ57a73wzRxnkQngHEbOh9MAC0LcgY6OSnwau5xy7kpvcS4RlA7IbcBwNAF2I9C0cvEZ4BAACG\nL8qzcPQNwRkAAGA8eh+gFx5+MzXAPnz9TGvbJjgDAACMT+8D9CyTQm6dYE1gBgAAwOAD9CSEYQAA\nANTBQYQAAABABQRoAAAAoAICNAAAAFABARoAAACogAANAAAAVECABgAAACogQAMAAAAVEKABAACA\nCoJcSMXM1pOb5939UsH0NUnbkpbcfaPTxgHAwNEHA8B8Oq9Am9mqpM2kU15K7menL0uSu29m7wMA\n5kcfDADzCzGEY0lS2mFvJfezLmq38pFOXxUAoCn0wQAwp86HcOR+DlyWdCU3y6KkO5n7p/LrSH5+\nXJekw88cb7qJADBYTffBz5w40XQTASB6wQ4iTH4WvOnuN6su6+4b7r7i7isHF4600DoAGLam+uBn\njn2nhdYBQNxaqUBnDlDJuuPuVzP3V4sOXtHuT4cnk9uLkm433T4AGDL6YABoVysBetZR22a27u7v\nJbdX3X3TzBbdfVu7PyeuJLMuSdpso40AMFT0wQDQrlBn4bhsZrfM7G5m0nVJSn9OTObbrvPzIgCg\nGH0wAMwvxEGEm5L2HXXi7hcytznvKAC0gD4YAObHlQgBAACACgjQAAAAQAUEaAAAAKCCzsdAAwAA\nAF068NKDZtfX6NoAAKNiz+6EbgIAFGo6NO9Zd2trBgCMwqQPqcdfHO24JUBcjn42faTsg7N8AW1a\nm8tqCYEAACAASURBVKF5z3Y62QoAYHSyH2SEaYzBrMA8a34CdX1dBecn2+t0awCAUUo/3AjSGJqq\nobnMugjS1XQdniUCNACgQwRpDEWTwXnSugnSs5UNz6+9+NXEaf9aZ7s1lgEAYC4HXnpAiEZvtRme\ni7ZDkC5WJjxPC85zbbuVtQIARuHws9/s+4D66MvTpZYlRKOPqobnY58Xh9/7Z8qv5+hnC4TonHnD\n8xsvfPzk9n+os/0aywAAMFH6oVUmSBOi0Sdlw/Ok0Fw0T9kgTYh+qm54zobmeXElQgBAK9r66RQI\nocnwnJ+/7DJdDR3pu7bDsxQoQJvZevLv8oTpl9P5um0ZAAxfk33wcwf+NHV6mRAd4gh6oA1Vw3NT\ny47JrP6iTHj+8bFf7/lXR+cB2sxWJW26+4akpeR+3rqZ3ZK01W3rAGDY2uiDm67sALEpU/mdFICP\nf/qw8F+VdSA+ISrQS5LSDnsruZ/3trufd/fN7poFAKPQSh9MiAb2mxSUZ01Dc4qqz00wd29kRbU2\nbnZN0iV3v5l7fF27Hfuyu79XsNy6pPSnxVckfTJHM16Q9PUcy/fd2J+/xGsw9ucvhX8N/sLdv9v1\nRumDozD25y/xGoz9+UvhX4PKfXCwAG1my5IuuvulKfNclnStzUq0md1w95W21h+7sT9/iddg7M9f\nGudrQB8ch7E/f4nXYOzPX+rna9DKaewmHHhyx92vZu6vFnXcybLpvLdV/PMiAGAC+mAAaFcrATo5\nOGUiM1tPfxY0s1V33zSzRXff1u7PhjeSWU9JYhw0AFRAHwwA7Qp1Fo7LZnbLzO5mJl2XpOSnwlUz\nW5N0Oz82rwVTP2hGYOzPX+I1GPvzl0b0GtAHR2fsz1/iNRj785d6+BoEPYgQQHhmdnnaOFgAQHvo\ng/uJS9pgj0kXVhgiM1szs9UxX7Anee5rodsRyqwLigBdG9O+SB9MH9znPpgAnejzH7EpY3ojJ2cg\nSH+ufnJ/bJKxsqO8YFHJC4qgI/TB9MFhWxQGfXB/+2ACtPr/R2zKyN7IFyVtJ7e39PTCEhiPMhcU\nQQfog3fRB2Nket0Ht3IWjh5aSv6lnVev/oioZVHSncz9U6EagjByZ6pYlnQlVFtAHzxC9MEj1/c+\nmACt/v8RAdSX/HR8s4OzTWAC+mBgvPraBxOgM/r6Ryyr5MUVxmJb0snk9qJ2LxiBcSq8oAi6Rx88\nKvTBSPWyDx5NgJ7nylxDMeviCiNzRVJ62dAljfRiEcm5fleSC2uMbv8ouqBI6DYNFX0wfXAOfbDo\ng/vcB3Me6ER25+3bH7EpyRv5Z5IujeGNnHygb0laGsPzxV7JgWo/1+44zJOS3hrj+z4W9MH0waHb\ng271vQ8mQKv/f0QA6DP6YAB9Q4AGAAAAKuA80AAAAEAFBGgAAACgAgI0AAAAUAEBGgAAAKiAAA0A\nAABUMJoLqQBlJediXdLu+Ul/KOmn7r4dtlUAMA70wegDKtBAhpktJVdGSzvrK2nHbWZLyflqAQAt\noA9GXxCggQx330puXpC06e43M5PXJC123yoAGAf6YPQFARrIMLPl5OaSu29n7ksSV0YDgBbRB6Mv\nGAMN7LVqZkuSriU/Fd4J3SAAGBH6YPQCl/IGSjKzdyWdEge0AEDn6IMREwI0AAAAUAFjoAEAAIAK\nCNAAAABABQRoAAAAoAICNAAAAFABARoAAACogAANAAAAVECABgAAACogQAMAAAAVEKABAACACgjQ\nAAAAQAUEaECSma2Z2XLodgAAgPgdCN0AoClmdlfSndzD77j7ZjLtnLtvFyy3JGlJ0raZLbn71pT1\nbyXbOClp090v1WjThaJ2AACAfjB3D90GoBEzQnLhNDNblHRduyF3W9Kiu79ZZv1mdlmSpoXoaW0C\nAAD9xBAOjFoSbN+WtCJpWdJbFRa/kiwDAABGhACNUUsq0D+TtCHpfe1Wo8v6iaRryXrWzezD5N9q\n8y0FAACxYAgHBiMzRvkJd7+QmTZxKIWZrUvacvfNCuu/4u7vJQcfXk6HfpjZh5J+5O7bBcvMHDcN\nAADiRoDGYNQZA93E+pOx0P/Z3a8m99+VtO3uG4yBBgBgeBjCAcxvMXf/VJBWAACAThCggfn9XNI7\n0pMx1WuSJg4FAQAA/cZ5oDE0n5pZ9v4ld99oc4PJeaaXzeyWdk+Fd2nSuaQBAED/MQYaAAAAqIAh\nHAAwQumFgCZMWzOz1eTsNACAHAI0AIxMEozXJkxblnaHJmXvAwCeIkADwMgkxwVMGqd/Ubtj+ZXM\nw4WBACCHAA0AyFqUdCdzn9MyAkAOARoAAACooPensTu4cMSPPHM8dDMAoLY/PP7d1+7+3dDtSGxL\nOpncXpR0Oz9DMoZ6XZK+853vXHj11Ve7ax0ANOzDDz+s3Af3PkAfeea4/uqFt0I3AwBq+49f/m//\nFroNZraYXHL+iqSV5OElFVwUKBlDvSFJKysrfuPGjc7aCQBNM7PKfTBDOABgZMxsTdJK7jR11yXJ\n3W8m86xK2k7vAwCe6n0FGgBQjbtflXQ199iFzO1Wr94JAH1HBRoAAACogAANAAAAVBAsQM+4jOzl\n5H8uIwsAAICoBAnQ0y4jm1g3s1uafKUsAAAAIIggAXrGZWQl6W13P+/u+06fBAAAAIQU6xjok2a2\nambvhm4IAAAAkBVlgHb3jaT6fCo5F+keZrZuZjfM7MajnYcBWggAAICxii5AJ+E4HR99W7tXwtoj\nCdgr7r5ycOFItw0EAADAqEUToM1sMbm5paeXjj0liWvEAgAAIBqhzsIx7TKym5JWk3lucxlZAAAA\nxCTIpbxLXEb26r6FAAAAgAhEM4QDAAAA6AMCNAAAAFABARoAAACogAANAAAAVECABgAAACogQAMA\nAAAVEKABAACACgjQAAAAQAUEaAAAAKACAjQAAABQAQEaAAAAqIAADQAAAFQQLECb2eUp09bMbNXM\n1rtsEwAAADBLkACdBOO1CdOWJcndN7P3AQAAgBgECdDuviFpa8Lki5K2k9tbklY7aRQAAABQQoxj\noBcl3cncPxWqIQAAAEBejAEaAAAAiNaB0A0osC3pZHJ7UdLt/AzJGOp1STq8cKy7lgHAAJjZmnb7\n2qVkSF2l6QAwdtFUoM1sMbl5RdJScntJ0mZ+XnffcPcVd185uHCkqyYCQO/NOlA7ub+VTN/iQG4A\n2C/UWTjWJK3kTlN3XZLc/WYyz6qk7fQ+AKARZQ7UTk8zukQfDAD7BRnC4e5XJV3NPXYhc5ufDAGg\nHVMP1Hb3m2a2ZWZ3Jb3dacsAoCdiHAMNAK179OrZuZY/+PFnDbUkLslwulvaDc8/M7Ob7r6Vm+fJ\ncSgvv/xy940EgMAI0ABGY97QPG1dPQrUsw7UXpe04e7bZrat3YtevZedIfmVcEOSVlZWvN3mAkB8\nCNAABq3J0Fx2O5GH6SuSVpLbTw7UNrNFd9+WpMz/m2a2VLgWABgxAjSAQeoqOM/admxhOhnjvFJw\noPZ1SRfc/T0ze1e7Bxie5JgUANiPAA1gcEKG57y0LTEF6aJQnDuQ+738dADAUwRoAIMRU3DOi7kq\nDQCohgANYBDmDc+/P3+o0vzP3/pz7W3ta+uXtVcFAAiAAA2g9+qE56qBedby8wRqAEC/EKAB9FaI\n4DxrvQRpABg+AjSAXqoantsKzpO2Q5AGgOFaCN0AAKiqSnj+/flDnYXn/HYBAMNEBRrAYFUJsff+\nwmpt4/i/Tb4QH9VoABgmAjSAXilTfS4bnOuG5qJ1zArShGgAGI4gQzjMbM3MVs1sfcL0y8n/hdMB\njFNT4fneX1gj4bntdQIA4tR5gDazZUly983s/Zx1M7ul3UvJAkAj4bmLkDtp/YyJBoDhCDGE46Kk\na8ntLUmrkm7m5nnb3a922ioA0WoqPM/y55cflW6TJB367cGJ25o2pAMA0G8hAvSipDuZ+6cK5jlp\nZquSlt39vW6aBaCvpoXnWcG5amietGw+TBeFaMZCA8AwRHkaO3ffSIZ4nEqC9B5mtm5mN8zsxqOd\nhwFaCKArs6rPdcPzn19+NFd4Llpfle0DAPorRIDelnQyub0o6XZ2YhKO15K7tyUt5VeQBOwVd185\nuHCk1cYCCKfN8FzW98/+Tt8/+7tS8zYZyAEA8QoxhOOKpJXk9pKk9GDCRXff1u646BvJ9FPpdADI\nqhOeywTcSWG56PHffPbdwm1kh3MwHhoAhqfzAO3uN81sJRmase3u6QGE1yVdcPfN5DR3knQ7Mx3A\niFS9VHeqTnguW2GetFw+SOdDNABgWIJcSMXdNwoeu5C5zRk4AEw0qfpcNTzXDc6T1pMN0tkQna1C\n5w8kPPjxZ420AQDQHa5ECCA606rPXYfnN1/8eM/9a1++OnHe75/93cQQDQAYDgI0gKjUGbrRZHjO\nB+ZZ0/OBOh+is21kLDQADEOUp7EDgCJVruZXFJ6nnVHjzRc/nhmeJy1XtJ2idqRBn6sSAkC/EaAB\nRKOpoRtVTydXJzjnl8+vY1KITqXPp+7BkgCAcAjQAHqrSnguqjzXrTpPUiZEZ9v8+/OHqEYDQA8x\nBhpAFOpUn8uaFJ5nWXtu+lk0r/5huXC92XHRRWOiGQ8NAP1GBRpAL5WtPtc5Vd3aczdnhud0viKT\nwnl+PHT6DwDQLwRoAFGbp/o87YDBImWDc5llstvID+Xgkt8A0G8EaADBVT2Qbp4DB6eF53lUCdES\nQRoA+owADSBaTVefmzxgsEjVEC1VP2MIACA8AjSAoLqsPk8yb/W5yrqmnYsaANAPBGgAUWpj7HOR\nJsPzJLMuttI1M1szs1UzW58wfTmZp3A6AIwdARpAMDGMfW5D2VAeohptZsuS5O6b2fs5P3H3q5IW\nJ0wHgFELEqBLVD+mTgeAvukywM9wUdJ2cntL0mp2opmtSfrPkuTu77l7+yV6AOiZzgP0rOpHyeoI\ngAHj6nytWpR0J3P/VG76DyWdSoZxvNtdswCgP0JUoKdWP0pMBzAAbQ3f4AC9RtxOK89JRXoPM1s3\nsxtmduN3v+P1BjA+IQL0rOrHrOkAELWiS3xHZFvSyeT2oqTbuem3tVu8SOf9YX4F7r7h7ivuvvLd\n7343PxkABu9A6AbUkYyNXpekwwvHAremulmVt4Mff9ZRS4D49H34xqTwfO3LVztuyURXJK0kt5ck\npcPlFt19W9JVSWnVeVHJeGgAwFMhKtCzqh+zpu+pfhxcONJaQ5vy6NWze/5VmR8Yopj27cirxY3L\nDM1YlbSdOUjwejJ9S9J2MnTjVHI2DgBARogK9KzqR+H0vmkqIDx69SwVaaAHpgXxiKrPknaLEAWP\nXSiYTngGgAKdB2h3v2lmKxOqHxemTI9eW1W1dL0EaaAdafitc1GVWRXsaeH5N58xfhgA+ijIGOgK\n1Y9e6OrnaKrRwGy/+ey7hWfiuPblqzPPxVw2SJcZ9kFwBoDh6uVBhLGIaRwn0Bd13jdFp7Br07zj\noieFZ4IzAAwDAbqGkMGZKjSGrO0zcJSpQs+7/iIEZwAYFgJ0SVSbgXCO/5tXqkJPGsYhPQ25TQXp\neYZqHPrtwUbaAADoFgF6hhiDM1VoYNeh3x4svBphGXWDdJkzakwLzoRmAOg/AvQUMYZnAOVMq0Jn\nNXmKOYIzAIwDAboAwRnoj2lV6LIheh4M0wCA8SkVoM3snKS/lXTV3X9jZn/j7v/UbtO69/+3dz+7\ncWTnGYffb+xICmx52mSU1XjGaC2SrUnODSStG4gp6AqmfQcSfAWBdAFZcK6AMJ11gKF8A5Y5XmQx\nWYj2GFoYmJBuYxJkNLDxZdGnxFKpuqtOd1XXn/N7AGO6WdXdp6iqrteH55yvieC8bhLU+y9fb/3+\nAOJkAbfJIL1taL77pTfWFgDA7tXtgf6ploVOHprZZ5IeSBpNgN40OMeuGJDfnzANvOv9l69XXlfr\nJhLWGQudD72xYbruKhrrgjOhGQDGo26A/lzSS3f/3Mz+WcsS24O3SXBuapmt7H02DdJMJATeFjOh\nsMll5ehtBoD01A3Ql5Jmkv7d3Z+b7baoQdNig3Oba9P++f5teqOBmqqWs9tmVY4YTYRmrnsAGK7K\nAG1mP3D330n6XfYzd3/eaqtaEhOc2y7oUPwsbqZAM7Jw23SQrjMZcF1w5hoHgPGo0wP9czM7dfff\nmtlPJLm7/7bthjWpr8F5GwzfwFitGwct1S+qkg+8m4TpuqtnVPU2E5wBYHzqBOgXkqZmdhnGQP9T\n241qypCCM73QQH2xlQmbXkqO0AwAaasToKeSFpKeheXsPpP0q20+1MyOw3tO3f2kZPtTd39iZvOy\n7VWGFJwBvKuqF1qKD9HbanJcM39BAoBhqxOgL939l5I+lSQz+5dtPtDMDiTJ3c/NbG5mB+5+Udht\nHkL2z+q85y5W04i5UTPrHljt1hevGitWlF1rbQXputdyneBMaAaA8agM0O7+SzP7cSig8hNJ97f8\nzEda9mJLN6t7FAP0J+5+VufN/E7cn2ZjgvOmN+X862LCNMM4gBt1eqEzTQXpmOuV0AwA6aq1jJ27\n/z7893Mt14TexkTSde75fsk+e2Y2k3Tg7s+KG81sLmkuSbdvT2p9aN0bcdM9WV9/ZPRIAxuKCdFS\neQAuu6a3uSYJzgCAuutA71Q27tnMHpjZzN3PS7afSNIP7n6w9k5Y5+bb9jhKQjTwtphhHLEhuqiJ\na4+xzQCAvFYCdOghLroOwzIWkvbCzyaSrkpem+17pS2qHlbddHc5AalJ3KSRmm1D9DafWwfXJACk\npZUAXbFyxqmko/B4Kulcksxs4u4LLcdFvwjb97PtMfrQ61z2efRCAzdiJxPuKkTHzEMgOANAmnY+\nhMPdL8zsKIxxXuRW4Hgu6TCsznEcyoVflazQsVIfgzOA5mThtskgvcnEXYIzAKStkzHQZT3U7n6Y\ne1xrBQ5J+usda3yCYEzVsqYLNFThxo0x2XRJu3zojQ3T26x0w/UHAJB6OomwaXWC8yalfrPXNRWi\nWcIOKdp2XehdXDcEZwBA3ugDdFV43jQ4F99j1z3RANpHcAYAlHmv6wa0aRfhuY33AlLTt6B664tX\nvWtTk8I8k9mKFZPy+z3eVZsAYEhGHaDXGWLgHfMNHejD+T324CxJZnYgSdn6+tnzkv1mkh7ssGkA\nMBijHcKxrvd5iOEZSEEWXrcZE73pZybkkaTPwuNLSTNJtVc7AgAk2APdZnje5r2rJkIleJNHwnZx\nvqfQ27zCRNJ17vl+cQczOyhWgAUA3BhtD3QZep6B4WijNzrRwLyJvepdACBdSQVoAMOTD72xYZrA\nXGqhm4A8kXSV31in9zlMPpxL0ocffthGGwGg1wjQA0AIAJa4FhpxKukoPJ5KyiYTTtx9IWlqZlMt\nQ/ZeCNRvjZEOxbBOJOno6Mh31nIA6InkxkD3EQVUAOxKFobDKhuLXDh+Hraf5arBTjpoIgD0Hj3Q\nBT/+4Ks3j3//6l7Ua9cVU7n7JZ00APoh9CAXf3ZYss87+wEA6IHuPf5kDQAA0C+dBWgze7pmW60q\nWZL0nW980EMghtx2AACAFHUSoEMwPl6xrVaVrKJiEO3TkIk+tQUAAADb6WQMtLufmNnDFZs7rZIV\nO+45s27886YYvgEAANA/fRwDXVkla5WhDYcYWnsBAAAw0FU48ov433nv+x23pp3VN+h9BgAA6KdW\nAvSKyX/XubVF11lbJUt6e3ml9//m7/3WF6/eVCh7/+Vr/fn+7Y3avUv0PgMAAAxTKwG6bI3RKrkq\nWKVVsqrkQ3Tm7peurz+yN89v/+GWXn/4bWzT1mLtZwAAgLR0tQrHsaSjQk91VgVrVZWs0ajqfWb4\nBgAAQH91tQrHmaSzws8Oc48HUf2K3mcAAID09HEVjo0Uh2+s0tRyc20sWyfR+wwAANB3owjQdcNz\nU6rC87reZyYPAgAADNvgA7Tfie8J3qb3mPAMAACQtsEH6KI2l7Bra9hGhuEbAAAA/Te6AF20qkf4\n9h9u1Q7Edfel9xkAAGD8BlmJcJVNep/zwTi/RnRsb/O2q27Q+wwAADAMowrQ29p0iEZVeKb3GQAA\nYDwGH6D/esc6Ld3dRHim9xkAAGA4Rj8GWmqvqAnFUgAAANKTRIBuQ53wTO8zAADA+CQToJvqLb77\npTcWngEAADA8gx8DHSMLvl9/ZBu/to664ZneZwAAgOHprAfazJ5WbTOzeRufHROG6/Y4Z+h5BgAA\nGLdOeqBDMD6W9GTFLnMzO5b0s7ba0MYEwJjwTO8zAADAMHXSA+3uJ5Iu1+zyibvfd/fzXbVpW4Rn\nAACANPR1EuGemc3M7HHXDamDYRsAAADp6GWAdveT0Pu8b2az4nYzm5vZCzN78Zf/+98OWngjNjzT\n+wwAADBsrYyBXjH579rdz2q+Ntv3StK0uE8YAnIiSd+79yMvC7FtVyfcpNeZ8AwAADB8rQToEHCj\nmNnE3Rdajo1+EX68L2ntOOjvfFM+GTAfcJsM05sO1yA8AwAAjENXq3AcSzoys3kubD+XdOju52Z2\nbGaSdOXuF9t+XjH0xgbqbcc4E54B9En4Dl5ImpZ1eOT+injf3VetlgQAyeokQIfhGWeFnx0Wtrdm\nl5P+CM8A+sTMDiQpdFbMzewg31ER5p2cu/ulmf3CzGZDWhEJAHahl5MIx4LwDKCHHmnZ+ywth8wV\nJ2pPcz+7VMk8FABIXVKlvHeJ8AygpyaSrnPP9/MbC0M6DiSd7qJRADAk9EC3gPAMYOjCUI+Lsnko\n+aVEv/rqqw5aBwDdIkA36NYXrwjPAPpuIWkvPJ5ouVxomdmqCYRhrf4jdz+6d+9eG20EgF4bVYDu\nKsASnAEMyKluxjVPFZYKNbNJtkNYIelZePxOMSsASN3gx0DbN9++E17Xhdlv//GDxj6b0AxgaNz9\nwsyOQjBe5IZoPJd0GH7+1MyeaNlT/bCrtgJAXw0+QMdaFXqrgjVhGcBYlK39nC0lGpas++HOGwUA\nA5JcgF6FgAwAAIA6RjUGGgAAAGgbARoAAACIQIAGAAAAIhCgAQAAgAgEaAAAACACARoAAACI0Mky\ndmY2Dw/vl5WKNbNjLcvNTsvWKwUAAAC6svMe6FDl6jwE42mxTKyZHUhvFvN/8xwAAADogy6GcEwl\nZaH5MjzPe6Rl73O2fSYAAACgJ3Y+hKMwJONA0mlhl4mk69zz/eJ7hCEgc0m68973m24iAAAAsFJn\nkwjD0IwLd7+Ifa27n7j7kbsf3Xrvb1toHQAAAFCulR7o3CTBvGt3P8s9n5VNINRy+MZeeDyRdNV0\n+wAAAIBNtRKgq1bOMLO5uz8Lj2fufm5mE3dfaDmk4yjsOpV03kYbAQAAgE10tQrHUzN7aWZ/ym16\nLknZkI6w32KTIR4AAABAW7qYRHgu6YclPz/MPWbtZwAAAPQSlQgBAACACARoAAAAIAIBGgAAAIhA\ngAYAAAAiEKABAACACARoAAAAIAIBGgAAAIhAgAYAAAAiEKABAACACARoAAAAIAIBGgAAAIjw3a4b\nAADYLTM7lrSQNHX3k9jtAJC6TnqgzWwe/vd0xfan2X67bRkAjJuZHUiSu5/nn9fdDgDoIECb2UzS\neejVmIbnRXMzeynpcretA4DRe6Rl77K0/I4tfgdXbQeA5HXRAz3VzRfyZXhe9Im73896QAAAjZlI\nus4934/cDgDJ2/kY6MJ4ugNJpyW77YWe6QN3f1bcGIZ2ZMM7/uc//vhv/7VFk/5O0n9v8fqhS/34\nJX4HqR+/1P3v4KMOPzta4Tv4tZn9Z5ft6UDX50sXOOZ0pHjc/xD7gs4mEYZxdRfuflHcloVsM3tg\nZrNiT3TY3sjEFjN74e5HTbzXEKV+/BK/g9SPX0rud7CQtBceTyRdRW5/6zs4sd+dJI45FSkes5Tm\ncZvZi9jXtBKgV0z+u3b3s9zzmbs/WfHabN8rlQ/xAABs5lRSdnOcSsomC07cfbFqOwDgRisBumrZ\nIzObZ0Mzsh7m3Jf3paTs/wnsiy9vAGiMu1+Y2VEYJrfI/RXwuaTDNdsBAMHOh3CEL+WnZvZEyz8T\nPgybsi/vczM7NjNJutrBl3fqa5ymfvwSv4PUj19K7HdQ1snh7ofrtq+R1O8u4JjTkOIxS2ked/Qx\nm7u30RAAAABglCjlDSRuVUEjAGkKfwWeVRUzM7PHu2oT0KR197265z8BGm9JKUzVvUjGLBz7cdft\n6EpVVVTcqLpexno91TjuUZ1DdStRhuGYD3bYtFbV+Hc+CPuM5vxO+Jpeed+LqcRKgA7G9iW4iZTC\nFOWKl8JY1yQrftasigqlW/67xnGP8RxKrhJlzfP352F1sMkYzu+a1/Rl2H45hmPOVNz3ap//BGiN\n9kswWmJhKrmbBN5RpyoqllIt/111XGM8hyorUZrZwcgqBa/9dzazY0m/liR3fzaSlWnqXLNZh+J0\nJMdcR+1KrATopTF+CWI9yhUnzt1PcqtNHOhm+Uy8K9Xy32uPK+FzaK96l0GpOn8/lrQfhnGMZdx3\n1bl9oWXP858K+yEgQCvpL0EgeeuqogJ1jOwcWluJcoS9z3W9WVY39EiPmplNJL2U9ImkT80slY7F\nykqsmc5KeffRyL4E31GzQmQqal8kGL3Sqqh4y9blvweq7nGN6RyqqlQ5DWFqT9JeCNRDv2dW/Ttf\n6WZ440LLHumh3zerjnku6cTdF2a20HJ+1LMdtm+nNqnEmkyA3qa8+FhEFkcYO8oV601PylGoDprc\n+VFWFbXrNvVUquW/q457dOdQjUqVZ9Kbe+qkq3Y2rOrf+Uw3E+wnCuOhB67y3M7993xMPdAr7nvR\nlVgppBLkf5Fj+BLcRDipPpX0JIUwFW4Al1pOkBj98eJt4QvyF1qO79uT9DDF676usuvFzH6TVTAc\n6/W07rg5h8aj5vl9LenjsXS01Tjmx2H73piu6aYQoMWNFAAAAPURoAEAAIAIrMIBAAAARCBAFTUP\nHgAAANJJREFUAwAAABEI0AAAAEAEAjQAAAAQgQANAAAAREimkApQV1gPe6rl+pcfS/rXbEF5AAAA\neqCBHDObhkpbWWA+zVUcm4Y1wwEAQMII0ECOu1+Gh4eSzgtlPI81ntK1AABgQwRoIMfMDsLDqbsv\ncs8lieqUAACAMdBAwczMppI+C8M1rrtuEAAA6BdKeQM1mdljSftiUiEAAEkjQAMAAAARGAMNAAAA\nRCBAAwAAABEI0AAAAEAEAjQAAAAQgQANAAAARCBAAwAAABEI0AAAAECE/wcyvQpnfH3ZcgAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9ee81aca50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot():\n",
    "    Xeval = OptDemo.design.FactorialDesign(101, domain).generate()\n",
    "    Yevala,_ = objective_model.predict_f(Xeval)\n",
    "    Yevalb,_ = constraint_model.predict_f(Xeval)\n",
    "    Yevalc = ei.evaluate(Xeval)\n",
    "    Yevald = pof.evaluate(Xeval)\n",
    "    Yevale = joint.evaluate(Xeval)\n",
    "    shape = (101, 101)\n",
    "    plots = [('Objective model', Yevala), ('Constraint model', Yevalb), \n",
    "             ('EI', Yevalc), ('PoF', Yevald), \n",
    "             ('EI * PoF', Yevale)]\n",
    "\n",
    "    f, axes = plt.subplots(3, 2, figsize=(10, 10))\n",
    "    for ax, plot in zip(axes.ravel(), plots):\n",
    "        ax.contourf(Xeval[:,0].reshape(shape), Xeval[:,1].reshape(shape), \n",
    "                    plot[1].reshape(shape))\n",
    "        ax.set_title(plot[0])\n",
    "        ax.set_xlabel('$x_1$')\n",
    "        ax.set_ylabel('$x_2$')\n",
    "    plt.tight_layout()\n",
    "    \n",
    "# Plot representing the model belief, and the belief mapped to EI and PoF\n",
    "plot()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the Bayesian optimization is the next step. For this, we must set up an appropriate strategy to optimize the joint acquisition function. Sometimes this can be a bit challenging as often large non-varying areas may occur. A typical strategy is to apply a Monte Carlo optimization step first, then optimize the point with the best value (several variations exist). This approach is followed here. We then run the Bayesian Optimization and allow it to select up to 50 additional decisions. \n",
    "\n",
    "The joint acquisition function assures the feasibility (w.r.t the constraint) is taken into account while selecting decisions for optimality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'list' object is not callable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-e12ce444d84d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;31m# Then run the BayesianOptimizer for 50 iterations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0moptimizer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mOptDemo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBayesianOptimizer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdomain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjoint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moptimizer\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0macquisition_opt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0moptimizer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtownsend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconstraint\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_iter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m50\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/bo.pyc\u001b[0m in \u001b[0;36moptimize\u001b[0;34m(self, objectivefx, **kwargs)\u001b[0m\n\u001b[1;32m     86\u001b[0m         \u001b[0;31m# Bayesian optimization is Gradient-free: incorporate.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     87\u001b[0m         \u001b[0mfx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mobjectivefx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mBayesianOptimizer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     90\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_optimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfxs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_iter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/optim.pyc\u001b[0m in \u001b[0;36moptimize\u001b[0;34m(self, objectivefx, **kwargs)\u001b[0m\n\u001b[1;32m     36\u001b[0m         \u001b[0mobjective\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mObjectiveWrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobjectivefx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_wrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m             \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_optimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobjective\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mKeyboardInterrupt\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             result = OptimizeResult(x=objective._previous_x,\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/bo.pyc\u001b[0m in \u001b[0;36m_optimize\u001b[0;34m(self, fxs, n_iter)\u001b[0m\n\u001b[1;32m    104\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_iter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    105\u001b[0m             \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimizer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_acquisition_wrapper\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 106\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_update_model_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_evaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfxs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    107\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    108\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_create_result\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"OK\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/bo.pyc\u001b[0m in \u001b[0;36m_evaluate\u001b[0;34m(self, X, fxs)\u001b[0m\n\u001b[1;32m     46\u001b[0m         \u001b[0mfxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0matleast_1d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfxs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 48\u001b[0;31m             \u001b[0mevaluations\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfxs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     49\u001b[0m             \u001b[0;32massert\u001b[0m \u001b[0mevaluations\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquisition\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mevaluations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/bo.pyc\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f)\u001b[0m\n\u001b[1;32m     46\u001b[0m         \u001b[0mfxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0matleast_1d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfxs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 48\u001b[0;31m             \u001b[0mevaluations\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfxs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     49\u001b[0m             \u001b[0;32massert\u001b[0m \u001b[0mevaluations\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquisition\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mevaluations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/optim.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0matleast_2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m         \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mObjectiveWrapper\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcounter\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_no_gradient\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflow-distrib/GPflow/model.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m     40\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     41\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 42\u001b[0;31m         \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_objective\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     43\u001b[0m         \u001b[0mg_is_fin\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg_is_fin\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/joachim/PycharmProjects/GPflowOpt/OptDemo/bo.pyc\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(X)\u001b[0m\n\u001b[1;32m     85\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0moptimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobjectivefx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     86\u001b[0m         \u001b[0;31m# Bayesian optimization is Gradient-free: incorporate.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 87\u001b[0;31m         \u001b[0mfx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mobjectivefx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     88\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mBayesianOptimizer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: 'list' object is not callable"
     ]
    }
   ],
   "source": [
    "# Run the Bayesian Optimization\n",
    "\n",
    "# First setup the optimization strategy for the acquisition function\n",
    "# Combining MC step followed by L-BFGS-B\n",
    "acquisition_opt = OptDemo.optim.StagedOptimizer([OptDemo.optim.MCOptimizer(domain, 200), \n",
    "                                                   OptDemo.optim.SciPyOptimizer(domain)])\n",
    "\n",
    "# Then run the BayesianOptimizer for 50 iterations\n",
    "optimizer = OptDemo.BayesianOptimizer(domain, joint, optimizer=acquisition_opt)\n",
    "optimizer.optimize([townsend, constraint], n_iter=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we now plot the belief, we clearly see the constraint model has improved significantly. More specifically, its PoF mapping is an accurate representation of the true constraint function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting belief again\n",
    "plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we inspect the sampling distribution, we can see that the amount of samples in the infeasible regions is limited. The optimization has focussed on the feasible areas. In addition, it has been active mostly in two optimal regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot function, overlayed by the constraint. Also plot the samples\n",
    "axes = plotfx()\n",
    "valid = joint.data[1][:,1] < 0\n",
    "axes.scatter(joint.data[0][valid,0], joint.data[0][valid,1])\n",
    "axes.scatter(joint.data[0][np.logical_not(valid),0], joint.data[0][np.logical_not(valid),1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}