function heMatrix = closestProjectionOntoConfinementPartition(...
    heMatrix, x0, A, projectionMatrix)
% find the closest point in some subdomain
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% This program 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.
%
% This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.


%% notes
% Given two points x and xOpt, there distance (without sqrt) is:
%     (xOpt - x)'*(xOpt - x)
%   giving a necessarily quadratic problem. The quadprog minimizes:
%     1/2*x'*H*x + f'*x
%   where H is a positive definite matrix. Reformulating the distance,
%   gives:
%     xOpt'*xOpt - 2*xOpt'*x + x'*x
%   for which we can ignore the last term and divide by two to obtain:
%     H = 1
%     f = -1*x

nE = size(heMatrix, 1);
nVec = size(heMatrix, 2);

problem.solver = 'quadprog';
problem.options = optimset('quadprog');
problem.options.Algorithm = 'active-set';
problem.options.Display = 'off';
problem.H = eye(nE);
problem.Aineq = A;
problem.bineq = zeros(1, size(A, 1));
problem.Aeq = eye(nE) - projectionMatrix;
problem.beq = zeros(nE, 1);
problem.x0 = x0;

problem.f = -1*heMatrix(:, 1);
objective = @(x)(1/2*x'*problem.H*x + problem.f'*x);
objective(x0)
objective(heMatrix(:,1))

for iVec = 1:nVec
    iVec/nVec*100
    problem.f = -1*heMatrix(:, iVec);
    heMatrix(:, iVec) = quadprog(problem);
%     [heMatrix(:, iVec),fval,exitflag,output,lambda] = quadprog(problem)
%     lambda.lower,lambda.upper,lambda.eqlin,lambda.ineqlin
end
    