generate the inequality constraints for valid polytopes in h-space THIS FUNCTION IS NO USER FUNCTION.
0001 function result = gatherValidityConeIneq(facetNormalMatrix, zerotol, verbose) 0002 % generate the inequality constraints for valid polytopes in h-space 0003 % 0004 % THIS FUNCTION IS NO USER FUNCTION. 0005 0006 % The elk-library: convex geometry applied to crystallization modeling. 0007 % Copyright (C) 2012 Alexander Reinhold 0008 % 0009 % This program is free software: you can redistribute it and/or modify it 0010 % under the terms of the GNU General Public License as published by the 0011 % Free Software Foundation, either version 3 of the License, or (at your 0012 % option) any later version. 0013 % 0014 % This program is distributed in the hope that it will be useful, but 0015 % WITHOUT ANY WARRANTY; without even the implied warranty of 0016 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0017 % General Public License for more details. 0018 % 0019 % You should have received a copy of the GNU General Public License along 0020 % with this program. If not, see <http://www.gnu.org/licenses/>. 0021 0022 %% info 0023 % This function generates inequality constraints in form of a 0024 % H-representation in the h-space that is associated with the provided 0025 % matrix of facet normals. For valid solutions the polytope is not empty 0026 % (at least contains a single point) and all support planes have a 0027 % nonempty intersection with the given polytope. The formal condition is 0028 % then: the polytope is valid for: 0029 % inequalityMatrix * h <= 0 0030 % This function also returns in cone.responsibleFacetMatrix a matrix that 0031 % stores pairs of inequality indices (first column) and positive indices 0032 % to the facet that grows out at this boundary. Negative indices are 0033 % entered if the dissolution of that facet triggers an empty polytope. 0034 0035 %% Input handling 0036 % some statistical numbers 0037 nDim = size(facetNormalMatrix, 2); 0038 nFacet = size(facetNormalMatrix, 1); 0039 0040 %% initialize matrices 0041 % generate tuples to check 0042 tupleMatrix = obtainOrderedSetMatrix(nFacet, nDim+1); 0043 % using a sparse matrix was tested, but does not pay off. 0044 inequalityMatrix = zeros(nDim, nFacet); 0045 responsibleFacetVector = nan(nDim, 1); 0046 thisInequalityIndex = 0; 0047 0048 %% start loop for all tuples 0049 if verbose,fprintf('......');end 0050 for iTuple = 1:size(tupleMatrix, 1) 0051 if verbose && mod(iTuple, 1000) == 1 0052 fprintf('\b\b\b\b\b\b%5.1f%%', iTuple / size(tupleMatrix, 1)*100); 0053 end 0054 0055 thisTuple = tupleMatrix(iTuple, :); 0056 0057 % find the n-tuple where the (n+1) element is in the positive OR 0058 % negative cone and the cone is full-dimensional. The last element is 0059 % contained in the positive span, if: 0060 % [a(1), .., a(n)] * lam = a(n+1) and lam >= 0 0061 % and in the negative span for lam <= 0. 0062 coneType = 0; % (not found) 0063 0064 for iFacet = thisTuple 0065 thisIndexVector = thisTuple(thisTuple ~= iFacet); 0066 [u,s,v] = svd(facetNormalMatrix(thisIndexVector, :)); 0067 s = diag(s); 0068 isNotZeroVector = (s <= -zerotol) | (s >= zerotol); 0069 if all(isNotZeroVector) 0070 s(isNotZeroVector) = 1./s(isNotZeroVector); 0071 lambda = u * diag(s) * v' * facetNormalMatrix(iFacet, :)'; 0072 0073 if all(lambda >= -zerotol) 0074 % tuplet combination found 0075 thisTuple = [thisIndexVector iFacet]; 0076 coneType = 1; % (positive cone) 0077 break 0078 elseif all(lambda <= zerotol) 0079 % tuplet combination found 0080 thisTuple = [thisIndexVector iFacet]; 0081 coneType = -1; % (negative cone) 0082 break 0083 end 0084 end 0085 end 0086 0087 %! Info: There is an alternative implementation, that uses fewer lines 0088 % of code and might be quicker to grasp, put it is not faster: 0089 0090 % % switch some warnings to errors to detect lower rank during mldivide 0091 % warningStateOne = warning('error', 'MATLAB:nearlySingularMatrix'); 0092 % warningStateTwo = warning('error', 'MATLAB:singularMatrix'); 0093 % 0094 % for iFacet = thisTuple 0095 % thisIndexVector = thisTuple(thisTuple ~= iFacet); 0096 % % use mldivide directly and catch low rank 0097 % try 0098 % lambda = facetNormalMatrix(thisIndexVector, :)' \ ... 0099 % facetNormalMatrix(iFacet, :)'; 0100 % catch 0101 % continue 0102 % end 0103 % %! Assignement as above! 0104 % end 0105 % % reset warning states 0106 % warning(warningStateOne); 0107 % warning(warningStateTwo); 0108 %/ 0109 0110 % result: we have evaluated the tuple to form a negative cone 0111 % (coneType=-1), positve cone (coneType=1) or none of these 0112 % (coneType=0). The array thisIndexVector also have the indces to the 0113 % cone-forming facet normals in the first three elements (that 0114 % defines the vertex). It contains the forth facet as the last 0115 % element. 0116 0117 % continue with next tuple (this tuple beeing invalid) or continue with 0118 % checking locality 0119 if coneType == 0 0120 %disp([num2str(iTuple) ': no full dimensional positive/negative span']) 0121 %thisNormalMatrix 0122 continue 0123 else 0124 % this tuple now contains the same (n+1) normal vectors that we 0125 % started with, but the first n generate the positive/negative 0126 % cone where the (n+1)st element is an element of. 0127 thisNormalMatrix = facetNormalMatrix(thisTuple, :); 0128 end 0129 0130 % the disabled code represents improvements by Karl (1995, "local 0131 % conditions") that were found to give no significant improvement. 0132 % % check if the n-tuple is local 0133 % % Another facet normal is part of the positive cone, whenever we have: 0134 % % [a(1), .., a(n)] * lam = a(i) and lam >= 0 0135 % % which can be rewritten to: 0136 % % inv([a(1), .., a(n)]) * a(i) >= 0 0137 % % the same facet normal is part of the negative cone, for: 0138 % % -inv([a(1), .., a(n)]) * a(i) >= 0 0139 % notInTuple = 1:nFacet; 0140 % notInTuple(thisTuple) = []; 0141 if coneType == 1 0142 testMatrix = inv(thisNormalMatrix(1:nDim, :)'); 0143 else 0144 testMatrix = inv(thisNormalMatrix(1:nDim, :)'); 0145 end 0146 % exitFlag = 0; 0147 % for iFacet = notInTuple 0148 % lambda = testMatrix * facetNormalMatrix(iFacet, :)'; 0149 % % as this is a "we save equations" condition, less samples should 0150 % % be valid to fulfill this condition. 0151 % if all(lambda >= zerotol) 0152 % % there is another facet in the positive/negative span. Hence, 0153 % % the tuple is not local. 0154 % exitFlag = 1; 0155 % break 0156 % end 0157 % end 0158 % if exitFlag 0159 % %disp([num2str(iTuple) ': not local positive span']) 0160 % %thisNormalMatrix 0161 % continue; 0162 % else 0163 % %disp([num2str(iTuple) ': add constraint']) 0164 % end 0165 0166 % now, add the condition 0167 % from Karl (1995) we have: 0168 % a(n+1)' * inv([a(1)'; ...; a(n)']) * [h(1); ..; h(n)] >= h(n+1) 0169 % which is equal to: 0170 % (inv([a(1), ..., a(n)]) * a(n+1))' * [h(1); ..; h(n)] >= h(n+1) 0171 % where the inverse of the second equation is the testMatrix from above 0172 % for positive cones. 0173 % if coneType == -1 0174 % testMatrix = -1*testMatrix; 0175 % end 0176 % Further more, for negative cones the inequality sign changes so that 0177 % we have for negative cones: 0178 % (inv([a(1), ..., a(n)]) * a(n+1))' * [h(1); ..; h(n)] - h(n+1) <= 0 0179 % 0180 % The common H-representation of the cone uses <= instead of >= so that 0181 % for the positive cone, we obtain: 0182 % -1*(inv([a(1), ..., a(n)]) * a(n+1))' * [h(1); ..; h(n)] + h(n+1) <= 0 0183 0184 thisInequalityIndex = thisInequalityIndex + 1; 0185 0186 % the following equals inv([a(1), ..., a(n)]) * a(n+1) 0187 coefficientVector = (testMatrix * thisNormalMatrix(end, :)')'; 0188 0189 if coneType == 1 0190 inequalityMatrix(thisInequalityIndex, thisTuple(end)) = 1; 0191 inequalityMatrix(thisInequalityIndex, thisTuple(1:(end-1))) = ... 0192 -1 * coefficientVector; 0193 responsibleFacetVector(thisInequalityIndex) = thisTuple(end); 0194 elseif coneType == -1 0195 inequalityMatrix(thisInequalityIndex, thisTuple(end)) = -1; 0196 inequalityMatrix(thisInequalityIndex, thisTuple(1:(end-1))) = ... 0197 coefficientVector; 0198 responsibleFacetVector(thisInequalityIndex) = -1*thisTuple(end); 0199 else 0200 error('elk:decomposition:numericalError', ... 0201 'cone is neither positive nor negative. This should not happen'); 0202 end 0203 0204 % reserve more space in sparce matrix, if required 0205 if thisInequalityIndex >= size(inequalityMatrix, 1) 0206 % allocate more space 0207 estimatedRowSize = thisInequalityIndex / iTuple * size(tupleMatrix, 1); 0208 % round and take some more 0209 estimatedRowSize = round(estimatedRowSize) + nDim; 0210 % enlarge matrix 0211 inequalityMatrix(estimatedRowSize, :) = 0; 0212 responsibleFacetVector((end+1):estimatedRowSize) = nan; 0213 end 0214 0215 end 0216 if verbose, fprintf('\b\b\b\b\b\b');end 0217 0218 %% remove the reservated rows 0219 inequalityMatrix((thisInequalityIndex+1):end, :) = []; 0220 responsibleFacetVector((thisInequalityIndex+1):end, :) = []; 0221 0222 result.A = inequalityMatrix; 0223 result.responsibleFacetVector = responsibleFacetVector;