code: def kernel(x1, x2, method='gaussian', sigma=1, **kwargs):
      
      if method.lower() in ['gaussian', 'gauss', 'rbf']:
          K = np.exp(-dist(x1, x2) / (2 * sigma**2))
      return K
docs: Compute kernel matrix