function dec = constructProperRepresentation(dec, validityDataOriginal, ...
    zerotol, verbose, allowImproper)
% check (and construct) proper representation
%
% error/warnings are thrown according to allowImproper. If the input is
%   improper, a temporary field dec.improperData is added.
%
% 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/>.

%% reject ill-conditioned input
if rank(dec.mappingReducedToFull) < dec.nC
    error('elk:decomposition:wrongInput', ['The input representation is ' ...
        'ill-conditioned, its dimension nC can be reduced without loss ' ...
        'of shape information']);
end

%% get validity cone data in hC-space
% transfer cone to hC-space
validityCone.A = validityDataOriginal.A * dec.mappingReducedToFull;
% normalize and trim (no removal of zero-rows)
validityCone.A = trimConstraintMatrix(validityCone.A, zerotol, 0);

%% decide properness
% we assume properness until something else is shown
isProper = 1;
% this facet filter, we try to sucessively reduce by handling facet by
%   facet
facetFilter = true(1, dec.nH);
% in this vector, we record the group index to which each facet belongs
facetGroupIndexVector = zeros(1, dec.nH);
% this matrix records the mapping from an improper hC vecctor to a proper
%   hE vector
mappingReducedToFull = zeros(dec.nH, 0);
mappingImproperToProper = zeros(0, dec.nC);
% we cycle through all facets and handle groups of coherent facets
while any(facetFilter)
    % init
    remainingFacetVector = find(facetFilter);
    thisFacet = remainingFacetVector(1);
    
    % find coherent facets
    referenceConstraintMatrix = validityCone.A(...
        validityDataOriginal.responsibleFacetVector == thisFacet, :);
    referenceConstraintMatrix = sortrows(referenceConstraintMatrix);
    coherentFacetFilter = facetFilter;
	remainingIndex = 2;
    while remainingIndex <= length(remainingFacetVector)
        constraintMatrix = validityCone.A(...
            validityDataOriginal.responsibleFacetVector == ...
            remainingFacetVector(remainingIndex), :);
        constraintMatrix = sortrows(constraintMatrix);
        if size(constraintMatrix, 1) == size(referenceConstraintMatrix, 1) && ...
                all(referenceConstraintMatrix(:) == constraintMatrix(:))
            % thats fine, coherent facet filter can remain as is
        else
            % this facet is not coherent, change coherentFacetIndex
            coherentFacetFilter(remainingFacetVector(remainingIndex)) = false;
        end
        
        remainingIndex = remainingIndex + 1;
    end
    
    % go through properness conditions
    matrixMg = dec.mappingReducedToFull;
    matrixMg(coherentFacetFilter, :) = 0;
    matrixA = dec.A;
    matrixA(coherentFacetFilter, :) = 0;
    if size(referenceConstraintMatrix, 1) == 0
        % if no constraint is selected, specifically for this facet, then
        %   we currently deal with the non-empty polytope cone, which means
        %   that the facet can only disappear together with all others.
%         disp('empty')
    elseif rank(matrixMg) < dec.nC
        % thats the usual properness where the correction in done according
        %   to the nullspace of matrixMg
%         disp('normal')
    elseif rank(dec.constraintMatrix * matrixA) < dec.nDim
        % in this case, validity mapping would use the translation of the
        %   polytope.
%         disp('translation')
    else
        % if none of the above applies, this set of coherent facets cannot
        %   be corrected, if they disappear and we have an improper
        %   representation
%         disp('improper')
        isProper = 0;
    end
    
    % remove checked facets
    facetFilter(coherentFacetFilter) = false;
    
    % assign facet groups which also considers unique rows of the group
    %   mapping matrix to construct the proper embedding
    while any(coherentFacetFilter)
        thisIndex = find(coherentFacetFilter, 1, 'first');
        subFacetFilter = coherentFacetFilter & all(bsxfun(@minus, ...
            dec.mappingReducedToFull, ...
            dec.mappingReducedToFull(thisIndex, :)) == 0, 2)';
        facetGroupIndexVector(subFacetFilter) = max(facetGroupIndexVector) + 1;
        % the facets in this subset will have equal distances given by the
        %   hE component that is added in this line:
        mappingReducedToFull(:, end+1) = double(subFacetFilter)';
        % mapping from hC-space to hE-space means to apply the row of the
        %   group mapping matrix to the hC-vector
        mappingImproperToProper(end+1, :) = dec.mappingReducedToFull(thisIndex, :);
        % the subset of facets must not be considered further
        coherentFacetFilter(subFacetFilter) = false;
    end
end

% Everything is done, if representation is proper
if isProper
    dec.isProper = 1;
    % nothing left to do
    return
else
    % error/warning handling
    if allowImproper == 0
        error('elk:decomposition:improper', ['The input representation is ' ...
            'improper. To allow constructing a proper representation, ' ...
            'change parameter allowImproper to 1']);
    elseif allowImproper == 1
        warning('elk:decomposition:improper', ['The input representation ' ...
            'is improper. The mapping/constraints are extended to a proper ' ...
            'representation. To handle improper representations, look for ' ...
            'approximation of decomposition data. To supress this warning, '...
            'use allowImproper with (-1)']);
    end   
    
    
    % this is temporary data for later projection..
    dec.isProper = 0;
    % ..copy current dec information
    dec.improperData = dec;
    % ..add mapping matrices
    dec.improperData.mappingImproperToProper = mappingImproperToProper;
    dec.improperData.mappingProperToImproper = pinv(mappingImproperToProper);
    
    
    dec.isProper = 1;
    
    % rewritte dec mapping
    dec.constraintMatrix = [];
    dec.mappingReducedToFull = mappingReducedToFull;
    % double-check position fix (this should not be an issue), but this 
    %   also regenerates the constraint matrix and reverse mapping.
    try
        dec = computeMappingAndConstraint(dec, zerotol, 0);
    catch me
        error('elk:decomposition:numericalError', ['Constructed proper ' ...
            'decomposition does not fulfill the fixed-position constraint. ' ...
            'This should not happen.']);
    end
    dec.nC = size(mappingImproperToProper, 1);
