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, validateDec0001 function [hcMatrixOut filterToZero] = mapToValidityCone(hcMatrix, dec, ... 0002 zerotol, hcMatrixOut, applyConstraintFilter) 0003 % map a hC vector to the cone of valid polytopes 0004 % 0005 % Syntax: hcMapped = mapToValidConeDec(hcMatrix, dec) 0006 % 0007 % Find the hC-vectors that represent the given polytope(s), but reside in 0008 % the cone of valid polytopes with: dec.validityCone.A * hcMapped <= 0. 0009 % 0010 % mapToValidConeDec(hcMatrix, dec, zerotol) sets the zerotol to determine 0011 % the validity and to select the projections from precomputed conditions 0012 % in dec.confinementMappingData.projectionList(). 0013 % 0014 % [hcMatrix, filterToZero] = mapToValidCone(..) returns a logical row 0015 % vector filterToZero that contains true for empty polytopes. Note that 0016 % empty polytopes, empty sets {}, cannot be represented by a valid vector 0017 % so that they must be represented by the origin, the set {0}. 0018 % 0019 % mapToValidityCone(hcMatrix, dec, zerotol, gcMatrix) explicitely 0020 % realises growth rate mapping. The violated constraints are determined 0021 % based on hcMatrix for which zerotol is used such that the validity 0022 % projection is applied gradually from points residing a distance of 0023 % zerotol inside the cone (0%) to a distance of 0 (100%). For a 0024 % projection to be applied, the vectors in gcMatrix must also indicate 0025 % growing out facets. The projection is finally applied to gcMatrix, not 0026 % to hcMatrix. 0027 % 0028 % See also: obtainDec, validateDec 0029 0030 % The elk-library: convex geometry applied to crystallization modeling. 0031 % Copyright (C) 2012 Alexander Reinhold 0032 % 0033 % This program is free software: you can redistribute it and/or modify it 0034 % under the terms of the GNU General Public License as published by the 0035 % Free Software Foundation, either version 3 of the License, or (at your 0036 % option) any later version. 0037 % 0038 % This program is distributed in the hope that it will be useful, but 0039 % WITHOUT ANY WARRANTY; without even the implied warranty of 0040 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0041 % General Public License for more details. 0042 % 0043 % You should have received a copy of the GNU General Public License along 0044 % with this program. If not, see <http://www.gnu.org/licenses/>. 0045 0046 %! Info: The underlying principle is described in the thesis. 0047 %/ 0048 0049 %% Comments on zertol=0 0050 % If zerotol=0 is applied, the user obviously wants to ensure that the 0051 % vectors reside inside the validity cone. According to numerical issues 0052 % that where observed when a set of vectors where used (on of the 0053 % projections was not possible) and when the same single vectors where 0054 % used (every vector performed right), I concluded that matrix 0055 % multiplication does not yield the same results if performed in subsets. 0056 % The following line should demonstrate it instantly: 0057 % A= rand(5),B=rand(5,10);C=A*B;A*B(:,1:5)-C(:,1:5) 0058 % Therefore when zerotol=0 is used the underlying projectFromData must use 0059 % a zerotol that is not zero. Otherwise, it will not be able to find a 0060 % projection. 0061 0062 %% Input handling 0063 % define zerotol, if not provided 0064 if ~exist('zerotol', 'var') || isempty(zerotol) 0065 zerotol = elkZerotol; 0066 end 0067 0068 if ~exist('applyConstraintFilter', 'var') 0069 applyConstraintFilter = true(1, size(dec.confinementMappingData.A, 1)); 0070 end 0071 0072 % the following switches between normal validity mapping and growth rate 0073 % mapping. The remainder of the code is oriented on growth rate mapping 0074 % while it is fully compatible to normal validity mapping 0075 if exist('hcMatrixOut', 'var') 0076 mapGrowth = 1; 0077 else 0078 hcMatrixOut = hcMatrix; 0079 mapGrowth = 0; 0080 end 0081 0082 % try to correct size of input 0083 if size(hcMatrix, 1) ~= dec.nC 0084 hcMatrix = hcMatrix'; 0085 end 0086 if size(hcMatrixOut, 1) ~= dec.nC 0087 hcMatrixOut = hcMatrixOut'; 0088 end 0089 0090 % handle incomplete representations (cannot be mapped) 0091 if isfield(dec, 'isComplete') && ~dec.isComplete 0092 if all(max(dec.validityConePart.A * hcMatrixForConstraints, [], 1) < zerotol) 0093 % no vector is invalid, everything is OK, nothing to do 0094 return 0095 else 0096 % this is not (yet) supported 0097 error('elk:decomposition:inputError', ['Validity mapping for ' ... 0098 'incomplete decompositions is not possible']); 0099 end 0100 end 0101 % return if there is nothing to do 0102 if isempty(hcMatrixOut) 0103 filterToZero = false(1, size(hcMatrixOut, 2)); 0104 return 0105 end 0106 0107 %% disappeared facets / disappearing facets 0108 % flagVerbose = 0; 0109 if mapGrowth 0110 % check the constraints for hcMatrix 0111 distanceMatrix = dec.confinementMappingData.A * hcMatrix; 0112 0113 % for boundaries that shall not be considered, we set the distance from 0114 % this boundary to a sufficienty high value 0115 distanceMatrix(~applyConstraintFilter, :) = -2*zerotol; 0116 0117 % get almost-violated boundaries 0118 [boundaryIndexVector hcIndexVector] = find(max(distanceMatrix + zerotol, 0)); 0119 % compute correction magnitude 0120 correctionMagnitudeVector = (distanceMatrix(sub2ind(size(distanceMatrix), ... 0121 boundaryIndexVector, hcIndexVector)) + zerotol) ./ zerotol; 0122 0123 % ensure proper limits 0124 correctionMagnitudeVector = min(correctionMagnitudeVector, 1); 0125 correctionMagnitudeVector = max(correctionMagnitudeVector, 0); 0126 % correct dimensions if distanceMatrix is a row vector. Usually this would 0127 % have been a matrix and the resulting vectors column vectors. 0128 if size(distanceMatrix, 1) == 1 0129 boundaryIndexVector = boundaryIndexVector'; 0130 hcIndexVector = hcIndexVector'; 0131 correctionMagnitudeVector = correctionMagnitudeVector'; 0132 end 0133 % make matrix 0134 invalidConstraintMatrix = [boundaryIndexVector(:), hcIndexVector(:), correctionMagnitudeVector(:)]; 0135 0136 % check the constraints for hcMatrixOut 0137 if ~isempty(invalidConstraintMatrix) 0138 valueVector = sum(... 0139 dec.confinementMappingData.A(invalidConstraintMatrix(:, 1), :)'.*... 0140 hcMatrixOut(:, invalidConstraintMatrix(:,2)), 1); 0141 % delete rows that are now fulfilled 0142 invalidConstraintMatrix(valueVector <= 0, :) = []; 0143 end 0144 else 0145 % create filterConstraintMatrix 0146 valueMatrix = dec.confinementMappingData.A * hcMatrix; 0147 [boundaryIndexVector hcIndexVector] = find(valueMatrix > zerotol); 0148 % build constraint matrix 0149 invalidConstraintMatrix = [boundaryIndexVector(:), hcIndexVector(:)]; 0150 % valueVector = valueMatrix(sub2ind(size(valueMatrix), ... 0151 % invalidConstraintMatrix(:,1), invalidConstraintMatrix(:,2)... 0152 % ))'; 0153 % correction magnitude is 1 0154 if ~isempty(invalidConstraintMatrix) 0155 invalidConstraintMatrix(:, 3) = 1; 0156 end 0157 end 0158 0159 % return, if there remains nothing to do 0160 if isempty(invalidConstraintMatrix) 0161 filterToZero = false(1, size(hcMatrixOut, 2)); 0162 return 0163 end 0164 0165 %% substitude boundaries indices by responsible facetgroup indices 0166 % the projection lists are merged for each facet group and the algorithm 0167 % continues with correcting each facet group seperately 0168 invalidConstraintMatrix(:, 1) = dec.confinementMappingData.responsibleVector(... 0169 invalidConstraintMatrix(:, 1)); 0170 0171 % Result and inputs of the for loop: 0172 % hcMatrixOut - to be corrected 0173 % hcMatrix - basis for constraints to be checked 0174 % invalidConstraintMatrix - set of contraints to consider (they are 0175 % violated), also contains correction magnitude 0176 % valueVector - accompanies invalidConstraintMatrix and must be updated 0177 % after each correction 0178 0179 %% handle zero mappings 0180 %! DEV: Zero mapping was treated differently in a prior version. The 0181 % improper version only ensured that the zero-mapping is done only 0182 % once in the following step. I've now merged it so that it directly 0183 % executes the zero mapping. Despite of that the confinement mapping 0184 % might be reconsidered to be FULLY executed in the embedding 0185 %/ 0186 0187 % initialize reminder vector 0188 filterToZero = false(1, size(hcMatrixOut, 2)); 0189 0190 if dec.isProper 0191 indexZero = dec.nC+1; 0192 else 0193 indexZero = 1; 0194 end 0195 0196 % affected vectors for zero-mapping 0197 vectorIndexVecor = invalidConstraintMatrix(... 0198 invalidConstraintMatrix(:, 1) == indexZero, 2); 0199 % correction 0200 hcMatrixOut(:, vectorIndexVecor) = 0*hcMatrixOut(:, vectorIndexVecor); 0201 % save which vectors where affected 0202 filterToZero(vectorIndexVecor) = 1; 0203 %! Dev: the following lines were problematic - when a lot of vectors are 0204 % handled, there are also a lot of constraints that might be triggered and 0205 % the bsxfun then combines a lot of rows with a lot of columns - creating 0206 % an incredibly large matrix. Not reasonable: 0207 % % % avoid correcting the same vectors again 0208 % % deleteFilter = any(bsxfun(@eq, invalidConstraintMatrix(:, 2), ... 0209 % % vectorIndexVecor(:)'), 2); 0210 % % invalidConstraintMatrix(deleteFilter, :) = []; 0211 % instead we cycle through the remaining vector indices 0212 if ~isempty(vectorIndexVecor) 0213 vectorIndexVecor = unique(vectorIndexVecor); 0214 for vectorIndex = vectorIndexVecor(:)' 0215 invalidConstraintMatrix(... 0216 invalidConstraintMatrix(:, 2) == vectorIndex, :) = []; 0217 end 0218 end 0219 0220 0221 % return, if there remains nothing to do 0222 if isempty(invalidConstraintMatrix);return;end 0223 0224 %% cycle through groups 0225 for iGroup = 1:dec.nC 0226 % these vectors are affected: 0227 thisVectorIndexVector = invalidConstraintMatrix(... 0228 invalidConstraintMatrix(:, 1) == iGroup, 2); 0229 % skip, if there is nothing to do 0230 if isempty(thisVectorIndexVector);continue;end 0231 0232 % these vectors should use this correction magnitude: 0233 correctionMagnitudeVector = invalidConstraintMatrix(... 0234 invalidConstraintMatrix(:, 1) == iGroup, 3)'; 0235 0236 % the vectors might be chosen several times at once (multiple entries 0237 % in responsibleMatrix after trafo to directions. This is corrected, 0238 % here: 0239 [thisVectorIndexVector originalIndexVector] = unique(thisVectorIndexVector); 0240 %! DEV: reconsider - this following line means, we choose a arbitrary 0241 % correction magnitude, here. 0242 correctionMagnitudeVector = correctionMagnitudeVector(originalIndexVector); 0243 %/ 0244 0245 % project from data (for increase of zerotol, see comments above) 0246 hcMatrixOut(:, thisVectorIndexVector) = ... 0247 bsxfun(@times, correctionMagnitudeVector, ... 0248 projectFromData(hcMatrix(:, thisVectorIndexVector), ... 0249 dec.confinementMappingData.projectionData(iGroup), ... 0250 zerotol + 10*eps, 2, hcMatrixOut(:, thisVectorIndexVector))) + ... 0251 bsxfun(@times, 1 - correctionMagnitudeVector, ... 0252 hcMatrixOut(:, thisVectorIndexVector)); 0253 end 0254 0255 %% handle output 0256 if nargout < 2 0257 clear filterToZero 0258 end