% function [Xsub,idx]=licols(X,tol,mode)
% Matt J (2021). Extract linearly independent subset of matrix columns 
%(https://www.mathworks.com/matlabcentral/fileexchange/77437-extract-linearly-independent-subset-of-matrix-columns), MATLAB Central File Exchange. Retrieved June 11, 2021.
%
% Significant modifications by Ethan Liang
%
%Extract a linearly independent set of columns of a given matrix X
%
%    [Xsub,idx]=licols(X)
%
%Inputs:
%
%  X: The given input matrix
%  tol: A rank estimation tolerance. Default=1e-10
%  mode: Rule for choice of columns
%          'desc': choose columns based on descending values of R
%          'LR': choose columns from left to right
%
%out:
%
% Xsub: The extracted columns of X
% idx:  The indices (into X) of the extracted columns
function [Xsub,idx]=licols(X,tol,mode)
     if ~nnz(X) %X has no non-zeros and hence no independent columns
         Xsub=[]; idx=[];
         return
     end
     
     if (nargin<3)
         mode = 'LR';
         if (nargin<2)
             tol=1e-10; 
         end
     end
     
     
     
     if (strcmpi(mode, 'desc'))
         [Q, R, E] = qr(X,'vector'); 
         if ~isvector(R)
             diagr = abs(diag(R));
         else
             diagr = abs(R(1));   
         end

           %Rank estimation
           r = find(diagr >= tol*diagr(1), 1, 'last'); %rank estimation
           idx=sort(E(1:r));
           Xsub=X(:,idx);
     elseif (strcmpi(mode, 'LR'))
         
         % Although the above code works, this selection prioritizes
         % selecting columns from left to right, which is more commonly
         % used in the examples in the textbook.
         [Q, R] = qr(X); 
         if ~isvector(R)
             diagr = abs(diag(R));
         else
             diagr = abs(R(1));   
         end
           
         idz = 1;
         idx = [];
         for idy = 1:size(X, 2)
             if (abs(R(idy, idz)) >= tol*diagr(1))
                 idx = [idx, idy];
                 idz = idz + 1;
             end
         end
         Xsub=X(:,idx);
         
     else
         error("Incorrect mode selected: " + string(mode));
     end
end
       
       