content: from __future__ import print_function
  from __future__ import division
  import numpy as np
  from datetime import datetime
  from scipy.stats import norm
  from scipy.optimize import minimize
  
  
  def acq_max(ac, gp, y_max, bounds, random_state):
      """
      A function to find the maximum of the acquisition function
  
      It uses a combination of random sampling (cheap) and the 'L-BFGS-B'
      optimization method. First by sampling 1e5 points at random, and then
      running L-BFGS-B from 250 random starting points.
  
      Parameters
      ----------
      :param ac:
          The acquisition function object that return its point-wise value.
  
      :param gp:
          A gaussian process fitted to the relevant data.
  
      :param y_max:
          The current maximum known value of the target function.
  
      :param bounds:
          The variables bounds to limit the search of the acq max.
  
  
      Returns
      -------
      :return: x_max, The arg max of the acquisition function.
      """
  
      # Warm up with random points
      x_tries = random_state.uniform(bounds[:, 0], bounds[:, 1],
                                   size=(100000, bounds.shape[0]))
      ys = ac(x_tries, gp=gp, y_max=y_max)
      x_max = x_tries[ys.argmax()]
      max_acq = ys.max()
  
      # Explore the parameter space more throughly
      x_seeds = random_state.uniform(bounds[:, 0], bounds[:, 1],
                                  size=(250, bounds.shape[0]))
      for x_try in x_seeds:
          # Find the minimum of minus the acquisition function
          res = minimize(lambda x: -ac(x.reshape(1, -1), gp=gp, y_max=y_max),
                         x_try.reshape(1, -1),
                         bounds=bounds,
                         method="L-BFGS-B")
  
          # Store it if better than previous minimum(maximum).
          if max_acq is None or -res.fun[0] >= max_acq:
              x_max = res.x
              max_acq = -res.fun[0]
  
      # Clip output to make sure it lies within the bounds. Due to floating
      # point technicalities this is not always the case.
      return np.clip(x_max, bounds[:, 0], bounds[:, 1])
  
  
  class UtilityFunction(object):
      """
      An object to compute the acquisition functions.
      """
  
      def __init__(self, kind, kappa, xi):
          """
          If UCB is to be used, a constant kappa is needed.
          """
          self.kappa = kappa
  
          self.xi = xi
  
          if kind not in ['ucb', 'ei', 'poi']:
              err = "The utility function " \
                    "{} has not been implemented, " \
                    "please choose one of ucb, ei, or poi.".format(kind)
              raise NotImplementedError(err)
          else:
              self.kind = kind
  
      def utility(self, x, gp, y_max):
          if self.kind == 'ucb':
              return self._ucb(x, gp, self.kappa)
          if self.kind == 'ei':
              return self._ei(x, gp, y_max, self.xi)
          if self.kind == 'poi':
              return self._poi(x, gp, y_max, self.xi)
  
      @staticmethod
      def _ucb(x, gp, kappa):
          mean, std = gp.predict(x, return_std=True)
          return mean + kappa * std
  
      @staticmethod
      def _ei(x, gp, y_max, xi):
          mean, std = gp.predict(x, return_std=True)
          z = (mean - y_max - xi)/std
          return (mean - y_max - xi) * norm.cdf(z) + std * norm.pdf(z)
  
      @staticmethod
      def _poi(x, gp, y_max, xi):
          mean, std = gp.predict(x, return_std=True)
          z = (mean - y_max - xi)/std
          return norm.cdf(z)
  
  
  def unique_rows(a):
      """
      A functions to trim repeated rows that may appear when optimizing.
      This is necessary to avoid the sklearn GP object from breaking
  
      :param a: array to trim repeated rows from
  
      :return: mask of unique rows
      """
  
      # Sort array and kep track of where things should go back to
      order = np.lexsort(a.T)
      reorder = np.argsort(order)
  
      a = a[order]
      diff = np.diff(a, axis=0)
      ui = np.ones(len(a), 'bool')
      ui[1:] = (diff != 0).any(axis=1)
  
      return ui[reorder]
  
  
  class BColours(object):
      BLUE = '\033[94m'
      CYAN = '\033[36m'
      GREEN = '\033[32m'
      MAGENTA = '\033[35m'
      RED = '\033[31m'
      ENDC = '\033[0m'
  
  
  class PrintLog(object):
  
      def __init__(self, params):
  
          self.ymax = None
          self.xmax = None
          self.params = params
          self.ite = 1
  
          self.start_time = datetime.now()
          self.last_round = datetime.now()
  
          # sizes of parameters name and all
          self.sizes = [max(len(ps), 7) for ps in params]
  
          # Sorted indexes to access parameters
          self.sorti = sorted(range(len(self.params)),
                              key=self.params.__getitem__)
  
      def reset_timer(self):
          self.start_time = datetime.now()
          self.last_round = datetime.now()
  
      def print_header(self, initialization=True):
  
          if initialization:
              print("{}Initialization{}".format(BColours.RED,
                                                BColours.ENDC))
          else:
              print("{}Bayesian Optimization{}".format(BColours.RED,
                                                       BColours.ENDC))
  
          print(BColours.BLUE + "-" * (29 + sum([s + 5 for s in self.sizes])) +
              BColours.ENDC)
  
          print("{0:>{1}}".format("Step", 5), end=" | ")
          print("{0:>{1}}".format("Time", 6), end=" | ")
          print("{0:>{1}}".format("Value", 10), end=" | ")
  
          for index in self.sorti:
              print("{0:>{1}}".format(self.params[index],
                                      self.sizes[index] + 2),
                    end=" | ")
          print('')
  
      def print_step(self, x, y, warning=False):
  
          print("{:>5d}".format(self.ite), end=" | ")
  
          m, s = divmod((datetime.now() - self.last_round).total_seconds(), 60)
          print("{:>02d}m{:>02d}s".format(int(m), int(s)), end=" | ")
  
          if self.ymax is None or self.ymax < y:
              self.ymax = y
              self.xmax = x
              print("{0}{2: >10.5f}{1}".format(BColours.MAGENTA,
                                               BColours.ENDC,
                                               y),
                    end=" | ")
  
              for index in self.sorti:
                  print("{0}{2: >{3}.{4}f}{1}".format(
                              BColours.GREEN, BColours.ENDC,
                              x[index],
                              self.sizes[index] + 2,
                              min(self.sizes[index] - 3, 6 - 2)
                          ),
                        end=" | ")
          else:
              print("{: >10.5f}".format(y), end=" | ")
              for index in self.sorti:
                  print("{0: >{1}.{2}f}".format(x[index],
                                                self.sizes[index] + 2,
                                                min(self.sizes[index] - 3, 6 - 2)),
                        end=" | ")
  
          if warning:
              print("{}Warning: Test point chose at "
                    "random due to repeated sample.{}".format(BColours.RED,
                                                              BColours.ENDC))
  
          print()
  
          self.last_round = datetime.now()
          self.ite += 1
  
      def print_summary(self):
          pass
  