function path = follow_entire_path(Ks,y,loss,ds,path_params,Ks_test,ytest)

% FOLLOW THE ENTIRE PATH OF REGULARIZATION
%
% INPUT
% Ks : kernel matrices
% y : response
% loss : loss.type= 'regression' or 'logistic'
% ds : weights of the block 1-norm
% path_params : parameters of the path (defined below)
% Ks_test, ytest (optional) : testing data
%
% OUTPUT
% path.alphas                   dual variables
% path.etas                     weights of kernels
% path.sigmas                   log(lambda)
% path.netas                    number of "nonzero" weights
% path.ks                       number of corrector steps
% path.training_errors          
% path.testing_errors
% path.training_predictions     values of K * alpha
% path.testing_predictions      values of K * alpha


% INITIALIZATION OF CONSTANTS
mu                  = path_params.mu;                   % parameter of log-barrier
EPS1                = path_params.EPS1;                 % precision parameters of Newton steps (very small and fixed)
EPS2                = path_params.EPS2;                 % precision parameter of tube around path (adaptive)
predictor_type      = path_params.predictor_type;       % 1 : first order predictor, 2 : second order
efficient_predictor = path_params.efficient_predictor;  % 1 : efficient predictor steps, 0 : full steps
efficient_eta       = path_params.efficient_eta;        % real value : threshold eta for corrector steps, 0 : no threshold
% i.e., only consider the significant ones, but it is buggy because, a zero
% eta might become significant and negative (when a constraint is broken)
maxsigma            = path_params.maxsigma;             % maximum value of sigma = -log(lambda);
newton_iter1        = path_params.newton_iter1 ;         % number of iterations with no modification
newton_iter2        = path_params.newton_iter2  ;         % delta number of iterations under which EPS2 is divided/multiplied by 2
maxdsigma           = path_params.maxdsigma;            % maximal step
mindsigma           = path_params.mindsigma;            % minimal step
maxvalue_EPS2       = path_params.maxvalue_EPS2;        % maximum value of tolerance for predictor steps

problemtype     = problem_type(loss); % classification or regression
m               = Ks.m;

test = 0;       % test = 1 if performs testing online
if nargin>=6,
    if ~isempty(Ks_test), test = 1; end
end

 MIN_DSIGMA_OCURRED = 0;
 
