% APRENDIZAJE AUTOMÁTICO (G. Pajares y J. M. de la Cruz, Eds.)
% Editorial: RA-MA
% Capítulo 5: Máquinas de Vectores Soporte
% Autores: Jesús Antonio Vega Sánchez y Sebastián Dormido Canto

function [svm,psv,b0,nsv]=wsvm(XX,clase,kernel,c,param);

% Función que realiza el entrenamiento de la máquina k.
%
% Los parámetros de entrada son:
%                 - XX: es la matriz de los vectores de entrenamiento (se genera con el programa norm_svm). Ver 
%                       norm_svm para observar el formato.
%
%                 - clase: (vector columna) se almacenan las etiquetas {1,-1} que identifican si los vectores de 
%                          entrenamiento pertenecen a la clase considerada {1} ó no (es decir, pertenecen al resto) {-1}
%
%                 - kernel: es un string donde se introduce la función de transformación a utilizar.   
%
%                 - c: almacena el valor del límite de los alphas de la máquina k 
%
%                 - param: es un vector de dimensión 2 que almacena los parámetros de las funciones kernel que
%                          se pueden utilizar. Por ejemplo param=[2 1]; 
%
% Los parámetros de salida son:
%                - svm: almacena los vectores soporte de la máquina k (se almacenan en filas), la penúltima
%                       columna almacena el valor de la etiqueta {1,-1} para saber si pertenece a la clase k ó
%                       al resto y en la última columna el valor de alpha
%
%                - psv: es una matriz de dos columnas, en la primera se almacenan los índices absolutos de los
%                       vectores soporte y en la segunda se almacena una etiqueta {1,-1} que identifica si el 
%                       vector soporte corresponde a la clase k {1} ó al resto {-1}
%
%                - b0: es el término bias de la máquina k
%
%                - nsv: número de vectores soporte de la máquina k


% Se inicializan variables
svm=[];
psv=[];


n = size(XX,1);

% tolerance for Support Vector Detection
%epsilon = svtol2(C);
epsilon = 1e-8;   %OJO con este valor 
    
% Construct the Kernel matrix
%fprintf('Constructing ...\n');
H = zeros(n,n);  
for i=1:n
    for j=1:n
       H(i,j) = clase(i)*clase(j)*svkernel2(kernel,XX(i,:),XX(j,:),param);
    end
end
cc = -ones(n,1);  

% Add small amount of zero order regularisation to 
% avoid problems when Hessian is badly conditioned. 
H = H+1e-10*eye(size(H));
    
% Set up the parameters for the Optimisation problem
vlb = zeros(n,1);      % Set the bounds: alphas >= 0
vub = c*ones(n,1);     %                 alphas <= C
x0 = zeros(n,1);       % The starting point is [0 0 0   0]
neqcstr = nobias2(kernel); % Set the number of equality constraints (1 or 0)  
if neqcstr
    A = clase';, b = 0;     % Set the constraint Ax = b
   else
    A = [];, b = [];
end

% Solve the Optimisation Problem
%fprintf('Optimising ...\n');
st = cputime;
    
%AJUSTE PARA EL MATLAB 2008b
%[alpha lambda how] = qp(H, cc, A, b, vlb, vub, x0, neqcstr);
OPTIONS = optimset('LargeScale','off','Display','off');
[alpha lambda how] = quadprog(H, cc, A, b, A, b, vlb, vub, x0, OPTIONS);
    
%fprintf('Matriz de alphas\n')
   
%fprintf('Execution time: %4.1f seconds\n',cputime - st);
%fprintf('Status : %s\n',how);
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
alphamax=max(alpha);

j=0;
for i=1:length(clase)
    if alpha(i)/alphamax > epsilon
       j=j+1;
       svm(j,:)=[XX(i,:),clase(i),alpha(i)];
       psv(j,:)=[i clase(i)];
    end
end
svi=psv(:,1);

 
%svi = find( alpha > epsilon);
nsv = length(svi);
%fprintf('Support Vectors : %d (%3.1f%%)\n',nsv,100*nsv/n);

% Implicit bias, b0
b0 = 0;

% Explicit bias, b0 
if nobias2(kernel) ~= 0
   % find b0 from average of support vectors on margin
   % SVs on margin have alphas: 0 < alpha < C
   svii = find( alpha > epsilon & alpha < (c - epsilon));
   if length(svii) > 0
     b0 =  (1/length(svii))*sum(clase(svii) - H(svii,svi)*alpha(svi).*clase(svii));
    else 
    % fprintf('No support vectors on margin - cannot compute bias.\n');
    end
end

%Cálculo del w0 de arriba utilizando la fórmula 9.22 de las fotocopias
%El cálculo está bien realizado y coincide con el valor de b0
 %I=find(alpha > epsilon);
 %nalpha=alpha(I) 
 
 nX=svm(:,1:end-2);
 nY=svm(:,end-1);


 nalpha=svm(:,end);
 
 nalphanY=nalpha.*nY;
 
 nnX = size(nX,1);
 nH = zeros(nnX,1);  
 
    for i=1:nnX
       
          nH(i,1) = svkernel2(kernel,nX(i,:),nX(1,:),param);
    end
   
   nHaux=nH.*nalphanY;
   b0=nY(1)-sum(nHaux);



return
 















