function [nsv, alpha, b0] = svc(X,Y,ker,C)
%SVC Support Vector Classification
%
%  Usage: [nsv alpha bias] = svc(X,Y,ker,C)
%
%  Parameters: X      - Training inputs
%              Y      - Training targets
%              ker    - kernel function
%              C      - upper bound (non-separable case)
%              nsv    - number of support vectors
%              alpha  - Lagrange Multipliers
%              b0     - bias term
%
%  Author: Steve Gunn (srg@ecs.soton.ac.uk)

  if (nargin <2 | nargin>4) % check correct number of arguments
    help svc
  else

    n = size(X,1);
    if (nargin<4) C=Inf;, end
    if (nargin<3) ker='linear';, end
    epsilon = 1e-10;
    
    % Construct the H matrix and c vector
    
    H = zeros(n,n);  
    for i=1:n
       for j=1:n
          H(i,j) = Y(i)*Y(j)*svkernel(ker,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. 

    if (abs(cond(H)) > 1e+10)
      fprintf('Hessian badly conditioned, regularising ....\n');
      fprintf('    Old condition number: %4.2g\n',cond(H));
      H = H+0.00000001*eye(size(H));
      fprintf('    New condition number: %4.2g\n',cond(H));
    end
    
    % 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 = [ ];              % The starting point is [0 0 0   0]
    neqcstr = nobias(ker); % Set the number of equality constraints (1 or 0)  
    if neqcstr
       A = Y';, b = 0;     % Set the constraint Ax = b
    else
       A = [];, b = [];
    end

    % Solve the Optimisation Problem
    
    st = cputime;
    
    if ( vlb == zeros(size(vlb)) & min(vub) == Inf & neqcstr == 0 )
       % Separable problem with Implicit Bias term
       % Use Non Negative Least Squares
       alpha = fnnls(H,-c);
    else
       % Otherwise
       % Use Quadratic Programming
       alpha = qp(H, c, A, b, vlb, vub, x0, neqcstr, -1);
    end

    fprintf('Execution time: %4.1f seconds\n',cputime - st);
    fprintf('|w0|^2    : %f\n',alpha'*H*alpha);  
    fprintf('Sum alpha : %f\n',sum(alpha));
    
        
    % Compute the number of Support Vectors
    svi = find( abs(alpha) > epsilon);
    nsv = length(svi);

    if neqcstr == 0
       % Implicit bias, b0
       b0 = 0;
    else
       % Explicit bias, b0; 
       % find b0 from pair of support vectors, one from each class
       classAsvi = find( abs(alpha) > epsilon & Y == 1);
       classBsvi = find( abs(alpha) > epsilon & Y == -1);
       nAsv = length( classAsvi );
       nBsv = length( classBsvi );
       if ( nAsv > 0 & nBsv > 0 )
          svpair = [classAsvi(1) classBsvi(1)];
          b0 = -(1/2)*sum(Y(svpair)'*H(svpair,svi)*alpha(svi));
       else 
          b0 = 0;
       end
    end
    
  end