function [betahat, lambdahat, Hisbeta, HisAclen,  Bic, Lam] = glmnet_path(X,y,opts)
% ========================================================================
% Approximate the  solution path of Lasso/Enet panalize linear regerssion 
% in High dimentional with n << p by Glmnet  to find the minimizer of 
%             1/2||X*beta-y||^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.mu :     stop parameter when ||beta^k||_0 >mu*n   (defaut: 0.5)
% opts.tol:     stop when the relative chage on the path smaller than 
%               opts.tol (default: 1e-4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
% 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.N = 100;
    opts.tau = 1e-8;  % if we set the defaut value of Glmnet (1e-4), it will not work well for high dinamic range signal
    opts.mu = 0.5;
    opts.tol = 1e-4;
elseif nargin > 3
    disp(' too much  inputment ')
end
% define the set Lambda and 
Xty = X'*y;          
norminf = norm(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;
Bic = [];
%% main loop for approximate the whole solution path 
disp(' run Glmnet on Lambda ...')
%% set options parameter default exlude 
% options.dfmax = opts.mu*n;
% options.lambda         =       Lambda;
options.weights        =           [];
options.alpha          =          1.0;
options.nlambda        =          100;
options.lambda_min     =            0;
options.lambda         =       Lambda;
options.standardize    =         true;
options.thresh         =     opts.tol;
options.dfmax          =    opts.mu*n;
options.pmax           =            0;
options.exclude        =           [];
options.penalty_factor =           [];
options.maxit          =          100;
options.HessianExact   =        false;
options.type           = 'covariance';
%run glmnet by Hui Jiang, jiangh@stanford.edu 
Hist = glmnet(X,y,'gaussian',options);
% Hist = glmnet(X,y);  % The default setting works not well  for high dinamic range signal
Hisbeta = Hist.beta;
HisAclen = Hist.df;
Lam = Hist.lambda;
for k = 1:size(Hisbeta,2)
bick = 0.5*norm(X*Hisbeta(:,k) - y)^2  + log(p)*log(n)*HisAclen(k);% Bic
Bic = [Bic,bick];  
end
% select the optimal betahat and lambdahat  
id  = find(Bic == min(Bic));
lambdahat = Lam(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


