'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):

    X = np.atleast_2d(X)
    length_scale = _check_length_scale(X, self.length_scale)
    if (Y is None):
        dists = pdist((X / length_scale), metric='euclidean')
    else:
        if eval_gradient:
            raise ValueError('Gradient    can    only    be    evaluated    when    Y    is    None.')
        dists = cdist((X / length_scale), (Y / length_scale), metric='euclidean')
    if (self.nu == 0.5):
        K = np.exp((- dists))
    elif (self.nu == 1.5):
        K = (dists * math.sqrt(3))
        K = ((1.0 + K) * np.exp((- K)))
    elif (self.nu == 2.5):
        K = (dists * math.sqrt(5))
        K = (((1.0 + K) + ((K ** 2) / 3.0)) * np.exp((- K)))
    else:
        K = dists
        K[(K == 0.0)] += np.finfo(float).eps
        tmp = (math.sqrt((2 * self.nu)) * K)
        K.fill(((2 ** (1.0 - self.nu)) / gamma(self.nu)))
        K *= (tmp ** self.nu)
        K *= kv(self.nu, tmp)
    if (Y is None):
        K = squareform(K)
        np.fill_diagonal(K, 1)
    if eval_gradient:
        if self.hyperparameter_length_scale.fixed:
            K_gradient = np.empty((X.shape[0], X.shape[0], 0))
            return (K, K_gradient)
        if self.anisotropic:
            D = (((X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2) / (length_scale ** 2))
        else:
            D = squareform((dists ** 2))[:, :, np.newaxis]
        if (self.nu == 0.5):
            K_gradient = ((K[..., np.newaxis] * D) / np.sqrt(D.sum(2))[:, :, np.newaxis])
            K_gradient[(~ np.isfinite(K_gradient))] = 0
        elif (self.nu == 1.5):
            K_gradient = ((3 * D) * np.exp((- np.sqrt((3 * D.sum((-1))))))[..., np.newaxis])
        elif (self.nu == 2.5):
            tmp = np.sqrt((5 * D.sum((-1))))[..., np.newaxis]
            K_gradient = ((((5.0 / 3.0) * D) * (tmp + 1)) * np.exp((- tmp)))
        else:
            def f(theta):
                return self.clone_with_theta(theta)(X, Y)
            return (K, _approx_fprime(self.theta, f, 1e-10))
        if (not self.anisotropic):
            return (K, K_gradient[:, :].sum((-1))[:, :, np.newaxis])
        else:
            return (K, K_gradient)
    else:
        return K