% FIRST STEP
beta = target(y,loss);
Kbeta = kernel_operation(Ks,1,beta);
start = sqrt( beta' * Kbeta ) ./ ds ;
lambda = 1.3 * max(start);
sigma = -log(lambda);
alpha =  beta / lambda;

% selecte subsets of eta
etas_comp = get_etas(Ks,y,loss,alpha,ds,lambda,mu);
indices = find( etas_comp' .* ds.^2 / mu  > efficient_eta );
 
% first corrector step
[alpha,lambda2,exit_status,k]= newton_method(Ks,y,loss,alpha,ds,lambda,mu,indices,0,EPS1,500);
    [f,grad,hessianinvchol0,eta,lambda2,da_dsigma,d2a_dsigma2] = all_derivatives(Ks,y,loss,alpha,ds,lambda,mu);

   

    da_dsigma_old = da_dsigma;
d2a_dsigma2_old = d2a_dsigma2;
alpha_old = alpha;
sigma_old = sigma;

% store first step
alphas = alpha;
sigmas = sigma;
etas = eta;
netas = length(find( eta' .* ds.^2 > mu * 4 ) );
ks = k;
lambda2s = lambda2;
normgrads = norm(grad)^2;
nevals = 0;
EPS2s = EPS2;

ypred = - kernel_operation(Ks,1,alpha) * eta;
switch problemtype,
    case 'regression'
        training_error = sum( (ypred-y).^2 ) / length(y);

    case 'classification'
        ypreds = sign( ypred );
        training_error = length(find( abs(ypreds-y)> 0 )) / length(y);

    case 'classification_unbalanced'
        ypreds = sign( ypred );
        training_error = sum( y(:,2) .* ( abs(ypreds-y(:,1))> 0 ) ) / length(y);


end
training_errors = training_error;
training_predictions = ypred;
path.ytrain = y;
path.problemtype = problemtype;

if test
    path.ytest = ytest;
    yhat = - kernel_operation_test(Ks_test,1,alpha) * eta;
    ypred = yhat;
    switch problemtype,
        case 'regression'
            testing_error = sum( (ypred-ytest).^2 ) / length(ytest);
        case 'classification'
            ypreds = sign( ypred );
            testing_error = length(find( abs(ypreds-ytest)> 0 )) / length(ytest);
        case 'classification_unbalanced'
            ypreds = sign( ypred );
            testing_error = sum( ytest(:,2) .* ( abs(ypreds-ytest(:,1))> 0 ) ) / length(ytest);

    end
    testing_errors = testing_error;
    testing_predictions = ypred;
end


fprintf('iter %d - sigma=%f - n_corr=%d - lambda2=%1.2e - EPS2=%1.2e - n_pred=%d - neta=%d\n',length(sigmas),sigmas(end),ks(end),lambda2s(end),EPS2,nevals(end),netas(end));
max_steps_no_progress = 6;

% PREDICTOR-CORRECTOR STEPS
% exits if upperbound on sigma is reached or Newton steps didnot converge
delta_sigma = 1;
try_all_predictors=0;
while sigma < maxsigma  & lambda2 < EPS2 & delta_sigma > 1e-6 & length(sigmas)<=2000

    % every 100 moves, reset to large tube around data
    if mod(length(sigmas),100)==1, EPS2                = path_params.EPS2;  end

    % choosing dsigma : predictor steps
    neval = 0;

    if length(sigmas)==1,
        boosted = 0;
        % first predictor step: try out several values
        dsigmas = 10.^[-5:.5:0];
        for idsigma=1:length(dsigmas);
            dsigma=dsigmas(idsigma);
            compute_newlambda2;
            neval = neval + 1;
        if isinf(lambda2) || isnan(lambda2) || (lambda2>EPS2),  break;  
            
            end;
        end
        dsigma = dsigmas(idsigma-1);
        switch predictor_type
            case 0, newalpha = alpha;
            case 1, newalpha = alpha+dsigma*da_dsigma;
            case 2, newalpha = alpha+dsigma*da_dsigma+ dsigma^2 * .5 * d2a_dsigma2;
        end
        newlambda = exp(-sigma-dsigma);

    else
        % for regression, if long stagnation with the same number of
        % kernels, the path is likely to be linear in 1/lambda
        if length(netas)>=10 & isequal(problemtype,'regression')
            if all(netas(end-9:end)==netas(end)), boosted = 1; fprintf(' boosted '); else boosted = 0; end
        else
            boosted = 0;
        end
        predictor_step;
    end

    predicted_type=predictor_type;
    if dsigma <= 1e-5
        switch predictor_type
            case 2,
                predictor_type = 1;
                predictor_step
                predicted_type = 1;
                if dsigma <= 1e-5
                    predictor_type = 0;
                    predictor_step
                    predicted_type=0;
                end
                predictor_type = 2;
            case 1,
                predictor_type = 0;
                predictor_step
                predictor_type = 1;
                predicted_type =0;
        end
    end

    if dsigma <= 1e-5,
        % not event the simpler predictor steps work, get out!
        break;
    end

    if try_all_predictors
        alpha_predicted = newalpha;
        dsigma_predicted = dsigma;
        predictor_type = 0;
        predictor_step;
        predictor_type = 1;
        if dsigma>dsigma_predicted,
            predicted_type=0;
            fprintf(' trivial is better ');
        else
            newalpha = alpha_predicted ;
            dsigma =  dsigma_predicted;
        end
    end


    go_on=1;
    counter = 1;

    % select etas
    etas_comp = get_etas(Ks,y,loss,newalpha,ds,lambda,mu);
    indices = find( etas_comp' .* ds.^2 / mu  > efficient_eta );

    % store for later use
    sigma0 = sigma;
    newalpha0 = newalpha;
    dsigma0=dsigma;
    alpha0=alpha;

    while go_on & counter <= max_steps_no_progress
        % try to perform Newton steps from predicted sigma
        % if does not converge, diminishes sigma and divides EPS2 by 8
        % if never converges after 4 times, exit
        counter = counter + 1;
        sigmapot = sigma + dsigma;
        try
            [alphapot,lambda2,exit_status,k,nevals_newton,exit_params]= newton_method(Ks,y,loss,newalpha,ds,exp(-sigmapot),mu,indices,0,EPS1,30);
        catch
            q=round(10*sum(clock));
            q
            save(sprintf('error_%d',q));
            error('error in newton method');
        end

        if strcmp(exit_status,'max_iterations') || lambda2 >  1e-5 || ...
                strcmp(exit_status,'infinite_lambda_0') || ...
                strcmp(exit_status,'no_inverse_hessian_0') || ...
                strcmp(exit_status,'nan_lambda_0')
            fprintf('Newton''s method takes too long!! dsigma=%e lambda2=%e\n',dsigma,lambda);
            dsigma = dsigma / 8;
            EPS2 = EPS2 / 8;
            sigmapot = sigma + dsigma;
            if boosted

                switch predicted_type
                    case 0, newalpha = alpha ;
                    case 1, newalpha = alpha + ( exp(dsigma) - 1 ) * da_dsigma;
                    case 2, newalpha = alpha + ( exp(dsigma) - 1 ) * da_dsigma + ...
                            +  .5 * ( exp(dsigma) - 1 )^2 * (  d2a_dsigma2 - da_dsigma );
                end

            else
                switch predicted_type
                    case 0, newalpha = alpha;
                    case 1, newalpha = alpha+dsigma*da_dsigma;
                    case 2, newalpha = alpha+dsigma*da_dsigma+ dsigma^2 * .5 * d2a_dsigma2;
                end
            end
        else go_on=0;
        end
    end
    if counter == max_steps_no_progress & go_on, break; end

    alpha_old = alpha;
    sigma_old = sigma;

    alpha = alphapot;
    delta_sigma = abs( sigma - sigmapot);
    sigma = sigmapot;
    lambda = exp(-sigma);


    % store step
    try
        da_dsigma_old = da_dsigma;
        d2a_dsigma2_old = d2a_dsigma2;
        [f,grad,hessianinvchol,eta,lambda2,da_dsigma,d2a_dsigma2] = all_derivatives(Ks,y,loss,alpha,ds,lambda,mu);


        if isinf(f)
            % if the value of f is infinite, then it means that the previous
            % corrector steps with less etas made one of the forgotten eta
            % active -> redo the Newton steps with all etas


            counter = 1;
            sigma = sigma0;
            dsigma = dsigma0;
            newalpha = newalpha0;
            alpha =alpha0;

            while go_on & counter <= max_steps_no_progress
                % try to perform Newton steps from predicted sigma
                % if does not converge, diminishes sigma and divides EPS2 by 8
                % if never converges after 4 times, exit
                counter = counter + 1;
                sigmapot = sigma + dsigma;
                try
                    [alphapot,lambda2,exit_status,k,nevals_newton,exit_params]= newton_method(Ks,y,loss,newalpha,ds,exp(-sigmapot),mu,1:m,0,EPS1,30);
                catch
                    q=round(10*sum(clock));
                    q
                    save(sprintf('error_%d',q));
                    error('error in newton method');
                end

                if strcmp(exit_status,'max_iterations') || lambda2 >  1e-5 || ...
                        strcmp(exit_status,'infinite_lambda_0') || ...
                        strcmp(exit_status,'no_inverse_hessian_0') || ...
                        strcmp(exit_status,'nan_lambda_0')
                    fprintf('Newton''s method takes too long!! dsigma=%e lambda2=%e\n',dsigma,lambda);
                    dsigma = dsigma / 8;
                    EPS2 = EPS2 / 8;
                    sigmapot = sigma + dsigma;
                    if boosted

                        switch predicted_type
                            case 0, newalpha = alpha ;
                            case 1, newalpha = alpha + ( exp(dsigma) - 1 ) * da_dsigma;
                            case 2, newalpha = alpha + ( exp(dsigma) - 1 ) * da_dsigma + ...
                                    +  .5 * ( exp(dsigma) - 1 )^2 * (  d2a_dsigma2 - da_dsigma );
                        end

                    else
                        switch predicted_type
                            case 0, newalpha = alpha;
                            case 1, newalpha = alpha+dsigma*da_dsigma;
                            case 2, newalpha = alpha+dsigma*da_dsigma+ dsigma^2 * .5 * d2a_dsigma2;
                        end
                    end
                else go_on=0;
                end
            end
            if counter == max_steps_no_progress & go_on, break; end

            alpha_old = alpha;
            sigma_old = sigma;

            alpha = alphapot;
            delta_sigma = abs( sigma - sigmapot);
            sigma = sigmapot;
            lambda = exp(-sigma);
            [f,grad,hessianinvchol,eta,lambda2,da_dsigma,d2a_dsigma2] = all_derivatives(Ks,y,loss,alpha,ds,lambda,mu);



        end

        % if the values of the derivatives with respect to the path are infinite or NAN, put them to zero
        if sum( isnan(da_dsigma) | isinf(da_dsigma) ),
            fprintf('derivative of path is undefined\n');
            fprintf('only doing zero order now (trivial predictors)\n');
            da_dsigma=zeros(size(da_dsigma));
            d2a_dsigma2=zeros(size(d2a_dsigma2));
            predictor_type = 0;
            try_all_predictors=0;

        end
        if sum( isnan(d2a_dsigma2) | isinf(d2a_dsigma2) ) & predictor_type ==2,
            fprintf('second derivative of path is undefined\n');
            fprintf('only doing first order now\n');
            d2a_dsigma2=zeros(size(d2a_dsigma2));
            predictor_type = 1;
            try_all_predictors=1;
        end


    catch
        q=round(10*sum(clock));
        q
        save(sprintf('error_%d',q));
        error('error in all derivatives - stack saved for examination');
    end
    lambda2s = [ lambda2s lambda2];
    normgrads = [ normgrads norm(grad)^2 ];
    alphas = [ alphas alpha];
    sigmas = [ sigmas sigma];
    etas = [ etas eta];
    netas = [ netas length(find( eta' .* ds.^2 > mu * 4 ) ) ];
    ks = [ ks k];
    nevals = [ nevals neval];

    % update EPS2
    if ~strcmp(exit_status,'df_small')
        EPS2=EPS2 * 2^ ( (newton_iter1 -k)/newton_iter2 );
    else
        EPS2=EPS2 * 2^ ( (newton_iter1 -k)/newton_iter2/2 );
    end
    if lambda2 >= EPS2/10, EPS2 = lambda2 * 10; end
    EPS2 = max(EPS2, 1e-5);
    EPS2 = min(EPS2, maxvalue_EPS2);
    EPS2s = [ EPS2s EPS2 ];
    ypred = - kernel_operation(Ks,1,alpha) * eta;
    switch problemtype,
        case 'regression'
            training_error = sum( (ypred-y).^2 ) / length(y);
        case 'classification'
            ypreds = sign( ypred );
            training_error = length(find( abs(ypreds-y)> 0 )) / length(y);
        case 'classification_unbalanced'
            ypreds = sign( ypred );
            training_error = sum( y(:,2) .* ( abs(ypreds-y(:,1))> 0 ) ) / length(y);
    end
    training_errors = [ training_errors training_error];
    training_predictions = [ training_predictions ypred];
    
    % if the training errors is zero, stop the path following soon
    if training_errors(end)<training_errors(1)*1e-6
       maxsigma = min(maxsigma,sigma+1); 
    end
    
    
    % perform testing if applicable
    if test
        yhat = - kernel_operation_test(Ks_test,1,alpha) * eta;
        ypred = yhat;
        switch problemtype,
            case 'regression'
                testing_error = sum( (ypred-ytest).^2 ) / length(ytest);
            case 'classification'
                ypreds = sign( ypred );
                testing_error = length(find( abs(ypreds-ytest)> 0 )) / length(ytest);
            case 'classification_unbalanced'
                ypreds = sign( ypred );
                testing_error = sum( ytest(:,2) .* ( abs(ypreds-ytest(:,1))> 0 ) ) / length(ytest);


        end
        testing_errors = [ testing_errors testing_error ];
        testing_predictions = [ testing_predictions ypred];

    end



    fprintf('iter %d - sigma=%f - n_corr=%d - lambda2=%1.2e - EPS2=%1.2e - n_pred=%d - neta=%d - dsigma=%f - status=%s\n',length(sigmas),sigmas(end),ks(end),lambda2s(end),EPS2,nevals(end),netas(end),dsigma,exit_status);
    if dsigma < mindsigma,
        % after 6 times in a row with a small progress, exit
        MIN_DSIGMA_OCURRED = MIN_DSIGMA_OCURRED + 1;
        if MIN_DSIGMA_OCURRED == 6
            break;
        end
    else
        MIN_DSIGMA_OCURRED = 0;
    end
end

path.alphas=alphas;                                 
path.etas=etas;
path.sigmas=sigmas;
path.netas=netas;
path.ks=ks;
% path.nevals=nevals;
% path.normgrads = normgrads;
% path.lambda2s = lambda2s;
% path.EPS2s = EPS2s;
path.training_errors = training_errors;
path.training_predictions = training_predictions;
if test, path.testing_errors = testing_errors;
    path.testing_predictions = testing_predictions;
end

if 0,
    % debugging
    q=round(10*sum(clock));
    q
    save(sprintf('debug_%d',q));
    fprintf('results of path following saved');
end
