function [betahat, lambdahat, Hisbeta, HisAclen, HisIter, Bic,Lam] = cdj_path(X,y,opts)
% ========================================================================
% Approximate the  solution path of Lasso panalize linear regerssion  in
% High dimentional with n << p by minmizing 
%                 1/2||X*beta-y||^2  + lambda||beta||_{1} 
% with Coodinate decesent (Jacobi version, Iterative thresholsing) algorithm 
% 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 cd_j for safe guard   (defaut:50)
% opts.mu :     stop parameter when ||beta^k||_0 >mu*n   (defaut: 0.5)
% opts.step:    step size for convergence (defaut: norm(X))
% opts.tol:     stop when the relative chage on the path smaller than 
%               opts.tol (default: 1e-8)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output: 
% betahat:      an estimate of regression coefficient selected by BIC
% lambdahat:    regularization parameter selected by BIC
% Hisbeta:      a matix  storing the history of 
%               beta_{lambda_{k}},k =1,...N 
% HisAClen:     a vector storing the history of acive set size of
%               beta_{lambda_{k}},k =1,...N 
% HisIter:      a vector storing  the iteration number when cd_j stop at
%               each  lambda_{k},k =1,...N
% Bic:          BIC value 
% ========================================================================
[n,p] = size(X);
if nargin < 2
    disp('error, not enough inputment ')
elseif  nargin > 3
    disp('error, too much inputment ')
elseif nargin == 2 || isempty(opts)
    opts.N = 100;
    opts.tau = 1e-8;
    opts.alpha = 0;
    opts.initial = zeros(p,1);
    opts.K = 50;
    opts.mu = 0.5;
    opts.step = norm(X);
    opts.tol = 1e-8;
end

% define the path Lambda
Xty = X'*y;
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 = norm(Xty,inf)*Lambda;
lenLam = length(Lambda);
Hisbeta = zeros(p,1);
HisAclen = [];
HisIter =[];
Bic = [];
Lam = [];
%% main loop for choosing lamda 
disp(' run CD-J on Lambda ...')
df = 0;
k = 0;
mun = opts.mu*n;
while df <mun && k < lenLam
         k = k + 1;
         lambda = Lambda(k);
          Lam = [Lam;lambda];
         [ebeta, AClen, AC, nIter] = cdj(X,y,lambda,opts);
         Hisbeta = [Hisbeta,ebeta];
         HisIter = [HisIter;nIter];
         df = AClen;
         HisAclen = [HisAclen;df];
         opts.initial = ebeta; % warm starting 
         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(1));
betahat = Hisbeta(:,id(1));
% % 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, AC, nIter] = cdj(X,y,lambda,opts)
% ========================================================================
% Coordinate decesent (Jacobi version, Iterative thresholding)  Algorithm 
% for min_{beta}{1/2||X beta-y||^2 + lambda||beta||_{1}
% Copyright (c), Des, 2013
% Yuling Jiao  Dept. Math.,Whuhan Univiversity
% Email: yulingjiaomath@whu.edu.cn
% ========================================================================
% Input:  
% X:            the covariate matrix (nomalized) in R^{n*p}
% y:            response vector in R^{n}
% lambda:       tuing parameter 
% opts.K:       maxiteration number in SSN for safe guard   (defaut: 50)
% opts.mu:      stop parameter when |df|>mu*n  (defaut: 0.5)
% opts.step:    step length that make convergence (default: norm(X))
% opts.initial: initial value, (default:0)
% opts.tol:     stop when norm(beta^k+1)/norm(beta^k+1-beta^k)<opts.tol 
%               (default,1e-8)
% Output: 
% ebeta:    an estimate of regression coefficient 
% AClen:    acive set size  of ebeta 
% AC:       acive set of ebeta 
% nIter:    the iteration number when CD-J stop
% ========================================================================
[n,p] = size(X);
if nargin < 3
    disp('error, not enough inputment')
elseif nargin > 4
    disp('error, too much inputment')
elseif nargin == 3 || isempty(opts)
    K = 50;
    mu = 0.5;
    betainitial = zeros(p,1);
    stoprule = 1e-8;
else
    K = opts.K;
    mu = opts.mu;
    betainitial = opts.initial;
    stepsize = opts.step;
    stoprule = opts.tol;
end
Xty = X'*y;
nIter = 0;
ebeta = zeros(p,1);
AC = [];
AClen = 0;
stop  = stoprule+1;
mun = mu*n;
Eps = 1e-8;
while stop > stoprule && nIter < K && AClen < mun
    nIter = nIter + 1;
    tempbeta = betainitial + (Xty - X'*(X*betainitial))/stepsize;
    AC = find(abs(tempbeta)>(lambda/stepsize));
    AClen = length(AC);
    ebeta = softth (tempbeta,lambda/stepsize);
    stop = norm(betainitial-ebeta)/(norm(betainitial)+ Eps);
    betainitial = ebeta;
end
end
%% subfunction 
function thz = softth(z,th)
    thz = max(abs(z)-th,0).*sign(z);
end

