function betas=lars(X,y, max_steps, tracep) 
% LARS Least Angle Regression 
% lars(X,y, max_steps) returns the coefficients at each step in LARS regression 
% INPUT: X is the design matrix of n by m, 
%          where n is the number of examples and m is the number of features 
%        y is the response vector of n by 1 
%        max_steps is max number of steps for stopping early
%        tracep is true/false for verbose display
% OUTPUT:
%        betas contains the regression coefficient at each step as 
%          columns. It is usually an m by min(m,max_steps) matrix.
% 
% Reference:
% Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani.
% Least angle regression
% Source: Ann. Statist.  32 (2004), no. 2, 407--499
% http://www-stat.stanford.edu/~hastie/pub.htm

  [n,m]=size(X); 
  assert(n>m, 'more observations than covariates');
  y=y(:); assert(length(y)==n, 'X is n x m and Y is of length n');
  
  eps_tol=1E-10;
  assert(max(abs(sum(X,1)))<eps_tol, 'zero mean covariates')
  assert(max(abs(std(X,0,1)-1))<eps_tol, 'covariates must have std dev 1');
  assert(abs(sum(y))<eps_tol, 'zero mean response');
  assert(abs(std(y,0)-1)<eps_tol, 'for stability the response better be normalizedto have std dev 1');
  
  if (~exist('max_steps','var') | isempty(max_steps)), max_steps=m; end;
  if (~exist('tracep', 'var') | isempty(tracep)), tracep=false; end;

  assert(max_steps<=m, 'max steps <= m');
  assert(islogical(tracep), 'tracep is a predicate true/false');
  
  mu_hat=zeros(n,1); % the current response on a subset of features
  beta=zeros(m,1); % the initial beta is all zero

  %%% start of main loop
  k=0;
  hit_all=false;  % whether we've hit all the covariates
  R=[];
  Acurl_old=[];
  while (k<max_steps & ~hit_all)
    k=k+1;
    C=X'*(y-mu_hat); % C is the current correlation
    C_hat=max(abs(C));  % the max of the correlations
    Acurl=find(abs(C)>C_hat-eps_tol); % the set of covariates that are maximally correlated
    cov_added = setdiff(Acurl, Acurl_old); cov_added=cov_added(:)';
    Acurl=[Acurl_old cov_added]; % so new covariates would always be appending, necessary for update Cholesky

    s=sign(C(Acurl)); s=s(:);% get the sign of the maximal correlations
    if (isempty(R)) % first time Cholesky decomposition
      XA=X(:,Acurl);
      try
        R=chol(XA'*XA);
      catch
        error(['LARS Variable ' sprintf(' %i', Acurl) ' are collinear ']);
      end
    else % otherwise build Cholesky incrementally
      Rout=updateCholXtX(R, X(:,Acurl_old), X(:,cov_added));
      R=Rout;
    end
% check for linear redundancy in data
    if size(R,2)==1, diagR=R(1,1); else diagR=diag(R); end;
    if (min(abs(diagR))<eps_tol)
      error(['LARS Variable ' sprintf(' %i', Acurl) ' are collinear ']);
    end
    
    Acurl_old=Acurl;
    if (tracep), 
      fprintf('LARS Step %i\tVariable', k);
      fprintf(' %i', cov_added);
      fprintf(' added\n');
    end
%   if (isempty(cov_added)), keyboard; end
    
    beta_inc=R\(R'\s);  % the increment in beta
    u=X(:,Acurl)*beta_inc; % the vector that makes equal angle to
% the columns of XA
    a=X'*u;  % the vector of inner products
    if (length(Acurl)==m) % all covariates are active
      gamma=C_hat;
      hit_all=true;
    else
      Acurl_c=setdiff([1:m], Acurl);  % set of non-maximally
% correlated covariates
      C=C(Acurl_c); a=a(Acurl_c); 
      gammas=[(C_hat-C)./(1-a); (C_hat+C)./(1+a)]; % compute the values where a new covariate would join the set of maximally correlated covariates
      gammas=gammas(gammas>eps_tol); % take only positive gammas
      gamma=min(gammas); % move to the soonest event
    end
    mu_hat=mu_hat+gamma*u; % update the current response
    beta(Acurl)=beta(Acurl)+gamma*beta_inc; % update beta
    
    betas(:,k)=beta; % record beta
    
  end

  % make sure the last beta is the Ordinary Least Squares solution
  if hit_all,
    assert(abs(betas(:,end)-X\y)<eps_tol*m, 'last beta is the OLS solution'); 
  end

      
      
% called as
%      R=updateCholXtX(R, X(:,Acurl_old), X(:,cov_added));
function Rout=updateCholXtX(R, Xold, xnew)
% Let Rout=[R v; 0 alpha], X = [Xold xnew], xnew could be multi-cols
% equate the blocks of Rout'*Rout=X'*X and it gives
%  assert(norm(R'*R-Xold'*Xold)<1E-5, 'precondition of updateChol')
  v=R'\(Xold'*xnew);
  alpha=chol(xnew'*xnew-v'*v);
  n1=size(R,1); n2=size(alpha,1);
  Rout=[R v;  zeros(n2,n1) alpha];
%  X=[Xold xnew];
%  assert(norm(Rout'*Rout-X'*X)<1E-5, 'postcondition of updateChol');