function dec = addValidityMappingData(dec, zerotol)
% add data to identify vanished facets and map to validity cone
%
% THIS FUNCTION IS NO USER FUNCTION.

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

%% data structure
% Draft of the gathered information in boundaryData:
%     .A       - boundary facet normals in nC-dimensional space
%     .nBoundary - total number of rows in A
%     .responsibleVector - a vector that indicates the facet (group) that
%                disappears given corresponding to the constraint in A.
%     .directionMatrix - the collection of directions (column vectors) in 
%                hC-space which is used to correct a specific boundary 
%                (column). It contains nan, when mapping is not possible
%                without changing a valid facet and inf when mapping
%                to the origin would be appropriate. In the latter case,
%                the boundary indicates the disappearance of all facets at
%                once.

%% prepare setup
% copy constraints
dec.confinementMappingData.A = dec.facetValidityData.A;
% predefine assigned directions (per constraint)
dec.confinementMappingData.responsibleVector = ...
    nan(1, size(dec.confinementMappingData.A, 1));
% predefine possible projection sets (per facet)
dec.confinementMappingData.projectionData = struct([]);
% predefine the directions
dec.confinementMappingData.directionMatrix = nan([dec.nC dec.nC+1]);
% local copy
responsibleMatrix = dec.facetValidityData.responsibleMatrix;

%% generate zero-projection data
dec.confinementMappingData.responsibleVector(...
    responsibleMatrix(responsibleMatrix(:, 2) == (dec.nH + 1), 1)) = dec.nC+1;
dec.confinementMappingData.projectionData(dec.nC+1).A = zeros(0, dec.nC);
dec.confinementMappingData.projectionData(dec.nC+1).responsibleMatrix = ...
    zeros(0, 3);
dec.confinementMappingData.projectionData(dec.nC+1).projectionList = ...
    {zeros(dec.nC)};
% delete all entries for these constraints, they appear for all facets
for iConstraint = find(dec.confinementMappingData.responsibleVector == (dec.nC+1))
    responsibleMatrix(responsibleMatrix(:, 1) == iConstraint, :) = [];
end

%% generate facet group projections
% construct unique non-zero rows
[directionMatrix facetIndexVector, directionIndexVector] = ...
    unique(dec.mappingReducedToFull, 'rows');
% Zero rows in the direction matrix occure when the position was fixed
%   by the decomposition algorithm. These direction must be considered
%   with some precaution. They must be treated separately. Therefore, we 
%   identify the original zero rows and multiply them appropriately.
zeroRowFilter = all(dec.mappingReducedToFull < zerotol, 2);
zeroRowIndexFilter = find(zeroRowFilter);
% if there is only one zero row then the above is already correct. If there
%   are more zero rows, we must append some zero-row directions.
if sum(zeroRowFilter) > 1
    % find the current zero-row index
    thisZeroIndex = find(all(directionMatrix < zerotol, 2));
    % add zero rows
    directionMatrix = [directionMatrix zeros(sum(zeroRowFilter)-1, dec.nC)];
    % substitute direction indices
    zeroDirIndexVector = find(directionIndexVector == thisZeroIndex);
    directionIndexVector(zeroDirIndexVector) = zeroRowIndexFilter;
end

nDir = size(directionMatrix, 1);

for iDir = 1:nDir

    % ---
    % identify facet(s) for direction
    thisFacetIndex = find(directionIndexVector == iDir, 1, 'first');
    isZeroDir = any(zeroRowIndexFilter == thisFacetIndex);
    
    % ---
    % identify mapping direction (null space)
    thisDirectionFilter = true(1, nDir);
    thisDirectionFilter(iDir) = false;
    % Things are different when iDir indicates a zero row, therefore the
    % separate treatment:
    if ~isZeroDir
        % The projection direction equals the null space of the oblique
        %   projection that is to be constructed. While this direction is
        %   not important for the handling of the proper representation, it
        %   becomes valuable as soon as an improper representation is
        %   constructed.
        dec.confinementMappingData.directionMatrix(:, iDir) = null(...
            directionMatrix(thisDirectionFilter, :));
        
        % For the construction of the obligue projection, we require the
        %   orthogonal complement of the null space (see wikipedia).
        %! DEV: this line was used beforehand.. ..and can't be right in
        %       general. It is only a proper result if the selected
        %       direction is already orthogonal to all other directions.
%         thisOrthNull = orth(directionMatrix(thisDirectionFilter, :)');
        % we use this new line:
        thisOrthNull = null(dec.confinementMappingData.directionMatrix(:, iDir)');
    else
        % Things are different for a zero row that was introduced to fix the
        %   position. The disappearance of this facet also implies that the
        %   position is not fixed anymore. The nullspace is, hence, computed
        %   from this degree of freedom. Assuming this facet has 
        %   disappeared, we temporarly ignore the facet
        tmpA = dec.A; tmpA(thisFacetIndex, :) = [];
        tmpMc = dec.constraintMatrix; tmpMc(:, thisFacetIndex) = [];
        tmpMapToH = dec.mappingReducedToFull; tmpMapToH(thisFacetIndex, :) = [];
        tmpMapFromH = pinv(tmpMapToH);
        % obtain free direction via singular value decomposition
        [U S V] = svd(tmpMc * tmpA);
        % The direction in which we can move freely is:
        %     x0 = V(:, end)
        %   which is applied according to:
        %     h' = h + A*x0.
        %   Scince h fulfills tmpMc*h=0 and h' is computed so that it
        %   fulfills tmpMc*h', its difference also fulfills this
        %   constraint so that hC is:
        thisDirHc = tmpMapFromH*tmpA*V(:, end);
        % This is the null space and we need the orthogonal complement.
        thisOrthNull = null(thisDirHc');
    end
    
    % ---
    % determine range and construct projection(s)
    
    % identify assigned constraints (range)
    dec.confinementMappingData.responsibleVector(responsibleMatrix(...
          responsibleMatrix(:, 2) == thisFacetIndex, 1)) = iDir;
    thisConstraintFilter = dec.confinementMappingData.responsibleVector == iDir;
    
    % fill projection data (one for each constraint)
    thisProjectionData = appendProjectionData([], dec.nC);
    for iConstraint = find(thisConstraintFilter)
        % the facet is spanned by the following base vectors that are 
        %   constructed from the orthogonal null space, spanned by the
        %   facet normal. Hence, we take the null space of the null space
        %   for the range of the projection onto the facet.
        thisRange = null(dec.confinementMappingData.A(iConstraint, :));
        
        % the obligue projection is now:
        thisProjection = thisRange * inv(...
            thisOrthNull' * thisRange) * thisOrthNull';
        % it must map to the group validity cone:
        %   Av * P * hC <= 0
        thisConstraint = ...
            dec.confinementMappingData.A(thisConstraintFilter, :) * ...
            thisProjection;
        % .. and - the projection must actually be valid. This is redundant
        %   to the constraints that are checked in confinementMappingData.A
        %   but this check is required for completeness when dealing with
        %   projections of proper decompositions. There, the whole validity
        %   of this projection is required.
        thisConstraint = [thisConstraint; ...
                          -1*dec.confinementMappingData.A(iConstraint, :)];
        % assign results
        thisProjectionData = appendProjectionData(thisProjectionData, ...
            thisConstraint, thisProjection, zerotol);
    end
    dec.confinementMappingData.projectionData(iDir) = thisProjectionData;  
end

% add something that is used only for improper representations
dec.confinementMappingData.ignoreFadeOut = [];
