%NU_SVRO Support Vector Optimizer
%
%   [V,J] = NU_SVRO(K,Y,C)
%
% INPUT
%   K     Similarity matrix
%   NLAB  Label list consisting of -1/+1
%   C     Scalar for weighting the errors (optional; default: 10)
%
% OUTPUT
%   V     Vector of weights for the support vectors
%   J     Index vector pointing to the support vectors
%
% DESCRIPTION
% A low level routine that optimizes the set of support vectors for a 2-class
% classification problem based on the similarity matrix K computed from the
% training set. SVO is called directly from SVC. The labels NLAB should indicate 
% the two classes by +1 and -1. Optimization is done by a quadratic programming. 
% If available, the QLD function is used, otherwise an appropriate Matlab routine.
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% NU_SVR, SVO, SVC

% Revisions:
% DR1, 07-05-2003
%      Sign error in calculation of offset

% Copyright: D.M.J. Tax, D. de Ridder, R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: nu_svro.m,v 1.2 2010/02/08 15:29:48 duin Exp $

function [v,J,epsilon_or_nu] = nu_svro(K,y,C,svr_type,nu_or_epsilon,pd,abort_on_error)

		if (nargin < 7) | isempty(abort_on_error)
    abort_on_error = 0;
  end

	if (nargin < 6) | isempty(pd)
    pd = 1;
  end

	if (nargin < 5) 
    nu_eps = [];
  end

	if (nargin < 4)  | isempty(svr_type)
    svr_type = 'nu';
  end

  switch svr_type
  case 'nu'
    if isempty(nu_or_epsilon)
      prwarning(3,'nu is not specified, assuming 0.25.');
      nu_or_epsilon = 0.25;
    end
    nu = nu_or_epsilon;
  case {'eps', 'epsilon'}
    svr_type = 'epsilon';
    if isempty(nu_or_epsilon)
      prwarning(3,'epsilon is not specified, assuming 1e-2.');
      nu_or_epsilon = 1e-2;
    end
    epsilon = nu_or_epsilon;
  end

	if (nargin < 3)
		prwarning(3,'C is not specified, assuming 1.');
    C = 1;
	end
	vmin = C*1e-9;		% Accuracy to determine when an object becomes the support object.

	% Set up the variables for the optimization.
	n  = size(K,1);
  D  = K;	
  switch svr_type
  case 'nu'
    f  = [-y', y'];  
    A  = [[ones(1,n), -ones(1,n)]; ones(1,2*n)];
    b  = [0; C*nu*n];
  case 'epsilon'
    f   = epsilon + [-y', y'];
  	A   = [ones(1,n), -ones(1,n)];
    b   = 0;
  end  
  
  lb = zeros(2*n,1);
	ub = repmat(C,2*n,1);
	p  = rand(2*n,1);

	if pd
    % Make the kernel matrix K positive definite.
	  i = -30;
  	while (pd_check (D + (10.0^i) * eye(n)) == 0)
	  	i = i + 1;
  	end
  	if (i > -30),
	  	prwarning(2,'K is not positive definite. The diagonal is regularized by 10.0^(%d)*I',i);
  	end
  	i = i+2;
  	D = D + (10.0^(i)) * eye(n);
  end
  D  = [[D, -D]; [-D, D]];	

	% Minimization procedure:
	% minimizes:  0.5 x' D x + f' x
	% subject to: Ax = b
	%
  if (exist('qld') == 3)
		v = qld (D, f, -A, b, lb, ub, p, length(b));
	elseif (exist('quadprog') == 2)
		prwarning(1,'QLD not found, the Matlab routine QUADPROG is used instead.')
		opt = optimset; opt.LargeScale='off'; opt.Display='off';
		v = quadprog(D, f, [], [], A, b, lb, ub,[],opt);
	else 
		prwarning(1,'QLD not found, the Matlab routine QP is used instead.')
		verbos = 0;
		negdef = 0;
		normalize = 1;
		v = qp(D, f, A, b, lb, ub, p, length(b), verbos, negdef, normalize);
	end
  % Find all the support vectors.
  
  if isempty(v) 
    ErrMsg = 'Quadratic Optimization failed.';
    [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
    return;
  end

  v  = v(1:n)-v((n+1):end);
  av = abs(v);
  J  = find(av > vmin);	    % sv's
  I  = J(av(J) < (C-vmin)); % on-tube sv's
  
  %plot(v,y-K(:,J)*v(J),'.')

  switch svr_type
  case 'nu'
    Ip = I(v(I) > 0);
    Im = I(v(I) < 0);  
    if isempty(Ip) | isempty(Im);
      prwarning(2,'epsilon and b are not unique: values from the admissible range will be taken');
      J0 = find(av <= vmin); % non-sv's
      if ~isempty(J0)
        y_minus_wx_0 = y(J0)-K(J0,J)*v(J); 
      else
        prwarning(2,'There are no non-sv''s: admissible (eps,b) range is partially unbounded');
      end
    end
    
    if ~isempty(Ip)
      b_plus_epsilon = mean(y(Ip)-K(Ip,J)*v(J),1); 
    else
      lp_bound = -inf;
      up_bound =  inf;

      if ~isempty(J0)
        lp_bound = max(y_minus_wx_0,[],1); 
      end
    
      Ipo = J((av(J) >= C-vmin) & v(J) > 0); % positive out of tube sv's
      if ~isempty(Ipo)
        up_bound = min(y(Ipo)-K(Ipo,J)*v(J),[],1); 
      end  
    
      if isinf(up_bound)
        Msg = 'Impossible situation: there are no positive sv''s.';    
        [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
        return;       
      elseif isinf(lp_bound)
        b_plus_epsilon = up_bound;
      else
        if lp_bound > up_bound
          ErrMsg = 'Impossible situation: admissible (eps,b) region is empty.';    
          [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
          return;       
        end
        b_plus_epsilon = 0.5*(lp_bound+up_bound);
      end  
    end
        
    if ~isempty(Im)
      b_minus_epsilon = mean(y(Im)-K(Im,J)*v(J),1); 
    else
      lm_bound = -inf;
      um_bound =  inf;
      
      Imo = J((av(J) >= C-vmin) & v(J) < 0); % positive out of tube sv's
      if ~isempty(Imo)
        lm_bound = max(y(Imo)-K(Imo,J)*v(J),[],1); 
      end  
       
      if ~isempty(J0)
        um_bound = min(y_minus_wx_0,[],1); 
      end
      
      if isinf(lm_bound)
        ErrMsg = 'Impossible situation: there are no negative sv''s.';    
        [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
        return;
      elseif isinf(um_bound)
        b_minus_epsilon = lm_bound;
      else
        if lm_bound > um_bound
          ErrMsg = 'Impossible situation: admissible (eps,b) range is empty.';    
          [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
          return;       
        end
        b_minus_epsilon = 0.5*(lm_bound+um_bound);
      end  
    end
    
    % one more paranoic check
    if exist('J0') == 1 & ~isempty(J0)
      ok = 1;
      if isempty(Ip) & ~isempty(Im) 
        ok = b_minus_epsilon <= min(y_minus_wx_0,[],1);
      elseif ~isempty(Ip) & isempty(Im)
        ok = b_plus_epsilon >= max(y_minus_wx_0,[],1);      
      end
      
      if ~ok
        ErrMsg = 'Impossible situation: incosistance in admissible (eps,b) region.';    
        [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
      end  
    end
    
    epsilon = 0.5*(b_plus_epsilon-b_minus_epsilon);
    b = 0.5*(b_plus_epsilon+b_minus_epsilon);

    epsilon_or_nu = epsilon;

  case 'epsilon'
    if ~isempty(I)      
      b = mean(y(I)-K(I,J)*v(J)-epsilon*sign(v(I)));
    else 
      prwarning(2,'b is not unique: value from the admissible range will be taken');
      
      lp_bound = -inf;
      up_bound =  inf;

      lm_bound = -inf;
      um_bound =  inf;

      J0 = find(av <= vmin); % non-sv's
      if ~isempty(J0)
        y_minus_wx_0 = y(J0)-K(J0,J)*v(J); 
        lp_bound = max(y_minus_wx_0,[],1)-epsilon; 
        um_bound = min(y_minus_wx_0,[],1)+epsilon; 
      else
        prwarning(2,'Thers are no non-sv''s');
      end
    
      Ipo = J((av(J) >= C-vmin) & v(J) > 0); % positive out of tube sv's
      if ~isempty(Ipo)
        up_bound = min(y(Ipo)-K(Ipo,J)*v(J),[],1)-epsilon; 
      end  
    
      Imo = J((av(J) >= C-vmin) & v(J) < 0); % negative out of tube sv's
      if ~isempty(Imo)
        lm_bound = max(y(Imo)-K(Imo,J)*v(J),[],1)+epsilon; 
      end  

      l_bound = max(lm_bound,lp_bound);
      u_bound = min(um_bound,up_bound);
      
      ErrMsg = '';
      if isinf(up_bound)
        ErrMsg = 'Impossible situation: there are no positive sv''s.';    
      elseif isinf(lm_bound)
        ErrMsg = 'Impossible situation: there are no negative sv''s.';    
      elseif l_bound > u_bound
        keyboard
        ErrMsg = 'Impossible situation: admissible b region is empty.';    
      end
      
      if ~isempty(ErrMsg)  
        [v,J,epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error);
        return;
      end         
      
      b = 0.5*(l_bound+u_bound); 
    end

    nu = sum(av(J))/(C*n);
    epsilon_or_nu = nu;
    
  end  
	
  v = [v(J); b];

return;

function [v,J, epsilon_or_nu] = ErrHandler(K,y,ErrMsg,abort_on_error)
  if abort_on_error
    error(ErrMsg);
  else
    prwarning(1,[ErrMsg 'Pseudoinverse Regression is computed instead.']);
    n = size(K,1);
    v = prpinv([K ones(n,1)])*y;
    J = [1:n]';
	  epsilon_or_nu = nan;		  
  end  
return