tutorials/turbo_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BO with TuRBO-1 and TS/qEI\n",
    "\n",
    "In this tutorial, we show how to implement Trust Region Bayesian Optimization (TuRBO) [1] in a closed loop in BoTorch.\n",
    "\n",
    "This implementation uses one trust region (TuRBO-1) and supports either parallel expected improvement (qEI) or Thompson sampling (TS). We optimize the $20D$ Ackley function on the domain $[-5, 10]^{20}$ and show that TuRBO-1 outperforms qEI as well as Sobol.\n",
    "\n",
    "Since botorch assumes a maximization problem, we will attempt to maximize $-f(x)$ to achieve $\\max_x -f(x)=0$.\n",
    "\n",
    "[1]: [Eriksson, David, et al. Scalable global optimization via local Bayesian optimization. Advances in Neural Information Processing Systems. 2019](https://proceedings.neurips.cc/paper/2019/file/6c990b7aca7bc7058f5e98ea909e924b-Paper.pdf)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import math\n",
    "from dataclasses import dataclass\n",
    "\n",
    "import torch\n",
    "from botorch.acquisition import qExpectedImprovement\n",
    "from botorch.fit import fit_gpytorch_model\n",
    "from botorch.generation import MaxPosteriorSampling\n",
    "from botorch.models import SingleTaskGP\n",
    "from botorch.optim import optimize_acqf\n",
    "from botorch.test_functions import Ackley\n",
    "from botorch.utils.transforms import unnormalize\n",
    "from torch.quasirandom import SobolEngine\n",
    "\n",
    "import gpytorch\n",
    "from gpytorch.constraints import Interval\n",
    "from gpytorch.kernels import MaternKernel, ScaleKernel\n",
    "from gpytorch.likelihoods import GaussianLikelihood\n",
    "from gpytorch.mlls import ExactMarginalLogLikelihood\n",
    "from gpytorch.priors import HorseshoePrior\n",
    "\n",
    "\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "dtype = torch.double\n",
    "SMOKE_TEST = os.environ.get(\"SMOKE_TEST\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimize the 20-dimensional Ackley function\n",
    "\n",
    "The goal is to minimize the popular Ackley function:\n",
    "\n",
    "$f(x_1,\\ldots,x_d) = -20\\exp\\left(-0.2 \\sqrt{\\frac{1}{d} \\sum_{j=1}^d x_j^2} \\right) -\\exp \\left( \\frac{1}{d} \\sum_{j=1}^d \\cos(2 \\pi x_j) \\right) + 20 + e$\n",
    "\n",
    "over the domain  $[-5, 10]^{20}$.  The global optimal value of $0$ is attained at $x_1 = \\ldots = x_d = 0$.\n",
    "\n",
    "As mentioned above, since botorch assumes a maximization problem, we instead maximize $-f(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fun = Ackley(dim=20, negate=True).to(dtype=dtype, device=device)\n",
    "fun.bounds[0, :].fill_(-5)\n",
    "fun.bounds[1, :].fill_(10)\n",
    "dim = fun.dim\n",
    "lb, ub = fun.bounds\n",
    "\n",
    "batch_size = 4\n",
    "n_init = 2 * dim\n",
    "max_cholesky_size = float(\"inf\")  # Always use Cholesky\n",
    "\n",
    "\n",
    "def eval_objective(x):\n",
    "    \"\"\"This is a helper function we use to unnormalize and evalaute a point\"\"\"\n",
    "    return fun(unnormalize(x, fun.bounds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maintain the TuRBO state\n",
    "TuRBO needs to maintain a state, which includes the length of the trust region, success and failure counters, success and failure tolerance, etc. \n",
    "\n",
    "In this tutorial we store the state in a dataclass and update the state of TuRBO after each batch evaluation. \n",
    "\n",
    "**Note**: These settings assume that the domain has been scaled to $[0, 1]^d$ and that the same batch size is used for each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class TurboState:\n",
    "    dim: int\n",
    "    batch_size: int\n",
    "    length: float = 0.8\n",
    "    length_min: float = 0.5 ** 7\n",
    "    length_max: float = 1.6\n",
    "    failure_counter: int = 0\n",
    "    failure_tolerance: int = float(\"nan\")  # Note: Post-initialized\n",
    "    success_counter: int = 0\n",
    "    success_tolerance: int = 10  # Note: The original paper uses 3\n",
    "    best_value: float = -float(\"inf\")\n",
    "    restart_triggered: bool = False\n",
    "\n",
    "    def __post_init__(self):\n",
    "        self.failure_tolerance = math.ceil(\n",
    "            max([4.0 / self.batch_size, float(self.dim) / self.batch_size])\n",
    "        )\n",
    "\n",
    "\n",
    "def update_state(state, Y_next):\n",
    "    if max(Y_next) > state.best_value + 1e-3 * math.fabs(state.best_value):\n",
    "        state.success_counter += 1\n",
    "        state.failure_counter = 0\n",
    "    else:\n",
    "        state.success_counter = 0\n",
    "        state.failure_counter += 1\n",
    "\n",
    "    if state.success_counter == state.success_tolerance:  # Expand trust region\n",
    "        state.length = min(2.0 * state.length, state.length_max)\n",
    "        state.success_counter = 0\n",
    "    elif state.failure_counter == state.failure_tolerance:  # Shrink trust region\n",
    "        state.length /= 2.0\n",
    "        state.failure_counter = 0\n",
    "\n",
    "    state.best_value = max(state.best_value, max(Y_next).item())\n",
    "    if state.length < state.length_min:\n",
    "        state.restart_triggered = True\n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take a look at the state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TurboState(dim=20, batch_size=4, length=0.8, length_min=0.0078125, length_max=1.6, failure_counter=0, failure_tolerance=5, success_counter=0, success_tolerance=10, best_value=-inf, restart_triggered=False)\n"
     ]
    }
   ],
   "source": [
    "state = TurboState(dim=dim, batch_size=batch_size)\n",
    "print(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate initial points\n",
    "This generates an initial set of Sobol points that we use to start of the BO loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_initial_points(dim, n_pts, seed=0):\n",
    "    sobol = SobolEngine(dimension=dim, scramble=True, seed=seed)\n",
    "    X_init = sobol.draw(n=n_pts).to(dtype=dtype, device=device)\n",
    "    return X_init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate new batch\n",
    "Given the current `state` and a probabilistic (GP) `model` built from observations `X` and `Y`, we generate a new batch of points.  \n",
    "\n",
    "This method works on the domain $[0, 1]^d$, so make sure to not pass in observations from the true domain.  `unnormalize` is called before the true function is evaluated which will first map the points back to the original domain.\n",
    "\n",
    "We support either TS and qEI which can be specified via the `acqf` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_batch(\n",
    "    state,\n",
    "    model,  # GP model\n",
    "    X,  # Evaluated points on the domain [0, 1]^d\n",
    "    Y,  # Function values\n",
    "    batch_size,\n",
    "    n_candidates=None,  # Number of candidates for Thompson sampling\n",
    "    num_restarts=10,\n",
    "    raw_samples=512,\n",
    "    acqf=\"ts\",  # \"ei\" or \"ts\"\n",
    "):\n",
    "    assert acqf in (\"ts\", \"ei\")\n",
    "    assert X.min() >= 0.0 and X.max() <= 1.0 and torch.all(torch.isfinite(Y))\n",
    "    if n_candidates is None:\n",
    "        n_candidates = min(5000, max(2000, 200 * X.shape[-1]))\n",
    "\n",
    "    # Scale the TR to be proportional to the lengthscales\n",
    "    x_center = X[Y.argmax(), :].clone()\n",
    "    weights = model.covar_module.base_kernel.lengthscale.squeeze().detach()\n",
    "    weights = weights / weights.mean()\n",
    "    weights = weights / torch.prod(weights.pow(1.0 / len(weights)))\n",
    "    tr_lb = torch.clamp(x_center - weights * state.length / 2.0, 0.0, 1.0)\n",
    "    tr_ub = torch.clamp(x_center + weights * state.length / 2.0, 0.0, 1.0)\n",
    "\n",
    "    if acqf == \"ts\":\n",
    "        dim = X.shape[-1]\n",
    "        sobol = SobolEngine(dim, scramble=True)\n",
    "        pert = sobol.draw(n_candidates).to(dtype=dtype, device=device)\n",
    "        pert = tr_lb + (tr_ub - tr_lb) * pert\n",
    "\n",
    "        # Create a perturbation mask\n",
    "        prob_perturb = min(20.0 / dim, 1.0)\n",
    "        mask = (\n",
    "            torch.rand(n_candidates, dim, dtype=dtype, device=device)\n",
    "            <= prob_perturb\n",
    "        )\n",
    "        ind = torch.where(mask.sum(dim=1) == 0)[0]\n",
    "        mask[ind, torch.randint(0, dim - 1, size=(len(ind),), device=device)] = 1\n",
    "\n",
    "        # Create candidate points from the perturbations and the mask        \n",
    "        X_cand = x_center.expand(n_candidates, dim).clone()\n",
    "        X_cand[mask] = pert[mask]\n",
    "\n",
    "        # Sample on the candidate points\n",
    "        thompson_sampling = MaxPosteriorSampling(model=model, replacement=False)\n",
    "        with torch.no_grad():  # We don't need gradients when using TS\n",
    "            X_next = thompson_sampling(X_cand, num_samples=batch_size)\n",
    "\n",
    "    elif acqf == \"ei\":\n",
    "        ei = qExpectedImprovement(model, train_Y.max(), maximize=True)\n",
    "        X_next, acq_value = optimize_acqf(\n",
    "            ei,\n",
    "            bounds=torch.stack([tr_lb, tr_ub]),\n",
    "            q=batch_size,\n",
    "            num_restarts=num_restarts,\n",
    "            raw_samples=raw_samples,\n",
    "        )\n",
    "\n",
    "    return X_next"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization loop\n",
    "This simple loop runs one instance of TuRBO-1 with Thompson sampling until convergence.\n",
    "\n",
    "TuRBO-1 is a local optimizer that can be used for a fixed evaluation budget in a multi-start fashion.  Once TuRBO converges, `state[\"restart_triggered\"]` will be set to true and the run should be aborted.  If you want to run more evaluations with TuRBO, you simply generate a new set of initial points and then keep generating batches until convergence or when the evaluation budget has been exceeded.  It's important to note that evaluations from previous instances are discarded when TuRBO restarts.\n",
    "\n",
    "NOTE: We use a `SingleTaskGP` with a noise constraint to keep the noise from getting too large as the problem is noise-free. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44) Best value: -1.11e+01, TR length: 8.00e-01\n",
      "48) Best value: -1.08e+01, TR length: 8.00e-01\n",
      "52) Best value: -1.01e+01, TR length: 8.00e-01\n",
      "56) Best value: -1.01e+01, TR length: 8.00e-01\n",
      "60) Best value: -1.01e+01, TR length: 8.00e-01\n",
      "64) Best value: -9.05e+00, TR length: 8.00e-01\n",
      "68) Best value: -9.05e+00, TR length: 8.00e-01\n",
      "72) Best value: -9.05e+00, TR length: 8.00e-01\n",
      "76) Best value: -9.05e+00, TR length: 8.00e-01\n",
      "80) Best value: -9.05e+00, TR length: 8.00e-01\n",
      "84) Best value: -8.11e+00, TR length: 8.00e-01\n",
      "88) Best value: -7.41e+00, TR length: 8.00e-01\n",
      "92) Best value: -7.41e+00, TR length: 8.00e-01\n",
      "96) Best value: -7.41e+00, TR length: 8.00e-01\n",
      "100) Best value: -7.41e+00, TR length: 8.00e-01\n",
      "104) Best value: -7.41e+00, TR length: 8.00e-01\n",
      "108) Best value: -7.41e+00, TR length: 4.00e-01\n",
      "112) Best value: -6.57e+00, TR length: 4.00e-01\n",
      "116) Best value: -6.27e+00, TR length: 4.00e-01\n",
      "120) Best value: -6.24e+00, TR length: 4.00e-01\n",
      "124) Best value: -5.58e+00, TR length: 4.00e-01\n",
      "128) Best value: -5.58e+00, TR length: 4.00e-01\n",
      "132) Best value: -5.57e+00, TR length: 4.00e-01\n",
      "136) Best value: -5.57e+00, TR length: 4.00e-01\n",
      "140) Best value: -5.57e+00, TR length: 4.00e-01\n",
      "144) Best value: -5.57e+00, TR length: 4.00e-01\n",
      "148) Best value: -5.57e+00, TR length: 4.00e-01\n",
      "152) Best value: -5.33e+00, TR length: 4.00e-01\n",
      "156) Best value: -5.32e+00, TR length: 4.00e-01\n",
      "160) Best value: -5.32e+00, TR length: 4.00e-01\n",
      "164) Best value: -5.32e+00, TR length: 4.00e-01\n",
      "168) Best value: -5.32e+00, TR length: 4.00e-01\n",
      "172) Best value: -5.32e+00, TR length: 2.00e-01\n",
      "176) Best value: -4.98e+00, TR length: 2.00e-01\n",
      "180) Best value: -4.27e+00, TR length: 2.00e-01\n",
      "184) Best value: -4.04e+00, TR length: 2.00e-01\n",
      "188) Best value: -4.04e+00, TR length: 2.00e-01\n",
      "192) Best value: -4.04e+00, TR length: 2.00e-01\n",
      "196) Best value: -4.04e+00, TR length: 2.00e-01\n",
      "200) Best value: -4.04e+00, TR length: 2.00e-01\n",
      "204) Best value: -4.02e+00, TR length: 2.00e-01\n",
      "208) Best value: -4.02e+00, TR length: 2.00e-01\n",
      "212) Best value: -3.90e+00, TR length: 2.00e-01\n",
      "216) Best value: -3.90e+00, TR length: 2.00e-01\n",
      "220) Best value: -3.84e+00, TR length: 2.00e-01\n",
      "224) Best value: -3.84e+00, TR length: 2.00e-01\n",
      "228) Best value: -3.84e+00, TR length: 2.00e-01\n",
      "232) Best value: -3.84e+00, TR length: 2.00e-01\n",
      "236) Best value: -3.84e+00, TR length: 2.00e-01\n",
      "240) Best value: -3.84e+00, TR length: 1.00e-01\n",
      "244) Best value: -3.65e+00, TR length: 1.00e-01\n",
      "248) Best value: -3.35e+00, TR length: 1.00e-01\n",
      "252) Best value: -3.35e+00, TR length: 1.00e-01\n",
      "256) Best value: -3.03e+00, TR length: 1.00e-01\n",
      "260) Best value: -3.03e+00, TR length: 1.00e-01\n",
      "264) Best value: -3.03e+00, TR length: 1.00e-01\n",
      "268) Best value: -2.74e+00, TR length: 1.00e-01\n",
      "272) Best value: -2.74e+00, TR length: 1.00e-01\n",
      "276) Best value: -2.74e+00, TR length: 1.00e-01\n",
      "280) Best value: -2.74e+00, TR length: 1.00e-01\n",
      "284) Best value: -2.52e+00, TR length: 1.00e-01\n",
      "288) Best value: -2.52e+00, TR length: 1.00e-01\n",
      "292) Best value: -2.49e+00, TR length: 1.00e-01\n",
      "296) Best value: -2.49e+00, TR length: 1.00e-01\n",
      "300) Best value: -2.49e+00, TR length: 1.00e-01\n",
      "304) Best value: -2.49e+00, TR length: 1.00e-01\n",
      "308) Best value: -2.49e+00, TR length: 1.00e-01\n",
      "312) Best value: -2.49e+00, TR length: 5.00e-02\n",
      "316) Best value: -2.09e+00, TR length: 5.00e-02\n",
      "320) Best value: -2.09e+00, TR length: 5.00e-02\n",
      "324) Best value: -2.09e+00, TR length: 5.00e-02\n",
      "328) Best value: -1.83e+00, TR length: 5.00e-02\n",
      "332) Best value: -1.83e+00, TR length: 5.00e-02\n",
      "336) Best value: -1.83e+00, TR length: 5.00e-02\n",
      "340) Best value: -1.81e+00, TR length: 5.00e-02\n",
      "344) Best value: -1.81e+00, TR length: 5.00e-02\n",
      "348) Best value: -1.81e+00, TR length: 5.00e-02\n",
      "352) Best value: -1.81e+00, TR length: 5.00e-02\n",
      "356) Best value: -1.81e+00, TR length: 5.00e-02\n",
      "360) Best value: -1.81e+00, TR length: 2.50e-02\n",
      "364) Best value: -1.38e+00, TR length: 2.50e-02\n",
      "368) Best value: -1.38e+00, TR length: 2.50e-02\n",
      "372) Best value: -1.38e+00, TR length: 2.50e-02\n",
      "376) Best value: -1.38e+00, TR length: 2.50e-02\n",
      "380) Best value: -1.11e+00, TR length: 2.50e-02\n",
      "384) Best value: -1.11e+00, TR length: 2.50e-02\n",
      "388) Best value: -1.11e+00, TR length: 2.50e-02\n",
      "392) Best value: -1.11e+00, TR length: 2.50e-02\n",
      "396) Best value: -1.11e+00, TR length: 2.50e-02\n",
      "400) Best value: -1.11e+00, TR length: 1.25e-02\n",
      "404) Best value: -9.81e-01, TR length: 1.25e-02\n",
      "408) Best value: -8.63e-01, TR length: 1.25e-02\n",
      "412) Best value: -8.63e-01, TR length: 1.25e-02\n",
      "416) Best value: -8.63e-01, TR length: 1.25e-02\n",
      "420) Best value: -8.63e-01, TR length: 1.25e-02\n",
      "424) Best value: -8.63e-01, TR length: 1.25e-02\n",
      "428) Best value: -8.63e-01, TR length: 6.25e-03\n"
     ]
    }
   ],
   "source": [
    "X_turbo = get_initial_points(dim, n_init)\n",
    "Y_turbo = torch.tensor(\n",
    "    [eval_objective(x) for x in X_turbo], dtype=dtype, device=device\n",
    ").unsqueeze(-1)\n",
    "\n",
    "state = TurboState(dim, batch_size=batch_size)\n",
    "\n",
    "NUM_RESTARTS = 10 if not SMOKE_TEST else 2\n",
    "RAW_SAMPLES = 512 if not SMOKE_TEST else 4\n",
    "N_CANDIDATES = min(5000, max(2000, 200 * dim)) if not SMOKE_TEST else 4\n",
    "\n",
    "\n",
    "while not state.restart_triggered:  # Run until TuRBO converges\n",
    "    # Fit a GP model\n",
    "    train_Y = (Y_turbo - Y_turbo.mean()) / Y_turbo.std()\n",
    "    likelihood = GaussianLikelihood(noise_constraint=Interval(1e-8, 1e-3))\n",
    "    covar_module = ScaleKernel(  # Use the same lengthscale prior as in the TuRBO paper\n",
    "        MaternKernel(nu=2.5, ard_num_dims=dim, lengthscale_constraint=Interval(0.005, 4.0))\n",
    "    )\n",
    "    model = SingleTaskGP(X_turbo, train_Y, covar_module=covar_module, likelihood=likelihood)\n",
    "    mll = ExactMarginalLogLikelihood(model.likelihood, model)\n",
    "\n",
    "    # Do the fitting and acquisition function optimization inside the Cholesky context\n",
    "    with gpytorch.settings.max_cholesky_size(max_cholesky_size):\n",
    "        # Fit the model\n",
    "        fit_gpytorch_model(mll)\n",
    "    \n",
    "        # Create a batch\n",
    "        X_next = generate_batch(\n",
    "            state=state,\n",
    "            model=model,\n",
    "            X=X_turbo,\n",
    "            Y=train_Y,\n",
    "            batch_size=batch_size,\n",
    "            n_candidates=N_CANDIDATES,\n",
    "            num_restarts=NUM_RESTARTS,\n",
    "            raw_samples=RAW_SAMPLES,\n",
    "            acqf=\"ts\",\n",
    "        )\n",
    "\n",
    "    Y_next = torch.tensor(\n",
    "        [eval_objective(x) for x in X_next], dtype=dtype, device=device\n",
    "    ).unsqueeze(-1)\n",
    "\n",
    "    # Update state\n",
    "    state = update_state(state=state, Y_next=Y_next)\n",
    "\n",
    "    # Append data\n",
    "    X_turbo = torch.cat((X_turbo, X_next), dim=0)\n",
    "    Y_turbo = torch.cat((Y_turbo, Y_next), dim=0)\n",
    "\n",
    "    # Print current status\n",
    "    print(\n",
    "        f\"{len(X_turbo)}) Best value: {state.best_value:.2e}, TR length: {state.length:.2e}\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GP-EI\n",
    "As a baseline, we compare TuRBO to qEI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44) Best value: -1.12e+01\n",
      "48) Best value: -1.06e+01\n",
      "52) Best value: -9.38e+00\n",
      "56) Best value: -9.38e+00\n",
      "60) Best value: -8.60e+00\n",
      "64) Best value: -8.60e+00\n",
      "68) Best value: -8.60e+00\n",
      "72) Best value: -8.60e+00\n",
      "76) Best value: -8.60e+00\n",
      "80) Best value: -8.60e+00\n",
      "84) Best value: -8.60e+00\n",
      "88) Best value: -8.60e+00\n",
      "92) Best value: -8.60e+00\n",
      "96) Best value: -8.60e+00\n",
      "100) Best value: -8.60e+00\n",
      "104) Best value: -8.60e+00\n",
      "108) Best value: -8.60e+00\n",
      "112) Best value: -8.60e+00\n",
      "116) Best value: -8.60e+00\n",
      "120) Best value: -8.60e+00\n",
      "124) Best value: -8.60e+00\n",
      "128) Best value: -8.60e+00\n",
      "132) Best value: -8.60e+00\n",
      "136) Best value: -8.60e+00\n",
      "140) Best value: -8.60e+00\n",
      "144) Best value: -8.60e+00\n",
      "148) Best value: -8.60e+00\n",
      "152) Best value: -8.60e+00\n",
      "156) Best value: -8.60e+00\n",
      "160) Best value: -8.60e+00\n",
      "164) Best value: -8.60e+00\n",
      "168) Best value: -8.60e+00\n",
      "172) Best value: -8.60e+00\n",
      "176) Best value: -8.60e+00\n",
      "180) Best value: -8.60e+00\n",
      "184) Best value: -8.60e+00\n",
      "188) Best value: -8.60e+00\n",
      "192) Best value: -8.60e+00\n",
      "196) Best value: -8.60e+00\n",
      "200) Best value: -8.60e+00\n",
      "204) Best value: -8.60e+00\n",
      "208) Best value: -8.60e+00\n",
      "212) Best value: -8.60e+00\n",
      "216) Best value: -8.60e+00\n",
      "220) Best value: -8.60e+00\n",
      "224) Best value: -8.60e+00\n",
      "228) Best value: -8.60e+00\n",
      "232) Best value: -8.60e+00\n",
      "236) Best value: -8.60e+00\n",
      "240) Best value: -8.60e+00\n",
      "244) Best value: -8.60e+00\n",
      "248) Best value: -8.60e+00\n",
      "252) Best value: -8.60e+00\n",
      "256) Best value: -8.60e+00\n",
      "260) Best value: -8.60e+00\n",
      "264) Best value: -8.60e+00\n",
      "268) Best value: -8.60e+00\n",
      "272) Best value: -8.60e+00\n",
      "276) Best value: -8.60e+00\n",
      "280) Best value: -8.60e+00\n",
      "284) Best value: -8.60e+00\n",
      "288) Best value: -8.60e+00\n",
      "292) Best value: -8.60e+00\n",
      "296) Best value: -8.60e+00\n",
      "300) Best value: -8.60e+00\n",
      "304) Best value: -8.60e+00\n",
      "308) Best value: -8.60e+00\n",
      "312) Best value: -8.60e+00\n",
      "316) Best value: -8.60e+00\n",
      "320) Best value: -8.60e+00\n",
      "324) Best value: -8.60e+00\n",
      "328) Best value: -8.60e+00\n",
      "332) Best value: -8.60e+00\n",
      "336) Best value: -8.60e+00\n",
      "340) Best value: -8.60e+00\n",
      "344) Best value: -8.60e+00\n",
      "348) Best value: -8.60e+00\n",
      "352) Best value: -8.60e+00\n",
      "356) Best value: -8.60e+00\n",
      "360) Best value: -8.60e+00\n",
      "364) Best value: -8.60e+00\n",
      "368) Best value: -8.60e+00\n",
      "372) Best value: -8.60e+00\n",
      "376) Best value: -8.60e+00\n",
      "380) Best value: -8.60e+00\n",
      "384) Best value: -8.60e+00\n",
      "388) Best value: -8.60e+00\n",
      "392) Best value: -8.60e+00\n",
      "396) Best value: -8.60e+00\n",
      "400) Best value: -8.60e+00\n",
      "404) Best value: -8.60e+00\n",
      "408) Best value: -8.60e+00\n",
      "412) Best value: -8.60e+00\n",
      "416) Best value: -8.60e+00\n",
      "420) Best value: -8.60e+00\n",
      "424) Best value: -8.60e+00\n",
      "428) Best value: -8.60e+00\n"
     ]
    }
   ],
   "source": [
    "X_ei = get_initial_points(dim, n_init)\n",
    "Y_ei = torch.tensor(\n",
    "    [eval_objective(x) for x in X_ei], dtype=dtype, device=device\n",
    ").unsqueeze(-1)\n",
    "\n",
    "while len(Y_ei) < len(Y_turbo):\n",
    "    train_Y = (Y_ei - Y_ei.mean()) / Y_ei.std()\n",
    "    likelihood = GaussianLikelihood(noise_constraint=Interval(1e-8, 1e-3))\n",
    "    model = SingleTaskGP(X_ei, train_Y, likelihood=likelihood)\n",
    "    mll = ExactMarginalLogLikelihood(model.likelihood, model)\n",
    "    fit_gpytorch_model(mll)\n",
    "\n",
    "    # Create a batch\n",
    "    ei = qExpectedImprovement(model, train_Y.max(), maximize=True)\n",
    "    candidate, acq_value = optimize_acqf(\n",
    "        ei,\n",
    "        bounds=torch.stack(\n",
    "            [\n",
    "                torch.zeros(dim, dtype=dtype, device=device),\n",
    "                torch.ones(dim, dtype=dtype, device=device),\n",
    "            ]\n",
    "        ),\n",
    "        q=batch_size,\n",
    "        num_restarts=NUM_RESTARTS,\n",
    "        raw_samples=RAW_SAMPLES,\n",
    "    )\n",
    "    Y_next = torch.tensor(\n",
    "        [eval_objective(x) for x in candidate], dtype=dtype, device=device\n",
    "    ).unsqueeze(-1)\n",
    "\n",
    "    # Append data\n",
    "    X_ei = torch.cat((X_ei, candidate), axis=0)\n",
    "    Y_ei = torch.cat((Y_ei, Y_next), axis=0)\n",
    "\n",
    "    # Print current status\n",
    "    print(f\"{len(X_ei)}) Best value: {Y_ei.max().item():.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sobol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_Sobol = SobolEngine(dim, scramble=True, seed=0).draw(len(X_turbo)).to(dtype=dtype, device=device)\n",
    "Y_Sobol = torch.tensor([eval_objective(x) for x in X_Sobol], dtype=dtype, device=device).unsqueeze(-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAHOCAYAAACVVnNvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABDX0lEQVR4nO3df3Qc5WHv/8/o9w9L3l3LwsYGBzkQAyGJ7aEG53dihwvpNbe3wkluA/qGEpum4XzjnBQnTs/XNT0F7ORC+4XcLzb5UUPbpGBuT3wbuFybkpymJpSxSZo4KAbLGGwwsrxaS9Zvreb7x86uRqvd1Ur7e/b9OsfHs7Mz8zzP7kr70fM8M2PYti0AAAAvqSh0BQAAALKNgAMAADyHgAMAADyHgAMAADyHgAMAADyHgAMgKcMw7jYM44BhGLbz70nDMHY7/550njtgGMa6LJfbZhjGk1k81mHDMHoNw9iZ5vYHnO13Z6MOAPLP4DRxwDsMw7hb0gJJq5xVu23b3pdk202SfJJCkpZLeinFtr2SLNu21yd4rk3Sbkkh27ZvybQNzjF3Srpbkt+27VCWjnlYkTZsTnP7A8pimwDkV1WhKwAgO5xQcF80EBiGsUrSc4ZhfCb+S9rZ9pxt27vc6wzDaHOvcwkqEoSmsW27S9J6p9fjySwFgmhA2yhpTxaOJ0ldc9g+kKWyAeQZQ1SABzi9MbvdvR22bR+RdIukdsMw2l3btknaFB9kbNveKukbhmH45liNaFkZDVc5wWynIgEjrd4WAIhHwAG8YbnTkzKFbdsHFel5cQeFzZIOJjlOl6RNc6mAE672KRJOMrHOqfduSasyCFwAyhgBB/CGu1NMyrUkma7H65R8uKZL0rR5NrPwj4qEkmxMOo7OB5pT4AJQ3gg4gDccVGSeTCKBuOdWSTqXZNsuTQ1Ds3XE+X9OIckJRgel2NyeI0pzmMowjE3OWV93O/OJZgxZrm3vdob50q3nTqe8Tc4ZZatcz91tGMZx56yz4+7nnOcPO88djn8OQPYwyRjwgERnN7ms0uwm6voyqEo0SLXNcf9VcXODdkva7Ux+TjpJ2Om96nLmEUXXrTIMoz3ZmWEumyTdp+TDdvFlHZb0RWeOk5whtMOGYdxi2/YR27Z3GYaxR1KvpK3R7Vw+KenJGd4zABmiBwfwMFevxE7nsc95HMpFea5Jzr4Um6USinv8hPN/0l4cZwL1One4cZip9nMskLTatu1dqQKUq6ydigSpWGhx2rxb0qNx6/ZJ+kyCw5iS4usKIMsIOIBHOWFmp6Rb0vnyzmKZc923XZOBRlIsKByU1J5oH8fO+P0cq2co725FTqufzWtzt6QDCdYf1OSp7VG7FTmrzBe3flWCXh0AWcYQFeBdTyoyRBIborFtO2QYhpTZMFQq0evGzCVQfUaR6+kkOmabYRjJgkGbpOPxK1Nd0M8ZZto6m4sIuubLLHefdu+qo9xDabZtHzQMI6Tp1/JJu0wAc0fAATzIGUp50rbt2V4kz6fMvoCjc28S9XIk5fRyHEhWX8MwbEWGmzbHrY+Wl26gCjivTfQ09OWzqGa0rH9MErQS1X2PInXeI8WGDBP1NgHIMoaoAI9xhl6Opwg3B5X8iz2gyGnlc3WLIrc3mGlib7yNSv3Fv8/ZZgrX8FK6k5p9tm1vjc7XSefeVC5zGVaKXssnWj9ftm49ASA1Ag7gIc7QSVd8uIk7ZfqIkgeCNkWGtuZqo6QvzmG/5TN88f+jJF+SU7+7lH5PjDuk3KLI9YPSOlXbFaYSnn6eaP6R61T3rU7IYe4NkCcEHMAjol/+SXpP3F/iu5XgWjfOF/QqzXEIxTlV25pt743zxT9tDo2b65iJ5tXsVILeHefYd6c45hFFho5mE+i2KvGZUVLyCxLudurX7lyhGUAeEHAAD3B6ITY7y+2uf5ucL/lYD4fTq7AnwfDMTiWfeBtQkonJhmG0OXfe9s3x2i47lfwihW4Jz6ZyeqsOGoaxO75emj43xxe372ZFJjCnNVTlXKMnGB+cnNc/4Twgp34+RU5JB5Anhm3bha4DgAwZhtGr1GdGbY2/uaYz4dWnyKRinyJDW/vitrlbkasSR4dl3ENf7tCzew49N+sU6d1oc+rwxWTHcHqH1jnlHVGkpyh+wvEmRYLccTl3P3fOZGpTJERF2/CEnCAXvUmps/6IIsNrXYpc0ya6/UFn+1iAcV6XBfFlpWjrbkk783W6PgACDgDknGEYd8cHTAC5xRAVAOSQM7eJnhsgzwg4AJBFzj2w3GdabZzDafMAMsSF/gAgu76hyFyhg07vTToTqAFkGXNwACCLnEnNsR6cOVxNGkAWEHAAAIDnMAcHAAB4TkHn4Hz961+n+wgAAGTk/vvvN+LXFXyS8Y4dOwpdhbzr7u5Wa2troauRd7S7vNDu8kK7y0sxtXv79u0J1zNEBQAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPIeAAwAAPKcqFwc1TbNdUkhSm2VZe3JRBgAAQDJZ78Fxwo0syzroPF6X7TIAAABSycUQ1bWSupzlLkmrclAGAABAcrZtZ/Xf6tWrd69evXqVs7xu9erVO5Ntu3XrVltSWv9uv/12e3h4eMq/22+/Pe39v/nNb07b/6abbkp7/4cffnja/itXrkx7/6eeeiq23xtvvGEPDw/bixcvTnv/Q4cOTSs/3X0l2V1dXVP27erqmtX+8WUfOnQo7X0XL148pd3Dw8P2U089lfb+K1eunFb+ww8/nPb+N91007T9v/nNb+bts/eVr3ylaD570X/l9tlz/yunz14x/d7js8dnL1efva1bt9qJMkYu5uCEJAWyfdChoSF1d3dPW5eugYGBafuPjIykvX9/f/+0/cfGxtLePxQKxfYPhUKSpHA4nPb+wWBwWvmz0dPTo6qqqimPZyO+7GAwmPa+4XBY3d3dsXZLmrI8k7GxsWnl9/f3p73/yMjItP0HBgbS3j/Tz97w8HDRfPaiyu2z51ZOn71i+r0XxWcvPXz2MvvsSbmZZPySJJ+z3CbpQDYOWl9fr9bW1mnr0tXY2Dht/9ra2rT3b2pqmrZ/dXV12vv7fL4p+7e2tqqysjLt/QOBwLTyZ6OlpWXK/uPj47PaP77sQCD9DFtZWRnbP/q/z+dLe//q6upp5Tc1NaW9f21t7bT9Gxsb094/089eXV1dUX32JJXlZy+qnD57xfZ7T+Kzly4+e5l99iTJsG17VjukwzTNuyUdkbTKsqxdybb7+te/bu/YsSPr5Re77u7ujH5oSxXtLi+0u7zQ7vJSTO3evn277r//fiN+fU5OE3eFmoO5OD4AAEAqXOgPAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4DgEHAAB4Tk7uJg4AAMrHWHhCDxx8Tb95qz+v5bY21WpekucIOAAAICN/+8Ib+v6hN/Je7qWBeq1N8hwBBwAApO3CyLhOhUY0VDkoSRqfsPW3L+Q/3MyEgAMAANLyr6+d05/8wy81FrYTPr+ouVb3/8HVMoz81KeuqlL/9N3/k/A5Ag4AAJiRbdva9eyxpOFGkm5fu0zXtwXyWCvpn5KsJ+AAAFCiegdH9UJXr8YnJnJe1uneYR3rHpAkVVUYuthXN+V5c5lPn7t2ac7rkS4CDgAAJWhkLKz2PS/pVO9Q3sv+L+8N6K/+cGXey50NroMDAEAJ+vVbfQUJN1UVhtrf15L3cmeLHhwAAErQb073xZYva2nQ1Yubc15mdaWhG997kZbMz/2QWKYIOAAAlKBfuy6qd+uaS/RHv3dJ3sru7u7OW1lzxRAVAAAl6DdvTfbgXHNx7ntvSg09OAAA5NHEhK2x8IROnBvU+aGxOR1jNGzr9XORC+1VVRh6z0XJblhQvgg4AADkyff+7aQeev64hsayN4fliovmqba6MmvH8woCDgAAeTAentDf/MtxjYxnd4KuucyX1eN5BQEHAIA8ePv88JRws8RXp4vn12V0W4NL/A36k49cloXaeQ8BBwCAPHjTdc2aVZfO1w//+NoC1sb7OIsKAIA8eCM4GXAu9TcUsCblgYADAEAevOHqwbkkUF/AmpQHAg4AAHnwZnAwtnwpASfnCDgAAOSBuweHIarcI+AAAJBjtm1PnYNDD07OcRYVAKAkjY5P6J9/fWZKcEhXd/+IfnEiqAsj4xnVYWLCVkXFb2fczralwdGwJKmxtlL+huqMysXMCDgAgJL0+Itvatf/ebXQ1Zi1dy1okJHJxW+QlqwHHNM0fZLanH/XWpa1NdtlAADwL787W+gqzFpDTaXu+OC7Cl2NspCLHpyNkmRZ1h7TNK81TXOTZVl7clAOAKBMjYcndNR1N+0/+chlqq5Mv1ekqsLQBy6Zr8tb56kig96Unp6zamlZmPb2DTWVqqli+ms+ZD3gxIWZNkkHsl0GAKC8vdo9ELth5eL5tfrKJ5cXpB6jdVXyMZ+mKOUsRpqm2SYpaFnWwVyVAQAoT786dT62/L4l8wtYExSrOfXgmKbZLikQt7orLsy0W5a1eaZjdXd3z6UKJS0UChW6CgVBu8sL7S4vc233WHhCj774jk6cG57VfidDI7HltvkVBfsu4f0uXnMKOJZl7Uv1vGma7ZZl7XKWV1mWdSTZtq2trXOpQsmj3eWFdpcX2p2+7/3bSf3w5cwmC39wxRK1tvozOkYmeL+LUy7Oolonaadpmt9wVnEWFQAgoZ/85kxG+1/W0qD3L2WICtPlYpLxQUmFme0FACgZb4WGdfStfklSdaWhv954jaor0p8aWlVpaNUlPs5KQkJc6A8AUBAHXpmcN7PmsoDWrSjuIQ+UFmIvAKAgnj/WE1v+1JXpX0sGSAcBBwCQd6PjE3r5zVDs8UcubylcZeBJBBwAQN795q0+DTsX6lvqr9fi+XUFrhG8hoADAMg762Qotvx7y3wFqwe8i0nGAFAmbNtW55kLOtOX+KJ6E7b0ytv9eulkr8Yn7Fkde2x0TNU1b6S9/Ymewdjyte8q3DVs4F0EHAAoE3//76f0l0//rtDVmOZaenCQAwxRAUCZeOLw6UJXYZoPvXuBlvrrC10NeBA9OABQBoIDo/rdOxckSZUVhj787gUJt5tXW6kPX96ii+fXyZjF8XtDvfL7ZjfUVFtdoasXN8swZlMSkB4CDgCUgZdO9saW37ekWbv/6ANZPX5391hB7wcFxCPgAEARGhkL66ev9qinfzQrx/uX303e0HINk3pRBgg4AFCE/vvB17T3F2/m5NhrLiPgwPuYZAwARejAK2dn3mgOAo3VWnmJLyfHBooJPTgAUGRCg2N663zkWjXVlYZuWbUkK8etrarQf37fItXXVGbleEAxI+AAQJF55Ux/bPmKi+Zp+++vKGBtgNLEEBUAFJlX3p4MOFcuaipgTYDSRcABgCLj7sG5ioADzAlDVACQZT/59Rn9z5ff0tDYhHoujGh4fGJW+wcHJk8Nv3IxAQeYCwIOAGRRcHBMW/+pU2Ph2d2sMhHDkN5z0bws1AooPwxRAUAWdXYPZSXcSFLHdZeqsZa/Q4G54CcHALLotZ7h2PLlrY3661uu0bw5hJS66kr5GqqzWTWgrBBwACCLjp8bii3ftuZSvbuVISagEBiiAoAscvfgrFhEuAEKhYADAFkyNBrWqfMjkqQKQ7qc3hugYBiiAoA0PfT8cX3v304mPe3bds0tXraggVsiAAVEwAGANPQOjuo7PzsxJcSkwgX6gMJiiAoA0tB55kLa4WaJr053fGhZbisEICV6cAAgDZ2u2yf815WL9Vcbrkq4XffZbl3U2irDMPJVNQAJEHAAIA2dZy7Elq9a1KSKisQBpsIwCDdAEWCICgDS0PnOZA/OCubXAEWPHhwASGAsPKH/ON2n4bGwbFs6fnYg9hz3hwKKHwEHAOLYtq3Nf/9L/dvx4LTnlvjq1FzPLRSAYscQFQDE6eoZTBhuJGn1pb78VgbAnNCDAwBx/uV3Z2PLi5pr1dbSGFmeX6u7Pra8UNUCMAsEHAAlY3gsrNOhYdnpXpBmjg680h1b/r8/sVz/deXFOS0PQPYRcACUhNfPDarjbw/rTN9I3so0DOmjl7fkrTwA2UPAAVAQr58b1F/+pFPdF0bT2v6dvmGdHxrPca2mWnWJTwvm1eS1TADZkdOAY5rmTsuytuayDACl6c9//Fu9dDI06/1qqip0ib8++xWKs6CxRtv+0xU5LwdAbuQs4JimuU5SW66OD6B0hAbHNDIejj1+tXtgTuGmvrpC/+9n3qePMGwEYAY5CTimabZJ6srFsQFk7qXXe/Xsb7s1PpG7ybpDQ0OqrTunzjP9+tWpvqTb3XBVq7700cvSOuai5jr5GrgGDYCZ5aoHp82yrIOmac64YXd394zbeE0oFCp0FQqCdheH4OCY7ni8U8PjuT0TKV2fvWa+AhVDaW07emFI3Rdm3q6Qiu39zhfaXV5Kod1zCjimabZLCsSt7nJCzTrLsg6me6zW1ta5VKHk0e7yUkztfvE/zuQ93BiG1DKvRoYmb0JZXWnos+ZSrb1qWV7rkg/F9H7nE+0uL8Xe7jkFHMuy9qV4OujMv/FJajNNc5VlWUfmUg6A7HvpZG9s+VNXLtR1bfF/q2THhf5+zWtqUn11pT787gVa2FSbk3IAIJGsD1FFw4xpmpsUCTkAiojlmtz7+TWXaM1luQk43d3dRf8XHgDvytlZVJZl7ZG0J1fHBzCzw2+E9O8nJntswhN27K7Y1ZWG3r90fqGqBgA5xYX+AI96tfuC/tv3rKTPX7OkWXXVlXmsEQDkD3cTBzzKmuE6M5+6kuEjAN5FDw7gUb2Dk7dAWH2pT9cu88UeX7qgQRvet6gAtQKA/CDgAB7VOzgWW15/5UJ9Ya33TscGgGQYogI8KjgwGXD8jdwwEkB5IeAAHhVyDVH5ub0BgDJDwAE8yj1ERcABUG4IOIBHBaf04DBEBaC8EHAAj3L34ATowQFQZgg4gAcNjYY1PDYhSaqpqlBDDRf0A1BeCDiABwXjJhgbhpFiawDwHgIO4EFMMAZQ7gg4gAdNnX/DBGMA5YeAA3hQL9fAAVDmCDiAB025ijEBB0AZIuAAHtTLNXAAlDkCDuBBJ3oGY8v+RnpwAJQfAg7gMb2Do/qX352NPV59qa9wlQGAAqkqdAUAZMePXjqlF04E9U7fiMbCtiTpmiXNWrGoqcA1A4D8I+AAHvCLE0Ft/+fOaetvWXVxAWoDAIXHEBXgAQd+2z1t3cXz6/TpaxYVoDYAUHj04AAlzrZt/ezVc7HHWz65XMsXNmr1pT7Nq+VHHEB54rcfkGO9Q+N65dUehSfs2D87i8fvGxrTm71DkqSGmkrdvnaZaqronAVQ3gg4QA6d6h3Sxsde0fB4NiNNctdd5ifcAICYgwPk1NO/eSdv4UaSPnVVa97KAoBiRg8OkEOH3wjFlq9Z0qzWplpVGpJhGFkv670XN+vm9y3O+nEBoBQRcIAcmZiwdcQVcL79h+/VuxY0FK5CAFBGCDgoaxMTuRs+OtZ9QX3D45Kklnk1Whaoz1lZAICpCDjwhDd7h/Tdn7+ungujM28sKTxh6zdv9elsmttnavWlvpwMSwEAEiPgwBPu/9/HdLDz7MwbFgj3gwKA/CLgwBN+986FOe+b646V9y5q0B9yywQAyCsCDjzh3MDkUNO3/vBq1VdXzrjP4vl1unJRkyorcptwuru7uaIwAOQZv3VR8oZGwxocDUuSqioM/edrFjHfBQDKHBf6Q8kLDk723gQaawg3AAACDkpfcGAstrygsbqANQEAFAsCDkpe0DX/ZkFjTQFrAgAoFjmZg2Oa5ipJbZJkWda+XJQBRLknGAcIOAAA5a4H5xtOsAmYptmWozIASVMDDj04AAApBz04pmlukvSSaZptlmXtyfbxgXhBenAAAHFyMUS13Pk/aJrmbklbLcsKJdu4u7s7B1UobqFQqNBVKIhctfutc/2x5erwcNF9pni/ywvtLi+0u3jNKeCYptkuKRC3usuyrIPO8nHLskKmaR6WtEnSrmTHam1tnUsVSh7tzp6B8OnY8rsWL1Br68Ksl5Ep3u/yQrvLC+0uTnMKODNMHH5Jk+HHJyk0lzKAdHEWFQAgXtYnGTvhx2ea5jrnMfNwkFOcRQUAiJeT08Qty4oOSR1MuSEwSz0XRnTwlbMaGovcmmFoLKwzfSOx5wk4AACJe1GhhNi2rU1/90sdfbs/4fO+hmo11Mx8k00AgPdxJWOUjDN9I0nDTX11hb5+w+V5rhEAoFjRg4OS8cqZyXCzxFen9VdGZvD7G6r1Bx9YrIua6wpVNQBAkSHgoGT81tV784n3LNQ3/tMVBawNAKCYMUSFktHp6sG5anFTAWsCACh2BByUDHcPzpWLCDgAgOQYokJR+u8HXtP+/3hb4Qk7tu7shcj1bqorDS1f2FioqgEASgABB0XnRM+A9vz89aTPv3tho2qq6HwEACTHtwSKjvvCffHqqiu06cPvyl9lAAAliR4cFJ3Q0Fhs+SOXL9Bf3XxV7HFTbZXquZgfAGAGBBwUnT5XwGltqlVrU20BawMAKEUMUaHonB8ajy0311UXsCYAgFJFwEHROe/qwZlfTycjAGD2CDgoOn3Dkz048+vpwQEAzB4BB0WHHhwAQKYIOCg6UwMOPTgAgNkj4KDouIeomGQMAJgLAg6KDkNUAIBMEXBQdBiiAgBkioCDojIentCFkbAkyTAiVy4GAGC2CDgoKu75N021VaqoMApYGwBAqSLgoKhwDRwAQDYQcFBU3PehamaCMQBgjgg4KCruO4n76MEBAMwRfyJjRuEJW4e6gno7NJzRcfr6+9T8xljKbf7Xr9+OLXMNHADAXBFwMKNHf/66HnzueJaOdjrtLRmiAgDMFUNUmNHPXztXkHJ/713+gpQLACh9/ImMGQUHJ4eVbriqVc11c/vYDA0Pq76ubsbtFsyr0Zp3+XV9W2BO5QAAQMDBjHoHR2PLf37Te9TaVDun43R3d6u1tTVb1QIAICmGqJDSxISt0CBnNgEASgsBBymdHx7ThB1ZnldbqZoqPjIAgOLHtxVS6nX13gQaawpYEwAA0kfAQUrugONvYHgKAFAaCDhIqXdgcoJxoIEeHABAaSDgICV6cAAApYiAg5SCrh4cP3NwAAAlIifXwTFNs11SSFKbZVl7clEG8oMeHABAKcp6D45pmuskdVmWdVBSl2maq7JdBvKHs6gAAKUoF0NUlqQnnWDTZlnWkRyUgTwJuq5iTA8OAKBUZH2IyrKskGmauyU9KengTNt3d3dnuwpFLxQKFboKaTt7fjC2bIxcUHe3PedjlVK7s4l2lxfaXV5od/GaU8Bx5tjE3wmxy7Ksg84Q1UHLsnaZprnTNM12y7L2JTtWud6bqBjb/fPXzunhn3apb3g8tu6N4HBsuW3JRWpd0JBRGcXY7nyg3eWFdpcX2l2c5hRwUgUWSassy9rlLN8naeNcykD+/eXTv9Pr5waTPs8cHABAqcjFWVR7TNPcJKlLnEVVMsbDEzoZTB5ubn7/IjXVcfN5AEBpyMkcHEmEmhITHByT7Uyv8dVX6+9vXx17rq66Ukv99QWqGQAAs8ef5JAkne0fiS23Ntfq3a3zClgbAAAyw5WMIUnquTB5OvjCecy1AQCUNgIOJEndrh6chfNqC1gTAAAyR8CBpKk9OC304AAAShwBB5KksxdcPThN9OAAAEobAQeSmIMDAPAWzqIqViP9Mi7k7zYWVb1duszolyRdMuGXcW4g+2WcPyejsj/rxy12tLu80O7yQrsLy65MHmMIOEWo4rWDqvqnP5YxPpS3Mr8jSdGRqadzU8bi3By26NHu8kK7ywvtLizb9y5JNyZ8jiGqIlTxH3+f13ADAIDX0INThIyBntiy3XSx7Kq6nJZn24rdpqHCkC4NZHZDzWTC4bAqKytzcuxiRrvLC+0uL7S7wOYvkc4kfoqAU4yGQ7HFsY3/ILv1qpwWd/zsgG56+AVJ0lJ/vZ6784M5Kae7u7vo7z6bC7S7vNDu8kK7i8D27QlXM0RVhIyhYGzZrg/kvLyunskJxZdyzykAgAcQcIqNbUtDocnH9b6cF3n87GTAWb6wMeflAQCQawScYjM6IGNiTJJkVzdIOZ5/IxFwAADeQ8ApNsO9k8t1vrwUebyHgAMA8BYCTpExhiYDjl3vz3l5ExO2TvQMxh6/m4ADAPAAzqIqNq6Ao7iA0zc0pnMDo8qmcwOjGhwNS5J8DdUKNHKbBgBA6SPgFJlkPTjP/vYdfe2poxodn8hZ2ctb6L0BAHgDQ1TFJskcnH1H3sppuJGkqxY35fT4AADkCz04RSZZD875obHY8uL5taqtyu4VJJcF6nXHB5dl9ZgAABQKAafYJJmDc2EkHFt+9PMrdXnrvHzWCgCAksIQVZGZ0oNT5w4447HlebXkUgAAUiHgFBvXfaim9uAQcAAASBcBp8hMvQ9VJOBMTNgacA1RNdQUwR1cAQAoYnQF5EnFawdUdWCbNNijJbYtwzASbzg6edG9aA9O9Do1UiTcVFYk2RcAAEgi4ORN5b/ulBE6KUlKN57YDS2SGJ4CAGC2+LbMk2i4SYdtVCq86gtSwwJJ8QGH4SkAAGZCwMmHsSEZw+clSXZFlU7/0b9q4cKLkm9fUSlV18ceuk8RpwcHAICZ8W2ZDwPdk8uNC2XXNEm16V/Hpt/Vg9NUx1sGAMBMOIsqD4wLkwHHbkzRc5MEc3AAAJgdAk4eGAPvxJbteXMIOMOTAaeROTgAAMyIgJMPFyYDjua1znr3AXpwAACYFQJOHkwdopp9wGGSMQAAs0PAyQPjQoZDVPTgAAAwKwScfHCfRZVxwGEODgAAMyHg5EHmPTgMUQEAMBsZf1uaptkuabNlWevj1oUktVmWtSfTMkrdlIDT2CqNzG7/KT04XAcHAIAZZdyDY1nWPvdjJ9zIsqyDzuN1mZZR0ibC0mDP5OM5nEXFHBwAAGYnF0NU10rqcpa7JK3KQRmlY+icDHtCkmTX+aXKmlkfgtPEAQCYnVx8W/riHi9ItXF3d3eqp0tedfCYFjnLZ8bn6S/+wdLo6Khqat5I+xhv9g7Flkf6Q+o2BrNcy/wIhUKFrkJB0O7yQrvLC+0uXjMGHGfIKRC3uis6BJVAKMH2SbW2zn7IppSE+ztjy2+MNurZ34WcR3MLKZdefJF8DdWZV6xAvP5+J0O7ywvtLi+0uzjNGHDi59ik4SVN9uK0STowy/09pT/4thqc5XN2c0bHumZJs+bXM0QFAMBMsnEW1brIf2a7ZVn7LMvaZ5rm3c56X4qenrIw0nd2crnGr52/f5X6+vrV3Nw0q+PUVVfqw+9eIMMwsl1FAAA8J+OA4wQYf9y6Xc5iWYcbSQr3TwacifoW/ZcPXKzu7u6i79oDAKCUcaG/HLMHJk8RH6tLe2oSAADIAAEnx4yhc7HlifqUJ5QBAIAsIeDkWNVwMLZsNLYUsCYAAJQPAk6O1Y1OBpzKeQQcAADygYCTYw1jvbHlmvlMLAYAIB8IOLk0Ma7Gif7Iom2oYf7CAlcIAIDywFXjcmlwcoJxSI3yNdZldLjwRFjPnHxGR4NHM61ZQQwNDqn+zfpCVyPvaHd5od3lhXYXlq/Gl/Q5Ak62TYRlvGVJY0Myzp+KrT5nz5evPrNbLPz87Z/rL178iwwrCACANyydt1Tv1/sTPkfAybLqJ/6bKk48P219UE26LMN7SB3rPZbR/gAAlAsCTjYNhRKGG0l61V6q99Vm9nKPhkdjyx+++MO6btF1GR0v3y5cuKB58+YVuhp5R7vLC+0uL7S7sOZVz9O/v/TvCZ8j4GSRceFMbNmuadTIRav00uu9etNeqMeq23VLRWb3kRqZGIktr25drc9c8ZmMjpdv5XqLCtpdXmh3eaHdhffvIuDknjvgLPqAXlu/V7d+5xeSpDZfQ7K90ubuwamtrM34eAAAeBWniWeRceGd2PI/n7C14X/8IvbYl+H8G0kaCU/24NRU1mR8PAAAvIqAk0UjvW/Fls/Yfk3Yk8+1NmXe4+LuwampIOAAAJAMASeLQt1vxpbfsX2x5cXza9Vx3aUZH390whVw6MEBACAp5uBk0Ujv27Hldy27TEc//wlJUmWFIcPIbIKxxBwcAADSRcDJIvdZVAsuulRVldntIGOICgCA9DBElUX1I2djyxdfsizrx3efJs4QFQAAyRFwsmAsPKG/PXRSgYnJO4cvW9aW9XIYogIAID0EnCz430e79T+efVm1xpgkaVB1apg3P+vlTBmiogcHAICkCDhZcOLNU9pf8+exx0N1ubm645Tr4DAHBwCApAg4WfDe7h9rWUV37PH81sxPCU+E08QBAEgPAScLmodOT3lsm7fnpBzm4AAAkB4CThY0jE1OLj7ygXs08Z7fz0k5U3pwGKICACApAk4WNI5PBpxwc/ZPD49ikjEAAOkh4GRBUzgUWzbmLcxJGeMT4wrbYUlSpVGpqgqu0QgAQDIEnCxodgWcyhwFHHpvAABIHwEnU+PDatSgJGnMrlRtUyAnxTD/BgCA9BFwMjV4LrYYVJPqa6pzUsyUa+DQgwMAQEoEnAwZAz2x5XP2fDXUVOaknLHwWGyZHhwAAFIj4GTIGJwMOD12c84CjvtGm1wDBwCA1Ag4GRrvn7yCcVDNqq40clIOk4wBAEgfASdD4f6zseXzFT4ZRm4CDnNwAABIHwEnQ+ELkwGnv9KXs3Km3KahgiEqAABSIeBkov+M5v9yd+zhQJU/Z0W5TxOvrszNmVoAAHgFAScD1ftunfJ4sDp3Acc9RMUkYwAAUsv4ev+mabZL2mxZ1nrnsU9Sm/PvWsuytmZaRlEauaCKM7+afGhX61TdFTkrjtPEAQBIX8Y9OJZl7YtbtVGSGV1vmuamTMsoRsaFM1Me3zL6/2ikLje3aZCmnibOJGMAAFLL+h0bLcva43rYJulAtssoCq6A89LEFfoPe7nW5+gaOFLcJGOGqAAASClnt6Q2TbNNUtCyrIOptuvu7k71dNFqeOtVLXCW37Ejc28qJsbSak8oFJp1eefOT94SYnxkvCRft7m02wtod3mh3eWFdhevGQOOM8cm/g6SXTMFF0ntlmVtnun4ra2tM21SlCpPDMaWu52A429qTLs9s213zdnJYSn/PH/Jvm6lWu9M0e7yQrvLC+0uTjMGnARzbGZkmma7ZVm7nOVVlmUdmUvlilr/5BBVtAcnV7dpkDhNHACA2ch4krFpmusi/5ntrsc7TdM8bJrmYU3v/fEE48I7seVu2ydJqq/OXcDhSsYAAKQv4zk4zlCVP+7x8kyPW+zcZ1G94zS/Poc9OGMTk6eJcyVjAABS40J/c+XqwcnHEBU9OAAApI+AMxe2HTdElYc5OJwmDgBA2nJ2mrgXvPnOcb3+9m+nra8ID+sSjUnV1Ro1qjVgn1eF+nQ+3KTXQoMJjjTVuQvndL7m/KzqEhwOxpbpwQEAIDUCThJ/+/QOPXJ+v8YNI/EGSxfHFhv1N5Kkh1+N/Ms1btUAAEBqDFEl8Xz3M8nDTYFd1HBRoasAAEBRowcngfHxMb1eM6po/rtk1E6YBMOqVNDwa9So0by6Ki1oTN2zYshwjj+uqqrZv/QVRoU+suQjWuFfMet9AQAoJwScBA53/lQDFZFI4wtP6Kk/ekkVldmbQNzd3V30V4AEAKCUMUSVwMvHJ+8PumysLqvhBgAA5F7R9+D0nj+rt3pO5LXMY6FfSs6Z2EurL85r2QAAIHNFHXD+7tmd+v/OPaGRijxP9nVdZubdC1blt2wAAJCxoh6i+ucz+/IfblwqbFvXXfn7BSsfAADMTdH24Pz61UN6rdaWFAkaC8J2Xsuvsg19qOb9unzZ+/JaLgAAyFxRBJzhkUH9+F93KzQweQPLztCR2FDRNSO1+t7/dahAtQMAAKWmKALON350s/61tnfqStc8mPc3X5vfCgEAgJJWFHNwXqk8l/S5+okJ3Xz9l/NYGwAAUOqKogdn0BWzrh1uUIWTu6qMSn300j/QssWXF6hmAACgFBU84IyPj2mwYjLh/M0fHVBNdW2KPQAAAFIr+BBVMPRObLlhYoJwAwAAMlbwgNPTN3nmVP1EASsCAAA8o+ABJ9Q32YNTP1G4i/oBAADvKHjAOT/QE1uuswteHQAA4AEFTxR9Q8HYcp1d8DnPAADAAwoecC4MT17gr67wJ3UBAAAPKHjAGRg9H1uuUU0BawIAALyi4AFneGwgtlxXwSniAAAgcwUPOEPjF2LLtRX1BawJAADwioIHnOHwYGy5rrKhgDUBAABeUfCAM2IPx5YbqpsKWBMAAOAVRRBwRmPLBBwAAJANBQ84wxqLLTfW+gpXEQAA4BkFDzgjxnhseV59oIA1AQAAXlHwgDNshGPLzQ3+AtYEAAB4RREEHDu27JvXWsCaAAAAryh4wBmqmAw4geZFBawJAADwioIHnIEKI7bc4ltcwJoAAACvKHjAGTciAafatjWvcX6BawMAALwg44Bjmma7aZoHkjy3M93jNEzYM28EAACQhowDjmVZ+xKtN01znaS2dI/TOGHMvBEAAEAacjJEZZpmm6Su2ezTPFGVi6oAAIAylKs5OG2WZc0q4Myza3NUFQAAUG5m7DYxTbNdUvwlhrssyzqYZPt1yZ5LpVH16u7unu1uJSkUChW6CgVBu8sL7S4vtLu8lEK7Zww4yebYpBB05t/4JLWZprnKsqwjM+3UXO1Xa2v5XOivnNrqRrvLC+0uL7S7vBR7u7NxFtW6yH9muyRZlnXE6cEJKBJy0uKrXZhpVQAAACSl0YMzEyfMTLuJlGVZeyTtSfc4/kauYgwAALKj4Bf6i1o4/5JCVwEAAHhE0QScpQvfXegqAAAAjyiagLNs8YpCVwEAAHhEUQScxokJNc+LPxMdAABgbooi4PjC3KYBAABkT1EEnCZu0wAAALKoKAIOt2kAAADZVBQBp86oKXQVAACAhxRFwKkm4AAAgCwqioBTYzBEBQAAsqc4Ak4FAQcAAGRPUQSc2sr6QlcBAAB4SHEEnCoCDgAAyJ4iCTgNha4CAADwkKIIOPXV8wpdBQAA4CHFEXBqmwtdBQAA4CFFEXAaagg4AAAge4oi4MyrJ+AAAIDsKfhdLj/9yqf1s1eO6mf/vL3QVQEAAB5h2LZd6DoAAABkVVEMUQEAAGQTAQcAAHgOAQcAAHgOAQcAAHhOQc6iMk2zXVJIUptlWXsKUYd8Mk2zV5Il6YBlWbucdZ58DZx2bbYsa33cupBcbfVa+5O029Pvu2maPkltzr9rLcva6qz39Pudot2efr8lyTTNdc7i+nJ5v6Wk7fb8+x1lmubOUny/896D47wQsizroPN4Xeo9POEWy7LWx/0QePI1sCxrn/txorZ6sf3x7XZ4/X3fKMmMtt00zU1l8n5Pa7ez3tPvt2maqyStctqzyjTNtnJ4vxO123nK0+93lNOONme5pN7vQgxRXSupy1nukrSqAHXIN5/rh0Iqr9cgUVvLpf2eft8ty9rj+mutTZE2ef79TtJuyfvv9xHLsnY5PVhdlmWVy/udqN2Sx99vSXLa1+VaVVLvdyECji/u8YIC1CHfApKCpmnudh774p738mvgi3u8IMk6LyqL9935JRh0/oLzxT3t2fc7rt1SmbzfkkxJx51lX9xznn2/NbXdUnm8322uQCeV2PtdiIATUuSDUTacv/hCkkKuscpyeQ1Cmt7WROs8p4ze93bLsjY7yyGVz/vtbnfZvN9OoFueoo2J1pW8uHZ7/v02TXOdK7xHhVRC73chJhm/pMnE1ybpQAHqkDfO+LxlWdYR1+pyeg0StdWXYJ2nlMv7bppmu2sOwiqVyfudoN2mPP5+m6a5U9JxZ3gupMiXmuff70TtLpOf76Azn8Ynqa0Uf74LcqsG0zTvlnREkYlbu/JegTxKccaFJ18D5wfiSUlfdE3CnNZWr7U/vt3l8L47bd6tyC99SdpqWdZBr7/fidqtyNk0Xn+/o+2TIhNsNzvrvf5+T2t3Ofx8RzlhbqsibT9SSu8396ICAACew4X+AACA5xBwAACA5xBwAACA5xBwAACA5xTkXlRAuQgEAndL2iwpFAwGV8c9Fz07oUvS1mAweCTBIbJVj1WSdjr1uCVX5cxFIBBok9TuPDyYy9chQdnR10XBYHD9DJsDKCH04AA5FAwGdynyBepzAo37uT3Oc7fk+kvdOf7OXJaRgSclRW99kNf72Divy9ZcHDsQCPiSrN8ZCAR2J3oOQPbQgwPkXlCRXhz3F3lUVzAYDOWxHkXFCQE+5zUo1PUzQjk67kZNf78l6R9zVB4AF3pwgDwIBoMHJVmBQKBYe1EKJaCpN/PzkoRDXsFg8Eg+h+GAckXAAfJns6S7nTkn8DAnyPoKXQ+gnDFEBeRJMBjsCgQCexS5zP+0v+6dCa+PRjaNTHh1TUTeGQwG97gnC0u6T5EeEJ+ka4PB4NZAIBCdw7JKkeGvfUnKkZxLzTvzhOK3iV563Scp4Cr7UUVuS/Cks//6mSYtO22IDo8FnLlHcuq6XlKbU17C+s5Qp+jrE1JkLlOXs+0Bp36bg8HgwUAgEJ3E3CbpiNOjlqiMGd8D17YJj+ms97naJUl7gsFgyAm3u6Wpk5pTvEYp3++4um/SZG+Yz30coCzZts0//vEvh//8fn+7a9nn9/t7/X7/OufxurhtV/n9/gNx63b7/f5Nrsfr/H7/Yb/f3+Za96Tf798Zt19vgmP3+v1+3wzlPen3+1e5Hu+MtsFdttOWu2do+5Nx9fS5y3OOcyDVMdKo0ya/3787/jWPtjO+jn6//0Bcndri6pTOezDTMacdI9lzabxGyd7vda7HmxIcY2ei8vnHv3L5xxAVkEfOZNqtcv6KT1Mo/jCKTMx1z11JNI8lmOBMHss9qdmZCxKI9vw4PQyr4uaI/KMiw2tTyg4Gg6FEvT9RTu9Dm7ueTtld8WeUpTJTnZxeio3x+7naea2rt0WK9ALN9mytUNzjTI4ZO1aar1Gy9zt+qDP6HkWPwWRmlDWGqIA8c4ZWtsYNKcxWKMG6c3M8VpciQ1oHFfmSDrmGuiRnuCVu+3SYSbY9Lml1gvXJpFOnJwKBwCbntfXJ9fq4h9CcsORThvNjsnjMdF+j0Az12RMIBJ4MBAK2Iu/jkwxPodwRcIDCiJ42PtdrsCQ65Ts059pM8ikyFyZ+jop7bky65fhSPBdIv0pp1Wm3InNn9kha557L4wSQrZIOK/LlH5pF2QnN9piBQKAtrgcmypdqN9dyylP8A4GALxgM3uKEO1PS1kAgsDoYDG5OtR/gZQxRAQUQPW1c0wNOSNO//H05rk6bIl/SUmSoJVtneR1Mcqzlkg7M4jgz1skZvvIlOUPtsCJXit4THzKSXIwvpJnfg9kec1WCdVL2XqNvSJGhqWAweNCZwMzZeihrBBwg965Nsn6zpn8JBTX9y7QtwbpEPSDx2yTS5v4CdoZ9uqLzW5zgFXSdaRXdzj1nJp1yoqGjy32saA/DbIZP0qyTFDnb6ElNhrXYGWNxF1P0ueqS6HVM+R6kecz4OTLubScLSv81mvH9TvB6ePX6QkBaDNu2C10HwLOcS/JvUmToZGv8VYsDgcDOBKf73q3IF2KXJueabNbkfau+oci8lPuCweAuZ7Jr9FTirc5p0Xc76/Yocnqz+0s0oMiX+EyniYec7RQMBvc5+8fKlnPqcxqvwd2a/LJt09RTpneme7xEdYp73ue0dXPc+p2KzE864uzb5ZR7wFkXq0P0tUj1HjivRdJjRuvlrq9rnbvNW12ngyd7jaa85jO83+5J2D5FTl0n5KBsEXAAAIDnMEQFAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8h4ADAAA8pyrdDQ8fPnxDVVXVdkmLRDACAAD5NyHpzPj4+I7Vq1c/m2rDtO4mfvjw4RtqamoeXrp06WhdXd1wRUUFtyAHAAB5NTExYQwPD9edOnWqZnR09MupQk5aPTFVVVXbly5dOtrQ0DBEuAEAAIVQUVFhNzQ0DC1dunTUGVVKvm2ax1xUV1c3nIW6AQAAZMTJJItSbZNuwKmg5wYAABQDJ5OkzDBMFgYAAJ5DwAEAAJ6T9mniXtHZ2Vnz4Q9/+Kobb7yx97LLLhsOBoNVP/rRjxZef/31fStXrhwIBoNVTz/9tP/LX/7ymS1btvTMdKzvfOc7Cx999NFFV1555eCGDRuCkhQMBqskafv27WdaWlrCM2178uTJ2jvuuOPshg0b+uPLOHToUMOzzz7b1NbWNipJXV1dNatXrx5MtG0iPT09lQ899FBLMBiseuihh07P7tWCW7L3UJJOnDhR98wzz/i//e1vn+zo6OiNvm87d+5cGt12x44d7xSy/qWqp6encseOHYsCgcC4JPl8vvCaNWsGe3p6KtP5Odi/f3/T9u3bl15yySUj+/fv75pLHQ4dOtSwbdu2JZL005/+9NW5HANAfqV1mvivfvWr11esWJHyy75U7N+/v6m3t7eqo6OjN7qurq5u9RNPPHEs+suys7Oz5oc//KE/3S+k5cuXv/fzn//8Wff2Dz74YMvDDz+86MUXX3wlGnKi28aHp56ensqlS5d+wF2H6DGef/755vhfynfdddeSQCAwPlP9om197rnnmiXpscceO5lOe5Baovdbmgyj7vXJtkX6NmzY0Pb973//ZPTnqKenp3L9+vVXfOUrXznj/jlOJdnP0mzs37+/6YEHHlhEwMlcZ2dnzb333rv4sssuG/b5fOFAIBC+/PLLR6TI+7thw4b+2YZKrwTZ2267bZmUu9/X+Whjvl7Hzs7Olve///3vSvZ82Q1RxYebRFasWDHq8/nCqbaJF7/9li1bek6fPl37+OOP+93rm5ubpx23paUlfN111/V997vfXRhdd+jQoYZ777136fe///1pH/KHHnro9N/93d8t3L9/f1OqOm3YsKG/o6Ojd/78+eOzaQvmZu3atYOz/dxgZkePHm1w/5HQ0tIS3rFjx6nZHCMQCPC+FIm9e/f677zzzmUPPPDAqR07dryzZcuWno6Ojt6enp7KW2+9ta23t7dKivw8ffWrXz2T7nE3bNjQf9ttt53NpG6zLTNT27dvvyh+XXt7e7C9vT2YaPtsyEcb8/06JlN2ASf6V8JMli9fntZ2yXR2dtbM5jhHjx5tXLly5UD08bZt25bceOONve5f7G433XRT75YtW5ZlUkdkx969e2Mhds2aNYOFrIsXnT59uvbQoUMN7nUbNmzoT/dnGcWjs7Oz5mtf+9qyffv2dcX/btuwYUP/Bz/4wbSG3pMptSD78ssvN8av27BhQ3+6UxCQWtnNwVm7dm1aX0AtLS3hj33sY5dLkS62np6eyq9+9atLn3jiiZZTp079MlnwkCK9L4888sjC3bt3d830QY121X7pS1962z2M8Ytf/KL5ox/9aNK/Uk3THHz00UdTXgMAudfT01P51FNP+aO9gul+vpC+L37xi2c+8YlPXLlx48aeT37yk32f/vSn+1paWsLu1zo618w9Vy1+WLCvr68yGkaDwWClFOlpnc0xkJk777xzWao/3Nrb24PRHhyvu+uuu5YUug5eN+cP0nu2H1ydzYpk4nc71h3O9jGjXWwPPPDAIikSeB577LGTTzzxREui7V9++eXGvXv3hiXpqaee8n/84x/vSzYU5t5WkqLj0NHHPT09lZIU/UWbiN/vH5ciYarUvlRr72stms/OyDe65/TZ+dnPfta8fft2BYPBqhdeeKE52/XKNvNHZtG85tZnrVm95g899NDptra2kccee2xh9Odv48aNPe45CuvXr7/iwIEDx6JfnHv37vXfdttty9zbvPnmm7Xun8nbbrtt2YMPPhgLOekco1j82Z/92cUPPfTQ4nS2veWWW3oef/zxKW249dZblz355JMJf5fFu+uuu97+1re+9ZZ73c9//vOGD33oQ7P+vfOLX/yi+eabb076em7YsKE/+vsvkWwE2UOHDjX09PRU9vb2VlmW1bB+/fq+2faYpKqHey7Qxz/+8b5AIBC2LKuhra1tJFqP/fv3N508ebL26NGjDdu3b7/I5/OFt2zZ0hM/dyX6uLm5OXzHHXeclaTDhw83+Hy+8I033tj34osvNgaDwcqXX3650f05zUYb9+7d67/nnnuWSNLjjz/etXbt2sG9e/f6v/a1ry278cYbe++8886z6ZYR365kHQbREwpM0xyca73dyiIp58PKlSsHor88Ozo6eleuXHllKBSqTPQXoHvbqI997GOXRz+k0Te7q6urJll5x48fr5WkK664gm76AvjoRz/aF31vo2f3IHe2bNnSE/1y2Lt3r3/z5s1tK1euHNiyZUtP9MvM3SvQ0dHRu3nz5rYHHnjgVHT91VdfPeULub29PfjHf/zHy2dzDMxdNLjMNIyU6rXORpD90z/902XRCeodHR29Cxcu/MDZs2d/OZu2pKrHhg0b+nt7e8/cc889S6KTnTs6Ono3bNjQJkU+y842VZLk/o6I/8N67dq1g1/4whd67rnnniVXXHHFyIoVK0bXrl07uHTp0g8sX778WLSdH/vYx/x79+6N9SRno40dHR29fr9/fPv27Uujf0R3dHT0BoPByi1btvSsXLnyynTLSLfDYM2aNVdGT8yJfo+uXbv22Fx//spuDk6+XHPNNYP79+8PpLv9F77whR73m33dddf1JRqfjerq6qpdsmTJSHQobfny5e+N/ptp8jGy63Of+1zsl6l7Pg4y19PTUxmdzxbV0dHRe99995388Y9/7Jcifwgkmkjf1NQUPnbsWG2yY19xxRUj/f39lZkcA+lraWkJNzU1haO9Km4PPvhgi/tf/JwrafJnKz6EPvHEEy3uXp9EQfbee+9dGn38wx/+8Hj8H5ipeo2yVY877rjjrLse6Yr21q9YsWLUXa77j9vm5uYpr2umbYzasGFD/6lTp2Jz4Pbv399066239mazjKi9e/f6m5ubw+7X9ZprrhmMP1FnNubcg5OLYSEvmT9//vipU6fS/sUY/0N/7733nv7EJz5xZU9PT2Wi9PqjH/1o4fe+973jUuldl2Ouw0LFKvqLR0rd61ZIsx0WKhY9PT2VzzzzTHP8ZSrWrFkz+PzzzzdLkaHc/fv3T/td1t/fX5mqhzMYDFYtWbJkJJNjFMq3vvWtt+KHjWbj8ccfPxk/bDUbcxmekqSrr7564Pnnn2+Ov8bYli1beqKXy9i6deupRMPuM4XQlpaWhHVyB1kpEhCiw0LRk0CS/Z5NJFv1mI1LLrlk2mfQXd/4szczbaPbl770pbcfeeSRhWvXrj15/Pjx2uiQUTbLkCZ/d7r/QG9vbw9m8vNHD04KV1xxxcibb74ZCymJ/qqICoVCUz64pmkO9vf3V0YTbfRN6+vrS/gBf+yxxxZu3Lgx9kO/du3awfvuu+/k7bffPu1Mqdtuu23ZZz/72YQXBkTh9PT0VKbqdcPc3HvvvUvje3F++MMf+qNzEjo6Onr7+vqm9PTs3bvX/6lPfWrKZFb3z7IkPfLIIwu//OUvn5nNMZCZRx555OQLL7zQHP9+SpM9PMnmHra1tY2eP38+oyArRYZBbrjhhv7oUFF0fbq9D3Otx7Fjx2rd9XDr7OysSfSazFWmbXS76667ep555hl/Z2dnjfus4GyWIUVe176+vsroWWTRf+4/IGerbOfgRCeJnThxok6SHnjggUWHDx9u+NznPtcbfUFXrFgxetNNN/U++OCDLdELUTU1NYVvv/32Zbt27TolRX7RXn311YMnTpyo2759+0XR8dSOjo7e5557ruerX/3q0pUrVw5cf/31A9u3b7/o9OnTtT/+8Y/90UAUDAarfv3rXzd86EMf6ou/0vCWLVt61qxZM3jXXXctMU1zUIqk3Pb29mA64SZ64bmnn37aL0WuuXDDDTf0l9qk5GIRvZLx6dOna+OHH6NXMr7xxht7pcnX3r0tZ+TMXktLS3jbtm2nXnzxxcZnnnmmORQKVQaDwSrTNKdczfvFF198JTo5MRgMVoZCoUr3xd78fv/4d77znZPR4YWurq6aT37yk1NOBEh1jEOHDjV897vfXXj06NHGBx98sGWmq5wjsRUrVox++9vfPvm5z31uuXsOi5T6D0gp8jv1nnvuWdLZ2VkT/R092yB76NChhr6+vqro70D3l/FPfvKT5nQuHJluPY4ePTqlPQ888MCiaD2kyCVLon/wHjt2rDbV7/REgcrN/Qd2Ntro1tLSEr766qsH7rzzzmXR0YK5lDFTh0FHR0fvX//1Xy9yv66dnZ01M702qZTdlYwBAIUVfyVjKXLNsLVr1w4eO3asdu3atYOHDh1quP/++xe98MILzdu2bTsVHcaKD6HuPxz279/f1NLSEn711VdrpUiQbWtrG42fdHzZZZcNr169elCKXPz1ueeea77zzjvPSlJ8mYnqP1M99u7d6//BD37QcvPNN/cGAoFwonpE6xK9/ln0LCp3+WvWrBnctm3bkqNHjzZu27bt1K233tq7Y8eORY8++uiijRs39mzbtu3tZ555pvnhhx9e1NzcHN6xY8cp54KHGbfRbf/+/U2HDx9ucLdxLmXcddddS9ra2kaiHQY333zzFddff33frl27Tq1YsWLU/bpKkT9MUoWbma5kTMABACCL9u7d63/qqaf8mdwaBDPjVg0AAORZsvmWyB8CDgAAWbJ///6mH/zgBy3R+VqFrk85K9tJxgAAZBv3kioe9OAAAADPIeAAAADPSTfgTExMTBg5rQkAAEAanEwykWqbdAPOmeHh4brMqwQAAJAZJ5OcSbVNWgFnfHx8x6lTp2oGBwfr6ckBAACFMDExYQwODtafOnWqZnx8fEeqbdO60J8kHT58+IaqqqrtkhaJuTsAACD/JiSdGR8f37F69epnU22YdsABAAAoFfTEAAAAzyHgAAAAzyHgAAAAzyHgAAAAzyHgAAAAz/n/AdKJN87EX0clAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib import rc\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "names = [\"TuRBO-1\", \"EI\", \"Sobol\"]\n",
    "runs = [Y_turbo, Y_ei, Y_Sobol]\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "for name, run in zip(names, runs):\n",
    "    fx = np.maximum.accumulate(run.cpu())\n",
    "    plt.plot(fx, marker=\"\", lw=3)\n",
    "\n",
    "plt.plot([0, len(Y_turbo)], [fun.optimal_value, fun.optimal_value], \"k--\", lw=3)\n",
    "plt.xlabel(\"Function value\", fontsize=18)\n",
    "plt.xlabel(\"Number of evaluations\", fontsize=18)\n",
    "plt.title(\"20D Ackley\", fontsize=24)\n",
    "plt.xlim([0, len(Y_turbo)])\n",
    "plt.ylim([-15, 1])\n",
    "\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.legend(\n",
    "    names + [\"Global optimal value\"],\n",
    "    loc=\"lower center\",\n",
    "    bbox_to_anchor=(0, -0.08, 1, 1),\n",
    "    bbox_transform=plt.gcf().transFigure,\n",
    "    ncol=4,\n",
    "    fontsize=16,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "bento_stylesheets": {
   "bento/extensions/flow/main.css": true,
   "bento/extensions/kernel_selector/main.css": true,
   "bento/extensions/kernel_ui/main.css": true,
   "bento/extensions/new_kernel/main.css": true,
   "bento/extensions/system_usage/main.css": true,
   "bento/extensions/theme/main.css": true
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}