code: def predict_wishart_embedding(self, Xnew, kern=None, mean=True, covariance=True):
          
          if kern is None:
              kern = self.kern
  
          mu_jac, var_jac = self.predict_jacobian(Xnew, kern, full_cov=False)
          mumuT = np.einsum('iqd,ipd->iqp', mu_jac, mu_jac)
          Sigma = np.zeros(mumuT.shape)
          if var_jac.ndim == 4: # Missing data
              Sigma = var_jac.sum(-1)
          else:
              Sigma = self.output_dim*var_jac
  
          G = 0.
          if mean:
              G += mumuT
          if covariance:
              G += Sigma
          return G
docs: Predict the wishart embedding G of the GP. This is the density of the
          input of the GP defined by the probabilistic function mapping f.
          G = J_mean.T*J_mean + output_dim*J_cov.
  
          :param array-like Xnew: The points at which to evaluate the magnification.
          :param :py:class:`~GPy.kern.Kern` kern: The kernel to use for the magnification.
  
          Supplying only a part of the learning kernel gives insights into the density
          of the specific kernel part of the input function. E.g. one can see how dense the
          linear part of a kernel is compared to the non-linear part etc.