docs: def fully_correlated_conditional_repeat(Kmn, Kmm, Knn, f, *, full_cov=False, full_output_cov=False, q_sqrt=None,
                                          white=False):
      """
      This function handles conditioning of multi-output GPs in the case where the conditioning
      points are all fully correlated, in both the prior and posterior.
  
      Note: This conditional can handle 'repetitions' R, given in `f` and `q_sqrt`.
  
      :param Kmn: LM x N x P
      :param Kmm: LM x LM
      :param Knn: N x P or N x P x N x P
      :param f: data matrix, LM x R
      :param q_sqrt: R x LM x LM  or R x ML
      :param full_cov: calculate covariance between inputs
      :param full_output_cov: calculate covariance between outputs
      :param white: use whitened representation
      :return:
          - mean: R x N x P
          - variance: R x N x P, R x N x P x P, R x P x N x N, R x N x P x N x P
      """
      logger.debug("fully correlated conditional")
      R = tf.shape(f)[1]
      M, N, K = [tf.shape(Kmn)[i] for i in range(Kmn.shape.ndims)]
      Lm = tf.cholesky(Kmm)
  
      # Compute the projection matrix A
      # Lm: M x M    Kmn: M x NK
      Kmn = tf.reshape(Kmn, (M, N * K))  # M x NK
      A = tf.matrix_triangular_solve(Lm, Kmn, lower=True)  # M x NK
      Ar = tf.reshape(A, (M, N, K))
  
      # compute the covariance due to the conditioning
      if full_cov and full_output_cov:
          # fvar = Knn - tf.matmul(Ar, Ar, transpose_a=True)  # NK x NK, then reshape?
          fvar = Knn - tf.tensordot(Ar, Ar, [[0], [0]])  # N x K x N x K
      elif full_cov and not full_output_cov:
          At = tf.transpose(Ar)  # K x N x M
          fvar = Knn - tf.matmul(At, At, transpose_b=True)  # K x N x N
      elif not full_cov and full_output_cov:
          # This transpose is annoying
          At = tf.transpose(Ar, [1, 0, 2])  # N x M x K
          # fvar = Knn - tf.einsum('mnk,mnl->nkl', Ar, Ar)
          fvar = Knn - tf.matmul(At, At, transpose_a=True)  # N x K x K
      elif not full_cov and not full_output_cov:
          # Knn: N x K
          fvar = Knn - tf.reshape(tf.reduce_sum(tf.square(A), [0]), (N, K))  # Can also do this with a matmul
  
      # another backsubstitution in the unwhitened case
      if not white:
          # A = tf.matrix_triangular_solve(tf.matrix_transpose(Lm), A, lower=False)  # M x NK
          raise NotImplementedError("Need to verify this.")  # pragma: no cover
  
      # f: M x R
      fmean = tf.matmul(f, A, transpose_a=True)  # R x M  *  M x NK  ->  R x NK
      fmean = tf.reshape(fmean, (R, N, K))  # R x N x K
  
      if q_sqrt is not None:
          Lf = tf.matrix_band_part(q_sqrt, -1, 0)  # R x M x M
          if q_sqrt.get_shape().ndims == 3:
              A_tiled = tf.tile(A[None, :, :], tf.stack([R, 1, 1]))  # R x M x NK
              LTA = tf.matmul(Lf, A_tiled, transpose_a=True)  # R x M x NK
          elif q_sqrt.get_shape().ndims == 2:  # pragma: no cover
              raise NotImplementedError("Does not support diagonal q_sqrt yet...")
          else:  # pragma: no cover
              raise ValueError("Bad dimension for q_sqrt: %s" %
                               str(q_sqrt.get_shape().ndims))
  
          if full_cov and full_output_cov:
              addvar = tf.matmul(LTA, LTA, transpose_a=True)  # R x NK x NK
              fvar = fvar[None, :, :, :, :] + tf.reshape(addvar, (R, N, K, N, K))
          elif full_cov and not full_output_cov:
              LTAr = tf.transpose(tf.reshape(LTA, [R, M, N, K]), [0, 3, 1, 2])  # R x K x M x N
              addvar = tf.matmul(LTAr, LTAr, transpose_a=True)  # R x K x N x N
              fvar = fvar[None, ...] + addvar  # R x K x N x N
          elif not full_cov and full_output_cov:
              LTAr = tf.transpose(tf.reshape(LTA, (R, M, N, K)), [0, 2, 3, 1])  # R x N x K x M
              fvar = fvar[None, ...] + tf.matmul(LTAr, LTAr, transpose_b=True)  # R x N x K x K
          elif not full_cov and not full_output_cov:
              addvar = tf.reshape(tf.reduce_sum(tf.square(LTA), axis=1), (R, N, K))  # R x N x K
              fvar = fvar[None, ...] + addvar  # R x N x K
      else:
          fvar = tf.broadcast_to(fvar[None], tf.shape(fmean))
      return fmean, fvar