code: def apply_inverse(self, y):
          
          self.recompute(quiet=False)
          r = np.array(y, dtype=np.float64, order="F")
          r = self._check_dimensions(r, check_dim=False)
  
          # Broadcast the mean function
          m = [slice(None)] + [np.newaxis for _ in range(len(r.shape) - 1)]
          r -= self._call_mean(self._x)[m]
  
          # Do the solve
          if len(r.shape) == 1:
              b = self.solver.apply_inverse(r, in_place=True).flatten()
          else:
              b = self.solver.apply_inverse(r, in_place=True)
          return b
docs: Self-consistently apply the inverse of the computed kernel matrix to
          some vector or matrix of samples. This method subtracts the mean,
          sorts the samples, then returns the samples in the correct (unsorted)
          order.
  
          :param y: ``(nsamples, )`` or ``(nsamples, K)``
              The vector (or matrix) of sample values.