code: def predict(self, X1new, X2new):
          
          k1xf = self.kern1.K(X1new, self.X1)
          k2xf = self.kern2.K(X2new, self.X2)
          A = k1xf.dot(self.U1)
          B = k2xf.dot(self.U2)
          mu = A.dot(self.Ytilde.reshape(self.num_data1, self.num_data2, order='F')).dot(B.T).flatten(order='F')
          k1xx = self.kern1.Kdiag(X1new)
          k2xx = self.kern2.Kdiag(X2new)
          BA = np.kron(B, A)
          var = np.kron(k2xx, k1xx) - np.sum(BA**2*self.Wi, 1) + self.likelihood.variance
  
          return mu[:, None], var[:, None]
docs: Return the predictive mean and variance at a series of new points X1new, X2new
          Only returns the diagonal of the predictive variance, for now.
  
          :param X1new: The points at which to make a prediction
          :type X1new: np.ndarray, Nnew x self.input_dim1
          :param X2new: The points at which to make a prediction
          :type X2new: np.ndarray, Nnew x self.input_dim2