calculates the set difference for a list of hrep's versus a single hrep THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE OF APPLICATIONS (decomposition case study; technical reports)
0001 function resultHrepList = subtractPolytopeHrep(... 0002 hrepMinuendList, hrepSubtrahend, zerotol) 0003 % calculates the set difference for a list of hrep's versus a single hrep 0004 % 0005 % THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 0006 % OF APPLICATIONS (decomposition case study; technical reports) 0007 % 0008 0009 %! ToDo: redo documentation 0010 % Syntax: vrepList = subtractPolytope(polyMinuend, polySubtrahend) 0011 % 0012 % The input minuend can be a V-representation, H-representation or a cell 0013 % array of these representations. While in the latter case it is expected 0014 % that this list contains only simplices. The subtrahend should be a 0015 % H-representation of V-representation only. 0016 % 0017 % subtractPolytope(.., zerotol) lets you specify the zerotol that is used 0018 % to identify values that are close to zero by zero. Default: elkZerotol. 0019 % 0020 % subtractPolytope(.., zerotol, isCone) is a specific logical parameter that 0021 % identifies the input as a convex cone that contains the origin as one 0022 % point. All other points must be at a hyperplane. This parameter is 0023 % passed to triangulateVrep. 0024 % 0025 % See also: triangulateVrep, identifyPolytope 0026 0027 % The elk-library: convex geometry applied to crystallization modeling. 0028 % Copyright (C) 2012 Alexander Reinhold 0029 % 0030 % This program is free software: you can redistribute it and/or modify it 0031 % under the terms of the GNU General Public License as published by the 0032 % Free Software Foundation, either version 3 of the License, or (at your 0033 % option) any later version. 0034 % 0035 % This program is distributed in the hope that it will be useful, but 0036 % WITHOUT ANY WARRANTY; without even the implied warranty of 0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0038 % General Public License for more details. 0039 % 0040 % You should have received a copy of the GNU General Public License along 0041 % with this program. If not, see <http://www.gnu.org/licenses/>. 0042 0043 %% Input handling 0044 % zerotol 0045 if ~exist('zerotol', 'var') 0046 zerotol = elkZerotol; 0047 end 0048 % sort the constraints of the subtrahend 0049 [hrepSubtrahend.A idxVector] = sortrows(hrepSubtrahend.A); 0050 hrepSubtrahend.h = hrepSubtrahend.h(idxVector); 0051 0052 % The principle is that every part of the minuend is treated seperately and 0053 % then we go through the disjoint convex compartements for the complement 0054 % of the subtrahend. 0055 % some numbers: 0056 nSubtrahendFacet = length(hrepSubtrahend.h); 0057 nMinuend = length(hrepMinuendList); 0058 nResult = 0; 0059 % common options for Chebychev ball computation: 0060 % problem.solver = 'linprog'; 0061 % problem.options = optimset('Display', 'off', 'largeScale', 'off', ... 0062 % 'TolFun', 0.01*zerotol, 'simplex', 'off'); 0063 % problem.f = [-1 zeros(1, size(hrepSubtrahend.A, 2))]; 0064 % problem.lb = -1; 0065 % problem.ub = 1; 0066 % problem.x0 = zeros(1 + size(hrepSubtrahend.A, 2), 1); 0067 problem.obj = [-1 zeros(1, size(hrepSubtrahend.A, 2))]; 0068 problem.lin = []; 0069 0070 0071 resultHrepList = cell(1, nSubtrahendFacet*nMinuend); 0072 0073 for iSubtrahendFacet = 1:length(hrepSubtrahend.h) 0074 for iMinuend = 1:length(hrepMinuendList) 0075 % construct hrep 0076 thisHrep.A = [hrepMinuendList{iMinuend}.A; ... 0077 -1*hrepSubtrahend.A(iSubtrahendFacet, :); ... 0078 hrepSubtrahend.A(1:(iSubtrahendFacet-1), :)]; 0079 thisHrep.h = [hrepMinuendList{iMinuend}.h; ... 0080 -1*hrepSubtrahend.h(iSubtrahendFacet, :); ... 0081 hrepSubtrahend.h(1:(iSubtrahendFacet-1), :)]; 0082 % thisHrep = reduceHrep(thisHrep); 0083 % find chebychev radius 0084 0085 % problem.Aineq = [ones(size(thisHrep.A, 1), 1) thisHrep.A]; 0086 % problem.bineq = thisHrep.h; 0087 % [x, ~, exitFlag] = linprog(problem); 0088 % if (exitFlag == 1) && (x(1) < zerotol) 0089 % % empty polytope, skip 0090 % continue 0091 % elseif exitFlag ~= 1 0092 % % no solution 0093 % error('elk:decomposition:numericalError', ['Chebychev ball ' ... 0094 % 'radius computation failed. This should not happen.']); 0095 % end 0096 problem.A = [ones(size(thisHrep.A, 1), 1) thisHrep.A]; 0097 problem.B = thisHrep.h; 0098 OPT = callCddLib('solve_lp', problem); 0099 if OPT.how == 1 && OPT.xopt(1) < 0.5*min(zerotol, zerotolCddlib) 0100 continue 0101 elseif OPT.how ~= 1 0102 error('elk:decomposition:numericalError', ['Chebychev ball ' ... 0103 'radius computation failed. This should not happen.']); 0104 end 0105 0106 nResult = nResult + 1; 0107 resultHrepList{nResult} = thisHrep; 0108 end 0109 end 0110 0111 % cleanup hrep 0112 resultHrepList((nResult+1):end) = [];