function [mu alpha_sv b pos history obj status] = RSVMPlusl1_ls(X, Y,componentDistances, C, C_r,option)
%% RSVM^+_{\mu} - presented in our paper.
%% large scale implementation, kernel is computed on the fly.
%  RSVMPlusl1 - the plus form of R-SVM with an explicit scaling factor Mu (with constrained mu \sum mu_k = 1),  use L_1 norm on slack variable (L_1 SVM)
% The radius is approximated by linear programming, i.e the maximum pair distance between instances. 
%% Params:
% Inputs: 
% X : data, rows are instances, columns are features, size [NxD]
% Y : label, size [Nx1]
% basicKernels: NxNxD matrix of  basic kernel, last dimension is the kernel index
% C         : SVM hyperparameter
% option    : RSVMPlusl1 algorithm hyperparameter
%       option.nbitermax : maximal number of iterations (default 1000)
%       option.seuil     : threshold for zeroing kernel coefficients
%                          (default 1e-12) during the first
%                          option.seuilitermax
%       option.seuilitermax : threshold weigth when the nb of iteration is
%                           lower than this value
%       option.sigmainit  : initial kernel coefficient (default average)
%       option.alphainit : initial Lagrangian coefficient
%       option.lambdareg : ridge regularization added to kernel's diagonal
%                          for SVM (default 1e-10)
%       option.numericalprecision    : force to 0 weigth lower than this
%                                     value (default = eps)
%
% Outputs:
% mu            : the weights of kernel.
% alpha_sv      : the weigthed lagrangian of the support vectors: alpha_sv = w =(alpha(supportvector).*ysupport);
% b             : the bias (b)
% pos           : the indices of Support  vectors
%
%% Copyright@Huyen Do 2012
% Some parts are inherited from SimpleMKL toolbox.

%% Initialize params
[n] = length(Y);

numFeatures = size(X,2);
numKernels = numFeatures;

if ~isfield(option,'nbitermax');
    nloopmax = 1200;
else
    nloopmax = option.nbitermax;
end;
if ~isfield(option,'seuil');
    seuil = 1e-12;
else
    seuil = option.seuil;
end
if ~isfield(option,'seuilitermax')
    option.seuilitermax = 20;
end;
if ~isfield(option,'seuildiffsigma');
    option.seuildiffsigma = 1e-5;
end
if ~isfield(option,'seuildiffconstraint');
    option.seuildiffconstraint = 0.05;
end

if ~isfield(option,'lambdareg');
    lambdareg = 1e-10;
    option.lambdareg = 1e-10;
else
    lambdareg = option.lambdareg;
end


if ~isfield(option,'numericalprecision');
    option.numericalprecision = 1e-5;
end;


%------------------------------------------------------
% choosing the stopping criterion
%------------------------------------------------------
option.stopvariation    = 1; % use variation of weights for stopping criterion
option.stopKKT          = 0; % set to 1 if you use KKTcondition for stopping criterion [not yet implemented]
option.stopdualitygap   = 1; % set to 1 for using duality gap for stopping criterion   [not yet implemented]

%------------------------------------------------------
% choosing the stopping criterion value
%------------------------------------------------------
option.seuildiffsigma       = 1e-2/numKernels;  % stopping criterion for weight variation
option.seuildiffconstraint  = 0.1;              % stopping criterion for KKT
option.seuildualitygap      = 0.01;              % stopping criterion for duality gap
%------------------------------------------------------
% some algorithms paramaters
%------------------------------------------------------
option.firstbasevariable =  'first'; % tie breaking method for choosing the base
% variable in the reduced gradient method
option.nbitermax    = 500;             % maximal number of iteration
option.seuil        = 0;               % forcing to zero weights lower than this
option.seuilitermax = 10;              % value, for iterations lower than this one

option.miniter      = 0;                 % minimal number of iterations
option.verbosesvm   = 0;              % verbosity of inner svm algorithm
option.efficientkernel      = 1;         % use efficient storage of kernels

verbose = 0;

%---------------------------------------------------------
%Initialize mu

if ~isfield(option,'muInit');
   % mu = 1:numKernels;
   % mu = mu'/sum(mu);
   %mu = ones(numKernels,1)/numKernels;
   mu = rand(numKernels,1);
   mu = mu/sum(mu);
else 
    mu = option.muInit ;
end;
%End Initialize mu
% -----------------------------------------------------------

if isfield(option,'alphainit');
    alphainit = option.alphainit;
else
    alphainit = [];
end;


%--------------------------------------------------------------------------------
% Options used in subroutines
%--------------------------------------------------------------------------------
if ~isfield(option,'goldensearch_deltmax');
    option.goldensearch_deltmax = 1e-1;
end
if ~isfield(option,'goldensearchmax');
    option.goldensearchmax = 1e-8;
end;
if ~isfield(option,'firstbasevariable');
    option.firstbasevariable = 'first';
end;

%------------------------------------------------------------------------------%
% Initialize
%--------------------------------------------------------------------------

