function [x,w] = nnls(E,f,tol)
%NNLS   Non-negative least-squares.
%   NNLS has been replaced with LSQNONNEG.  NNLS currently works but  
%   will be removed in the future.  Use LSQNONNEG instead.
%
%   X = NNLS(A,b) returns the vector X that minimizes NORM(A*X - b)
%   subject to X >= 0.   A and b must be real.
%
%   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
%   when x(i) = 0 and w(i) is approximately 0 when x(i) > 0.
%
%   See also LSCOV, SLASH.

%   L. Shure 5-8-87
%   Revised, 12-15-88,8-31-89 LS.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.7 $  $Date: 1998/08/31 23:09:41 $

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

if nargin<2, error('Not enough input arguments.'); end
if ~isreal(E) | ~isreal(f), error('A and b must be real.'); end

% initialize variables
if nargin < 3
    tol = 10*eps*norm(E,1)*length(E);
end
[m,n] = size(E);
P = zeros(1,n);
Z = 1:n;
x = P';
ZZ=Z;
w = E'*(f-E*x);

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

% outer loop to put variables into set to hold positive coefficients
while any(Z) & any(w(ZZ) > tol)
    [wt,t] = max(w(ZZ));
    t = ZZ(t);
    P(1,t) = t;
    Z(t) = 0;
    PP = find(P);
    ZZ = find(Z);
    nzz = size(ZZ);
    EP(1:m,PP) = E(:,PP);
    EP(:,ZZ) = zeros(m,nzz(2));
    z = pinv(EP)*f;
    z(ZZ) = zeros(nzz(2),nzz(1));
% inner loop to remove elements from the positive set which no longer belong
    while any((z(PP) <= tol))
        iter = iter + 1;
        if iter > itmax
            error(['Iteration count is exceeded.', ...
                   '  Try raising the tolerance.'])
        end
        QQ = find((z <= tol) & P');
        alpha = min(x(QQ)./(x(QQ) - z(QQ)));
        x = x + alpha*(z - x);
        ij = find(abs(x) < tol & P' ~= 0);
        Z(ij)=ij';
        P(ij)=zeros(1,length(ij));
        PP = find(P);
        ZZ = find(Z);
        nzz = size(ZZ);
        EP(1:m,PP) = E(:,PP);
        EP(:,ZZ) = zeros(m,nzz(2));
        z = pinv(EP)*f;
        z(ZZ) = zeros(nzz(2),nzz(1));
    end
    x = z;
    w = E'*(f-E*x);
end
