% svmdemo.m -- demonstration of svm algorithm for classification
% Utilizing libsvm C programs  svmtrain.exe
% call datasepf.m, datagen.m, mat2libsvm.m, readmodel.m
% (c) copyright 2000 by Yu Hen Hu

clear all; clf;
nfeature=2; nclass=2;  % 2D feature space, 2 classes
eymat=eye(nclass);

disp('Enter 0 (default) to generate linearly separable data, or ...');
disp('enter 1 to generate non-separable data');
chos=input('enter 2 to load from a data file: ')
if isempty(chos), chos=0; end

if chos==0,
   n=input('Total number of samples to be generated = ');
   n1=ceil(n/2); n0=n-n1;
   % generating linearly separable or non-separable 2D data files
   disp('generate linearly separable data files');
   pn=1; % generate +1 -1 label
   % data is n x 3 with last col being target +1, -1
   [data,slope]=datasepf(n,n1,pn);  
   % convert into 1 0, 0 1 target format
   data=[data(:,1:2) [data(:,3)==-1] [data(:,3)==1]]; 
   C=1e9;
elseif chos==1,
   n=input('Total number of samples to be generated = ');
   n1=ceil(n/2); n0=n-n1;
   Nvec=[n0 n1];
   mean_var=[-0.5 0 0.2; 0.5 0 0.2]';
   data=datagen(Nvec,mean_var,1); % fr is n x 5, training data
   % convert into 1 0, 0 1 target format
   data=[data(:,1:2) [data(:,3)==0] [data(:,3)==1]]; 
   C=1000;
elseif chos==3,
   fname=input('Enter filename, default extension .txt, within single quotes: ');
   eval(['load ' fname '.txt']); % variable name is fname
   data=eval(fname); eval(['clear ' fname]);
   [n,n2]=size(data);
   C=1000;
end
mat2libsvm(data,2,'data')

disp('Choose kernel type:')
disp('0. Linear kernel;')
disp('1. Polynomial kernel (defult);')
disp('2. RBF kernel;')
disp('3. Sigmoid kernel;')
ktype=input('Enter kernel type: ');
if isempty(ktype), ktype=1; end

switch ktype
case 0, % kernel='linear';
   para=['-t 0 '];
case 1, %kernel='poly';
   degree=input('Enter order of polynomial (default = 2) = ');
   if isempty(degree), degree=2; end
   % set parameters, K = (u'v + 1)^degree
   para=['-t 1 -d ' int2str(degree) ' -g 1 -r 1 '];
case 2, % kernel='rbf';
   para=['-t 2 -g 1 '];
case 3, % kernel = 'sigmoid';
   para=['-t 3 -g 1 -r 0 '];
end

% call svmtrain.exe in dos using shell escape!
% the model file is data.model
eval(['!svmtrain -s 0 ' para ' -c ' int2str(C) ' data']);
[sv,model]=readmodel('data.model',2); % read support vectors from data.model
alpha=sv(:,1); % sv(:,1): alpha value model.nsv X 1
svs=sv(:,2:nfeature+1); % support vectors  model.nsv X 2
label= sv(:,nfeature+2:nfeature+1+nclass)*[-1;+1]; %target values, model.nsv X 1

if ktype==0,
   disp('the weights are:'); % w = sum(alpha(i)*d(i)*x(i))
   w=svs'*(alpha.*label)  % w is 2 by 1
   % b=sum(1-w'*(d(i).*x(i)))/sum(d(i)) over all support vectors
   % b=(model.nsv-sum(diag(label)*svs*w))/sum(label)
   b=(sum(label)-sum(svs*w))/model.nsv
   disp(['the margin (2*rho) = ' num2str(2/norm(w))]);
end

figure(1)
plot(data(find(data(:,3)==1),1),data(find(data(:,3)==1),2),'g.',...
     data(find(data(:,4)==1),1),data(find(data(:,4)==1),2),'b.',...
     svs(find(alpha==C),1),svs(find(alpha==C),2),'ro',...
     svs(find(alpha<C),1),svs(find(alpha<C),2),'rs')
legend('Class 1 data','Class 2 data','SV, alpha=C','SV, on margin')
axis equal
title(['SVM results, C = ' int2str(C)])