function [x,   		  ...
          SquaredNormOfResidualVector,    ...
          ResidualVectorbMinusMx,   ...
          exitflag] = ...
          BoundedVariablesLeastSquares0to1GJK(M,b)	;
% Purpose		Use the GJK algorithm to solve the bounded variables least squares problem 
%				Mx=b, where each entry of x is greater than or equal to 0 and less than or equal to 1.
%
% Description	The bounded variables least squares (BVLS) problem is to find a vector x that 
%				comes closest to satisfying Mx=b, where M is a known matrix of m rows and
%				n columns, and b is a known column vector with n entries.  BVLS constrains
%				x so that each entry of x has an upper and lower bound.  By a suitable scaling
%				and translation, one can always find an equivalent BVLS problem such that 
%				each entry of x is greater than or equal to 0 and less than or equal to 1.
%
%				The algorithm behind this routine is described in [Centore2015].  Rather than
%				more common approaches involving quadratic programming (QP), the routine
%				implements the Gilbert-Johnson-Keerthi (GJK) algorithm ([Gilbert1988],
%				[Rabbitz1994]), which finds the minimal
%				distance between two convex polytopes.  The routine uses a special case of
%				the algorithm, in which the minimal distance is found between a point and a
%				convex polytope.  The approach in [Centore2015] considers Mx as a subset of
%				the vector space R^m.  That subset is determined by both the columns of M 
%				and the constraints on x.  With the given constraints, Mx is the zonotope
%				generated by the columns of M, where each column is treated as a vector in
%				R^m.  The target point b is also in R^m.  Since a zonotope is a convex
%				polytope, the GJK algorithm can find the minimum distance between b and the
%				zonotope, GJK also finds the point b_Z on the zonotope that is closest to b,
%				and an expression for b_Z as a convex combination of the columns of M.  This
%				information is sufficient to work backwards to a minimizing vector x.  See
%				[Centore2015] for further details.
%
%				The GJK approach is proposed here as a BVLS algorithm that has the potential
%				to improve on current algorithms, but whose effectiveness has not been
%				tested.  The GJK approach has the advantage of being easy to understand and
%				implement, but has the disadvantage of using 2^n points as a generating set
%				for the zonotope.  The GJK algorithm speed depends approximately linearly on
%				the size of the generating set, so the GJK approach could be exponential in
%				the number of entries in the vector of x.  On the other hand, there is no
%				need to solve the zonotope (i.e. determine its vertices), not to express it
%				in terms of half-space inequalities.  If there are not too many variables, then
%				the GJK approach might be preferable.  In any event, the algorithm is presented
%				as open source so that other researchers may make tests and comparisons as they 
%				wish, and also add any desired improvements.
%
% 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.
%
%				M,b				Entries in the BVLS 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.
%
%				ResidualVectorbMinusMx		The residual vector, given by b-Mx, where x is the
%								vector which minimizes the norm of Mx-b.
%
%				SquaredNormOfResidualVector		The squared norm of the residual vector b-Mx.			
%
%				ErrorF	Each vertex is a linear combination of the generating vectors.
%									The ith row of the matrix VertexCoefficients gives the
%									coefficients of the generating vectors for the ith vertex, which
%									is given in the ith row of the matrix Vertices.
%
% Author		Paul Centore (January 25, 2015)
%
% Copyright 2014 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                     = 1	;

% Extract information about the number of rows and columns of M
[m,n]        = size(M)			;

% The index matrix IndexJ expresses a generating set for the zonotope P (a generating set
% for P is a set V of vectors whose convex hull is P) as a set of linear combinations of
% the columns of M.  A node of a zonotope is a linear combination of the vectors that
% generate the zonotope, where every coefficient in the combination is either 0 or 1. (The
% generators of the zonotope should not be confused with a generating set for the zonotope.
% Every point in the zonotope is a convex combination of vectors in the generating set, but
% many points are not convex combinations of the generators.)  The index matrix IndexJ will
% therefore consist of all 2^n possible sequences of 0's and 1's.  Each row gives one such
% sequence.
[IndexJ,ind] = combn([0 1],n)	;

% The generators V of P are linear combinations of the columns of M, with coefficients
% given by the index matrix J.  Calculate these combinations.  The result will be a
% matrix PolytopeGeneratorsV.  Each row of PolytopeGeneratorsV is one of the generators 
% of P. 
NumberOfGenerators  = rows(IndexJ)	;
PolytopeGeneratorsV = NaN * ones(NumberOfGenerators,m)	;
for ctr = 1:NumberOfGenerators
	% Calculate a linear combination of columns of M, term by term
	TempGenerator = zeros(1,m)	;
	for ColumnCtr = 1:n
		TempGenerator = TempGenerator + ...
						IndexJ(ctr,ColumnCtr) .* transpose(M(:,ColumnCtr))	;
	end
	PolytopeGeneratorsV(ctr,:) = TempGenerator	;
end

% The constrained least squares solution to Mx=b corresponds to the point on the 
% zonotope P that is nearest to the origin.  Call the GJK algorithm to calculate
% this point.  The return variable alphaCoefficients gives the coefficients in a convex
% combination of the generators of P; this convex combination is the closest point b_Z
% on P to b.
[ b_Z,   				...
  Distance,       		...
  PointInPolytope,		...
  alphaCoefficients,  	...
  ErrorFlag] = 			...
  ClosestPointInConvexPolytopeGJK(				...
						PolytopeGeneratorsV,	...
						reshape(b,1,length(b)))	;
					
% If the GJK algorithm did not produce a solution, then relay this failure by setting the
% exit flag of the current routine to 0, and return						
if ErrorFlag
	exitflag = 0	;
	return
end						
						
% [Centore2015] shows that the entries of x given by dot products of the coefficients alpha
% and the columns of the index matrix J give a BVLS solution that satisfies the constraints. 						
x = NaN * ones(n,1)	;
for ctr = 1:n
	x(ctr,1) = dot(alphaCoefficients, IndexJ(:,ctr))	;
end		

% The residual vector is the vector between the target point b and the point b_Z on the 
% zonotope which is closest to b.  To be consistent with the returned solution, b_Z is
% calculated as Mx rather than using the GJK output.  Following Octave/MATLAB conventions
% for lsqnonneg, this routine also returns the squared norm of the residual vector
ResidualVectorbMinusMx      = reshape(b,length(b),1) - (M * reshape(x,length(x),1)) ;
SquaredNormOfResidualVector = (norm(ResidualVectorbMinusMx))^2						;