%SVO_NU Support Vector Optimizer: NU algorithm
%
%   [V,J,C] = SVO(K,NLAB,NU,PD)
%
% INPUT
%   K     Similarity matrix
%   NLAB  Label list consisting of -1/+1
%   NU    Regularization parameter (0 < NU < 1): expected fraction of SV (optional; default: 0.25)
%
%   PD    Do or do not the check of the positive definiteness (optional; default: 1 (to do))
%
% OUTPUT
%   V     Vector of weights for the support vectors
%   J     Index vector pointing to the support vectors
%   C     Equivalent C regularization parameter of SVM-C algorithm
%
% 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.
%
% NU is bounded from above by NU_MAX = (1 - ABS(Lp-Lm)/(Lp+Lm)), where
% Lp (Lm) is the number of positive (negative) samples. If NU > NU_MAX is supplied 
% to the routine it will be changed to the NU_MAX.
%
% If NU is less than some NU_MIN which depends on the overlap between classes 
% algorithm will typically take long time to converge (if at all). 
% So, it is advisable to set NU larger than expected overlap.
%
% Weights V are rescaled in a such manner as if they were returned by SVO with the parameter C.
%
% SEE ALSO
% SVC_NU, SVO, SVC

% Copyright: S.Verzakov, s.verzakov@ewi.tudelft.nl 
% Based on SVO.M by D.M.J. Tax, D. de Ridder, R.P.W. Duin
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

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

function [v,J,C] = svo_nu(K,y,nu,pd)
	  if (nargin < 4) 
    pd = 1;
	end

	if (nargin < 3)
		prwarning(3,'Third parameter (nu) not specified, assuming 0.25.');
    nu = 0.25;
	end

  nu_max = 1 - abs(nnz(y == 1) - nnz(y == -1))/length(y);

  if nu > nu_max
    prwarning(3,['nu==' num2str(nu)  ' is not feasible; set to ' num2str(nu_max)]);
    nu = nu_max;
  end 

  vmin = 1e-9;		% Accuracy to determine when an object becomes the support object.

	% Set up the variables for the optimization.
	n  = size(K,1);
	D  = (y*y').*K;
	f  = zeros(1,n);
	A  = [y';ones(1,n)];
	b  = [0; nu*n];
	lb = zeros(n,1);
	ub = ones(n,1);
	p  = rand(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
   
  % Minimization procedure initialization:
	% 'qp' 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.')
		v = quadprog(D, f, [], [], A, b, lb, ub);
	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.
	J = find(v > vmin);

	% First find the SV on the boundary
	I  = J(v(J) < 1-vmin);
	Ip = I(y(I) ==  1);
	Im = I(y(I) == -1);
	if (isempty(v) | isempty(Ip) | isempty(Im))
		%error('Quadratic Optimization failed. Pseudo-Fisher is computed instead.');
		prwarning(1,'Quadratic Optimization failed. Pseudo-Fisher is computed instead.');
		v = prpinv([K ones(n,1)])*y;
		J = [1:n]';
		C = nan;
    return;
	end

  v = y.*v;

  %wxI = K(I,J)*v(J);
  wxIp = mean(K(Ip,J)*v(J),1); %  rho-b
  wxIm = mean(K(Im,J)*v(J),1); % -rho-b

  rho =  0.5*(wxIp-wxIm);
	b   = -0.5*(wxIp+wxIm);
  %b = mean(rho*y(I) - wxI);

	v = [v(J); b]/rho;
  C = 1/rho;   

  return;