'Return the kernel k(X, Y) and optionally its gradient.
Parameters
X : array, shape (n_samples_X, n_features)
Left argument of the returned kernel k(X, Y)
Y : array, shape (n_samples_Y, n_features), (optional, default=None)
Right argument of the returned kernel k(X, Y). If None, k(X, X)
if evaluated instead.
eval_gradient : bool (optional, default=False)
Determines whether the gradient with respect to the kernel
hyperparameter is determined. Only supported when Y is None.
Returns
K : array, shape (n_samples_X, n_samples_Y)
Kernel k(X, Y)
K_gradient : array (opt.), shape (n_samples_X, n_samples_X, n_dims)
The gradient of the kernel k(X, X) with respect to the
hyperparameter of the kernel. Only returned when eval_gradient
is True.'

def __call__(self, X, Y=None, eval_gradient=False):

    pairwise_kernels_kwargs = self.pairwise_kernels_kwargs
    if (self.pairwise_kernels_kwargs is None):
        pairwise_kernels_kwargs = {}
    X = np.atleast_2d(X)
    K = pairwise_kernels(X, Y, metric=self.metric, gamma=self.gamma, filter_params=True, **pairwise_kernels_kwargs)
    if eval_gradient:
        if self.hyperparameter_gamma.fixed:
            return (K, np.empty((X.shape[0], X.shape[0], 0)))
        else:
            def f(gamma):
                return pairwise_kernels(X, Y, metric=self.metric, gamma=np.exp(gamma), filter_params=True, **pairwise_kernels_kwargs)
            return (K, _approx_fprime(self.theta, f, 1e-10))
    else:
        return K
