%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
%
function [x,w] = nnls(A,b,tol)
% function [x,w] = nnls(A,b,tol)
% NNLS	Non-negative least-squares.
%
%	x = nnls(A,b) returns the vector X that solves x = pinv(A)*b
%	in a least squares sense, subject to x >= 0.
%	Differently stated it solves the problem min ||y - Xx|| if
%	A = X'*X and b = X'*y.
%
%	A default tolerance of TOL = MAX(SIZE(A)) * NORM(A,1) * EPS
%	is used for deciding when elements of x are less than zero.
%	This can be overridden with x = nnls(A,b,TOL).
%
%	[x,w] = nnls(A,b) also returns dual vector w where
%	w(i) < 0 where x(i) = 0 and w(i) = 0 where x(i) > 0.

% 	Adapted from NNLS of Mathworks, Inc.
%
%	L. Shure 5-8-87
%	Revised, 12-15-88,8-31-89 LS.
%	(Partly) Copyright (c) 1984-94 by The MathWorks, Inc.

%	Modified by R. Bro 5-7-96 according to
%       Bro R., de Jong S., Journal of Chemometrics, 1997, 11, 393-401
% 	Corresponds to the FNNLSa algorithm in the paper
%	
%	Modified by N. P. Pitsianis 9/23/04
%	to accomodate rectangular matrices,
%       cleanup the handling of the zero and non-zero index sets.

%  Reference:
%  Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974.

% initialize variables
if nargin < 3
  tol = 10*eps*norm(A,1)*max(size(A));
end

[m,n] = size(A);

% if m ~= n
  b = A'*b;
  A = A'*A;
% end

x = zeros(n,1);
r = b-A*x;

P = zeros(n,1);
z = zeros(n,1);
Z = 1:n;
ZZ = 1:n;

% set up iteration criterion
iter = 0;
itmax = 30*n;

% outer loop to put variables into set to hold positive coefficients
while any(Z) & any(r(ZZ) > tol)
  [rt,t] = max(r(ZZ));
  t = ZZ(t);
  P(t) = t;
  Z(t) = 0;
  PP = find(P);
  ZZ = find(Z);

  z(PP) = A(PP,PP) \ b(PP);  
  z(ZZ) = 0;

  % inner loop to remove elements from the positive set which no longer belong
  while any((z(PP) <= tol)) & iter < itmax
    
    iter = iter + 1;
    QQ = find((z <= tol) & P);
    alpha = min(x(QQ)./(x(QQ) - z(QQ)));
    x = x + alpha*(z - x);
    ij = find(abs(x) < tol & P);
    Z(ij) = ij';
    
    P(ij) = 0;
  
    PP = find(P);
    ZZ = find(Z);
    
    z(PP) = A(PP,PP) \ b(PP);
    z(ZZ) = 0;
  end
  x = z;
  r = b-A*x;
end