project lower dimensional (flat) polytope into embedding space
Syntax: hrep = reduceVrepDimension(hrep, zerotol)
or [hrep, rot, x0] = reduceVrepDimension(hrep, zerotol)
The input polytope is reduced in its dimension by first finding the
linear constraints that define the embedding linear subspace. In a
second step the polytope is rotated such that this subspace spans the
last coordinate axes. Note that herewith, the position of the polytope
in its original space is lost so that this function returns a rotatoin
matrix and an offset. Note that an empty polytope returns:
hrep = [], rot = eye(nDim), x0 = []
while a zero-dimensional polytope returns:
hrep.A = [], hrep.h = []
and a proper offset vector x0.
Use liftHrepDimension to restore the original polytope. Though, according
to the allowed tolerances, the result might have an position offset in
the range of zerotol.
See also: liftHrepDimension, reduceVrepDimension, obtainPolytope,
rotatePolytope, movePolytope0001 function [hrep rot x0] = reduceHrepDimension(hrep, zerotol) 0002 % project lower dimensional (flat) polytope into embedding space 0003 % 0004 % Syntax: hrep = reduceVrepDimension(hrep, zerotol) 0005 % or [hrep, rot, x0] = reduceVrepDimension(hrep, zerotol) 0006 % 0007 % The input polytope is reduced in its dimension by first finding the 0008 % linear constraints that define the embedding linear subspace. In a 0009 % second step the polytope is rotated such that this subspace spans the 0010 % last coordinate axes. Note that herewith, the position of the polytope 0011 % in its original space is lost so that this function returns a rotatoin 0012 % matrix and an offset. Note that an empty polytope returns: 0013 % hrep = [], rot = eye(nDim), x0 = [] 0014 % while a zero-dimensional polytope returns: 0015 % hrep.A = [], hrep.h = [] 0016 % and a proper offset vector x0. 0017 % 0018 % Use liftHrepDimension to restore the original polytope. Though, according 0019 % to the allowed tolerances, the result might have an position offset in 0020 % the range of zerotol. 0021 % 0022 % See also: liftHrepDimension, reduceVrepDimension, obtainPolytope, 0023 % rotatePolytope, movePolytope 0024 0025 % If the polytope vrep has a lower dimension r<n than the space it is 0026 % embedded in, this fucntion moves and rotates the polytope such that the 0027 % values for the last (n-r) coordinate axes are zero for all points in 0028 % vrep. 0029 % 0030 % The zerotol determines how flat a polytope must be to be recognized. 0031 % Points can be approximately a distance of zerotol away from the flat 0032 % polytope. Default: see elkZerotol. 0033 % 0034 % With the second syntax, the rotation matrix rot and the offset x0 is 0035 % returned so that the original polytope can be obtained by: 0036 % removedCoordinates = zeros(size(vrep, 1), size(V, 1) - size(vrep, 2)); 0037 % vrepOriginal = [vrep, removedCoordinates]*rot'; 0038 % vrepOriginal = moveVrep(vrepOriginal, x0); 0039 % 0040 % See also: reduceVrep, scaleVrep 0041 0042 % The elk-library: convex geometry applied to crystallization modeling. 0043 % Copyright (C) 2012 Alexander Reinhold 0044 % 0045 % This program is free software: you can redistribute it and/or modify it 0046 % under the terms of the GNU General Public License as published by the 0047 % Free Software Foundation, either version 3 of the License, or (at your 0048 % option) any later version. 0049 % 0050 % This program is distributed in the hope that it will be useful, but 0051 % WITHOUT ANY WARRANTY; without even the implied warranty of 0052 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0053 % General Public License for more details. 0054 % 0055 % You should have received a copy of the GNU General Public License along 0056 % with this program. If not, see <http://www.gnu.org/licenses/>. 0057 0058 %% input 0059 if ~exist('zerotol', 'var') 0060 zerotol = elkZerotol; 0061 end 0062 % reject empty polytopes 0063 if isempty(hrep) 0064 rot = []; 0065 x0 = []; 0066 return 0067 end 0068 % reject unbounded or zero-dim polytopes 0069 if isempty(hrep.A) 0070 rot = []; 0071 x0 = []; 0072 return 0073 end 0074 0075 %% algorithm 0076 % The algorithm idea is based on some lecture notes by Komei Fukuda, Zürich 0077 % 2011 obtained, here: 0078 % http://www.ifor.math.ethz.ch/teaching/lectures/poly_comp_ss11/lecture8 0079 % However, the algorithm is altered quite a bit, mainly to allow for some 0080 % tolerance that is required in the scope of this library. 0081 % First, we use the lagrangian multipliers that are also returned by the 0082 % constructed optimization problems. These lagrangian multipliers 0083 % indicate which constraints are limiting the result. How to obtain 0084 % these variables was not clear in the given source. 0085 % Secondly, we do not use equality constraints. Instead, we add the 0086 % identified equality constraints as inequalities. This is necessary if 0087 % you (for example) have two opposing facets that are some distance away 0088 % - even though this distance is smaller than zerotol. The algorithm 0089 % would then select ONE of these constraints while the two constraints 0090 % would not exactly give the same subdomain. 0091 % 0092 % To make the desciption general for all steps, assume that a set of 0093 % constraints j from J are already found that limit the dimension by: 0094 % a_j * x = h_j 0095 % while the remainder of constraints i still gives: 0096 % a_i * x <= h_i 0097 % The optimization problem shall maximize a radius or general (possible) 0098 % decrease for all h_i values r so that the constraints: 0099 % a_i * x + r <= h_i equaling a_i * x <= h_i - r 0100 % is fulfilled together with the equality constraints. 0101 % In the initial run where the set of equations is emtpy, r can be negative 0102 % which implies an empty polytope. In subsequent runs, this is not 0103 % possible. 0104 % If r is positive (here, larger than zerotol) than the polytope is 0105 % fulldimensional when it is embedded in the subspace indicated by the 0106 % equality constraints. 0107 % If r is zero (here, abs(r) < zerotol) then the polytope is degenerated 0108 % and the set of equations must be extended. Since only a minimal number 0109 % of inequalities will be marked limiting by the Lagrangian multipliers, 0110 % the full set should be identified by finding the linearly dependent 0111 % rows from [A h]. 0112 % The optimization problem might be unbounded, indicating a cone. In this 0113 % case, r is treated simply as positive. 0114 % The optimization problem cannot be infeasible since some negative value 0115 % for r always exists to fulfill Ax + r <= h. 0116 0117 %% additional scaling 0118 % the success of the linear program is dependent on scaling and centering. 0119 % While polytopes are, for this task, somehow necessarily badly scaled 0120 % (degenerate polytopes that are reduced in dimension), there position 0121 % might be adopted. 0122 extraX0 = 0.5*hrep.A'*(hrep.h./sum(hrep.A.^2, 2)); 0123 hrepCentered = moveHrep(hrep, -1*extraX0); 0124 % (hrep.A*[0 0 0]' - hrep.h)' 0125 % (hrepCentered.A*[0 0 0]' - hrepCentered.h)' 0126 0127 %% initialize loop 0128 nDim = size(hrep.A, 2); 0129 nFacet = size(hrep.A, 1); 0130 inequalityFilter = true(1, nFacet); 0131 isFirstRun = true; 0132 keepRunning = true; 0133 dimension = []; 0134 0135 problem.solver = 'linprog'; 0136 problem.options = optimset('Display', 'off', 'diagnostics', 'off', ... 0137 'largeScale', 'off', ... 0138 'TolFun', 0.01*zerotol, 'simplex', 'on'); 0139 % simplex 'on' is choosen to automatically find a feasible interior point 0140 % simplex 'on' produces problems for decomposition validation of generic 0141 % lactose 7 0142 problem.f = [-1; zeros(nDim,1)]; 0143 % problem.lb = [-1 -1*sqrt(1/zerotol)*ones(1, nDim)]; 0144 % problem.up = [1 sqrt(1/zerotol)*ones(1, nDim)]; 0145 problem.lb = [-1*sqrt(1/zerotol) -inf(1, nDim)]; 0146 problem.up = [sqrt(1/zerotol) inf(1, nDim)]; 0147 % the euqaion matrix reflects entries of the type [hrep.A hrep.h]: 0148 equationMatrix = zeros(0, nDim+1); 0149 nCycle = 0; 0150 0151 while keepRunning 0152 % count the cycles: we need to know the first clycle, at least 0153 nCycle = nCycle + 1; 0154 0155 %% solve optimization problem 0156 % inequalityFilter 0157 % figure,viewPolytope(struct('A', hrep.A(inequalityFilter, :), 'h', hrep.h(inequalityFilter))); 0158 % the ineqaulity constraints contain the original inequality 0159 % constrainty, but additionally the equations are also formulated as 0160 % (tight) inequality constraints to model a tolerance on the original 0161 % equality constraints. This is required since from two limiting and 0162 % opposing facets, only one is usually selected to reflect the 0163 % constraint while x-vectors living on this constraint might (for 0164 % numerical reasons) not intersect anymore with the original 0165 % polytope. 0166 nEq = size(equationMatrix, 1); 0167 if nCycle == 1 0168 problem.Aineq = [ones(nFacet, 1), hrep.A]; 0169 problem.bineq = [hrepCentered.h]; 0170 else 0171 problem.Aineq = [ones(sum(inequalityFilter), 1), ... 0172 hrep.A(inequalityFilter, :); ... 0173 zeros(nEq, 1) equationMatrix(:, 1:(end-1)); ... 0174 zeros(nEq, 1) -1*equationMatrix(:, 1:(end-1))]; 0175 problem.bineq = [hrep.h(inequalityFilter); ... 0176 equationMatrix(:, end) + zerotol; ... 0177 -1*equationMatrix(:, end) + zerotol]; 0178 end 0179 % max(problem.Aineq * [-0.1 0 0 0]' - problem.bineq) 0180 0181 % inequalityFilter 0182 if any(inequalityFilter) 0183 [x, ~, exitFlag, ~, lambda] = linprog(problem); 0184 % inequalityFilter 0185 % x' 0186 % exitFlag 0187 % lambda.ineqlin' 0188 % (problem.Aineq*x - problem.bineq)' 0189 lambdaIneq = lambda.ineqlin(1:(end-2*nEq)); 0190 else 0191 % there are no inequalities left - this happens when the polytope 0192 % is zero-dimensional. The safe result here is to state that the 0193 % result of the optimization problem is unbounded. 0194 x = []; 0195 exitFlag = -3; 0196 end 0197 % only the first run shall use the simplex approach, all subsequent 0198 % runs shall use active-set so that we can provide an initial 0199 % feasible point. 0200 % Note that the simplex algorithm might not find a feasible point 0201 % anymore when the feasible region is too narrow. 0202 % x is only empty for an initially infeasible problem 0203 if nCycle == 1 && ~isempty(x) 0204 problem.options.Simplex = 'off'; 0205 problem.x0 = [0; x(2:end)+extraX0]; 0206 x0 = x(2:end) + extraX0; 0207 elseif nCycle == 1 0208 x0 = []; 0209 end 0210 0211 % handle the exitFlag 0212 switch exitFlag 0213 case 1 0214 % solution found 0215 if x(1) < -zerotol 0216 % empty polytope 0217 reductionFlag = -1; 0218 keepRunning = false; 0219 elseif x(1) > zerotol 0220 % (now) full-dimensional polytope 0221 % disp('full') 0222 reductionFlag = false; 0223 keepRunning = false; 0224 else 0225 % not full-dimensional, more processing required 0226 % disp('reduce') 0227 reductionFlag = 1; 0228 end 0229 case 0 0230 % no solution / maximum iterations 0231 error('elk:polytope:numericalError', ['Dimension reduction ' ... 0232 'failed by reaching maximum iterations for linprog. '... 0233 'This should not happen.']); 0234 case -2 0235 % no feasible point this should not happen since x(1) can be used 0236 % to enforce the existence of a solution 0237 error('elk:polytope:numericalError', ['Dimension reduction ' ... 0238 'failed by finding no feasible solution. '... 0239 'This should not happen.']); 0240 case -3 0241 % problem is unbounded - well in this case, the polytope is 0242 % obviously full-dimensional even though we cannot return the 0243 % inner point 0244 reductionFlag = 0; 0245 keepRunning = 0; 0246 otherwise 0247 % other errors 0248 error('elk:polytope:numericalError', ['Dimension reduction ' ... 0249 'failed for an unknown reason. '... 0250 'This should not happen.']); 0251 end 0252 0253 if reductionFlag 0254 % disp('reducing') 0255 % lambda.ineqlin 0256 % index values of current inequality constraints 0257 inequalityIndexVector = find(inequalityFilter); 0258 % determine the rows in hrep.A that become equations 0259 % inequalityIndexVector 0260 % lambdaIneq 0261 newEqualityIndexVector = inequalityIndexVector(lambdaIneq > zerotol); 0262 % remove the same rows from the inequalities 0263 inequalityIndexVector(lambdaIneq > zerotol) = []; 0264 inequalityFilter(newEqualityIndexVector) = false; 0265 0266 % get (reduced) row echolon form of the new equation rows 0267 [R, jb] = rref([equationMatrix; ... 0268 hrep.A(newEqualityIndexVector, :) ... 0269 hrep.h(newEqualityIndexVector)], 10*zerotol); 0270 % remove the rows that contain zero-rows to obtain new equation 0271 % matrix 0272 equationMatrix = R(1:length(jb), :); 0273 0274 % find linearly dependent rows in the remaining [hrep.A hrep.h] 0275 for iRow = inequalityIndexVector 0276 thisRow = [hrep.A(iRow, :) hrep.h(iRow)]; 0277 thisRow = thisRow - thisRow(jb)*equationMatrix; 0278 if all(abs(thisRow) < 10*zerotol) 0279 inequalityFilter(iRow) = false; 0280 end 0281 end 0282 end 0283 end 0284 0285 %% compute output 0286 % this is based on the last reductionFlag 0287 switch reductionFlag 0288 case -1 0289 % polytope is empty 0290 hrep = []; 0291 rot = eye(nDim); 0292 case 0 0293 % disp('properly full in some dim') 0294 % polytope is properly full dimensional in some subspace 0295 if all(inequalityFilter) 0296 % polytope is full-dimensional as provided 0297 % hrep remains 0298 rot = eye(nDim); 0299 else 0300 % disp('full in subspace, obtaining rot') 0301 % polytope is full-dimensional in some subspace for which the 0302 % rotation matrix is obtained by singular value 0303 % decomposition: 0304 [u,s,v] = svd(hrep.A(~inequalityFilter, :)); 0305 dimensionMinus = sum(diag(s) > zerotol); 0306 rot = v'; 0307 % we flip the dimensions 0308 rot = rot(end:-1:1, :); 0309 % take rotated versions from A 0310 hrep.A = hrep.A(inequalityFilter, :)*rot'; 0311 % cut off dimensions 0312 hrep.A = hrep.A(:, 1:(end-dimensionMinus)); 0313 hrep.h = hrep.h(inequalityFilter); 0314 % remove zero rows 0315 eliminationFilter = all(abs(hrep.A) < zerotol, 2); 0316 hrep.A(eliminationFilter, :) = []; 0317 hrep.h(eliminationFilter) = []; 0318 end 0319 otherwise 0320 end