docs: def _conditional(Xnew, X, kern, f, *, full_cov=False, q_sqrt=None, white=False, full_output_cov=None):
      """
      Given f, representing the GP at the points X, produce the mean and
      (co-)variance of the GP at the points Xnew.
  
      Additionally, there may be Gaussian uncertainty about f as represented by
      q_sqrt. In this case `f` represents the mean of the distribution and
      q_sqrt the square-root of the covariance.
  
      Additionally, the GP may have been centered (whitened) so that
          p(v) = N(0, I)
          f = L v
      thus
          p(f) = N(0, LL^T) = N(0, K).
      In this case `f` represents the values taken by v.
  
      The method can either return the diagonals of the covariance matrix for
      each output (default) or the full covariance matrix (full_cov=True).
  
      We assume R independent GPs, represented by the columns of f (and the
      first dimension of q_sqrt).
  
      :param Xnew: data matrix, size N x D. Evaluate the GP at these new points
      :param X: data points, size M x D.
      :param kern: GPflow kernel.
      :param f: data matrix, M x R, representing the function values at X,
          for R functions.
      :param q_sqrt: matrix of standard-deviations or Cholesky matrices,
          size M x R or R x M x M.
      :param white: boolean of whether to use the whitened representation as
          described above.
      :return:
          - mean:     N x R
          - variance: N x R (full_cov = False), R x N x N (full_cov = True)
      """
      logger.debug("Conditional: Kernel")
      num_data = tf.shape(X)[-2]  # M
      Kmm = kern.K(X) + tf.eye(num_data, dtype=settings.float_type) * settings.numerics.jitter_level  #  [..., M, M]
      Kmn = kern.K(X, Xnew)  # [M, ..., N]
  
      if full_cov:
          Knn = kern.K(Xnew)  # [...,N,N]
      else:
          Knn = kern.Kdiag(Xnew)  # [...,N]
  
      mean, var = base_conditional(Kmn, Kmm, Knn, f, full_cov=full_cov, q_sqrt=q_sqrt, white=white)
  
      return mean, var