docs: def uncertain_conditional(Xnew_mu, Xnew_var, feat, kern, q_mu, q_sqrt, *,
                            mean_function=None, full_output_cov=False, full_cov=False, white=False):
      """
      Calculates the conditional for uncertain inputs Xnew, p(Xnew) = N(Xnew_mu, Xnew_var).
      See ``conditional`` documentation for further reference.
  
      :param Xnew_mu: mean of the inputs, size N x Din
      :param Xnew_var: covariance matrix of the inputs, size N x Din x Din
      :param feat: gpflow.InducingFeature object, only InducingPoints is supported
      :param kern: gpflow kernel object.
      :param q_mu: mean inducing points, size M x Dout
      :param q_sqrt: cholesky of the covariance matrix of the inducing points, size Dout x M x M
      :param full_output_cov: boolean wheter to compute covariance between output dimension.
                              Influences the shape of return value ``fvar``. Default is False
      :param white: boolean whether to use whitened representation. Default is False.
  
      :return fmean, fvar: mean and covariance of the conditional, size ``fmean`` is N x Dout,
              size ``fvar`` depends on ``full_output_cov``: if True ``f_var`` is N x Dout x Dout,
              if False then ``f_var`` is N x Dout
      """
  
      # TODO(VD): Tensorflow 1.7 doesn't support broadcasting in``tf.matmul`` and
      # ``tf.matrix_triangular_solve``. This is reported in issue 216.
      # As a temporary workaround, we are using ``tf.einsum`` for the matrix
      # multiplications and tiling in the triangular solves.
      # The code that should be used once the bug is resolved is added in comments.
  
      if not isinstance(feat, InducingPoints):
          raise NotImplementedError
  
      if full_cov:
          # TODO(VD): ``full_cov`` True would return a ``fvar`` of shape N x N x D x D,
          # encoding the covariance between input datapoints as well.
          # This is not implemented as this feature is only used for plotting purposes.
          raise NotImplementedError
  
      pXnew = Gaussian(Xnew_mu, Xnew_var)
  
      num_data = tf.shape(Xnew_mu)[0]  # number of new inputs (N)
      num_ind = tf.shape(q_mu)[0]  # number of inducing points (M)
      num_func = tf.shape(q_mu)[1]  # output dimension (D)
  
      q_sqrt_r = tf.matrix_band_part(q_sqrt, -1, 0)  # D x M x M
  
      eKuf = tf.transpose(expectation(pXnew, (kern, feat)))  # M x N (psi1)
      Kuu = features.Kuu(feat, kern, jitter=settings.jitter)  # M x M
      Luu = tf.cholesky(Kuu)  # M x M
  
      if not white:
          q_mu = tf.matrix_triangular_solve(Luu, q_mu, lower=True)
          Luu_tiled = tf.tile(Luu[None, :, :], [num_func, 1, 1])  # remove line once issue 216 is fixed
          q_sqrt_r = tf.matrix_triangular_solve(Luu_tiled, q_sqrt_r, lower=True)
  
      Li_eKuf = tf.matrix_triangular_solve(Luu, eKuf, lower=True)  # M x N
      fmean = tf.matmul(Li_eKuf, q_mu, transpose_a=True)
  
      eKff = expectation(pXnew, kern)  # N (psi0)
      eKuffu = expectation(pXnew, (kern, feat), (kern, feat))  # N x M x M (psi2)
      Luu_tiled = tf.tile(Luu[None, :, :], [num_data, 1, 1])  # remove this line, once issue 216 is fixed
      Li_eKuffu = tf.matrix_triangular_solve(Luu_tiled, eKuffu, lower=True)
      Li_eKuffu_Lit = tf.matrix_triangular_solve(Luu_tiled, tf.matrix_transpose(Li_eKuffu), lower=True)  # N x M x M
      cov = tf.matmul(q_sqrt_r, q_sqrt_r, transpose_b=True)  # D x M x M
  
      if mean_function is None or isinstance(mean_function, mean_functions.Zero):
          e_related_to_mean = tf.zeros((num_data, num_func, num_func), dtype=settings.float_type)
      else:
          # Update mean: \mu(x) + m(x)
          fmean = fmean + expectation(pXnew, mean_function)
  
          # Calculate: m(x) m(x)^T + m(x) \mu(x)^T + \mu(x) m(x)^T,
          # where m(x) is the mean_function and \mu(x) is fmean
          e_mean_mean = expectation(pXnew, mean_function, mean_function)  # N x D x D
          Lit_q_mu = tf.matrix_triangular_solve(Luu, q_mu, adjoint=True)
          e_mean_Kuf = expectation(pXnew, mean_function, (kern, feat))  # N x D x M
          # einsum isn't able to infer the rank of e_mean_Kuf, hence we explicitly set the rank of the tensor:
          e_mean_Kuf = tf.reshape(e_mean_Kuf, [num_data, num_func, num_ind])
          e_fmean_mean = tf.einsum("nqm,mz->nqz", e_mean_Kuf, Lit_q_mu)  # N x D x D
          e_related_to_mean = e_fmean_mean + tf.matrix_transpose(e_fmean_mean) + e_mean_mean
  
      if full_output_cov:
          fvar = (
                  tf.matrix_diag(tf.tile((eKff - tf.trace(Li_eKuffu_Lit))[:, None], [1, num_func])) +
                  tf.matrix_diag(tf.einsum("nij,dji->nd", Li_eKuffu_Lit, cov)) +
                  # tf.matrix_diag(tf.trace(tf.matmul(Li_eKuffu_Lit, cov))) +
                  tf.einsum("ig,nij,jh->ngh", q_mu, Li_eKuffu_Lit, q_mu) -
                  # tf.matmul(q_mu, tf.matmul(Li_eKuffu_Lit, q_mu), transpose_a=True) -
                  fmean[:, :, None] * fmean[:, None, :] +
                  e_related_to_mean
          )
      else:
          fvar = (
                  (eKff - tf.trace(Li_eKuffu_Lit))[:, None] +
                  tf.einsum("nij,dji->nd", Li_eKuffu_Lit, cov) +
                  tf.einsum("ig,nij,jg->ng", q_mu, Li_eKuffu_Lit, q_mu) -
                  fmean ** 2 +
                  tf.matrix_diag_part(e_related_to_mean)
          )
  
      return fmean, fvar