function [xMatrixOutT projectionIndexVector constraintMatrixList] = ...
    projectFromData(xMatrixT, projectionData, zerotol, checkCompletion, ...
    xMatrixOutT)
% project points based on precalculated projection data
%
% Syntax: xMatrix = projectOrthogonalFromData(xMatrix, projectionData)
%
% The input/output xMatrix contains vectors in each column. The
%   projectionData contains several possible projections as a cell array
%   .projectionList. It also contains a set of linear constraints in the
%   rows of projectionData.A for which .responsibleMatrix indicates which
%   constraint line (1st column) belongs to which projection (2nd column)
%   with a possible sign change (3rd column). Despite of the selected 
%   subset, the indicated projection is chosen by:
%      projectionData.A * xMatrix <= 0
%   with sign=1 and by:
%      projectionData.A * xMatrix > 0
%   with sign=-1.
%
% The additional argument is checkCompletion enable (=1) warnings when not
%   exactly one projection is selected. For checkCompletion=2, the warnings
%   are followed by an error.
%
% THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 
%   OF APPLICATIONS (decomposition case study; technical reports)
%
% See also: generateOrthogonalProjectionDataForCone

% 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/>.


%% Input handling
if ~exist('xMatrixOutT', 'var')
    xMatrixOutT = xMatrixT;
end
if size(xMatrixT, 1) ~= size(projectionData.projectionList{1}, 1) || ... 
   size(xMatrixOutT, 1) ~= size(projectionData.projectionList{1}, 1)
    error('elk:decomposition:wrongInput', ['Please, check the provided ' ...
        'xMatrix, the vectors are the columns of this matrix and must ' ...
        'equal the dimensions of the projection matrices']);
end
if ~exist('zerotol', 'var') || isempty(zerotol)
    zerotol = 0;
end
if ~exist('checkCompletion', 'var') || isempty(checkCompletion)
    checkCompletion = 0;
end

%% Initialize
% return vector for used projections
if nargout >= 2
    projectionIndexVector = zeros(1, size(xMatrixT, 2));
end
if nargout >= 3
    constraintMatrixList = cell(1, size(xMatrixT, 2));
end
% vector for checking whether it was projected properly
numberProjectionsVector = zeros(1, size(xMatrixT, 2));

%% Evaluate constraints
% determine which constraints are fulfilled, and which not - this CAN be
%   based on a second matrix, when provided. This behavior is required for
%   the grwoth rate limitations.
constraintMatrix = projectionData.A*xMatrixT;

% constraints are evaluated in favour of the negative halfspace; zerotol
%   must be applied because the demand of a projection already uses this
%   zerotol.
%! DEV: We change this decision such that a double-projection is possible -
%    but only in regions where the projection should be almost the same,
%    anyway.. ..in the scope of zerotol.
% old code:
% constraintMatrixNegative = constraintMatrix <= zerotol;
% constraintMatrixPositive = constraintMatrix > zerotol;
% new code:
constraintMatrixNegative = bsxfun(@le, constraintMatrix, zerotol);
constraintMatrixPositive = bsxfun(@ge, constraintMatrix, -1*zerotol);
%/

%% Find projection
% To execute the projection, it is necessary to cycle either through the
%   vectors or through the projections. The number of projections should
%   typically be less in applications.
for iProjection = 1:length(projectionData.projectionList)
    thisConstraintIndexVector = projectionData.responsibleMatrix(...
        projectionData.responsibleMatrix(:, 2) == iProjection, 1);
    thisConstraintSignVector = projectionData.responsibleMatrix(...
        projectionData.responsibleMatrix(:, 2) == iProjection, 3);
    % the projection is valid, when all constraints are fulfilled
    thisFilterMatrix = ...
        (bsxfun(@or, thisConstraintSignVector < 0, ...
                constraintMatrixNegative(thisConstraintIndexVector, :)) & ...
         bsxfun(@or, thisConstraintSignVector > 0, ...
                constraintMatrixPositive(thisConstraintIndexVector, :)));
    thisFilterVector = all(thisFilterMatrix, 1);
    
    if ~isempty(thisFilterVector)
        xMatrixOutT(:, thisFilterVector) = ...
            projectionData.projectionList{iProjection} * ...
            xMatrixOutT(:, thisFilterVector);
        if nargout >= 2
            projectionIndexVector(thisFilterVector) = iProjection;
        end
        if nargout >= 3
            [constraintMatrixList{thisFilterVector}] = deal(...
                bsxfun(@times, ...
                  projectionData.A(thisConstraintIndexVector, :), ...
                  thisConstraintSignVector(:)));
        end
        numberProjectionsVector(thisFilterVector) = ...
          numberProjectionsVector(thisFilterVector) + 1;
    end
end

%% error handling
flagError = 0;
if checkCompletion
    %! Info: any of the following warnings / errors is only thrown with
    %    respect to validity mapping. Otherwise, the "should not happen" is
    %    wrong.
%     if max(numberProjectionsVector) > 1
%         warning('elk:decomposition:general', ['A vector was ' ...
%             'projected more than once. This should not happen.']);
%         flagError = 1;
%     end
    if min(numberProjectionsVector) < 1
        warning('elk:decomposition:general', ['A vector was ' ...
            'not projected. This should not happen.']);
        flagError = 1;
    end
    if checkCompletion == 2 && flagError
        error('elk:decomposition:numericalError', 'See the warning above.');
    end
end