end

return

% -----------------------------------------
% This needs to be realized by the new code:
%   1) the mapping between proper embedding / improper rep
%   2) the indices to representative facets that constitute the new group
%      mapping matrix.


    
    
% --------------------------------------------------------------------


%% get unique direction matrix
[mappingImproperToProper facetIndexVector directionIndexVector] = ...
    unique(dec.mappingReducedToFull, 'rows');
nDir = size(mappingImproperToProper, 1);
directionMatrix = eye(size(mappingImproperToProper, 1));


%% pre-consideration
% the decomposition cannot be proper when more directions in the mapping
%   matrix exist than there are degrees of freedom. BUT, zero rows do not
%   count as a unique row.
if size(eliminateRows(mappingImproperToProper, ...
                      zeros(1, size(mappingImproperToProper, 2)), 0), ...
        1) > dec.nC
    isProper = 0;
else
    isProper = 1;
end
%% generate data for proper decomposition
% cycle over directions; then collect all facets that are related to the
%   first entry in thisFacetIndexVector. Repeat, unless
%   nextFacdetIndexVector is empty.
for iDir = 1:nDir
    thisFacetIndexVector = find(directionIndexVector == iDir);
    flagKeepChecking = 1;
    while flagKeepChecking
        % generate reference cone
        thisReferenceCone = validityCone.A(...
            validityDataOriginal.responsibleFacetVector == ...
            thisFacetIndexVector(1), :);
        thisReferenceCone = unique(thisReferenceCone, 'rows');
        % initialize
        nextFacetIndexVector = [];
        % check remaining cones
        for i = 2:length(thisFacetIndexVector)            
            thisCone = validityCone.A(...
            validityDataOriginal.responsibleFacetVector == ...
            thisFacetIndexVector(i), :);
            thisCone = unique(thisCone, 'rows');
            if size(thisCone, 1) == size(thisReferenceCone, 1) && ...
                    (any(thisCone(:) == thisReferenceCone(:)) || ...
                     isempty(thisCone))
                % facets are related, the information in
                %   directionIndexVector can remain as is. Nothing to do.
            else
                % facets are not related >> improper
                isProper = 0;
                nextFacetIndexVector(end+1) = thisFacetIndexVector(i);
                %! Debug:
%                 [iDir, i, thisFacetIndexVector(1), thisFacetIndexVector(i)]
                
            end
        end
        if isempty(nextFacetIndexVector)
            % everything proper, now, stop checking
            flagKeepChecking = 0;
        else
            % some facets were not related, we need to restart with the new
            % facet set and add a "new" direction:
            directionMatrix(end+1, end+1) = 1;
            mappingImproperToProper(end+1, :) = ...
                mappingImproperToProper(iDir , :);
            directionIndexVector(nextFacetIndexVector) = ...
                size(directionMatrix, 1);
            % report details
            if verbose
                disp(['Info: facet ' num2str(thisFacetIndexVector(1)) ...
                    ' and facet(s) ' num2str(nextFacetIndexVector) ' are '...
                    'not related, but share the same direction in '...
                    'hC-space.']);
            end
            thisFacetIndexVector = nextFacetIndexVector;
        end
    end
end
% Result: We know now, whether the input representation is proper, or not.
%   We also stored in directionIndexVector, which facets are assigned to
%   which direction

%% handle results
% resulting decomposition is always proper
dec.isProper = 1;
if isProper
    % nothing to do
else
    % error/warning handling
    if allowImproper == 0
        error('elk:decomposition:improper', ['The input representation is ' ...
            'improper. To allow constructing a proper representation, ' ...
            'change parameter allowImproper to 1']);
    elseif allowImproper == 1
        warning('elk:decomposition:improper', ['The input representation ' ...
            'is improper. The mapping/constraints are extended to a proper ' ...
            'representation. To handle improper representations, look for ' ...
            'approximation of decomposition data. To supress this warning, '...
            'use allowImproper with (-1)']);
    end
    % this is temporary data for later projection..
    % ..copy current dec information
    dec.improperData = dec;
    % ..add mapping matrices
    dec.improperData.mappingImproperToProper = mappingImproperToProper;
    dec.improperData.mappingProperToImproper = ...
        pinv(mappingImproperToProper);    
    
    % rewritte dec mapping
    dec.constraintMatrix = [];
    dec.mappingReducedToFull = directionMatrix(directionIndexVector, :);
    % double-check position fix (this should not be an issue), but this 
    %   also regenerates the constraint matrix and reverse mapping.
    try
        dec = computeMappingAndConstraint(dec, zerotol, 0);
    catch me
        error('elk:decomposition:numericalError', ['Constructed proper ' ...
            'decomposition does not fulfill the fixed-position constraint. ' ...
            'This should not happen.']);
    end
    dec.nC = max(directionIndexVector);
end

end