import numpy as np


def linear_kernel(X1, X2):
    return X1 @ X2.T


def rbf_kernel(X1, X2, gamma=1.0):
    sq_dists = (
        np.sum(X1 ** 2, axis=1)[:, None]
        + np.sum(X2 ** 2, axis=1)[None, :]
        - 2 * X1 @ X2.T
    )
    return np.exp(-gamma * sq_dists)


def simplified_smo(
    X, y, C=1.0, kernel=linear_kernel, tol=1e-3, max_passes=20
):
    n = len(y)
    alpha = np.zeros(n)
    b = 0.0
    K = kernel(X, X)

    def decision(i):
        return np.sum(alpha * y * K[:, i]) + b

    passes = 0
    while passes < max_passes:
        num_changed = 0
        for i in range(n):
            Ei = decision(i) - y[i]
            if (y[i] * Ei < -tol and alpha[i] < C) or (
                y[i] * Ei > tol and alpha[i] > 0
            ):
                j = np.random.randint(0, n - 1)
                j = j if j < i else j + 1
                Ej = decision(j) - y[j]

                alpha_i_old, alpha_j_old = alpha[i], alpha[j]
                if y[i] != y[j]:
                    L = max(0, alpha[j] - alpha[i])
                    H = min(C, C + alpha[j] - alpha[i])
                else:
                    L = max(0, alpha[i] + alpha[j] - C)
                    H = min(C, alpha[i] + alpha[j])
                if L == H:
                    continue

                eta = 2 * K[i, j] - K[i, i] - K[j, j]
                if eta >= 0:
                    continue

                alpha[j] -= y[j] * (Ei - Ej) / eta
                alpha[j] = np.clip(alpha[j], L, H)
                if abs(alpha[j] - alpha_j_old) < 1e-5:
                    continue

                alpha[i] += y[i] * y[j] * (alpha_j_old - alpha[j])

                b1 = (
                    b - Ei
                    - y[i] * (alpha[i] - alpha_i_old) * K[i, i]
                    - y[j] * (alpha[j] - alpha_j_old) * K[i, j]
                )
                b2 = (
                    b - Ej
                    - y[i] * (alpha[i] - alpha_i_old) * K[i, j]
                    - y[j] * (alpha[j] - alpha_j_old) * K[j, j]
                )
                if 0 < alpha[i] < C:
                    b = b1
                elif 0 < alpha[j] < C:
                    b = b2
                else:
                    b = (b1 + b2) / 2

                num_changed += 1
        passes = passes + 1 if num_changed == 0 else 0
    return alpha, b