nloop   = 0;
loop    = 1;
status  = 0;
numericalaccuracy        = 1e-9;
goldensearch_deltmaxinit = option.goldensearch_deltmax;


% End Initialize parameters
%% -----------------------------------------------------------


%-----------------------------------------
% Initializing SVM
%------------------------------------------
% kernel matrix with a fixed value of mu
kernel = combinedKernels_onfly(X, mu);  %% large scale implementation

% standard SVM (L1)  with that fixed value of mu
[alpha, alpha_sv, obj, w_sqd, b, pos, train_err] = SVM(X, Y, kernel, option.lambdareg, C, alphainit);
% the weighted lagrangian of the support vectors: alpha_sv =w =(alpha(sv).*ysupport);

% the real objective value should be added by C_r/2 * r
% radius
[r beta] = radius_by_LP(X, componentDistances, mu); 

obj = obj + C_r/2 *r;

[grad] = gradRSVMPlusl1_ls(X,pos,alpha_sv, beta, componentDistances, C_r);

muOld   = mu; 

history.obj = [];
history.mu  = [];

%------------------------------------------------------------------------------%
% Update Main loop
%------------------------------------------------------------------------------%

while loop & nloopmax >0 ;

    nloop = nloop+1;
    
    history.mu  = [history.mu;mu];
    history.obj = [history.obj obj];

    %-----------------------------------------
    % Update weigths mu
    %-----------------------------------------
    t = cputime ;
    [mu,alpha_sv,b,pos,obj] = RSVMPlusL1Update_ls(X,mu,pos,alpha_sv,b,C,C_r,Y,grad,obj, componentDistances, option) ;
    
    %-----------------------------------------
    % Thresholding
    %-----------------------------------------

    if seuil ~=0 & max(mu) > seuil & nloop < option.seuilitermax
        mu=(mu.*(mu > seuil))*sum(mu)/sum(mu.*(mu > seuil));
    end;

    %-------------------------------
    % Numerical cleaning
    %-------------------------------
    mu(find(abs(mu < option.numericalprecision))) = 0;
    mu = mu / sum(mu);
    %-----------------------------------------------------------
    % Enhance accuracy of line search if necessary
    %-----------------------------------------------------------
    if max(abs(mu-muOld))<option.numericalprecision & option.goldensearch_deltmax > option.goldensearchmax
        option.goldensearch_deltmax = option.goldensearch_deltmax/10;
    elseif option.goldensearch_deltmax ~= goldensearch_deltmaxinit
        option.goldensearch_deltmax*10;
    end;


    %----------------------------------------------------
    % process approximate KKT conditions
    %----------------------------------------------------
   
    %----------------------------------------------------
    % process duality gap
    %----------------------------------------------------
    
    %------------------------------------------
    %  verbosity
    %------------------------------------------
    if verbose
        if nloop == 1 || rem(nloop,10)==0
            fprintf('--------------------------------------------------\n');
            fprintf('Iter | Obj.    | DiffBetas  |\n');
            fprintf('--------------------------------------------------\n');
        end;
        fprintf('%d   | %8.4f | %6.4f   \n',[nloop obj   max(abs(mu-muOld)) ]);
    end


    %----------------------------------------------------
    % check variation of mu conditions
    %----------------------------------------------------
    if  option.stopvariation==1 & option.stopKKT== 0 & max(abs(mu - muOld))<option.seuildiffsigma;
        loop = 0;
     
        fprintf(1,'variation convergence criteria reached \n');
        history.mu  = [history.mu;mu];
        history.obj = [history.obj obj];
    end;
    %----------------------------------------------------
    % check approximate KKT conditions
    %----------------------------------------------------
   

    %----------------------------------------------------
    % check KKT and variation of mu conditions
    %----------------------------------------------------
    

    %----------------------------------------------------
    % check duality gap
    %----------------------------------------------------
    
    %----------------------------------------------------
    % check for premature convergence
    %----------------------------------------------------
% 
%     if  max(abs(mu-muOld))<option.numericalprecision & (length(indpos)==1 | abs(max(mu(indpos))-1) < option.numericalprecision) ...
%             & option.goldensearch_deltmax <=option.goldensearchmax
%         loop = 0;
%         history.mu  = [history.mu;mu];
%         history.obj = [history.obj obj];
%         status = 1;
%         fprintf(1,'Premature Convergence KKT- Duality Gap Not Satisfied\n');
%     end;

    %-----------------------------------------------------
    % check nbiteration conditions
    %----------------------------------------------------
    if nloop >= nloopmax ,
        loop = 0;
        history.mu  = [history.mu;mu];
        history.obj = [history.obj obj];
        status=2;
        fprintf(1,'maximum number of iterations reached\n')
    end;
    if nloop < option.miniter & loop==0
        loop = 1;
    end;

    %-----------------------------------------------------
    % Updating Variables
    %----------------------------------------------------

    muOld  = mu ;

end;

