function [history, stopReason] = lars(yin, xin, XTX, type, stopCriterion, regularization, trace, quiet)
% %{
% This program implements "LARS" algorithm and its lasso modification
% introduced by Efron et. al. 2003.  Read the paper to understand codes
% of this function.  Each line of this file has corresponding equation
% number in Efron et. al. 2003 for reader's convenience.
%
%
% *** CAUTION
% history(1).mu_OLS contains original 'yin' to provide convinience in
% writing user defined stop criterion function. Actual mu_OLS of the first
% step should be just mean of yin which is a simple array of copy of
% history(1).mu. So, if history(1).mu_OLS contains that information, it is
% redundant. In fact, this is also contained in history(1).b, which is a
% bias of the output. Therefore, to provide more information to user who
% want to write his/her own stop criterion function, history(1).mu_OLS
% contains yin.
% 
% 
% Example 1: moderate size x.
%     stopCrioterion = {};
%     stopCrioterion{1,1} = 'maxKernels';
%     stopCrioterion{1,2} = 100;
%
%     XTX = lars_getXTX(x_original);      % this takes long time.
%     sol = lars(y, x, 'lasso', XTX, stopCriterion);
% 
% Example 2: very small size x, or a really really big size x
%     stopCrioterion = {};
%     stopCrioterion{1,1} = 'maxKernels';
%     stopCrioterion{1,2} = 100;
%
%     sol = lars(y, x, 'lasso', XTX, stopCriterion);
% 
% Note:
%     Users can add any kind of stop criterion by editing
%     the corresponding portion of this file.  See the code
%     for existing examples.
% 
% Note 2:
%     This m-file does not implement routine for missing data.
% %}
% 

global USING_CLUSTER;
global RESOLUTION_OF_LARS;
global REGULARIZATION_FACTOR;
lars_init();

regularization_factor = REGULARIZATION_FACTOR; % Tikhonov regularization factor (or the ridge regression factor)
                              % -> This should be small enough in this case to get
                              %    a reasonable pseudoinverse.

stopReason = {};

