calculate the extreme points of degenerated cones
The extreme point matrix contains the extreme points in its rows. The
projection matrix can be used to project each extreme point back to the
lower-dimensional space that is determined by the equation matrix:
reducedMatrix = (projectionMatrix * extremePointMatrix')'
THIS FUNCTION IS NO USER FUNCTION.0001 function extremePointMatrix = ... 0002 computeExtremePointsOfCone(equationMatrix, inequalityMatrix, zerotol) 0003 % calculate the extreme points of degenerated cones 0004 % 0005 % The extreme point matrix contains the extreme points in its rows. The 0006 % projection matrix can be used to project each extreme point back to the 0007 % lower-dimensional space that is determined by the equation matrix: 0008 % reducedMatrix = (projectionMatrix * extremePointMatrix')' 0009 % 0010 % THIS FUNCTION IS NO USER FUNCTION. 0011 0012 % The elk-library: convex geometry applied to crystallization modeling. 0013 % Copyright (C) 2012 Alexander Reinhold 0014 % 0015 % This program is free software: you can redistribute it and/or modify it 0016 % under the terms of the GNU General Public License as published by the 0017 % Free Software Foundation, either version 3 of the License, or (at your 0018 % option) any later version. 0019 % 0020 % This program is distributed in the hope that it will be useful, but 0021 % WITHOUT ANY WARRANTY; without even the implied warranty of 0022 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0023 % General Public License for more details. 0024 % 0025 % You should have received a copy of the GNU General Public License along 0026 % with this program. If not, see <http://www.gnu.org/licenses/>. 0027 0028 %% Singular value decomposition 0029 % Given are homogenous equations and inequalities: 0030 % eq * h = 0 and ineq * h = 0 0031 % With a singular value decomposition, we obtain: 0032 % eq = U*S*V' with U'*U=E and V'*V=E. 0033 % If we transform the h-vector into a new space by using: 0034 % hDash = V'*h 0035 % the inequality conditions can be solved in this lower-dimensional space: 0036 % ineqDash = ineq*V (from ineq*h = 0 over ineq*V*V'*h = 0) 0037 % where only the last r coordinates are free that have eigenvalues of 0: 0038 % V = V(:, r). 0039 % When the inequalities are solved for this subspace that give summands 0040 % in hDash. These are finally transformed back to h-values by: 0041 % h = V*hDash 0042 0043 projectToReduced = computeProjectionFromConstraint(equationMatrix, zerotol); 0044 reducedDimension = size(projectToReduced, 2); 0045 0046 0047 % transfer inequalities 0048 reducedInequalityMatrix = inequalityMatrix * projectToReduced; 0049 0050 % remove zero rows and trim 0051 reducedInequalityMatrix = trimConstraintMatrix(reducedInequalityMatrix, ... 0052 zerotol, 1); 0053 % remove duplicate rows 0054 reducedInequalityMatrix = unique(reducedInequalityMatrix, 'rows'); 0055 0056 % find extreme points of polytope 0057 extremePointMatrixReduced = callCddLib('extreme', struct(... 0058 'A', reducedInequalityMatrix, ... 0059 'B', zeros(size(reducedInequalityMatrix, 1), 1)... 0060 )); 0061 if ~isempty(extremePointMatrixReduced.V) 0062 error('elk:decomposition:numericalProblem', ... 0063 ['Conversion of cone should only return rays, no extreme points. ', ... 0064 'This should not happen.']); 0065 end 0066 % these are the rays: 0067 extremePointMatrixReduced = extremePointMatrixReduced.R; 0068 0069 % normalize rays to points on plane 0070 planeNormal = mean(extremePointMatrixReduced, 1); 0071 for iSummand = 1:size(extremePointMatrixReduced, 1) 0072 extremePointMatrixReduced(iSummand, :) = ... 0073 extremePointMatrixReduced(iSummand, :) / ... 0074 (planeNormal * extremePointMatrixReduced(iSummand, :)'); 0075 end 0076 0077 extremePointMatrix = (projectToReduced * extremePointMatrixReduced')'; 0078 % projectionMatrix = projectToReduced';