% 4.5 separating hyperplanes 
% written by Hongjing Lu at UCLA, 2/8/2008

clear all;
close all;
% input data and scatter plot
reddata = [2 22; 6 10; 12 8; 26 4; 30 2; 30 8; 22 26; 22 21; 27 17; 30 19];
greendata = [9 49; 26 40; 26 35; 30 25; 37 27; 42 21; 46 27; 43 48; 45 58; 65 58];
plot(reddata(:,1),reddata(:,2),'+r'); hold on;
plot(greendata(:,1),greendata(:,2),'og'); hold on;
x1range = 0:70; x2range=0:60;
xlim([x1range(1) x1range(end)]); ylim([x2range(1) x2range(end)]);
axis square;

% Y matrix: binary group membership (-1: reddata or 1: greendata)
y=[-ones(size(reddata,1),1); ones(size(greendata,1),1)];
% X matrix including intercept term
x=[reddata; greendata];
n=size(x,1); % # of data points

% lease square solution
[beta,bint,r,rint,stats] = regress(y,[ones(n,1) x]);
plot(x1range,-(beta(1)+beta(2)*x1range)/beta(3),'Color',[1 0.5 0.5],'LineWidth',2); hold on;

% Perceptron learning
rho=1; 
stepsize = 50000;
for simi=1:2
	beta2 = rand(3,1); % [beta0 beta1 beta2]
	for step=1:stepsize
	for i=1:n
        if sign([1 x(i,:)]*beta2)~=sign(y(i))  % find misclassified points
            beta2 = beta2+rho*y(i)*[1 x(i,:)]';  % update betas
        end;
	end;
	end;
	plot(x1range,-(beta2(1)+beta2(2)*x1range)/beta2(3),'--b','LineWidth',2); hold on;
    pause;
end;

% optimal separating hyperplane
% Construct the Kernel matrix
fprintf('Constructing ...\n');
H = zeros(n,n);  
for i=1:n
   for j=1:n
      H(i,j) = y(i)*y(j)*(x(i,:)*x(j,:)');
   end
end
c = -ones(n,1);  

% Add small amount of zero order regularisation to 
% avoid problems when Hessian is badly conditioned. 
H = H+10^(-10)*eye(size(H));

% Set up the parameters for the Optimisation problem
vlb = zeros(n,1);      % Set the bounds: alphas >= 0
vub = Inf*ones(n,1);     %                 alphas <= C
x0 = zeros(n,1);       % The starting point is [0 0 0   0]

% Solve the Optimisation Problem
fprintf('Optimising ...\n');
st = cputime;

[alpha]=quadprog(H,c,-eye(n,n),c*0,y',0,vlb,vub,x0);
w2 = alpha'*H*alpha;
fprintf('|w0|^2    : %f\n',w2);
fprintf('Margin    : %f\n',2/sqrt(w2));
fprintf('Sum alpha : %f\n',sum(alpha));

% Compute the number of Support Vectors
epsilon = 10^(-5);  % tolerance for Support Vector Detection
svi = find( alpha > epsilon);
nsv = length(svi);
fprintf('# of Support Vectors : %d (%3.1f%%)\n',nsv,100*nsv/n);

% compute betas
for i=1:size(x,2)
    beta3(i,1) = sum(y.*alpha.*x(:,i));
end;
count=0;
for i=1:n
    if abs(alpha(i))>epsilon
        count=count+1;
        beta3_0(count) = 1/y(i)-x(i,:)*beta3;
    end;
end;

% optimal hyperplane plot
figure; 
plot(reddata(:,1),reddata(:,2),'+r'); hold on;
plot(greendata(:,1),greendata(:,2),'og'); hold on;
x1range = 0:70; x2range=0:60;
xlim([x1range(1) x1range(end)]); ylim([x2range(1) x2range(end)]);
axis square;
plot(x1range,-(beta3_0(1)+beta3(1)*x1range)/beta3(2),'-b','LineWidth',2); hold on;  
sangle=atan(beta3(1)/beta3(2));
plot(x1range,-(beta3_0(1)+beta3(1)*x1range)/beta3(2)+1/sqrt(w2)/cos(sangle),'--b','LineWidth',1); hold on;
plot(x1range,-(beta3_0(1)+beta3(1)*x1range)/beta3(2)-1/sqrt(w2)/cos(sangle),'--b','LineWidth',1); hold on;