function [x,   		  					...
          SquaredNormOfResidualVector,	...
          ResidualVectorbMinusMx,  		...
          exitflag] = 					...
          lsqnonnegGJK(M,b)	;
% Purpose		Solve the non-negative least squares problem Mx=b, where each entry of x is
%				zero or positive.  Use the GJK-based algorithm introduced in [Centore2015].
%
% Description	The non-negative least squares (NNLS) problem is to find a vector x, such
%				that the norm of Mx-b is minimized, where M is a known matrix, and b is a
% 				known vector; there is a requirement, however, that each entry of x be
%				zero or positive.  If the minimum norm is 0, then x satisfies Mx=b exactly,
%				if not, then x can be thought of as a vector which comes as close as possible
%				to satisfying Mx=b.
%
%				This routine uses an algorithm introduced in [Centore2015].  That algorithm
%				is based on the Gilbert-Johnson-Keerthi (GJK) algorithm that was introduced 
%				in 1988 ([Gilbert1988]).  Readable descriptions appear in [Rabbitz1994] (although 
%				the term GJK algorithm is not used there), [Ericson2004], and [McCutchan2006].
%				In its general form, the algorithm finds the minimum distance between two
%				convex polytopes.  We actually need only a special case in which one polytope
%				is a point.  
%
%				The algorithm in [Centore2015] considers the term Mx in the NNLS equation as
%				a set in a vector space.  Since each x is non-negative, Mx is the convex cone
%				generated by the columns of M.  b can be thought of
%				as a target vector.  The least squares problem finds the point on the convex
%				cone that is closest to b.
%
%				Rather than dealing with the cone as a whole, it is truncated at a certain
%				distance Gamma.  The truncated cone is the convex polytope P_Gamma generated by 
%				(adjusted) column vectors of M, and the zero vector.  The GJK algorithm
%				can then be used to find the nearest point on P_Gamma to b, and the shortest
%				distance.  As Gamma increases, the distance (as a function of Gamma) will
%				reach a plateau, and the nearest point will be identical for any larger
%				Gammas.  When this happens, that point will also be the minimizing point
%				for the entire convex cone.  This routine starts off with a small Gamma,
%				and increases it by a factor of 10 (up to some maximum), until the convergence
%				is noted, at which time it is concluded that a solution has been found.
%
%				This novel algorithm is presented here because it seems to have the 
%				potential to improve on current algorithms, both because the GJK algorithm is
%				very fast, and  because the calculation of the necessary convex polytopes is 
%				almost immediate from the matrix M.  No speed tests or formal comparisons
%				have been made with other algorithms, however, so any improvements are, for
%				now at least, just promising hypotheses.  The new algorithm does have the
%				advantage of being easy to understand and interpret.   Some issues with
%				numerical stability have also been noted, which might cause problems for any
%				algorithm.  Attempts were made to mitigate instability, but doubtless more
%				sophisticated approaches are available.  Other researchers with
%				more experience in NNLS problems are encourage to use, test, and modify this
%				open source code, in the process of which the algorithm s usefulness can be
%				assessed.
%
% References	[Centore2015] Paul Centore, A GJK-Based Algorithm for Constrained Least Squares, 
%								available at wwww.99main.com/~centore, 2015
%				[Gilbert1988] 	Elmer G. Gilbert, Daniel W. Johnson, & Sathiya Keerthi, A Fast
%								Procedure for Computing the Distance Between Complex Objects in
%								Three-Dimensional Space, IEEE Journal of Robotics and Automation,
%								Vol. 4, No. 2, April 1988, pp. 193-203.
%				[Rabbitz1994]	Rich Rabbitz, Fast Collision Detection of Moving Convex Polyhedra,
%								in Section I.8 of Graphics Gems IV (IBM Version), ed. Paul Heckbert,
%								Academic Press, 1994.
%				[Ericson2004]	Christer Ericson, The Gilbert-Johnson-Keerthi Algorithm, 
%								http://realtimecollisiondetection.net/pubs/SIGGRAPH04_Ericson_GJK_notes.pdf
%				[McCutchan2006]	John McCutchan, Introduction to GJK, 9 November 2006,
%								http://www.cas.mcmaster.ca/~carette/SE3GB3/2006/notes/gjk1_pres.pdf
%
% Variables		M,b				Entries in the NNLS problem Mx=b.  M is a matrix with
%								m rows and n columns, and b is a column vector with m entries.
%
%				x			    A column vector of n entries, that minimizes the norm of Mx-b.
%
%				SquaredNormOfResidualVector		The squared norm of the residual vector b-Mx.			
%
%				ResidualVectorbMinusMx		The residual vector, given by b-Mx, where x is the
%								non-negative vector which minimizes the norm of Mx-b.
%
%				exitflag		A performance flag that is 0 when the algorithm fails to 
%								produce a result, and a positive integer otherwise.
%
% Author		Paul Centore (January 8, 2015)
%
% Copyright 2015 Paul Centore
%
%    This file is part of MunsellAndKubelkaMunkToolbox.
%
%    MunsellAndKubelkaMunkToolbox is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    MunsellAndKubelkaMunkToolbox is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with MunsellAndKubelkaMunkToolbox.  If not, see <http://www.gnu.org/licenses/>.

