% Demonstrate SVMs on a 2D two-class toy dataset, using Matlab's quadprog.

% Dataset: linearly separable by a little bit
rng(1778);
x1 = linspace(-1,1,100); x2 = 1 - x1.^2;
X = [x1(1:50);x2(1:50)]; Y = -ones(1,size(X,2));
X = [X [x1+1;0.3-x2]]'; Y = [Y ones(1,length(x1))]'; [N,D] = size(X);
X = X + randn(N,D)/10;
ax = [-1.4 2.3 -1.2 1.4];		% axes for 2D plots

% $$$ % Dataset: Gaussian classes with varying overlap
% $$$ t = 0.5;				% take t in [0 pi] to change angle
% $$$ gm.p = [2 1]'; gm.c = [0 0;5 0];
% $$$ gm.S = zeros(2,2,2); gm.S(:,:,1) = [1 0;0 1];
% $$$ R = [cos(t) sin(t);-sin(t) cos(t)]; gm.S(:,:,2) = R'*[9 0;0 1]*R;
% $$$ gm.type = 'F'; gm.p = gm.p/sum(gm.p);
% $$$ N = 200; [X,Y] = GMsample(N,gm); Y = 2*(1.5 - Y);
% $$$ ax = [-4 10 -5 5];			% axes for 2D plots

% The code below verifies that the primal QP and dual QP give the same result.
% We compare the solution (weight vector and bias, Lagrange multipliers, SVs)
% for problems "1" and "2" by using max(abs("1"-"2")), if the result is a
% small number (< 1e-6, say) then the solutions are (presumably) the same.
% Also, observe the return code e from Matlab's quadprog: 1 means the solution
% was found; -2 means the QP was infeasible and there is no solution.

% - Optimal separating hyperplane:
C = NaN;
[f1,e1,sv1,l1] = linsvmtrain(X,Y,C,0); e1		% primal QP
[f2,e2,sv2,l2] = linsvmtrain(X,Y,C,1); e2		% dual QP
[f3,e3,sv3,l3] = linsvmtrain(X(sv1,:),Y(sv1),C,0); e3	% primal on SVs only
% primal = dual?
max(abs([f1.w;f1.w0]-[f2.w;f2.w0])), max(abs(l1-l2)), max(abs(sv1-sv2))
% primal on X = primal on SVs?
max(abs([f1.w;f1.w0]-[f3.w;f3.w0])), max(abs(l1(sv1)-l3))

% - Soft margin hyperplane (take any C > 0):
C = 1;
[f1,e1,sv1,l1,m1,xi1] = linsvmtrain(X,Y,C,0); e1	% primal QP
[f2,e2,sv2,l2,m2,xi2] = linsvmtrain(X,Y,C,1); e2	% dual QP
% primal = dual?
max(abs([f1.w;f1.w0]-[f2.w;f2.w0])), max(abs(l1-l2)), max(abs(sv1-sv2)), max(abs(xi1-xi2))


% ---------------------------------- Plots ----------------------------------

% Dataset
figure(1); clf; hold on;
set(gcf,'Name','1D dataset');
J = find(Y<0); plot(X(J,1),X(J,2),'m+');
J = find(Y>0); plot(X(J,1),X(J,2),'c+');
hold off; axis(ax); daspect([1 1 1]); box on; xlabel('x_1'); ylabel('x_2');

% Linear SVM
allC = [NaN 10.^(-4:2)];
figure(2);
set(gcf,'Name','1D dataset: linear SVMs for different hyperparameter values');
for k=1:length(allC)
  [f,e,sv,l,m] = linsvmtrain(X,Y,allC(k),0); w = [f.w;f.w0];
  fX = linsvm(X,f); fX = 2*(fX > 0) - 1;		% predicted labels
  err = sum(fX~=Y)/N;					% classification error
  set(0,'CurrentFigure',2); clf; hold on;
  fill([ax(1:2) fliplr(ax(1:2))],...
       [-(w(1)*ax(1:2)+w(3)+1)/w(2) fliplr(-(w(1)*ax(1:2)+w(3)-1)/w(2))],...
       [1 1 1]*.8);
  J = find(fX<0); plot(X(J,1),X(J,2),'r+');
  J = find(fX>0); plot(X(J,1),X(J,2),'b+');
  plot(ax(1:2),-(w(1)*ax(1:2)+w(3))/w(2),'k-','LineWidth',2);
  plot(X(sv,1),X(sv,2),'ko','MarkerSize',12);
  hold off; axis(ax); daspect([1 1 1]); box on; xlabel('x_1'); ylabel('x_2');
  if ~isnan(allC(k))
    msg = ['Linear SVM (soft margin hyperplane), C = ' num2str(allC(k)) ':'];
  else
    if e==-2 msg1 = 'NOT '; else msg1 = []; end
    msg = ['Linear SVM (optimal separating hyperplane): dataset ' msg1 ...
           'linearly separable'];
  end
  title({msg,...
         ['margin 1/|w| = ' num2str(m) ',   #SVs = ' num2str(length(sv)) ...
          ',   training error: ' num2str(100*err,2) '%']});
  pause
end


% Suggestions of things to try:
% - Different C.
% - Change overlap in the Gaussian classes dataset.
% - Different datasets.

