scale polytope to reasonable size and position Syntax: [hrep, scale, x0] = scaleHrep(hrep) or [hrep, scale, x0] = scaleHrep(hrep, zerotol) The Chebychev ball is the largest ball that fits into a given polytope. The input polytope is scaled so that the center of the Chebychev ball is the origin of the coordinate system and its radius is 1. The polytope remains unscaled whenever the radius of the Chebychev ball is below zerotol or above 1/zerotol. For abs(r)<=zerotol the output is assigned to scale=0 (flat polytope). For r=inf the output is assigned to scale=inf (unbounded polytope) and for r<(-1*zerotol) the output is assigned to scale=-inf (empty polytope). Default: elkZerotol. x0 is empty in unscaled cases and hrep is left unchanged. The polytope is scaled back to original by: originalHrep = stretchHrep(scaledHrep, scale); originalHrep = moveHrep(originalHrep, x0); See also: convertHrepToVrep, scaleHrep, reduceVrepDimension, reduceVrep
0001 function [hrep, scale, x0] = scaleHrep(hrep, zerotol) 0002 % scale polytope to reasonable size and position 0003 % 0004 % Syntax: [hrep, scale, x0] = scaleHrep(hrep) 0005 % or [hrep, scale, x0] = scaleHrep(hrep, zerotol) 0006 % 0007 % The Chebychev ball is the largest ball that fits into a given polytope. 0008 % The input polytope is scaled so that the center of the Chebychev ball 0009 % is the origin of the coordinate system and its radius is 1. 0010 % The polytope remains unscaled whenever the radius of the Chebychev ball 0011 % is below zerotol or above 1/zerotol. For abs(r)<=zerotol the output 0012 % is assigned to scale=0 (flat polytope). For r=inf the output is 0013 % assigned to scale=inf (unbounded polytope) and for r<(-1*zerotol) the 0014 % output is assigned to scale=-inf (empty polytope). Default: elkZerotol. 0015 % x0 is empty in unscaled cases and hrep is left unchanged. 0016 % 0017 % The polytope is scaled back to original by: 0018 % originalHrep = stretchHrep(scaledHrep, scale); 0019 % originalHrep = moveHrep(originalHrep, x0); 0020 % 0021 % See also: convertHrepToVrep, scaleHrep, reduceVrepDimension, reduceVrep 0022 0023 % The elk-library: convex geometry applied to crystallization modeling. 0024 % Copyright (C) 2012 Alexander Reinhold 0025 % 0026 % This program is free software: you can redistribute it and/or modify it 0027 % under the terms of the GNU General Public License as published by the 0028 % Free Software Foundation, either version 3 of the License, or (at your 0029 % option) any later version. 0030 % 0031 % This program is distributed in the hope that it will be useful, but 0032 % WITHOUT ANY WARRANTY; without even the implied warranty of 0033 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0034 % General Public License for more details. 0035 % 0036 % You should have received a copy of the GNU General Public License along 0037 % with this program. If not, see <http://www.gnu.org/licenses/>. 0038 0039 %! DEV: This function could in principle use a non-uniform scaling, given 0040 % the procedures in reduceHrepDimension. Analogue to this procedure, 0041 % we could identify a limiting set of constraints (and therefore a 0042 % subspace for scaling) at a time. 0043 % However, for most cases this simple size-scaling should be sufficient. 0044 %/ 0045 0046 %% Input 0047 if ~exist('zerotol', 'var') 0048 zerotol = elkZerotol; 0049 end 0050 nDim = size(hrep.A, 2); 0051 nFacet = size(hrep.A, 1); 0052 x0 = []; 0053 0054 %% additional scaling 0055 % the success of the linear program is dependent on scaling and centering. 0056 % While polytopes are, for this task, somehow necessarily badly scaled 0057 % (degenerate polytopes that are reduced in dimension), there position 0058 % might be adopted. 0059 extraX0 = 0.5*hrep.A'*(hrep.h./sum(hrep.A.^2, 2)); 0060 hrepCentered = moveHrep(hrep, -1*extraX0); 0061 0062 %% find Chebychev ball 0063 % f = [1 0 0 0] 0064 % b = [0 h] 0065 % A = [1 A] 0066 problem.solver = 'linprog'; 0067 problem.options = optimset('Display', 'off', 'largeScale', 'off', ... 0068 'TolFun', 0.01*zerotol, 'simplex', 'on'); 0069 problem.f = [-1 zeros(1, nDim)]; 0070 problem.Aineq = [ones(nFacet, 1) hrep.A]; 0071 problem.bineq = hrepCentered.h; 0072 % problem.lb = [-1 -1*sqrt(1/zerotol)*ones(1, nDim)]; 0073 % problem.ub = [1 sqrt(1/zerotol)*ones(1, nDim)]; 0074 % have to add -inf to match prior implementation but avoid errors in MATLAB 0075 % R2012b 0076 problem.lb = [-1*sqrt(1/zerotol) -inf(1, nDim)]; 0077 problem.up = [sqrt(1/zerotol) inf(1,nDim)]; 0078 [x, ~, exitFlag, ~, lambda] = linprog(problem); 0079 if lambda.upper(1) ~= 0 0080 % this also means unbounded 0081 exitFlag = -3; 0082 end 0083 % This optimization problem should always find a feasible point since the 0084 % Chebychev radius is allowed to be negative to ensure a feasible point. 0085 % Note also that the simplex algorithm should always find this feasible 0086 % point automatically. 0087 0088 switch exitFlag 0089 case 1 0090 % solution found we have to distinguish the different radii cases 0091 % while an internal point is provided: 0092 if x(1) < -zerotol 0093 % empty polytope 0094 hrep = []; 0095 scale = -inf; 0096 return 0097 elseif x(1) > 0.9/zerotol 0098 % we consider these cases as unbounded even if the probloem in 0099 % itself is bounded. But we still should not scale this 0100 % result. 0101 % Note that it happens that it exits with flag 1 even though 0102 % the problem is quite unbounded. In this case it has hit the 0103 % upper bounds. 0104 scale = inf; 0105 return 0106 elseif x(1) > zerotol 0107 % properly scalable case 0108 scale = x(1); 0109 x0 = x(2:end) + extraX0; 0110 else 0111 % this leaves the case of a flat polytope (abs(x) <= zerotol) 0112 % which should not be scaled 0113 scale = 0; 0114 return 0115 end 0116 case 0 0117 % no solution / maximum iterations 0118 error('elk:polytope:numericalError', ['Scaling failed by ' ... 0119 'reaching maximum iterations of linprog. This should not ' ... 0120 'happen.']); 0121 case -2 0122 % no feasible point / polytope empty - see above comments for the 0123 % formulated optimization problem 0124 error('elk:polytope:numericalError', ['Scaling failed by ' ... 0125 'not finding any feasible solution to the optimization ' ... 0126 'problem. This should not happen.']); 0127 case -3 0128 % polytope unbounded 0129 scale = inf; 0130 return 0131 otherwise 0132 % other errors 0133 error('elk:polytope:numericalError', ... 0134 'Scaling failed for an unknown reason. This should not happen.'); 0135 end 0136 0137 %% Scale polytope 0138 hrep = moveHrep(hrep, -1*x0); 0139 hrep = stretchHrep(hrep, 1/scale, zeros(1, nDim)); 0140 hrep = normalizeHrep(hrep);