'Return the kernel k(X, Y) and optionally its gradient.
Note that this compound kernel returns the results of all simple kernel
stacked along an additional axis.
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.
Returns
K : array, shape (n_samples_X, n_samples_Y, n_kernels)
Kernel k(X, Y)
K_gradient : array, shape (n_samples_X, n_samples_X, n_dims, n_kernels)
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):

    if eval_gradient:
        K = []
        K_grad = []
        for kernel in self.kernels:
            (K_single, K_grad_single) = kernel(X, Y, eval_gradient)
            K.append(K_single)
            K_grad.append(K_grad_single[..., np.newaxis])
        return (np.dstack(K), np.concatenate(K_grad, 3))
    else:
        return np.dstack([kernel(X, Y, eval_gradient) for kernel in self.kernels])
