function [hcMatrixOut filterToZero] = mapToValidityCone(hcMatrix, dec, ...
    zerotol, hcMatrixOut, applyConstraintFilter)
% map a hC vector to the cone of valid polytopes
%
% Syntax: hcMapped = mapToValidConeDec(hcMatrix, dec)
%
% Find the hC-vectors that represent the given polytope(s), but reside in
%   the cone of valid polytopes with: dec.validityCone.A * hcMapped <= 0.
%
% mapToValidConeDec(hcMatrix, dec, zerotol) sets the zerotol to determine
%   the validity and to select the projections from precomputed conditions
%   in dec.confinementMappingData.projectionList().
%
% [hcMatrix, filterToZero] = mapToValidCone(..) returns a logical row 
%   vector filterToZero that contains true for empty polytopes. Note that
%   empty polytopes, empty sets {}, cannot be represented by a valid vector
%   so that they must be represented by the origin, the set {0}.
%
% mapToValidityCone(hcMatrix, dec, zerotol, gcMatrix) explicitely
%   realises growth rate mapping. The violated constraints are determined
%   based on hcMatrix for which zerotol is used such that the validity
%   projection is applied gradually from points residing a distance of
%   zerotol inside the cone (0%) to a distance of 0 (100%). For a
%   projection to be applied, the vectors in gcMatrix must also indicate
%   growing out facets. The projection is finally applied to gcMatrix, not
%   to hcMatrix.
%
% See also: obtainDec, validateDec

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 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/>.

%! Info: The underlying principle is described in the thesis.
%/

%% Comments on zertol=0
% If zerotol=0 is applied, the user obviously wants to ensure that the
%   vectors reside inside the validity cone. According to numerical issues
%   that where observed when a set of vectors where used (on of the
%   projections was not possible) and when the same single vectors where
%   used (every vector performed right), I concluded that matrix
%   multiplication does not yield the same results if performed in subsets.
%   The following line should demonstrate it instantly:
%     A= rand(5),B=rand(5,10);C=A*B;A*B(:,1:5)-C(:,1:5)
% Therefore when zerotol=0 is used the underlying projectFromData must use
%   a zerotol that is not zero. Otherwise, it will not be able to find a
%   projection.

%% Input handling
% define zerotol, if not provided
if ~exist('zerotol', 'var') || isempty(zerotol)
    zerotol = elkZerotol;
end

if ~exist('applyConstraintFilter', 'var')
    applyConstraintFilter = true(1, size(dec.confinementMappingData.A, 1));
end

% the following switches between normal validity mapping and growth rate
%   mapping. The remainder of the code is oriented on growth rate mapping
%   while it is fully compatible to normal validity mapping
if exist('hcMatrixOut', 'var')
    mapGrowth = 1;
else
    hcMatrixOut = hcMatrix;
    mapGrowth = 0;
end

% try to correct size of input
if size(hcMatrix, 1) ~= dec.nC
    hcMatrix = hcMatrix';
end
if size(hcMatrixOut, 1) ~= dec.nC
    hcMatrixOut = hcMatrixOut';
end

% handle incomplete representations (cannot be mapped)
if isfield(dec, 'isComplete') && ~dec.isComplete
    if all(max(dec.validityConePart.A * hcMatrixForConstraints, [], 1) < zerotol)
        % no vector is invalid, everything is OK, nothing to do
        return
    else
        % this is not (yet) supported
        error('elk:decomposition:inputError', ['Validity mapping for ' ...
        'incomplete decompositions is not possible']);
    end
end
% return if there is nothing to do
if isempty(hcMatrixOut)
    filterToZero = false(1, size(hcMatrixOut, 2));
    return
end

%% disappeared facets / disappearing facets
% flagVerbose = 0;
if mapGrowth
    % check the constraints for hcMatrix
    distanceMatrix = dec.confinementMappingData.A * hcMatrix;
    
    % for boundaries that shall not be considered, we set the distance from
    % this boundary to a sufficienty high value
    distanceMatrix(~applyConstraintFilter, :) = -2*zerotol;
    
    % get almost-violated boundaries
    [boundaryIndexVector hcIndexVector] = find(max(distanceMatrix + zerotol, 0));
    % compute correction magnitude
    correctionMagnitudeVector = (distanceMatrix(sub2ind(size(distanceMatrix), ...
        boundaryIndexVector, hcIndexVector)) + zerotol) ./ zerotol;   
    
    % ensure proper limits
    correctionMagnitudeVector = min(correctionMagnitudeVector, 1);
    correctionMagnitudeVector = max(correctionMagnitudeVector, 0);
    % correct dimensions if distanceMatrix is a row vector. Usually this would
    %   have been a matrix and the resulting vectors column vectors.
    if size(distanceMatrix, 1) == 1
        boundaryIndexVector = boundaryIndexVector';
        hcIndexVector = hcIndexVector';
        correctionMagnitudeVector = correctionMagnitudeVector';
    end
    % make matrix
    invalidConstraintMatrix = [boundaryIndexVector(:), hcIndexVector(:), correctionMagnitudeVector(:)];
        
    % check the constraints for hcMatrixOut
    if ~isempty(invalidConstraintMatrix)
        valueVector = sum(...
            dec.confinementMappingData.A(invalidConstraintMatrix(:, 1), :)'.*...
            hcMatrixOut(:, invalidConstraintMatrix(:,2)), 1);
        % delete rows that are now fulfilled
        invalidConstraintMatrix(valueVector <= 0, :) = [];
    end
