function [x,stuff] = nnlsinv(b,TAij,sfactr)
% NNLS inverse converted to function call: 970314 -- cfg
%
% inputs:  b is vector of N Sv's
%          TAij is MxN matrix of TS's converted to intensities
%          sfactr is L-M parameter
%
% outputs: x is numerical abundance for M sizes
%          stuff is a M-vector containing
%              RNORM, XNORM, ANORM, and a lot of zeros


b_unscaled=b;

% Test b for valid data and to
% accordingly adjust the length of the b vector and the
% Aij matrix before doing the underdetermined inverse
% calculation

vpi = finite(b);               % Find indices for valid Sv values
[mfreqs,nsizes]=size(TAij);    % Size of input model matrix
stuff=zeros(1,nsizes);
if length(vpi)>2
   
   b = b(vpi);          % Include only valid Sv's
   TAijt = TAij(vpi,:); % Decrement TAij matrix (rows) to match length
   % of b vector
   
   b = (10 .^ (b/10))'; % Convert to linear form

%        Scale b-vector and A-matrix (row scaling)

[m n] = size(TAijt);

% Row scale b-vector
bscale = b;
b = b ./ bscale;

% Row scale TAijt
TAsij=[];
for im=1:m
  for jn=1:n
    TAsij(im,jn) = TAijt(im,jn) ./ bscale(im);
  end
end
%	Should this be sqrt(sum(sum(TAsij.*TAsij)))/(m*n);?

anorm = sqrt(sum(sum(TAsij .* TAsij)))/(m*n);
%anorm = sqrt(sum(sum(TAsij .* TAsij))/(m*n));

%          Append zeroes to b vector
b = [b;zeros(n,1)];

%          Compute and append minimization matrix
taij = [TAsij; eye(n,n)*anorm*sfactr];

tol=eps/2; % May be machine dependent

%          Compute inverse solution

x = nnls(taij,b,tol);

%     Compute Euclidean norms for residual and abundance vectors

rnorm = norm(taij(1:m,1:n) * x - b(1:m),'fro');
xnorm = norm(x,'fro');

stuff=zeros(1,n);stuff(1)=rnorm;stuff(2)=xnorm;stuff(3)=anorm;

else
   x=zeros(nsizes,1);
end
