function [b,idxfit]=NonNegativeLS(S,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General description:
% Multiple linear regression with nonnegative constraints. 
%
%       y = Sb,  sub: b>=0
%
% The algorithm is on paper:
% M. Tenenhaus, "Canonical analysis of two convex polyhedral cones and 
% applications," Psychometrika, vol. 53, no. 4, pp. 503--524.
%
% input:
%       S - subspace
%       y - target vector
%
% output: 
%       b - coefficients of contributions of all components. b_i=0 if no
%           contribution
%       idxfit - indices of those positive coefficents
%
% Written by: Wei Wang
% Machine Learning for Signal Processing Lab
% University of Maryland, Baltimore County
% 1/30/2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,p]=size(S);

%%%%% Step 1
ifit=zeros(p,1); % flag of all the elements in fitting model],
                 % e.g, [0 0 1 0] means the 3rd entered model
idxfit=find(ifit==1);
b=zeros(p,1);
a=zeros(p,1);
bt0=zeros(p,1);
at0=zeros(p,1);
count=0;
while 1
    count=count+1;
    if count>200
        disp('Not converge yet in NonNegativeLS!');
        break;
    end
    
    %%%%% Step 2
    R=S(:,idxfit);
    b=R\y;
    %% All elements enter model
    if length(idxfit)==p
        break;
    end
    
    %% (y-y_R)Rc are all nonpositive
    idxRc=find(ifit==0);
    Rc=S(:,idxRc);
    if length(idxfit)==0
        t=y'*Rc;
    else
        t=(y-R*b)'*Rc;
    end
    if sum(sign(-t))==length(t)
        break;
    end
    
    %%%% Step 3
    imax=idxRc(find(t==max(t))); % imax: the index of the maximum of t    
    idxfitpre=idxfit;
    a=b;
    ifit(imax)=1;     % imax enter the model
    idxfit=find(ifit==1);  % idxfit: the index of those element entered model
    if length(idxfit)==1
        continue;
    end
    
    %%%% Step 4
    y_R=R*b;
    count_inner=0;
    while 1
        if count_inner>200
            disp('Not converge yet in NonNegativeLS! (inner loop)');
            break;
        end
        count_inner=count_inner+1;        
        idxfit=find(ifit==1);  % idxfit: the index of those element entered model
        R=S(:,idxfit);        
        b=R\y;
        if sum(sign(b))==length(b)
            break;        
        end
    
    %%%% Step 5
        bt=bt0;
        at=at0;
        bt(idxfit)=b;
        at(idxfitpre)=a;
        idxnegb=idxfit(find(b<=0));    
        tnegb=bt(idxnegb)./(bt(idxnegb)-at(idxnegb));
        lambda=max(tnegb);
        
        gamma=bt(idxnegb)+lambda*(at(idxnegb)-bt(idxnegb));
        
        idxfitpre=idxfit;
        a=bt(idxfitpre)+lambda*(at(idxfitpre)-bt(idxfitpre));
        
        maxlpos=idxnegb(find(gamma==0));
        ifit(maxlpos)=0;      
        idxfit=find(ifit==1);
                
    end
end

bt=zeros(p,1);
bt(idxfit)=b;
b=bt;
