clc
clear all 
close all
addpath(genpath(fileparts(mfilename('fullpath'))));
warning off 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate data standanlized and normlized 
p = 1000;            % number of features 
n = 400;             % number of observation 
T = 10;              % number of non zero elemet in  real solution 
cor = 0.5;           % corelation number of columns in X 
ratio = 1;           % range of value in betatrue = 10^ratio;
isnorm = 1;          % 1 normlized and stantalized; 0 not
seednum = 0;         % the seed number 
sigma = 0.1;         % varince of  noise 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generate data 
[X,D,y,ye,betatrue,supp] = gendata(n,p,T,sigma,cor,ratio,seednum,isnorm);
%% get solution path and selecte the optimal solution 
tic,
[betahat, lambdahat, Hisbeta, HisAclen,  Bic, Lam] = glmnet_path(X,y);
if isnorm ==1 
   Hisbeta = D*Hisbeta;
   betahat = D*betahat;
end
id  = find(Bic == min(Bic));
lambdahat = D(id(end),id(end))\lambdahat;
toc
%% compute the error and  number of iterations used to get the whole path  
rel2error = norm(betatrue - betahat)/norm(betatrue);
abslinferror = norm(betatrue - betahat,inf);
esupp = find(betahat);
display(sprintf('regularization parameter = %g, nonzero number is = %g, ralative l_2 error = %g,absolute l_inf error = %g',lambdahat,length(esupp),rel2error, abslinferror))
%% plot result 
figure(1), plot(1:p,betatrue','ko',1:p,betahat','r*'), title('target (o) and estimated (*)'),
figure(2),
subplot(1,2,1), plot(HisAclen,'r+','LineWidth',1),ylabel(' ||\beta(\lambda_{k})||_{0} of Glmnet at each \lambda_{k}'),
xlabel('k,  the idex of   \lambda_{k} ')
subplot(1,2,2), semilogy(Bic,'go','LineWidth',1),ylabel(' BIC value of Glmnet at each \lambda_{k}'),
xlabel('k,  idex of   \lambda_{k} ')
figure(3)
for k = 1:p
         plot(Hisbeta(k,:),'.','LineWidth',6), hold on,
         plot(Hisbeta(k,:),'-','LineWidth',0.1), hold on 
end
% show the selected solution on the path by vetical line 
id =  find(Bic == min(Bic));
xid = id(end)*ones(p,1);
vmin = min(min(Hisbeta));
vmax = max(max(Hisbeta));
yid =  1.1*linspace(vmin,vmax,p);
plot(xid,yid,'r-','LineWidth',2)
hold off
xlabel('k,  the idex of   \lambda_{k} ')
ylabel( 'coefficient  ')
 title('Glmnet  Path');




