function [betahat, lambdahat, Hisbeta, HisAclen, HisIter, Bic, Lam] = snap_path(X,y,opts)
% ========================================================================
% Approximate the  solution path of Lasso/Enet panalize linear regerssion 
% in High dimentional with n << p by  using Semi-smooth Newton algorithm 
% (SNA) to find the minimizer of 
%     1/2||X*beta-y||^2 + alpha/2|beta||^2 + lambda||beta||_{1}, 
% on 
%                 Lambda ={lambda_{1},...,lambda_{N}}, 
% with warm-starting. 
% 
% Copyright (c), Dec, 2013
% Yuling Jiao  Dept. Math.,Whuhan Univiversity
% Email: yulingjiaomath@whu.edu.cn
% ========================================================================
% Input:  
% X:      the covariate matrix in R^{n*p} ( unite normalized) 
% y:      response vector in R^{n}
% opts is a structure contained the following:
% opts.N:       the size of Lambda  (defaut: 100)
% opts.tau:     Lambda = norm(X'y,inf)*[opts.tau,1] (defaut: 1e-8) 
% opts.alpha:   Enet parameter (default: 0)
% opts.initial: intial value (default: 0) 
% opts.K:       maxiteration number in SNA for safe guard   (defaut: 1)
% opts.mu :     stop parameter when ||beta^k||_0 >mu*n   (defaut: 0.5)
% opts.change:  stop when  changes on the path <= opts.change (defaut: 0.05)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output: 
% betahat:      an estimator  selected by BIC
% lambdahat:    regularization parameter selected by BIC
% Hisbeta:      a matrix  storing the history of 
%               beta_{lambda_{k}} at each lambda_{k},k =1,...N 
% HisAClen:     a vector storing the history of acive set size of
%               beta_{lambda_{k}}  at each lambda_{k},k =1,...N
% HisIter:      a vector storing  the iteration number used when SNA stop at
%               each  lambda_{k},k =1,...N
% Bic:          BIC value  at each  lambda_{k},k =1,...N
% Lam:          subset of Lambda actrully used for computing the whole path
% ========================================================================
[n,p] = size(X);
if nargin < 2
    disp(' error not enough inputment ')
elseif nargin == 2 || isempty(opts)
    opts.initial = zeros(p,1);
    opts.alpha = 0;
    opts.N = 100;
    opts.tau = 1e-8;
    opts.K = 1;
    opts.mu = 0.5;
    opts.change = 0.05;
elseif nargin > 3
    disp(' too much  inputment ')
end
% define the set Lambda and 
opts.Xty = X'*y;          % X'*y used in SNA   
norminf = norm(opts.Xty, inf);  
lambdaMin = opts.tau;
lambdaMax = 1;
loghi = log(lambdaMax);
loglo = log(lambdaMin);
logrange = loghi - loglo;
interval = -logrange/(opts.N-1);
Lambda = exp(loghi:interval:loglo)';
Lambda = Lambda(2:end);
Lambda = norminf*Lambda;
Hisbeta = zeros(p,1);
HisAclen = [];
HisIter =[];
Bic = [];
Lam = [];
%% main loop for approximate the whole solution path 
disp(' run SNA on Lambda ...')
df = 0;
k = 0;
lenLam = length(Lambda);
mun = opts.mu*n;
change = opts.change + 1;
while   df < mun && k < lenLam && change >opts.change
         k = k + 1;
         lambda = Lambda(k);
         Lam = [Lam;lambda];
         obeta = opts.initial;
         [ebeta, Aclen, nIter] = sna(X,lambda,opts);
         Hisbeta = [Hisbeta,ebeta];
         HisIter = [HisIter;nIter];
         df = Aclen;
         HisAclen = [HisAclen;df];
         change = norm(obeta-ebeta,inf);
         opts.initial = ebeta;  % warmstart 
         bick = 0.5*norm(X*ebeta - y)^2  + log(p)*log(n)*df;% Bic
         Bic = [Bic,bick];  
end
% select the optimal betahat and lambdahat  
id  = find(Bic == min(Bic));
lambdahat = Lambda(id(end));
betahat = Hisbeta(:,id(end));
% % debias
% supp = find(betahat);
% Xsupp =  X(:,supp);
% betahat(supp) = (Xsupp'*Xsupp)\(Xsupp'*y);
disp('The approximation is getton and optimal solution is selected')
end


function [ebeta,  Aclen,  nIter] = sna(X,lambda,opts)
% ========================================================================
% Semi-smooth Newton  Algorithm for solving 
% min_{beta}{1/2||X*beta-y||^2 + alpha/2|beta||^2 + lambda||beta||_{1}
% Copyright (c), Des, 2013
% Yuling Jiao  Dept. Math.,Whuhan Univiversity
% Email: yulingjiaomath@whu.edu.cn
% ========================================================================
% Input:  
% X:           the normalized covariate matrix in R^{n*p},  
% lambda:      regularization parameter 
% opts a structrue contained the following:
% opts.Xty:    X'*y;
% opts.alpha : Enet  parameter  (defaut: 0)
% opts.K:      maxiteration number of SNA   (defaut: 1)
% opts.mu:     stop parameter when |A_{k}|>mu*n  (defaut: 0.5)
% Output: 
% ebeta:    an estimate of regression coefficient 
% Aclen:    acive set size  of ebeta
% nIter:    iteration number when SNA stop
% ========================================================================
[n,p] = size(X);
if nargin < 3 || isempty(opts)
    disp('error, not enough inputment ')
elseif nargin > 3 
    disp('error, too much inputment ')
else
    Xty = opts.Xty;
    initial  = opts.initial;
    alpha = opts.alpha;
    K = opts.K;
    mu = opts.mu;
end
pd = initial + Xty  - X'*(X*initial);
Ac = find(abs(pd)>=lambda);
nIter = 0;
Aclen = length(Ac);
ebeta = initial;
stop  = 1;
mun = mu*n;
while stop ~= 0 && nIter < K && Aclen < mun
   ebeta = zeros(p,1);
   XtyAc = Xty(Ac);
   nIter = nIter + 1;
   ted = lambda*sign(pd(Ac));
   oAc = Ac;
   oAclen = Aclen;
   tempright  = (XtyAc - ted) ;
   Xac = X(:,Ac);
   G = (Xac'*Xac) + alpha*eye(Aclen);
   tbeta = G\tempright;
   ebeta(Ac) = tbeta;
   ed =  Xty  - X'*(Xac*tbeta);
   ed(Ac) = ted;
   pd = ebeta + ed;
   Ac = find(abs(pd)>=lambda);
   Aclen = length(Ac);
   if Aclen == oAclen  
        if oAc == Ac
          stop = 0;
        end
   end
end
end