% Initialize variables
x                            = []	;
SquaredNormOfResidualVector  = NaN	;
ResidualVectorbMinusMx       = []	;
exitflag                     = 0	;

% Minimum distance for some comparisons, to avoid numerical roundoff error
EpsDist = 1e-12	;	

[Rows, Cols] = size(M)	;

% Loop through a set of truncated cones, whose size increases by a factor of 10 at each
% iteration.  Apply the GJK algorithm to find shortest distance from the truncated cone to
% the target vector b.  When the distance stops decreasing, the constant distance is the
% distance to the convex cone Mx, and the GJK algorithm will have found the closest point
% on Mx to the target b; this is (or leads to) an NNLS solution.
SolutionFound = false	;
PowerOfTen    = 0		;
MaxPowerOfTen = 15		;
% The following matrices keep track of information over each iteration.  The ith row of
% each matrix corresponds to the ith truncated cone.
ClosestPoints    = []	;
Distances        = []	;
PointInPolytopes = []	;
CoefficientsList = []	;
NormsOfMxMinusb  = []	;
while not(SolutionFound) && PowerOfTen <= MaxPowerOfTen

	% A truncated cone is the convex hull of the columns of M and the zero vector; the
	% truncated cone is scaled by a factor which increases tenfold at each iteration.  
	PolytopeGenerators = [transpose(M); zeros(1,Rows)]			;
	PolytopeGenerators = (10^PowerOfTen) * PolytopeGenerators	;

	% Apply the GJK algorithm to a truncated cone, to find the nearest point to b on
	% that truncated cone, and related information
	[ClosestPoint,   ...
	 Distance,       ...
	 PointInPolytope,...
	 Coefficients] = ...
	 ClosestPointInConvexPolytopeGJK(...
							PolytopeGenerators,...
							reshape(b,1,length(b)) );
	% The variable Coefficients give the coefficients in a convex combination of the
	% polytope generators (the zero vector and the scaled columns of M) that produces the
	% nearest point to b on the truncated cone.  They must be rescaled to give similar
	% coefficients in terms of the original columns of M,							
    Coefficients = Coefficients*(10^PowerOfTen)	;		
    
    % Record the information for the new truncated cone as a new row in the information matrices
    MxMinusb         = M *reshape(Coefficients(1:(end-1)),length(Coefficients)-1,1) - reshape(b,length(b),1)	;
    NormsOfMxMinusb  = [NormsOfMxMinusb; norm(MxMinusb)]								;
	ClosestPoints    = [ClosestPoints; ClosestPoint]									;		
	Distances        = [Distances;	Distance ]											;					
	PointInPolytopes = [PointInPolytopes; PointInPolytope]								;		
	CoefficientsList = [CoefficientsList; reshape(Coefficients,1,length(Coefficients))]	;	
	
	% Check whether the point b is inside or on a truncated cone, in which case it is also 
	% in the convex cone Mx.  If so, then a solution has been found.
	if PointInPolytope
		SolutionFound = true	;
	% If not, then see whether the distance to the truncated cone has stopped decreasing.
	% If it has, then there will be no improvement as the truncated cones increase in
	% size, and we have found the nearest point
	elseif length(Distances) >= 2
		if abs(Distances(end) - Distances(end-1)) < EpsDist
			SolutionFound = true	;
		end
	end		
	
	% Make next truncated cone 10 times the size of the current one
	PowerOfTen = PowerOfTen + 1	;	
end

% If a solution has been found, then calculate the returned values, such as the residual
% vector and its norm.
if SolutionFound
	% The next sorting step was introduced because some numerical instabilities were
	% seen.  Since each truncated cone properly contains the previous truncated cone, the
	% minimum distance to the point b should constantly be decreasing.  The fact that it
	% is not indicates numerical problems.  To try to mitigate instability when it occurs,
	% check directly for the cone with the smallest distance, and use information from
	% that cone.
	[SortedNorms, SortedIndices] = sort(NormsOfMxMinusb)									;
	x 		                     = CoefficientsList(SortedIndices(1),1:(end-1))				;
	ResidualVectorbMinusMx       = reshape(b,length(b),1) - (M * reshape(x,length(x),1)) 	;
	SquaredNormOfResidualVector  = (Distances(SortedIndices(1)))^2							;
	exitflag                     = 1														;
end