else
    % create filterConstraintMatrix
    valueMatrix = dec.confinementMappingData.A * hcMatrix;
    [boundaryIndexVector hcIndexVector] = find(valueMatrix > zerotol);
    % build constraint matrix
    invalidConstraintMatrix = [boundaryIndexVector(:), hcIndexVector(:)];
%     valueVector = valueMatrix(sub2ind(size(valueMatrix), ...
%         invalidConstraintMatrix(:,1), invalidConstraintMatrix(:,2)...
%         ))';
    % correction magnitude is 1
    if ~isempty(invalidConstraintMatrix)
        invalidConstraintMatrix(:, 3) = 1;
    end
end

% return, if there remains nothing to do
if isempty(invalidConstraintMatrix)
    filterToZero = false(1, size(hcMatrixOut, 2));
    return
end

%% substitude boundaries indices by responsible facetgroup indices
% the projection lists are merged for each facet group and the algorithm
%   continues with correcting each facet group seperately
invalidConstraintMatrix(:, 1) = dec.confinementMappingData.responsibleVector(...
    invalidConstraintMatrix(:, 1));
    
% Result and inputs of the for loop:
%   hcMatrixOut - to be corrected
%   hcMatrix - basis for constraints to be checked
%   invalidConstraintMatrix - set of contraints to consider (they are
%              violated), also contains correction magnitude
%   valueVector - accompanies invalidConstraintMatrix and must be updated
%              after each correction

%% handle zero mappings
%! DEV: Zero mapping was treated differently in a prior version. The
%       improper version only ensured that the zero-mapping is done only
%       once in the following step. I've now merged it so that it directly
%       executes the zero mapping. Despite of that the confinement mapping
%       might be reconsidered to be FULLY executed in the embedding
%/

% initialize reminder vector
filterToZero = false(1, size(hcMatrixOut, 2));

if dec.isProper
    indexZero = dec.nC+1;
else
    indexZero = 1;
end

% affected vectors for zero-mapping
vectorIndexVecor = invalidConstraintMatrix(...
    invalidConstraintMatrix(:, 1) == indexZero, 2);
% correction
hcMatrixOut(:, vectorIndexVecor) = 0*hcMatrixOut(:, vectorIndexVecor);
% save which vectors where affected
filterToZero(vectorIndexVecor) = 1;
%! Dev: the following lines were problematic - when a lot of vectors are
%  handled, there are also a lot of constraints that might be triggered and
%  the bsxfun then combines a lot of rows with a lot of columns - creating
%  an incredibly large matrix. Not reasonable:
% % % avoid correcting the same vectors again
% % deleteFilter = any(bsxfun(@eq, invalidConstraintMatrix(:, 2), ...
% %     vectorIndexVecor(:)'), 2);
% % invalidConstraintMatrix(deleteFilter, :) = [];
% instead we cycle through the remaining vector indices
if ~isempty(vectorIndexVecor)
    vectorIndexVecor = unique(vectorIndexVecor);
    for vectorIndex = vectorIndexVecor(:)'
        invalidConstraintMatrix(...
            invalidConstraintMatrix(:, 2) == vectorIndex, :) = [];
    end
end


% return, if there remains nothing to do
if isempty(invalidConstraintMatrix);return;end

%% cycle through groups
for iGroup = 1:dec.nC        
    % these vectors are affected:
    thisVectorIndexVector = invalidConstraintMatrix(...
        invalidConstraintMatrix(:, 1) == iGroup, 2);
    % skip, if there is nothing to do
    if isempty(thisVectorIndexVector);continue;end
    
    % these vectors should use this correction magnitude:
    correctionMagnitudeVector = invalidConstraintMatrix(...
        invalidConstraintMatrix(:, 1) == iGroup, 3)';
    
    % the vectors might be chosen several times at once (multiple entries
    %   in responsibleMatrix after trafo to directions. This is corrected,
    %   here:
    [thisVectorIndexVector originalIndexVector] = unique(thisVectorIndexVector);
    %! DEV: reconsider - this following line means, we choose a arbitrary
    %    correction magnitude, here.
    correctionMagnitudeVector = correctionMagnitudeVector(originalIndexVector);
    %/
    
    % project from data (for increase of zerotol, see comments above)
    hcMatrixOut(:, thisVectorIndexVector) = ...
        bsxfun(@times, correctionMagnitudeVector, ...
          projectFromData(hcMatrix(:, thisVectorIndexVector), ...
          dec.confinementMappingData.projectionData(iGroup), ...
          zerotol + 10*eps, 2, hcMatrixOut(:, thisVectorIndexVector))) + ...
        bsxfun(@times, 1 - correctionMagnitudeVector, ...
          hcMatrixOut(:, thisVectorIndexVector));      
end

%% handle output
if nargout < 2
    clear filterToZero
end