% Copyright (c) 2016 by Miguel A. Carreira-Perpinan
% for use in CSE176 Introduction to Machine Learning at UC Merced

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

% Dataset: two classes, linearly separable
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);
v = randperm(N,50); Xv = X(v,:); Yv = Y(v,:);
X = X + randn(N,D)/10; Xv = Xv + randn(size(Xv))/10;

% Grid of points for contour plots
lx1 = 100; lx2 = 100; ax = [-1.4 2.3 -1.2 1.4];
x1 = linspace(ax(1),ax(2),lx1); x2 = linspace(ax(3),ax(4),lx2);
[XX,YY] = meshgrid(x1,x2);

% Linear SVM, optimal separating hyperplane:
C = Inf;
[f,~,e,sv,l,m] = svm2train(X,Y,{'lin'},C,[],0); e	% primal QP
w = [f.w;f.w0];
[f2,~,e,sv2,l2] = svm2train(X,Y,{'lin'},C,[],1); e	% dual QP
[f3,~,e,sv3,l3] = svm2train(X(sv,:),Y(sv),{'lin'},C,[],0); % primal on SVs only
max(abs(w-[f2.w;f2.w0])),max(abs(l-l2)),max(abs(sv-sv2)) % primal = dual?
max(abs(w-[f3.w;f3.w0])),max(abs(l(sv)-l3))		% primal on X = on SVs?

% Linear SVM, soft margin hyperplane:
C = 1;
[f,~,e,sv,l,m,xi] = svm2train(X,Y,{'lin'},C,[],0); e	% primal QP
w = [f.w;f.w0];
[f2,~,e,sv2,l2,m2,xi2] = svm2train(X,Y,{'lin'},C,[],1); e % dual QP
max(abs(w-[f2.w;f2.w0])),max(abs(l-l2)),max(abs(sv-sv2)),max(abs(xi-xi2)) % primal = dual?

% Kernel SVM, soft margin hyperplane:
C = 1;
s = 0.5; [f,~,e,sv,l,m,xi] = svm2train(X,Y,{'rbf',s},C,[]); e
q = 1; [f,~,e,sv,l,m,xi] = svm2train(X,Y,{'poly',q},C,[]); e


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

% Dataset
figure(1); clf; hold on;
J = find(Y<0); plot(X(J,1),X(J,2),'r+');
J = find(Y>0); plot(X(J,1),X(J,2),'b+');
hold off; axis(ax); daspect([1 1 1]); box on; xlabel('x_1'); ylabel('x_2');

% Linear SVM
allC = [10.^(-2:1) Inf];
figure(2);
for k=1:length(allC)
  [f,fX,e,sv,l,m] = svm2train(X,Y,{'lin'},allC(k)); w = [f.w;f.w0];
  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);
  %plot(ax(1:2),-(w(1)*ax(1:2)+w(3)+1)/w(2),'k:',...
  %     ax(1:2),-(w(1)*ax(1:2)+w(3)-1)/w(2),'k:');
  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');
  title(['ker = ' f.ker ', C = ' num2str(allC(k)) ...
         ', margin 1/|w| = ' num2str(m) ', #SVs = ' num2str(length(sv))]);
  pause
end

% Kernel SVM
allK = {'poly','rbf'};					% Kernels
allqs = {[1 2 3 5 10],...				% q
         [0.05 0.1 0.2 0.5 5]};				% s
allC = 10.^(-2:1);					% C
for i=1:length(allK)
  figure(i+2);
  for j=1:length(allqs{i})
    for k=1:length(allC)
      [f,fX,e,sv,l,m,xi] = svm2train(X,Y,{allK{i},allqs{i}(j)},allC(k));
      fXY = reshape(svm2([XX(:) YY(:)],f,0),lx1,lx2);
      set(0,'CurrentFigure',i+2); clf; hold on;
      contour(x1,x2,fXY);
      contour(x1,x2,fXY,[0 0],'LineColor','k','LineWidth',2);
      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(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');
      switch f.ker
       case 'rbf', kerp = ['(\sigma=' num2str(f.s) ')'];
       case 'poly', kerp = ['(q=' num2str(f.q) ')'];
      end
      title(['ker = ' f.ker kerp ', C = ' num2str(allC(k)) ...
             ', margin 1/|w| = ' num2str(m) ', #SVs = ' num2str(length(sv))]);
      pause
    end
  end
end


% Suggestions of things to try:
% - Different kernels: 'lin', 'poly', 'rbf'
% - Different hyperparameters for each kernel: q (poly), s (rbf)
% - Different C
% - Different datasets

