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

  if (nargin <3 | nargin>5) % check correct number of arguments
    help svr
  else

    n = size(X,1);
    if (nargin<5) C=Inf;, end
    if (nargin<4) e=0.05;, end
    if (nargin<3) ker='linear';, end  
    epsilon = 1e-10; % tolerance for Support Vector Detection
    
    % Construct the H matrix and c vector
    
    H = zeros(n,n);  
    for i=1:n
       for j=1:n
          H(i,j) = svkernel(ker,X(i,:),X(j,:));
       end
    end
    Hb = [H -H; -H H];
    c = [(e*ones(n,1) - Y); (e*ones(n,1) + Y)];  

    % Add small amount of zero order regularisation to 
    % avoid problems when Hessian is badly conditioned. 
    % Rank is always less than or equal to n.
    % Note that adding to much reg will peturb solution

    if (abs(cond(Hb)) > 1e+10)
      fprintf('Hessian badly conditioned, regularising ....\n');
      fprintf('    Old condition number: %4.2g\n',cond(Hb));
      Hb = Hb+0.0000000000001*eye(size(Hb));
      fprintf('    New condition number: %4.2g\n',cond(Hb));
    end

    % Set up the parameters for the Optimisation problem

    vlb = zeros(2*n,1);    % Set the bounds: alphas >= 0
    vub = C*ones(2*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 = [ones(1,n) -ones(1,n)];, 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(Hb,-c);
    else
       % Otherwise
       % Use Quadratic Programming
       alpha = qp(Hb, c, A, b, vlb, vub, x0, neqcstr, -1);
    end

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

    if neqcstr == 0
       % Implicit bias, b0
       b0 = 0;
    else
       % Explicit bias, b0; 
       % compute using robust method of Smola
       % find b0 from average of support vectors with interpolation error e
       svbi = find( abs(beta) > epsilon & abs(beta) < C );
       nsvb = length(svbi);
       if nsvb > 0
          b0 = (1/nsvb)*sum(Y(svbi) + e*sign(beta(svbi)) + H(svbi,svi)*beta(svi));
       else 
          b0 = (max(Y)+min(Y))/2;
       end
    end

  end