determine the partition of a given hC-vector
Syntax: p = computePartitionIndexDec(dec, hcMatrix)
The simplex partition (the index to dec.simplexPartition) is returned so
that dec.simplexPartition{p}.A*hC <= 0 is fulfilled.
getMeasure(..., zerotol, unified) lets you control the tolerance
(zerotol) that is used for algorithms (default elkZerotol) and lets the
function alternatively return the unified partition (the index to
dec.unifiedPartition) for unified=1 (defaulf unified=0).
See also: obtainDec, validateDec, addMeasureDataDec, validateMeasureDec,
computeMeasureDec0001 function p = computePartitionIndexDec(dec, hcMatrix, zerotol, unified) 0002 % determine the partition of a given hC-vector 0003 % 0004 % Syntax: p = computePartitionIndexDec(dec, hcMatrix) 0005 % 0006 % The simplex partition (the index to dec.simplexPartition) is returned so 0007 % that dec.simplexPartition{p}.A*hC <= 0 is fulfilled. 0008 % 0009 % getMeasure(..., zerotol, unified) lets you control the tolerance 0010 % (zerotol) that is used for algorithms (default elkZerotol) and lets the 0011 % function alternatively return the unified partition (the index to 0012 % dec.unifiedPartition) for unified=1 (defaulf unified=0). 0013 % 0014 % See also: obtainDec, validateDec, addMeasureDataDec, validateMeasureDec, 0015 % computeMeasureDec 0016 0017 % The elk-library: convex geometry applied to crystallization modeling. 0018 % Copyright (C) 2012 Alexander Reinhold 0019 % 0020 % This program is free software: you can redistribute it and/or modify it 0021 % under the terms of the GNU General Public License as published by the 0022 % Free Software Foundation, either version 3 of the License, or (at your 0023 % option) any later version. 0024 % 0025 % This program is distributed in the hope that it will be useful, but 0026 % WITHOUT ANY WARRANTY; without even the implied warranty of 0027 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0028 % General Public License for more details. 0029 % 0030 % You should have received a copy of the GNU General Public License along 0031 % with this program. If not, see <http://www.gnu.org/licenses/>. 0032 0033 %% Input handling 0034 % ensure hc being a column vector or matrix of hc-vectors 0035 if size(hcMatrix, 1) ~= dec.nC 0036 if size(hcMatrix, 2) == dec.nC 0037 hcMatrix = hcMatrix'; 0038 else 0039 error('elk:decomposition:wrongInput', ... 0040 ['the provided vector must be a constrained h-vector that does ' ... 0041 'have the dimension size(dec.se, 2)']); 0042 end 0043 end 0044 0045 % set zerotol 0046 if ~exist('zerotol', 'var') || isempty(zerotol) 0047 zerotol = elkZerotol*ones(1, size(hcMatrix,2)); 0048 elseif length(zerotol) < size(hcMatrix, 2) 0049 zerotol = zerotol*ones(1, size(hcMatrix,2)); 0050 end 0051 % set unified option 0052 if ~exist('unified', 'var') || isempty(unified) 0053 unified = 0; 0054 end 0055 0056 % catch incomplete parititions 0057 if ~dec.isComplete 0058 p = ones(1, size(hcMatrix, 2)); 0059 if ~unified 0060 p = dec.unifiedToSimplexVector*p; 0061 end 0062 return 0063 end 0064 0065 %% gather validity/confinement information 0066 % ensure that vector is in cone.. ..the confinement cone was 0067 % renamed according to the type of hC-representation to get 0068 % pin-pointed where the type makes a difference. 0069 if dec.isComplete && dec.isOrdinary 0070 isConfinedVector = max(dec.confinementCone.A * hcMatrix, [], 1) < zerotol; 0071 elseif dec.isComplete && ~dec.isOrdinary 0072 isConfinedVector = max(dec.confinementConeHull.A * hcMatrix, [], 1) < zerotol; 0073 else 0074 error('elk:decomposition:internalError', 'This should not happen'); 0075 end 0076 0077 %% Determine partition 0078 % the default is: "not in cone" 0079 pU = zeros(1, size(hcMatrix, 2)); 0080 pS = zeros(1, size(hcMatrix, 2)); 0081 for iVector = 1:size(hcMatrix, 2) 0082 %! DEV: the partition search is definitely a candidate for MEX 0083 % implementation. The numerical integrations in cs-solver with 0084 % Vegas took about 0.33 seconds for the tensor computations while 0085 % 158 seconds (almost all time) went into the identification of 0086 % the unified partition. It was a 7D paracetamol example that - 0087 % indeed had a lot of partitions to search. 0088 % if mod(iVector, 100) == 0, disp(iVector/size(hc, 2));end 0089 %/ 0090 0091 if unified && isConfinedVector(iVector) 0092 % partition invalid if not found: 0093 pU(iVector) = -1; 0094 for iPart = 1:length(dec.unifiedPartition) 0095 if all(dec.unifiedPartition{iPart}.A * hcMatrix(:, iVector) < zerotol(iVector)) 0096 pU(iVector) = iPart; 0097 break; 0098 end 0099 end 0100 elseif isConfinedVector(iVector) 0101 % partition invalid if not found: 0102 pS(iVector) = -1; 0103 % find simplex partition 0104 for iPart = 1:length(dec.simplexPartition) 0105 if prod(double((dec.simplexPartition{iPart}.A * hcMatrix(:, iVector)) < zerotol(iVector))) 0106 pS(iVector) = iPart; 0107 break; 0108 end 0109 end 0110 end 0111 0112 end 0113 0114 % the following seems to be odd, but adds some safety to the current 0115 % implementation 0116 if unified 0117 p = pU; 0118 else 0119 p = pS; 0120 end 0121 0122 % error handling 0123 if any(p == 0) 0124 error('elk:decomposition:wrongInput', ... 0125 ['the given vector is outside of the confinement cone: ' ... 0126 ' the vector must fulfill (dec.validityCone.A * vec <= 0)']); 0127 elseif any(p == -1) && dec.isOrdinary 0128 % this error could mean that the decomposition is not complete 0129 error('elk:decomposition:numericalProblem', ... 0130 ['the simplex partition of the given vector could not be found - ' ... 0131 'this should not happen.']); 0132 elseif any(p == -1) && ~dec.isOrdinary 0133 warning('elk:decomposition:generalProblem', ... 0134 ['the unified/simplex partition of the given vector(s) could ' ... 0135 'not be found - the input decomposition is singular (not ' ... 0136 'ordinary) so that this can happen.']); 0137 end 0138 0139 end