import numpy as np
from scipy.spatial.distance import cdist
from sklearn.datasets import load_iris
from sklearn.metrics.pairwise import rbf_kernel as skrbf_kernel

def rbf_kernel(X, Y=None, gamma=None):
    if Y is None:
        Y = X
    if gamma is None:
        gamma = 1 / X.shape[1]
    K = np.zeros((X.shape[0], Y.shape[0]))
    for i in range(X.shape[0]):
        for j in range(Y.shape[0]):
            K[i, j] = np.exp(-gamma * np.sum(np.square(X[i] - Y[j])))
    return K

X, _ = load_iris(return_X_y=True)
K1 = rbf_kernel(X)
K2 = skrbf_kernel(X)
assert np.allclose(K1, K2)

def rbf_kernel(X, Y=None, gamma=None):
    if Y is None:
        Y = X
    if gamma is None:
        gamma = 1 / X.shape[1]
    return np.exp(-gamma * cdist(X, Y, metric='sqeuclidean'))

X, _ = load_iris(return_X_y=True)
K1 = rbf_kernel(X)
K2 = skrbf_kernel(X)
assert np.allclose(K1, K2)