text: """
  Utilities for Gaussian process (GP) inference.
  """
  
  import numpy as np
  from scipy.linalg import solve_triangular
  from scipy.spatial.distance import cdist
  import itertools
  
  
  def kern_exp_quad(xmat1, xmat2, ls, alpha):
      """
      Exponentiated quadratic kernel function (aka squared exponential kernel aka
      RBF kernel).
      """
      return alpha ** 2 * kern_exp_quad_noscale(xmat1, xmat2, ls)
  
  
  def kern_exp_quad_noscale(xmat1, xmat2, ls):
      """
      Exponentiated quadratic kernel function (aka squared exponential kernel aka
      RBF kernel), without scale parameter.
      """
      distmat = squared_euc_distmat(xmat1, xmat2)
      sq_norm = -distmat / (2 * ls ** 2)
      return np.exp(sq_norm)
  
  
  def squared_euc_distmat(xmat1, xmat2, coef=1.0):
      """
      Distance matrix of squared euclidean distance (multiplied by coef) between
      points in xmat1 and xmat2.
      """
      return coef * cdist(xmat1, xmat2, 'sqeuclidean')
  
  
  def kern_distmat(xmat1, xmat2, ls, alpha, distfn):
      """
      Kernel for a given distmat, via passed in distfn (which is assumed to be fn
      of xmat1 and xmat2 only).
      """
      distmat = distfn(xmat1, xmat2)
      kernmat = alpha ** 2 * np.exp(-distmat / (2 * ls ** 2))
      return kernmat
  
  
  def kern_simple_list(xlist1, xlist2, ls, alpha, base_dist=5.0):
      """
      Kernel for two lists containing elements that can be compared for equality.
      K(a,b) = 1 + base_dist if a and b are equal and K(a,b) = base_dist otherwise.
      """
      distmat = simple_list_distmat(xlist1, xlist2)
      distmat = distmat + base_dist
      kernmat = alpha ** 2 * np.exp(-distmat / (2 * ls ** 2))
      return kernmat
  
  
  def simple_list_distmat(xlist1, xlist2, weight=1.0, additive=False):
      """
      Return distance matrix containing zeros when xlist1[i] == xlist2[j] and 0 otherwise.
      """
      prod_list = list(itertools.product(xlist1, xlist2))
      len1 = len(xlist1)
      len2 = len(xlist2)
      try:
          binary_mat = np.array([x[0] != x[1] for x in prod_list]).astype(int)
      except:
          # For cases where comparison returns iterable of bools
          binary_mat = np.array([all(x[0] != x[1]) for x in prod_list]).astype(int)
  
      binary_mat = binary_mat.reshape(len1, len2)
  
      if additive:
          distmat = weight + binary_mat
      else:
          distmat = weight * binary_mat
  
      return distmat
  
  
  def get_product_kernel(kernel_list, additive=False):
      """Given a list of kernel functions, return product kernel."""
  
      def product_kernel(x1, x2, ls, alpha):
          """Kernel returning elementwise-product of kernel matrices from kernel_list."""
          mat_prod = kernel_list[0](x1, x2, ls, 1.0)
          for kernel in kernel_list[1:]:
              if additive:
                  mat_prod = mat_prod + kernel(x1, x2, ls, 1.0)
              else:
                  mat_prod = mat_prod * kernel(x1, x2, ls, 1.0)
          mat_prod = alpha ** 2 * mat_prod
          return mat_prod
  
      return product_kernel
  
  
  def get_cholesky_decomp(k11_nonoise, sigma, psd_str):
      """Return cholesky decomposition."""
      if psd_str == 'try_first':
          k11 = k11_nonoise + sigma ** 2 * np.eye(k11_nonoise.shape[0])
          try:
              return stable_cholesky(k11, False)
          except np.linalg.linalg.LinAlgError:
              return get_cholesky_decomp(k11_nonoise, sigma, 'project_first')
      elif psd_str == 'project_first':
          k11_nonoise = project_symmetric_to_psd_cone(k11_nonoise)
          return get_cholesky_decomp(k11_nonoise, sigma, 'is_psd')
      elif psd_str == 'is_psd':
          k11 = k11_nonoise + sigma ** 2 * np.eye(k11_nonoise.shape[0])
          return stable_cholesky(k11)
  
  
  def stable_cholesky(mmat, make_psd=True, verbose=False):
      """Return a 'stable' cholesky decomposition of mmat."""
      if mmat.size == 0:
          return mmat
      try:
          lmat = np.linalg.cholesky(mmat)
      except np.linalg.linalg.LinAlgError as e:
          if not make_psd:
              raise e
          diag_noise_power = -11
          max_mmat = np.diag(mmat).max()
          diag_noise = np.diag(mmat).max() * 1e-11
          break_loop = False
          while not break_loop:
              try:
                  lmat = np.linalg.cholesky(
                      mmat + ((10 ** diag_noise_power) * max_mmat) * np.eye(mmat.shape[0])
                  )
                  break_loop = True
              except np.linalg.linalg.LinAlgError:
                  if diag_noise_power > -9:
                      if verbose:
                          print(
                              '\t*stable_cholesky failed with '
                              'diag_noise_power=%d.' % (diag_noise_power)
                          )
                  diag_noise_power += 1
              if diag_noise_power >= 5:
                  print('\t*stable_cholesky failed: added diag noise = %e' % (diag_noise))
      return lmat
  
  
  def project_symmetric_to_psd_cone(mmat, is_symmetric=True, epsilon=0):
      """Project symmetric matrix mmat to the PSD cone."""
      if is_symmetric:
          try:
              eigvals, eigvecs = np.linalg.eigh(mmat)
          except np.linalg.LinAlgError:
              print('\tLinAlgError encountered with np.eigh. Defaulting to eig.')
              eigvals, eigvecs = np.linalg.eig(mmat)
              eigvals = np.real(eigvals)
              eigvecs = np.real(eigvecs)
      else:
          eigvals, eigvecs = np.linalg.eig(mmat)
      clipped_eigvals = np.clip(eigvals, epsilon, np.inf)
      return (eigvecs * clipped_eigvals).dot(eigvecs.T)
  
  
  def solve_lower_triangular(amat, b):
      """Solves amat*x=b when amat is lower triangular."""
      return solve_triangular_base(amat, b, lower=True)
  
  
  def solve_upper_triangular(amat, b):
      """Solves amat*x=b when amat is upper triangular."""
      return solve_triangular_base(amat, b, lower=False)
  
  
  def solve_triangular_base(amat, b, lower):
      """Solves amat*x=b when amat is a triangular matrix."""
      if amat.size == 0 and b.shape[0] == 0:
          return np.zeros((b.shape))
      else:
          return solve_triangular(amat, b, lower=lower)
  
  
  def sample_mvn(mu, covmat, nsamp):
      """
      Sample from multivariate normal distribution with mean mu and covariance
      matrix covmat.
      """
      mu = mu.reshape(-1)
      ndim = len(mu)
      lmat = stable_cholesky(covmat)
      umat = np.random.normal(size=(ndim, nsamp))
      return lmat.dot(umat).T + mu
  
  
  def gp_post(x_train, y_train, x_pred, ls, alpha, sigma, kernel, full_cov=True):
      """Compute parameters of GP posterior"""
      k11_nonoise = kernel(x_train, x_train, ls, alpha)
      lmat = get_cholesky_decomp(k11_nonoise, sigma, 'try_first')
      smat = solve_upper_triangular(lmat.T, solve_lower_triangular(lmat, y_train))
      k21 = kernel(x_pred, x_train, ls, alpha)
      mu2 = k21.dot(smat)
      k22 = kernel(x_pred, x_pred, ls, alpha)
      vmat = solve_lower_triangular(lmat, k21.T)
      k2 = k22 - vmat.T.dot(vmat)
      if full_cov is False:
          k2 = np.sqrt(np.diag(k2))
      return mu2, k2
  