calculate measure from decomposition data
Syntax: measure = computeMeasureDec(dec, 'hc', hC, measure)
or measure = computeMeasureDec(dec, 'se', {p, lam}, measure)
computeMeasureDec(dec, 'hC', hC) calculates the measure based on the
hC-vector of the polytope (Hc-representation). This vector must have
the length dec.nC.
computeMeasureDec(dec, 'hC', {p, hC}) also provides the simplex
partitions for the vectors in hC.
computeMeasureDec(SR, 'se', {p, lam}) calculates the measure bases on the
vector lam that represents the linear combination of the structuring
elements in the simplex partition that is given by p. The vector lam
must have the length dec.nC in this case.
The simple forms above are sufficient and simply calculate all measures
that are defined in dec.measureData.
getMeasure(.., measure, ..) specifies the measure name (the fields
in dec.measureData) and eventually the parameter par that is required.
If this list starts with 'grad' or 'gradient' the gradient of the
subsequent measures are calculated. Calculating measures and measure
gradients cannot be combined.
Examples: {'volume'}, {'volume', 'feret', [1 0 0;1 1 0]}, {'grad',
'surface'}, {'grad', 'volume', 'support', [1 0 0], 'feret',
[1 1 1; 0 0 1], 'facet', [1 0 0]}
getMeasure(..., zerotol, simplex) lets you control the tolerance
(zerotol) that is used for algorithms (default elkZerotol) and enforces
the algorithm to search and use the proper simplex partition
(simplex=1) instead of the unified partition (simplex=0). This also
effects how p in {p, hc} is interpreted (default is simplex=0).
See also: obtainMeasureInfo, validateMeasureDataDec, obtainDec,
validateDec0001 function M = computeMeasureDec(dec, vecType, vec, measureList, ... 0002 zerotol, useSimplex, enforceValid) 0003 % calculate measure from decomposition data 0004 % 0005 % Syntax: measure = computeMeasureDec(dec, 'hc', hC, measure) 0006 % or measure = computeMeasureDec(dec, 'se', {p, lam}, measure) 0007 % 0008 % computeMeasureDec(dec, 'hC', hC) calculates the measure based on the 0009 % hC-vector of the polytope (Hc-representation). This vector must have 0010 % the length dec.nC. 0011 % 0012 % computeMeasureDec(dec, 'hC', {p, hC}) also provides the simplex 0013 % partitions for the vectors in hC. 0014 % 0015 % computeMeasureDec(SR, 'se', {p, lam}) calculates the measure bases on the 0016 % vector lam that represents the linear combination of the structuring 0017 % elements in the simplex partition that is given by p. The vector lam 0018 % must have the length dec.nC in this case. 0019 % 0020 % The simple forms above are sufficient and simply calculate all measures 0021 % that are defined in dec.measureData. 0022 % 0023 % getMeasure(.., measure, ..) specifies the measure name (the fields 0024 % in dec.measureData) and eventually the parameter par that is required. 0025 % If this list starts with 'grad' or 'gradient' the gradient of the 0026 % subsequent measures are calculated. Calculating measures and measure 0027 % gradients cannot be combined. 0028 % Examples: {'volume'}, {'volume', 'feret', [1 0 0;1 1 0]}, {'grad', 0029 % 'surface'}, {'grad', 'volume', 'support', [1 0 0], 'feret', 0030 % [1 1 1; 0 0 1], 'facet', [1 0 0]} 0031 % 0032 % getMeasure(..., zerotol, simplex) lets you control the tolerance 0033 % (zerotol) that is used for algorithms (default elkZerotol) and enforces 0034 % the algorithm to search and use the proper simplex partition 0035 % (simplex=1) instead of the unified partition (simplex=0). This also 0036 % effects how p in {p, hc} is interpreted (default is simplex=0). 0037 % 0038 % See also: obtainMeasureInfo, validateMeasureDataDec, obtainDec, 0039 % validateDec 0040 0041 % NOT YET SUPPORTED AND, HENCE, NOT DOCUMENTED: 0042 % 0043 % Symbolic vectors hc or se can be used to obtain the polynomials that 0044 % correspond to the specified partitions (columns) and the specified 0045 % measures (rows). The output is a symbolic matrix. For the input, it 0046 % should be a cell string {p, vec} where p is a vector of partition 0047 % indices and vec is a single symbolic vector. 0048 % Status: draft 0049 0050 % The elk-library: convex geometry applied to crystallization modeling. 0051 % Copyright (C) 2012 Alexander Reinhold 0052 % 0053 % This program is free software: you can redistribute it and/or modify it 0054 % under the terms of the GNU General Public License as published by the 0055 % Free Software Foundation, either version 3 of the License, or (at your 0056 % option) any later version. 0057 % 0058 % This program is distributed in the hope that it will be useful, but 0059 % WITHOUT ANY WARRANTY; without even the implied warranty of 0060 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0061 % General Public License for more details. 0062 % 0063 % You should have received a copy of the GNU General Public License along 0064 % with this program. If not, see <http://www.gnu.org/licenses/>. 0065 0066 %% make sure mex-file exists 0067 try 0068 evaluateMeasurePolynome([], [], 1); 0069 % if ~exist([mexPath 'evaluateMeasurePolynome.' mexext], 'file') 0070 catch noError 0071 mexPath = mfilename('fullpath'); 0072 mexPath = [mexPath(1:(end-17)) 'private' filesep]; 0073 try 0074 compileMessage = evalc(['mex ' mexPath 'evaluateMeasurePolynome.c -outdir ' mexPath]); 0075 catch thisError 0076 % disp(compileMessage) 0077 rethrow(thisError) 0078 end 0079 try 0080 rehashMessage = evalc('rehash'); 0081 catch thisError 0082 disp(rehashMessage) 0083 rethrow(thisError) 0084 end 0085 end 0086 0087 %% Input handling 0088 % set zerotol 0089 if ~exist('zerotol', 'var') || isempty(zerotol) 0090 zerotol = elkZerotol; 0091 end 0092 % set simplex option 0093 if ~exist('useSimplex', 'var') || isempty(useSimplex) 0094 useSimplex = 0; 0095 end 0096 % set enforceValid option 0097 if ~exist('enforceValid', 'var') || isempty(enforceValid) 0098 enforceValid = 0; 0099 end 0100 % enforce measureName to cell 0101 if ~iscell(measureList) 0102 measureList = {measureList}; 0103 end 0104 % if input decomposition is not proper, we fall back to the proper 0105 % embedding 0106 if isfield(dec, 'properData') && ~isfield(dec, 'measureData') 0107 if strcmpi(vecType, 'hc') 0108 if iscell(vec) 0109 warning('elk:decomposition:input', ['You have provided the ' ... 0110 'unified partition for an improper decomposition with no ' ... 0111 'measureData field. I compute the measure based on the ' ... 0112 'proper embedding, but ignore the provided partition data.']); 0113 vec = vec{2}; 0114 end 0115 M = computeMeasureDec(dec.properData, vecType, ... 0116 dec.properData.mappingNewToProper*vec, ... 0117 measureList, zerotol, useSimplex, enforceValid); 0118 return 0119 else 0120 error('elk:decomposition:wrongInput', ['The input decomposition ' ... 0121 'does not contain measure data. Utilizing the proper embedding ' ... 0122 'is only supported for hC-vectors.']); 0123 end 0124 end 0125 0126 0127 [vecType, vec, partitionIndexVector] = clarifyVectorInformation(... 0128 lower(vecType), vec, useSimplex, dec, ... 0129 zerotol, enforceValid); 0130 % the resulting partition index is always a simplex partition(!) 0131 nVec =length(partitionIndexVector); 0132 0133 % set symbolic handler 0134 isSymbolic = strcmp(class(vec), 'sym'); 0135 0136 %% strip gradient option 0137 % if the measureList starts with 'grad' or 'gradient', the partial 0138 % derivative of all following measures is calculated. 0139 if strcmpi(measureList{1}, 'grad') || strcmpi(measureList{1}, 'gradient') 0140 measureList(1) = []; 0141 gradientOption = 1; 0142 else 0143 gradientOption = 0; 0144 end 0145 0146 %% determine total number of measures 0147 % the created vector has the number of measures assigned to each 0148 % position where a string would be in measName. 0149 nMeasureVec = nan(1, length(measureList)); 0150 isParVec = nMeasureVec; 0151 for iMeasure = 1:length(measureList) 0152 if ischar(measureList{iMeasure}) 0153 nMeasureVec(iMeasure) = 1; 0154 isParVec(iMeasure) = 0; 0155 else 0156 nMeasureVec(iMeasure-1) = size(measureList{iMeasure}, 1); 0157 isParVec(iMeasure-1) = 1; 0158 nMeasureVec(iMeasure) = 0; 0159 isParVec(iMeasure) = nan; 0160 end 0161 end 0162 nMeasure = sum(nMeasureVec); 0163 0164 %% Initialize output 0165 caseList = {'measure', 'gradient', 'symbolic measure', 'symbolic gradient'}; 0166 outputCase = caseList{1.0 + logical(gradientOption)*1 + logical(isSymbolic)*2}; 0167 switch outputCase 0168 case 'measure' 0169 M = zeros(nMeasure, nVec); 0170 case 'gradient' 0171 M = zeros(nMeasure, nVec, dec.nC); 0172 case 'symbolic measure' 0173 M = sym('tmp', [nMeasure, nVec])*0; 0174 case 'symbolic gradient' 0175 error('elk:decomposition:wrongInput', ... 0176 'Symbolic gradients for measure calculation are not supported.'); 0177 otherwise 0178 error('elk:decomposition:internalError', ... 0179 'This should not happen at all.'); 0180 end 0181 0182 %% Prepare "easy-way-through (TM)" 0183 % In the following there are numurous decisions to make that include 0184 % optional extra cycles: 0185 % 0186 % (1) cycle: measure names 0187 % (2) cycle: parameters (if so) 0188 % (3) cycle: gradients (if so) 0189 % (4) cycle: vectors 0190 % 0191 % To simplify the loop and cross-decisions, we summarize the measure and 0192 % parameter loop. This needs some preparation: 0193 measureIndexVector = []; 0194 parameterIndexVector = []; 0195 for iMeasure = find(nMeasureVec) 0196 measureIndexVector = [measureIndexVector ... 0197 iMeasure*ones(1, nMeasureVec(iMeasure))]; 0198 parameterIndexVector = [parameterIndexVector 1:nMeasureVec(iMeasure)]; 0199 end 0200 0201 % cycle through measures/parameters AND simplex partition indices 0202 if useSimplex 0203 partitionIndexCycleVector = 1:dec.nPS; 0204 else 0205 partitionIndexCycleVector = dec.unifiedToSimplexVector(1:dec.nPU); 0206 % note that iPart is always a simplex partition (see below) while 0207 % partitionIndexVector denotes either simplex or unified partitions. 0208 % This is corrected, here. 0209 end 0210 0211 for iMeasPar = 1:length(measureIndexVector) 0212 for iPart = partitionIndexCycleVector 0213 vectorFilter = (partitionIndexVector == iPart); 0214 0215 % handle no vectors to deal with 0216 if ~any(vectorFilter) 0217 continue; 0218 end 0219 0220 % handle invalid vectors (they have partitionIndex 0) 0221 % if iPart == 0 0222 % % there is nothing to do since all fields are prefilled with zeros. 0223 % continue; 0224 % end 0225 0226 % determine tensor with/without parameter 0227 if isParVec(iMeasure) 0228 % obtain tensor (with parameter) 0229 [thisTensor tensorDim] = obtainTensorDec(dec, ... 0230 measureList{measureIndexVector(iMeasPar)}, ... 0231 measureList{measureIndexVector(iMeasPar)+1}(... 0232 parameterIndexVector(iMeasPar), : ... 0233 ), vecType, iPart); 0234 else 0235 % obtain tensor (without parameter) 0236 [thisTensor tensorDim] = obtainTensorDec(dec, ... 0237 measureList{measureIndexVector(iMeasPar)}, [], ... 0238 vecType, iPart); 0239 end 0240 0241 %! Debug: Here, I would need to distinguish symbolic and normal output 0242 % But "M" must be initialized properly in that case. .. enforce this 0243 % cast once, whenever a symbolic vector is met. 0244 % > if [vector is symbolic] and [M not symbolic] 0245 % > then [cast to symbolic] 0246 % > end 0247 % Hu?! Since vec is a matrix with all vectors, the symbolic thing is a 0248 % single decision based on the input-vec only. This can already be 0249 % considered with the initialization of M!! 0250 % 0251 % 1) Decide on symbolic/nonsymbolic and thisOutputCase in 0252 % initialization step 0253 % 2) Handle calculations case-based like below 0254 % 3) That's it - Isn't it?? 0255 %/ 0256 switch outputCase 0257 case 'measure' 0258 M(iMeasPar, vectorFilter) = evaluateMeasurePolynome(... 0259 thisTensor, vec(:, vectorFilter), tensorDim); 0260 case 'gradient' 0261 for iGradDim = 1:dec.nC 0262 switch tensorDim 0263 case 1 0264 M(iMeasPar, iVec, iGradDim) = thisTensor(iGradDim); 0265 case 2 0266 tmpTensor = permute(thisTensor(iGradDim, :), [2 1]) ... 0267 + thisTensor(:, iGradDim); 0268 M(iMeasPar, vectorFilter, iGradDim) = ... 0269 evaluateMeasurePolynome(... 0270 tmpTensor, vec(:, vectorFilter), 1); 0271 case 3 0272 tmpTensor = permute(thisTensor(iGradDim, :, :), [2 3 1]) ... 0273 + permute(thisTensor(:, iGradDim, :), [1 3 2]) ... 0274 + thisTensor(:, :, iGradDim); 0275 M(iMeasPar, iVec, iGradDim) = evaluateMeasurePolynome(... 0276 tmpTensor, vec(:, vectorFilter), 2); 0277 end 0278 end 0279 case 'symbolic measure' 0280 %! Dev: Obtaining symbolic expressions was implemented for 0281 % shape approximation algorithms that used polynomial 0282 % optimization. Since this approach did not work 0283 % sufficiently well, the following code was not adopted. 0284 % Minor changes would be required. The outer loop went over 0285 % all vectors (iVec) before and is now running over all 0286 % partitions (iPart). 0287 error('elk:decomposition:wrongInput', ['Extracting symbolic ' ... 0288 'expressions is currently not supported']); 0289 % zeroTensor = thisTensor*0; 0290 % thisTensor = sym('tmp', size(thisTensor(:)'))*0 + thisTensor(:)'; 0291 % switch tensorDim 0292 % case 1 0293 % for iDim = (1:size(dec.se, 2)) 0294 % thisTensor(iDim) = thisTensor(iDim) * vec(iDim); 0295 % end 0296 % case 2 0297 % for iDim = (1:size(dec.se, 2)) 0298 % thisFilter = zeroTensor;thisFilter(iDim,:) = 1; 0299 % thisFilter = thisFilter(:) == 1; 0300 % thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim); 0301 % 0302 % thisFilter = zeroTensor;thisFilter(:,iDim) = 1; 0303 % thisFilter = thisFilter(:) == 1; 0304 % thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim); 0305 % end 0306 % case 3 0307 % for iDim = (1:size(dec.se, 2)) 0308 % thisFilter = zeroTensor;thisFilter(iDim,:,:) = 1; 0309 % thisFilter = thisFilter(:) == 1; 0310 % thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim); 0311 % 0312 % thisFilter = zeroTensor;thisFilter(:,iDim,:) = 1; 0313 % thisFilter = thisFilter(:) == 1; 0314 % thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim); 0315 % 0316 % thisFilter = zeroTensor;thisFilter(:,:,iDim) = 1; 0317 % thisFilter = thisFilter(:) == 1; 0318 % thisTensor(thisFilter) = thisTensor(thisFilter) * vec(iDim); 0319 % end 0320 % end 0321 % % the meaning in this symbolic case is more iPartition than iVec 0322 % M(iMeasPar, iVec) = sum(thisTensor); 0323 otherwise 0324 error('elk:decomposition:internalError', 'This should not happen.'); 0325 end 0326 end 0327 end 0328 0329 end