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 
mu = 0.5;
disp('Lars is running')
tic
[Hisbeta, numIters, activationHist, duals,Lam] = SolveLasso(X, y, mu*n);
% select betahat from sols by BIC
[n1,n2] = size(Hisbeta);
Bic = [];
for k = 1:n2
    tempbeta = Hisbeta(:,k);
    id = find(tempbeta);
    Bick = 0.5*norm(X*tempbeta-y)^2 + log(p)*log(n)*length(id);
    Bic = [Bic;Bick];
end
id = find(Bic == min(Bic));
betahat = Hisbeta(:,id(1));
lambdahat = Lam(id(1));
% % debias
% supp = find(betahat);
% Xsupp =  X(:,supp);
% betahat(supp) = (Xsupp'*Xsupp)\(Xsupp'*y);
if isnorm ==1 
   Hisbeta = D*Hisbeta;
   betahat = D*betahat;
   id  = find(Bic == min(Bic));
   lambdahat = D(id(end),id(end))\lambdahat;
end
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),
semilogy(Bic,'go:','LineWidth',1),ylabel(' BIC value of Lars 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(1)*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('Lars  Path');















