validate decomposition of valid cone
Syntax: dec = validateDecomposition(dec)
This function checks whether dec is a proper decomposition of the set of
facet normals in dec.A. The structuring elements are double-checked for
(constrained) indecomposability and the unified partitions are checked
against the cone of valid polytopes. The simplex partitions are not
verified as they are only used internally.
The decomposition structure is extended with the fields:
dec.validation.tolerance - chosen zerotol
dec.validation.result - scalar logical 1 (passed) or 0 (failed)
dec.timing.validateDec - total time that was taken by validation
validateDec(.., 'zerotol', zerotol) sets the zerotol. Default: elkZerotol
validateDec(.., 'verbose', verbose) sets the verbosity. For verbose=0,
no messages are displayed, verbose=1 only states which tests went wrong
and verbose=2 also displays the succesful tests including some
statistical data.
validateDec(.., 'output', out) specifies the output that is returned.
With out='dec' the decomposition structure is returned. With
out='result' only the test result is returned. With out='both' both
results are returned using [dec result] = validateDec(..).
Quick checks that are redundant to the more time consuming ones:
(Q1) All structuring elements represent an extreme ray of at least one
unified partition. The algorithm stores indecomposable
polytopes, even if the partition is rejected. A failure here
may result from the Minkowski decomposition (the polytope is
inproper) or from an incomplete decomposition (Q4 should fail
as well).
(Q2) All structuring elements are part of the cone of valid polytopes.
(Q3) No volume of a unified partition is zero. Together with (L2), this
ensures that the intersection of two unified partitions has at
most dimension (nC-1).
(Q4) The sum of all volumes of the unified partitions equals the volume
of the valid cone. Together with (L2) this ensures that the
decomposition of the valid cone is complete.
Tests that take some more time:
(L1) All structuring elements are (constrained) indecomposable. This
double-checks the validity of the structuring elements. A
failure here maps to a failure in Minkowski decomposition.
(L2) All intersections of two unified partitions have the volume 0.
Together with (Q3), this ensures that the intersection of two
unified partitions has at most dimension (nC-1).
(L3) If the Minkowski decomposition for an inner point (not the
boundary) of a unified partition is computed. The resulting
summands equal the extreme rays of the unified partition. The
partition was found this way. A failure here maps to a failure
in Minkowski decomposition.
See also: obtainDec, validateMeasureDataDec0001 function varargout = validateDec(dec, varargin) 0002 % validate decomposition of valid cone 0003 % 0004 % Syntax: dec = validateDecomposition(dec) 0005 % 0006 % This function checks whether dec is a proper decomposition of the set of 0007 % facet normals in dec.A. The structuring elements are double-checked for 0008 % (constrained) indecomposability and the unified partitions are checked 0009 % against the cone of valid polytopes. The simplex partitions are not 0010 % verified as they are only used internally. 0011 % 0012 % The decomposition structure is extended with the fields: 0013 % dec.validation.tolerance - chosen zerotol 0014 % dec.validation.result - scalar logical 1 (passed) or 0 (failed) 0015 % dec.timing.validateDec - total time that was taken by validation 0016 % 0017 % validateDec(.., 'zerotol', zerotol) sets the zerotol. Default: elkZerotol 0018 % 0019 % validateDec(.., 'verbose', verbose) sets the verbosity. For verbose=0, 0020 % no messages are displayed, verbose=1 only states which tests went wrong 0021 % and verbose=2 also displays the succesful tests including some 0022 % statistical data. 0023 % 0024 % validateDec(.., 'output', out) specifies the output that is returned. 0025 % With out='dec' the decomposition structure is returned. With 0026 % out='result' only the test result is returned. With out='both' both 0027 % results are returned using [dec result] = validateDec(..). 0028 % 0029 % Quick checks that are redundant to the more time consuming ones: 0030 % (Q1) All structuring elements represent an extreme ray of at least one 0031 % unified partition. The algorithm stores indecomposable 0032 % polytopes, even if the partition is rejected. A failure here 0033 % may result from the Minkowski decomposition (the polytope is 0034 % inproper) or from an incomplete decomposition (Q4 should fail 0035 % as well). 0036 % (Q2) All structuring elements are part of the cone of valid polytopes. 0037 % (Q3) No volume of a unified partition is zero. Together with (L2), this 0038 % ensures that the intersection of two unified partitions has at 0039 % most dimension (nC-1). 0040 % (Q4) The sum of all volumes of the unified partitions equals the volume 0041 % of the valid cone. Together with (L2) this ensures that the 0042 % decomposition of the valid cone is complete. 0043 % 0044 % Tests that take some more time: 0045 % (L1) All structuring elements are (constrained) indecomposable. This 0046 % double-checks the validity of the structuring elements. A 0047 % failure here maps to a failure in Minkowski decomposition. 0048 % (L2) All intersections of two unified partitions have the volume 0. 0049 % Together with (Q3), this ensures that the intersection of two 0050 % unified partitions has at most dimension (nC-1). 0051 % (L3) If the Minkowski decomposition for an inner point (not the 0052 % boundary) of a unified partition is computed. The resulting 0053 % summands equal the extreme rays of the unified partition. The 0054 % partition was found this way. A failure here maps to a failure 0055 % in Minkowski decomposition. 0056 % 0057 % See also: obtainDec, validateMeasureDataDec 0058 0059 % The elk-library: convex geometry applied to crystallization modeling. 0060 % Copyright (C) 2012 Alexander Reinhold 0061 % 0062 % This program is free software: you can redistribute it and/or modify it 0063 % under the terms of the GNU General Public License as published by the 0064 % Free Software Foundation, either version 3 of the License, or (at your 0065 % option) any later version. 0066 % 0067 % This program is distributed in the hope that it will be useful, but 0068 % WITHOUT ANY WARRANTY; without even the implied warranty of 0069 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0070 % General Public License for more details. 0071 % 0072 % You should have received a copy of the GNU General Public License along 0073 % with this program. If not, see <http://www.gnu.org/licenses/>. 0074 0075 % timing 0076 tic 0077 0078 %% Input 0079 % options 0080 opt = mapOptionStruct(varargin, ... 0081 'verbose', 2, ... 0082 'zerotol', elkZerotol, ... 0083 'output', 'undefined', ... 0084 'L1', 1); 0085 0086 if strcmpi(opt.output, 'undefined') 0087 if nargout == 0 0088 opt.output = 'none'; 0089 elseif nargout == 1 0090 opt.output = 'dec'; 0091 elseif nargout >= 2 0092 opt.output = 'both'; 0093 end 0094 end 0095 0096 0097 % if no test fails, accept decomposition 0098 acceptDec = 1; 0099 0100 % Get volumes in a vector 0101 partitionVolumeVector = zeros(1, dec.nPU); 0102 unifiedSeIndexMatrix = zeros(dec.nPU, dec.nS); 0103 for iPartOne = 1:dec.nPU 0104 partitionVolumeVector(iPartOne) = dec.unifiedPartition{iPartOne}.volumeScaled; 0105 unifiedSeIndexMatrix(iPartOne, dec.unifiedPartition{iPartOne}.seIndexVector) = 1; 0106 end 0107 0108 %% Output volume 0109 if opt.verbose == 2 0110 disp(['Domain volume is: ' num2str(dec.confinementCone.volume)]); 0111 disp(['Sum of individual scaled partition volumes: ' ... 0112 num2str(sum(partitionVolumeVector))]); 0113 disp(['Minimal scaled partition volume is: ' ... 0114 num2str(min(partitionVolumeVector))]); 0115 disp(['Maximal scaled partition volume is: ' ... 0116 num2str(max(partitionVolumeVector))]); 0117 end 0118 0119 %% Q1 - all structuring elements used 0120 seUsedCount = sum(unifiedSeIndexMatrix, 1); 0121 if any(seUsedCount == 0) 0122 if opt.verbose 0123 disp(['[-] Q1: The structuring element(s) ' ... 0124 num2str(find(seUsedCount == 0)) ... 0125 ' are no extreme ray of any partition']); 0126 end 0127 acceptDec = 0; 0128 else 0129 if opt.verbose == 2 0130 disp('[+] Q1: All structuring elements are used for at least one partition'); 0131 end 0132 end 0133 0134 %% Q2 - all structuring elements are in the cone of valid polytopes 0135 thisTestAccepted = 1; 0136 for iSe = 1:dec.nS 0137 if all(dec.confinementCone.A * dec.seScaledToClose(iSe, :)' < opt.zerotol) 0138 % nothing to do, everything OK 0139 else 0140 if opt.verbose 0141 disp(['[-] Q2: The following structuring element is not in the valid cone:' ... 0142 num2str(iSe)]); 0143 end 0144 thisTestAccepted = 0; 0145 acceptDec = 0; 0146 end 0147 end 0148 % print something, if everything is alright 0149 if opt.verbose == 2 && thisTestAccepted 0150 disp('[+] Q2: Every structuring element is contained in the valid cone'); 0151 end 0152 0153 %% Q3 - no zero volume of partitions 0154 if any(isZero(partitionVolumeVector, opt.zerotol)) 0155 if opt.verbose 0156 disp(['[-] Q3: The volume of the following unified partition(s) is zero: '... 0157 num2str(isZero(partitionVolumeVector, opt.zerotol))]); 0158 end 0159 acceptDec = 0; 0160 else 0161 if opt.verbose == 2 0162 disp('[+] Q3: No partition volume equals zero'); 0163 end 0164 end 0165 0166 %% Q4 - sum of volumes 0167 if isZero(sum(partitionVolumeVector) - 1, opt.zerotol) 0168 if opt.verbose == 2 0169 disp('[+] Q4: The sum of partition volumes equals the domain volume'); 0170 end 0171 else 0172 if opt.verbose 0173 disp('[-] Q4: The sum of partition volumes does not equal the domain volume'); 0174 disp([' ..the scaled difference (sum(volumes) - 1) is: ' ... 0175 num2str(sum(partitionVolumeVector)-1)]); 0176 end 0177 acceptDec = 0; 0178 end 0179 0180 %% L1 - all structuring elements are indecomposable 0181 if opt.L1 0182 thisTestAccepted = 1; 0183 thisHrep.A = dec.A; 0184 for iSe = 1:dec.nS 0185 thisHrep.h = dec.mappingReducedToFull * dec.seScaledToClose(iSe, :)'; 0186 summandMatrix = decomposeMinkowskiHrep(thisHrep, ... 0187 'constraint', dec.constraintMatrix, ... 0188 'result', 'indecomposable', 'matrix', 1, ... 0189 'zerotol', opt.zerotol); 0190 summandMatrix = (dec.mappingFullToReduced * summandMatrix')'; 0191 summandMatrix = scaleToClose(summandMatrix, dec); 0192 if size(summandMatrix, 1) == 1 && ... 0193 all(isZero(summandMatrix - dec.seScaledToClose(iSe, :), opt.zerotol)) 0194 % nothing to do, everything OK 0195 else 0196 if opt.verbose 0197 disp(['[-] L1: The following structuring element is not (constrained) indecomposable: ' ... 0198 num2str(iSe)]); 0199 end 0200 thisTestAccepted = 0; 0201 acceptDec = 0; 0202 end 0203 end 0204 % print something, if everything is alright 0205 if opt.verbose == 2 && thisTestAccepted 0206 disp('[+] L1: Every structuring element is (constrained) indecomposable'); 0207 end 0208 end 0209 0210 %% L2 - All intersections of two unified partitions have the volume 0. 0211 % totalIntersectionVolume = 0; 0212 thisTestAccepted = 1; 0213 for iPartOne = 1:dec.nPU 0214 hrepOne = closeCone(dec.unifiedPartition{iPartOne}, dec.centerRay); 0215 0216 for iPartTwo = (iPartOne+1):dec.nPU 0217 0218 hrepTwo = closeCone(dec.unifiedPartition{iPartTwo}, dec.centerRay); 0219 0220 hrepIntersection = computeIntersectionHrep(hrepOne, hrepTwo, 0); 0221 0222 hrepIntersectionRed = reduceHrepDimension(hrepIntersection, opt.zerotol); 0223 0224 % volumeIntersection = computeVolumeHrep(hrepIntersection, zerotol); 0225 0226 % if volumeIntersection > zerotol * dec.confinementCone.volume 0227 if size(hrepIntersectionRed.A, 2) == size(hrepIntersection.A, 2) 0228 disp(['[-] L2: Parition ' num2str(iPartOne) ' and ' num2str(iPartTwo) ' have a full-dimensional intersection']) 0229 % disp([' ..the volume of this intersection is: ' ... 0230 % num2str(volumeIntersection/dec.confinementCone.volume)]) 0231 % disp([' ..there individual volumes are ' ... 0232 % num2str(partitionVolumeVector(iPartOne)) ... 0233 % ' and ' ... 0234 % num2str(partitionVolumeVector(iPartTwo)) ... 0235 % ', repsectively']); 0236 % totalIntersectionVolume = totalIntersectionVolume + volumeIntersection; 0237 thisTestAccepted = 0; 0238 acceptDec = 0; 0239 else 0240 % nothing to do, everything is fine 0241 end 0242 0243 end 0244 end 0245 if opt.verbose == 2 && thisTestAccepted 0246 disp('[+] L2: No two partitions have a full-dimensional intersection'); 0247 end 0248 0249 %% L3 - Minkowski decomposition for an inner point 0250 thisTestAccepted = 1; 0251 thisHrep.A = dec.A; 0252 for iPart = 1:dec.nPU 0253 thisSummandMatrix = dec.seScaledToClose(... 0254 dec.unifiedPartition{iPart}.seIndexVector, :); 0255 0256 thisHrep.h = dec.mappingReducedToFull * sum(thisSummandMatrix, 1)'; 0257 0258 summandMatrix = decomposeMinkowskiHrep(thisHrep, ... 0259 'constraint', dec.constraintMatrix, ... 0260 'result', 'indecomposable', 'matrix', 1, ... 0261 'zerotol', opt.zerotol); 0262 0263 summandMatrix = (dec.mappingFullToReduced * summandMatrix')'; 0264 0265 summandMatrix = scaleToClose(summandMatrix, dec); 0266 0267 if size(summandMatrix, 1) == size(thisSummandMatrix, 1) && ... 0268 isempty(eliminateRows(thisSummandMatrix, summandMatrix, opt.zerotol)) 0269 % everything OK 0270 else 0271 if opt.verbose 0272 disp(['[-] L3: The Minkowski decomposition of an inner point of ' ... 0273 'parition ' num2str(iPart) ' does not result in the ' ... 0274 'extreme rays of that partition']); 0275 end 0276 0277 thisTestAccepted = 0; 0278 acceptDec = 0; 0279 end 0280 end 0281 if opt.verbose == 2 && thisTestAccepted 0282 disp('[+] L3: Every decomposition of an inner point of a unified partition results in that partition'); 0283 end 0284 0285 %% Output 0286 dec.validationData.tolerance = opt.zerotol; 0287 dec.validationData.result = acceptDec; 0288 dec.timingData.validateDec = toc; 0289 switch lower(opt.output) 0290 case 'none' 0291 % nothing to do 0292 case 'dec' 0293 varargout{1} = dec; 0294 case 'result' 0295 varargout{1} = acceptDec; 0296 case 'both' 0297 varargout{1} = dec; 0298 varargout{2} = acceptDec; 0299 otherwise 0300 error('elk:decomposition:inputError', ... 0301 'The output option must be one of: ''dec'', ''result'', ''both''.'); 0302 end