function [w b alphas] = simpleSMO(y,X)

m = length(X(:,1)); % Number of examples

% Build up the kernel matrix
%------------------------------------------------------
K = zeros(m,m);
for i = 1:(m-1)
    for j = (i+1):m
        K(i,j) = kernel(X(i,:)',X(j,:)');
    end
end
K = K + K';
for i = 1:m
    K(i,i) = kernel(X(i,:)',X(i,:)');
end

% SMO and numerical parameters
%------------------------------------------------------
alphas = zeros(m,1); % Lagrange multiplers
C = 1; % Regularization parameter
b = 0;  % Intercept parameter

tol = 1e-4; % Numerical tolerance

max_passes = 100;
passes = 0;
while(passes<max_passes)
    num_changed_alphas = 0;
    for i = 1:m % Outer loop, looks for KKT violation
        % Ei = y(i)*[y(i)*f(x(i))-1] and
        % y(i)*Ei =  y(i)*f(x(i))-1
        Ei = alphas'*(y.*K(:,i))+b-y(i);
        if( ( (y(i)*Ei<-tol) && (alphas(i)<C) ) || ...
            ( (y(i)*Ei> tol) && (alphas(i)>0) ) )
            % 1. Pick a second index j (~=i) randomly
            j = randi(m-1,1);
            if (j>=i)
                j = j+1;
            end
            Ej = alphas'*(y.*K(:,j))+b-y(j);
            % 2. Save old alphas
            alphas_old = alphas;
            % 3. Compute bounds
            if(y(i)~=y(j))
                L = max([0 alphas(j)-alphas(i)]);
                H = min([C C+alphas(j)-alphas(i)]);
            else
                L = max([0 alphas(i)+alphas(j)-C]);
                H = min([C alphas(i)+alphas(j)]);
            end
            if(L==H)
                continue;
            end
            % 4. Compute eta
            eta = 2*K(i,j)-K(i,i)-K(j,j);
            if(abs(eta)<tol)
                continue;
            end
            % 5. Get new value for alphaj and clip
            alphas(j) = alphas_old(j) + -y(j)*(Ei-Ej)/eta;
            if(alphas(j)>H)
                alphas(j) = H;
            elseif (alphas(j)<L)
                alphas(j) = L;
            end
            if(abs(alphas(j)-alphas_old(j))<tol) % Did we actually move?
                continue;
            end
            % 6. Update value for alphas(i)
            alphas(i) = alphas(i) + y(i)*y(j)*(alphas_old(j)-alphas(j));
            % 7. Compute b1 and b2
            b1 = b-Ei-y(i)*(alphas(i)-alphas_old(i))*K(i,i)+...
                     -y(j)*(alphas(j)-alphas_old(j))*K(i,j);
            b2 = b-Ej-y(i)*(alphas(i)-alphas_old(i))*K(i,j)+...
                     -y(j)*(alphas(j)-alphas_old(j))*K(j,j);
            if( (tol<alphas(i)) && (alphas(i)<(C-tol)) )
                b = b1;
            elseif( (tol<alphas(j)) && (alphas(j)<(C-tol)) )
                b = b2;
            else
                b = 1/2*(b1+b2);
            end
            
            num_changed_alphas = num_changed_alphas + 1;
        end
    end
    if(num_changed_alphas==0)
        passes = passes + 1;
    else
        passes = 0;
    end
end

w = (y.*alphas)'*X; w = w';