%%% Check parameters
if length(yin)==0 | length(xin)==0
    warning('\nInput or Output has zero length.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
if size(yin,1) ~= size(xin,1)
    warning('\nSize of y does not match to that of x.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
if ~strcmp(type, 'lasso') & ~strcmp(type, 'lars') & ~strcmp(type, 'forward_stepwise')
    warning('\nUnknown type of regression.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
if strcmp(type, 'forward_stepwise')
    warning('\nForward_stepwise is not implemented.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
    

if exist('regularization','var') & ~isempty(regularization)
    regularization = 10;
else
    regularization  = 0;
end

if ~exist('trace','var') | isempty(trace)
    trace=0;
end

if ~exist('quiet','var') | isempty(quiet)
    quiet=0;
elseif quiet==1
    trace=0;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data preparation

% Program automatically centers and standardizes predictors.
if ~exist('XTX','var')
    XTX=[];
end
no_xtx = 0;
if ~isempty(XTX)
    if ~quiet & trace >=0
        fprintf('\nLars is using the provided xtx.\n');
    end
elseif size(xin,2)^2 > 10^6
    if ~quiet & trace >=0
        fprintf('Too large matrix (size(x,2)^2 > 10^6). lars will not pre-calculate xtx.\n');
    end
    no_xtx      = 1;
    XTX         = lars_getXTX(xin,no_xtx);
else
%     fprintf('\nCalculating xtx.\n');
    XTX         = lars_getXTX(xin);
end

x               = XTX.x;            % normalized xin
mx              = XTX.mx;           % mean xin
sx              = XTX.sx;           % length of each column of xin
ignores         = XTX.ignores;      % indices for constant terms
all_candidate   = XTX.all_candidate;% indices for all possible columns
if ~no_xtx
    xtx         = XTX.xtx;          % xtx matrix
    dup_columns = XTX.dup_columns;  % duplicated columns which will be automatically ignored.
end

my              = mean(yin);
y               = yin-my;

n               = size(x,1);        % # of samples
m               = size(x,2);        % # of predictors

% Now, we can determine the maximum number of kernels.
%maxKernels = min(maxKernels, min(size(xin,1)-1, length(all_candidate)));
%maxKernels = min(maxKernels, min(rank(xin), length(all_candidate)));
existMaxKernels = 0;
existMSE        = 0;
for is = 1:size(stopCriterion,1)
    if strcmp(stopCriterion{is,1},'maxKernels')
        existMaxKernels = 1;
%        stopCriterion{is,2} = min(stopCriterion{is,2}, min(size(xin,1)-1, length(all_candidate)));
        stopCriterion{is,2} = min(stopCriterion{is,2}, min(rank(xin), length(all_candidate)));
        if stopCriterion{is,2}<1
            warning('Max Kernel is less than 1. It must be larger than 0.\n');
            stopCriterion{is,2} = 1;
        end
    end
    if strcmp(stopCriterion{is,1},'MSE')
        existMSE        = 1;
        if stopCriterion{is,2}<1.0e-10
            warning('Maximum MSE is too small. Automatically set to 1.0e-10\n');
            stopCriterion{is,2} = 1.0e-10;
        end
    end
end
if ~existMaxKernels
    is = size(stopCriterion,1);
    stopCriterion{is+1,1} = 'maxKernels';
%    stopCriterion{is+1,2} = min(size(xin,1)-1, length(all_candidate));  % Stop when size of active set is data.maxKernels.
    stopCriterion{is+1,2} = min(rank(xin), length(all_candidate));  % Stop when size of active set is data.maxKernels.
end
if ~existMSE
    is = size(stopCriterion,1);
    stopCriterion{is+1,1} = 'MSE';
    stopCriterion{is+1,2} = 1.0e-10;
end
    




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization

active      = [];               % active set
inactive    = all_candidate;    % inactive set

mu_a        = zeros(n,1);       % current estimate (eq. 2.8)
mu_a_plus   = 0;                % next estimate (eq. 2.12)
mu_a_OLS    = 0;                % OLS estimate  (eq. 2.19)

beta        = zeros(1,size(x,2));
beta_new    = beta;
beta_OLS    = beta;

history.active_set          = [];
history.add                 = [];
history.drop                = [];
history.beta_norm           = [];
history.beta                = [];
history.b                   = my;
history.mu                  = my;
history.beta_OLS_norm       = [];
history.beta_OLS            = [];
history.b_OLS               = my;
history.mu_OLS              = my*ones(size(yin));
history.MSE                 = sum(y.^2)/length(y);
history.R_square            = 0;
history.resolution_warning  = [];


if var(yin)==0
    stopReason{1} = 'zeroVarY';
    stopReason{2} = var(yin);
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


c               = 0;                % correlation vector
C_max           = max(abs(c));
C_max_ind       = [];
C_max_ind_pl    = [];
drop            = [];               % used for 'lasso'
k               = 1;                % iteration index
if ~quiet & trace >= 0
    fprintf('Active predictors / total  :::::  Current iteration\n                            ');
end
while 1

    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%           Exit Criterions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if exist('stopCriterion','var')     % If there is any stop criterion...
        % Any of these is satisfied, algorithm stops.
        %
        % Other criterions: Make your own cost function.
        
        for is = 1:size(stopCriterion,1)  % there can be a number of stop criterions.

            % Default Criterions.
            % Maximum number of consecutive drops or maximum number of drops within a window.
            if strcmp(stopCriterion{is,1},'maxDrops')
                drop_window = stopCriterion{is,2}(1);   % number of backward history to be considerd. 0 means all history.
                if drop_window==0
                    drop_window=k;
                end
                drop_n      = min(drop_window,stopCriterion{is,2}(2));   % number of maximum drops within the window.
                drop_vector = [];
                for z = max(k-drop_window+1,1):k
                    drop_vector=[drop_vector, history(z).drop];
                end
                if length(drop_vector)>=drop_n
                    stopReason{1} = 'maxDrops';
                    stopReason{2} = drop_n;
                    break;
                end
            end
            % Maximum number of kernels.
            if strcmp(stopCriterion{is,1},'maxKernels')
                if length(active) >= min(stopCriterion{is,2}, min(size(xin,1)-1, length(all_candidate)))
                    stopReason{1} = 'maxKernels';
                    stopReason{2} = length(active);
                    break;
                end
            end
            % Maximum number of iterations.
            if strcmp(stopCriterion{is,1},'maxIterations')
                if k >= stopCriterion{is,2}
                    stopReason{1} = 'maxIterations';
                    stopReason{2} = k;
                    break;
                end
            end
            % MSE.
            if strcmp(stopCriterion{is,1},'MSE')
                if history(k).MSE <= stopCriterion{is,2}
                    stopReason{1} = 'MSE';
                    stopReason{2} = history(k).MSE;
                    break;
                end
            end
            
            % User defined stop criterion.
            if strcmp(stopCriterion{is,1},'userDefinedCriterion')
                fhandle = stopCriterion{is,2}.fhandle;
                r_fhandle = fhandle(history, stopCriterion{is,2}.data);
                if r_fhandle.stop
                    stopReason{1} = 'userDefinedCriterion';
                    stopReason{2} = r_fhandle;
                    break;
                end
            end
            
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    end % end of stop criterion checking :  if exist('stopCriterion','var')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    
    if length(stopReason)>0 % if there is any reason of stopping the algorithm, exit loop.
        break;
    end
    
    
    
    
    
    
    
    
    
    
    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%           LARS Algorithm
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %if ~USING_CLUSTER
    %fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%5d/%5d    :::::   %5d',length(active),length(all_candidate),k);
    %end
    % start of algorithm
    c                   = x'*(y-mu_a);                      % eq 2.8
    [C_max,C_max_ind]   = max(abs(c(inactive)));            % eq 2.9
    C_max_ind           = inactive(C_max_ind);
    % active          = sort(union(active,C_max_ind));      % If there is no machine limit, this can be used.
    % But because of machine limit, there can be multiple new predictors.
    % This dramatically improves the overall precision of the result,
    % and speeds up the whole process.
    C_max_ind_pl        = abs(c(inactive))>C_max-RESOLUTION_OF_LARS;
    C_max_ind_pl        = inactive(C_max_ind_pl);
    active          = sort(union(active,C_max_ind_pl));  
    inactive        = setdiff(all_candidate, active);   % eq 2.9
    if strcmp(type,'lasso')
        if ~isempty(drop) & length(find(drop==C_max_ind))==0           % If there is a drop, that must have the maximum correlation in inactive set.
            if ~quiet & trace >=0
                fprintf('\n');
                warning('Dropped item and index of maximum correlation is not the same. But it is being ignored here...');
                fprintf('\n                            ');
            end
            active(find(active==C_max_ind))=[];
        end
        if ~isempty(drop)
            C_max_ind   = [];
            C_max_ind_pl= [];
        end
        active          = setdiff(active,drop);             % eq 3.6
        inactive        = sort(union(inactive,drop));       % eq 3.6
    end

    s       = sign(c(active));              % eq 2.10
    xa      = x(:,active).*repmat(s',n,1);  % eq 2.4
    if ~no_xtx
        ga  = xtx(active,active).*(s*s');   % eq 2.5
    else
        ga  = xa'*xa;                       % eq 2.5
    end
    if regularization > 2
        ga  = ga+eye(length(ga))*regularization_factor; % This routine will make the test below
    end
    invga   = ga\eye(size(ga,1));           % eq 2.5
    aa      = sum(sum(invga))^(-1/2);       % eq 2.5
    wa      = aa*sum(invga,2);              % eq 2.6
    ua      = xa*wa;                        % eq 2.6
    
    % test using eq 2.7
    test_1      = xa'*ua;
    test_2      = aa*ones(size(test_1));
    test_1_2    = sum(sum(abs(test_1-test_2)));
    test_3      = norm(ua) - 1;
    
    history(k+1).resolution_warning=0;
    if test_1_2 > RESOLUTION_OF_LARS*100 | abs(test_3 ) > RESOLUTION_OF_LARS*100
        if regularization <=2
            if ~quiet & trace>0
                fprintf('\n');
                warning('Eq 2.7 test failure.');
                fprintf('\n                            ');
            end
            regularization = regularization + 1;
            if regularization > 2
                if ~quiet & trace>0
                    fprintf('\n');
                    warning('Lots of Eq 2.7 test failure. Regularization will be applied from now on.');
                    fprintf('\n                            ');
                end
            end
        end
        history(k+1).resolution_warning=1;
    end
    


    a               = x'*ua;                % eq 2.11
    tmp_1           = (C_max - c(inactive))./(aa - a(inactive));
    tmp_2           = (C_max + c(inactive))./(aa + a(inactive));
    tmp_3           = [tmp_1, tmp_2];
    tmp             = tmp_3(find(tmp_3>0));
    gamma = min(tmp);                       % eq 2.13
    if length(gamma)==0 % if this is the last step (i.e. length(active)==maxKernels)
        gamma = C_max/aa;                   % eq 2.19, eq 2.21 and 5 lines below eq 2.22
    end
    
    d               = zeros(1,m);
    d(active)       = s.*wa;

    if length(find(d(active)==0))
        fprintf('\n');
        warning('Something wrong with vector d: eq 3.4.');
        fprintf('\n                            ');
    end

    tmp             = zeros(1,m);
    tmp(active)     = -1*beta(active)./d(active);   % eq 3.4
    tmp2            = tmp(find(tmp>0));
    
    drop            = [];
    gamma_tilde = inf;                              % eq 3.5
    if ~isempty(tmp2) & gamma >= min(tmp2)
        gamma_tilde = min(tmp2);                    % eq 3.5
        drop        = find(tmp==gamma_tilde);       % eq 3.6
    end
    
    if strcmp(type, 'lars')
        mu_a_plus   = mu_a + gamma*ua;              % eq 2.12
        beta_new    = beta + gamma*d;               % eq 3.3
        drop        = [];
    elseif strcmp(type, 'lasso')
        mu_a_plus   = mu_a + min(gamma, gamma_tilde)*ua;    % eq 3.6
        beta_new    = beta + min(gamma, gamma_tilde)*d;     % eq 3.3
        active      = setdiff(active,drop);                 % eq 3.6
        inactive    = setdiff(all_candidate,active);
        beta_new(drop) = 0;
    elseif strcmp(type, 'forward_stepwise')
        drop        = [];
        error('forward.stepwise has not been implemented yet.');
        return;
    end
    
    mu_a_OLS    = mu_a + C_max/aa*ua;          % eq 2.19, 2.21
    beta_OLS    = beta + C_max/aa*d;           % eq 2.19, 2.21
    MSE         = sum((y - mu_a_OLS).^2)/length(y);
    
    
    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % update and save
    mu_a = mu_a_plus;
    beta = beta_new;
    
    % history with scale correction
    k = k+1;
    history(k).active_set   = active;
    history(k).drop         = drop;
    history(k).add          = C_max_ind_pl;
    history(k).beta_norm    = beta(active);
    history(k).beta         = beta(active)./sx(active);
    history(k).b            = my - sum(mx./sx.*beta);
    history(k).mu           = xin * (beta./sx)' + history(k).b;
    history(k).beta_OLS_norm= beta_OLS(active);
    history(k).beta_OLS     = beta_OLS(active)./sx(active);
    history(k).b_OLS        = my - sum(mx./sx.*beta_OLS);
    history(k).mu_OLS       = xin * (beta_OLS./sx)' + history(k).b_OLS;
    history(k).MSE          = MSE;
    history(k).R_square     = 1 - var(yin  - history(k).mu_OLS)/var(yin);
    
    
    % exit if exact mathing is achieved.
    if abs(C_max/aa - min(gamma,gamma_tilde)) < RESOLUTION_OF_LARS
        stopReason{1}   = 'ExactMatching';
        stopReason{2}   = 0;
        break;
    end

    
    
    
end % end of while loop
if ~quiet & trace >=0
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%5d/%5d    :::::   %5d\n',length(active),length(all_candidate),k);
end    


return;



