Description of constructProperRepresentation

0001 function dec = constructProperRepresentation(dec, validityDataOriginal, ...
0002     zerotol, verbose, allowImproper)
0003 % check (and construct) proper representation
0004 %
0005 % error/warnings are thrown according to allowImproper. If the input is
0006 %   improper, a temporary field dec.improperData is added.
0007 %
0008 % THIS IS NO USER FUNCTION
0009 
0010 % The elk-library: convex geometry applied to crystallization modeling.
0011 %   Copyright (C) 2013 Alexander Reinhold
0012 %
0013 % This program is free software: you can redistribute it and/or modify it
0014 %   under the terms of the GNU General Public License as published by the
0015 %   Free Software Foundation, either version 3 of the License, or (at your
0016 %   option) any later version.
0017 %
0018 % This program is distributed in the hope that it will be useful, but
0019 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0020 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0021 %   General Public License for more details.
0022 %
0023 % You should have received a copy of the GNU General Public License along
0024 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0025 
0026 %% reject ill-conditioned input
0027 if rank(dec.mappingReducedToFull) < dec.nC
0028     error('elk:decomposition:wrongInput', ['The input representation is ' ...
0029         'ill-conditioned, its dimension nC can be reduced without loss ' ...
0030         'of shape information']);
0031 end
0032 
0033 %% get validity cone data in hC-space
0034 % transfer cone to hC-space
0035 validityCone.A = validityDataOriginal.A * dec.mappingReducedToFull;
0036 % normalize and trim (no removal of zero-rows)
0037 validityCone.A = trimConstraintMatrix(validityCone.A, zerotol, 0);
0038 
0039 %% decide properness
0040 % we assume properness until something else is shown
0041 isProper = 1;
0042 % this facet filter, we try to sucessively reduce by handling facet by
0043 %   facet
0044 facetFilter = true(1, dec.nH);
0045 % in this vector, we record the group index to which each facet belongs
0046 facetGroupIndexVector = zeros(1, dec.nH);
0047 % this matrix records the mapping from an improper hC vecctor to a proper
0048 %   hE vector
0049 mappingReducedToFull = zeros(dec.nH, 0);
0050 mappingImproperToProper = zeros(0, dec.nC);
0051 % we cycle through all facets and handle groups of coherent facets
0052 while any(facetFilter)
0053     % init
0054     remainingFacetVector = find(facetFilter);
0055     thisFacet = remainingFacetVector(1);
0056     
0057     % find coherent facets
0058     referenceConstraintMatrix = validityCone.A(...
0059         validityDataOriginal.responsibleFacetVector == thisFacet, :);
0060     referenceConstraintMatrix = sortrows(referenceConstraintMatrix);
0061     coherentFacetFilter = facetFilter;
0062     remainingIndex = 2;
0063     while remainingIndex <= length(remainingFacetVector)
0064         constraintMatrix = validityCone.A(...
0065             validityDataOriginal.responsibleFacetVector == ...
0066             remainingFacetVector(remainingIndex), :);
0067         constraintMatrix = sortrows(constraintMatrix);
0068         if size(constraintMatrix, 1) == size(referenceConstraintMatrix, 1) && ...
0069                 all(referenceConstraintMatrix(:) == constraintMatrix(:))
0070             % thats fine, coherent facet filter can remain as is
0071         else
0072             % this facet is not coherent, change coherentFacetIndex
0073             coherentFacetFilter(remainingFacetVector(remainingIndex)) = false;
0074         end
0075         
0076         remainingIndex = remainingIndex + 1;
0077     end
0078     
0079     % go through properness conditions
0080     matrixMg = dec.mappingReducedToFull;
0081     matrixMg(coherentFacetFilter, :) = 0;
0082     matrixA = dec.A;
0083     matrixA(coherentFacetFilter, :) = 0;
0084     if size(referenceConstraintMatrix, 1) == 0
0085         % if no constraint is selected, specifically for this facet, then
0086         %   we currently deal with the non-empty polytope cone, which means
0087         %   that the facet can only disappear together with all others.
0088 %         disp('empty')
0089     elseif rank(matrixMg) < dec.nC
0090         % thats the usual properness where the correction in done according
0091         %   to the nullspace of matrixMg
0092 %         disp('normal')
0093     elseif rank(dec.constraintMatrix * matrixA) < dec.nDim
0094         % in this case, validity mapping would use the translation of the
0095         %   polytope.
0096 %         disp('translation')
0097     else
0098         % if none of the above applies, this set of coherent facets cannot
0099         %   be corrected, if they disappear and we have an improper
0100         %   representation
0101 %         disp('improper')
0102         isProper = 0;
0103     end
0104     
0105     % remove checked facets
0106     facetFilter(coherentFacetFilter) = false;
0107     
0108     % assign facet groups which also considers unique rows of the group
0109     %   mapping matrix to construct the proper embedding
0110     while any(coherentFacetFilter)
0111         thisIndex = find(coherentFacetFilter, 1, 'first');
0112         subFacetFilter = coherentFacetFilter & all(bsxfun(@minus, ...
0113             dec.mappingReducedToFull, ...
0114             dec.mappingReducedToFull(thisIndex, :)) == 0, 2)';
0115         facetGroupIndexVector(subFacetFilter) = max(facetGroupIndexVector) + 1;
0116         % the facets in this subset will have equal distances given by the
0117         %   hE component that is added in this line:
0118         mappingReducedToFull(:, end+1) = double(subFacetFilter)';
0119         % mapping from hC-space to hE-space means to apply the row of the
0120         %   group mapping matrix to the hC-vector
0121         mappingImproperToProper(end+1, :) = dec.mappingReducedToFull(thisIndex, :);
0122         % the subset of facets must not be considered further
0123         coherentFacetFilter(subFacetFilter) = false;
0124     end
0125 end
0126 
0127 % Everything is done, if representation is proper
0128 if isProper
0129     dec.isProper = 1;
0130     % nothing left to do
0131     return
0132 else
0133     % error/warning handling
0134     if allowImproper == 0
0135         error('elk:decomposition:improper', ['The input representation is ' ...
0136             'improper. To allow constructing a proper representation, ' ...
0137             'change parameter allowImproper to 1']);
0138     elseif allowImproper == 1
0139         warning('elk:decomposition:improper', ['The input representation ' ...
0140             'is improper. The mapping/constraints are extended to a proper ' ...
0141             'representation. To handle improper representations, look for ' ...
0142             'approximation of decomposition data. To supress this warning, '...
0143             'use allowImproper with (-1)']);
0144     end   
0145     
0146     
0147     % this is temporary data for later projection..
0148     dec.isProper = 0;
0149     % ..copy current dec information
0150     dec.improperData = dec;
0151     % ..add mapping matrices
0152     dec.improperData.mappingImproperToProper = mappingImproperToProper;
0153     dec.improperData.mappingProperToImproper = pinv(mappingImproperToProper);
0154     
0155     
0156     dec.isProper = 1;
0157     
0158     % rewritte dec mapping
0159     dec.constraintMatrix = [];
0160     dec.mappingReducedToFull = mappingReducedToFull;
0161     % double-check position fix (this should not be an issue), but this
0162     %   also regenerates the constraint matrix and reverse mapping.
0163     try
0164         dec = computeMappingAndConstraint(dec, zerotol, 0);
0165     catch me
0166         error('elk:decomposition:numericalError', ['Constructed proper ' ...
0167             'decomposition does not fulfill the fixed-position constraint. ' ...
0168             'This should not happen.']);
0169     end
0170     dec.nC = size(mappingImproperToProper, 1);
0171 end
0172 
0173 return
0174 
0175 % -----------------------------------------
0176 % This needs to be realized by the new code:
0177 %   1) the mapping between proper embedding / improper rep
0178 %   2) the indices to representative facets that constitute the new group
0179 %      mapping matrix.
0180 
0181 
0182     
0183     
0184 % --------------------------------------------------------------------
0185 
0186 
0187 %% get unique direction matrix
0188 [mappingImproperToProper facetIndexVector directionIndexVector] = ...
0189     unique(dec.mappingReducedToFull, 'rows');
0190 nDir = size(mappingImproperToProper, 1);
0191 directionMatrix = eye(size(mappingImproperToProper, 1));
0192 
0193 
0194 %% pre-consideration
0195 % the decomposition cannot be proper when more directions in the mapping
0196 %   matrix exist than there are degrees of freedom. BUT, zero rows do not
0197 %   count as a unique row.
0198 if size(eliminateRows(mappingImproperToProper, ...
0199                       zeros(1, size(mappingImproperToProper, 2)), 0), ...
0200         1) > dec.nC
0201     isProper = 0;
0202 else
0203     isProper = 1;
0204 end
0205 %% generate data for proper decomposition
0206 % cycle over directions; then collect all facets that are related to the
0207 %   first entry in thisFacetIndexVector. Repeat, unless
0208 %   nextFacdetIndexVector is empty.
0209 for iDir = 1:nDir
0210     thisFacetIndexVector = find(directionIndexVector == iDir);
0211     flagKeepChecking = 1;
0212     while flagKeepChecking
0213         % generate reference cone
0214         thisReferenceCone = validityCone.A(...
0215             validityDataOriginal.responsibleFacetVector == ...
0216             thisFacetIndexVector(1), :);
0217         thisReferenceCone = unique(thisReferenceCone, 'rows');
0218         % initialize
0219         nextFacetIndexVector = [];
0220         % check remaining cones
0221         for i = 2:length(thisFacetIndexVector)            
0222             thisCone = validityCone.A(...
0223             validityDataOriginal.responsibleFacetVector == ...
0224             thisFacetIndexVector(i), :);
0225             thisCone = unique(thisCone, 'rows');
0226             if size(thisCone, 1) == size(thisReferenceCone, 1) && ...
0227                     (any(thisCone(:) == thisReferenceCone(:)) || ...
0228                      isempty(thisCone))
0229                 % facets are related, the information in
0230                 %   directionIndexVector can remain as is. Nothing to do.
0231             else
0232                 % facets are not related >> improper
0233                 isProper = 0;
0234                 nextFacetIndexVector(end+1) = thisFacetIndexVector(i);
0235                 %! Debug:
0236 %                 [iDir, i, thisFacetIndexVector(1), thisFacetIndexVector(i)]
0237                 
0238             end
0239         end
0240         if isempty(nextFacetIndexVector)
0241             % everything proper, now, stop checking
0242             flagKeepChecking = 0;
0243         else
0244             % some facets were not related, we need to restart with the new
0245             % facet set and add a "new" direction:
0246             directionMatrix(end+1, end+1) = 1;
0247             mappingImproperToProper(end+1, :) = ...
0248                 mappingImproperToProper(iDir , :);
0249             directionIndexVector(nextFacetIndexVector) = ...
0250                 size(directionMatrix, 1);
0251             % report details
0252             if verbose
0253                 disp(['Info: facet ' num2str(thisFacetIndexVector(1)) ...
0254                     ' and facet(s) ' num2str(nextFacetIndexVector) ' are '...
0255                     'not related, but share the same direction in '...
0256                     'hC-space.']);
0257             end
0258             thisFacetIndexVector = nextFacetIndexVector;
0259         end
0260     end
0261 end
0262 % Result: We know now, whether the input representation is proper, or not.
0263 %   We also stored in directionIndexVector, which facets are assigned to
0264 %   which direction
0265 
0266 %% handle results
0267 % resulting decomposition is always proper
0268 dec.isProper = 1;
0269 if isProper
0270     % nothing to do
0271 else
0272     % error/warning handling
0273     if allowImproper == 0
0274         error('elk:decomposition:improper', ['The input representation is ' ...
0275             'improper. To allow constructing a proper representation, ' ...
0276             'change parameter allowImproper to 1']);
0277     elseif allowImproper == 1
0278         warning('elk:decomposition:improper', ['The input representation ' ...
0279             'is improper. The mapping/constraints are extended to a proper ' ...
0280             'representation. To handle improper representations, look for ' ...
0281             'approximation of decomposition data. To supress this warning, '...
0282             'use allowImproper with (-1)']);
0283     end
0284     % this is temporary data for later projection..
0285     % ..copy current dec information
0286     dec.improperData = dec;
0287     % ..add mapping matrices
0288     dec.improperData.mappingImproperToProper = mappingImproperToProper;
0289     dec.improperData.mappingProperToImproper = ...
0290         pinv(mappingImproperToProper);    
0291     
0292     % rewritte dec mapping
0293     dec.constraintMatrix = [];
0294     dec.mappingReducedToFull = directionMatrix(directionIndexVector, :);
0295     % double-check position fix (this should not be an issue), but this
0296     %   also regenerates the constraint matrix and reverse mapping.
0297     try
0298         dec = computeMappingAndConstraint(dec, zerotol, 0);
0299     catch me
0300         error('elk:decomposition:numericalError', ['Constructed proper ' ...
0301             'decomposition does not fulfill the fixed-position constraint. ' ...
0302             'This should not happen.']);
0303     end
0304     dec.nC = max(directionIndexVector);
0305 end
0306 
0307 end
constructProperRepresentation

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
