computeExtremeRaysFromCone

 compute extreme rays of a cone

 THIS IS NO USER FUNCTION

• convertHrepToVrep convert the H-representation of a polytope into V-representation

• addValidityCone data to identify valid hC-vectors
• projectCompleteDec project a complete proper representation

0001 function rayMatrix = computeExtremeRaysFromCone(A, zerotol)
0002 % compute extreme rays of a cone
0003 %
0004 % THIS IS NO USER FUNCTION
0005 
0006 % The elk-library: convex geometry applied to crystallization modeling.
0007 %   Copyright (C) 2013 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 %% input verification
0023 if rank(A, zerotol) < size(A, 2)
0024     error('elk:decomposition:numericalError', ['Rank of input should be ' ...
0025         'full column rank. This should not happen.']);
0026 end
0027 
0028 %% closing ray
0029 % according to Davis, 1954: "Theory of Positive Linear Dependence", a
0030 %   closing ray can be constructed from a linear basis of the Rn.
0031 %   Therefore, we first need to reduce the Matrix A
0032 tmpA = A;
0033 iRow = 1;
0034 while size(tmpA, 1) > size(tmpA, 2)
0035     thisFilter = (1:size(tmpA,1)) ~= iRow;
0036     if rank(tmpA(thisFilter, :), zerotol) == size(A, 2)
0037         tmpA = tmpA(thisFilter, :);
0038     else
0039         iRow = iRow + 1;
0040     end
0041     if iRow > size(tmpA, 1)
0042         error('elk:decomposition:numericalError', 'This should not happen.');
0043     end
0044 end
0045 % Now, Davis states that a positive basis (and, hence, the missing closing
0046 %   Ray) can be constructed by taking any negative combination of the input
0047 %   vectors:
0048 closeRay = -1*sum(tmpA, 1);
0049 closeRay = closeRay / norm(closeRay, 2);
0050 
0051 %% polytope operations
0052 closedConeHrep = struct('A', [A; closeRay], ...
0053                         'h', [zeros(size(A,1), 1); 1]);
0054 closedConeVrep = convertHrepToVrep(closedConeHrep, zerotol);
0055 
0056 %% ray extraction
0057 rayMatrix = eliminateRows(closedConeVrep.V, ...
0058                           zeros(1, size(closedConeVrep.V, 2)), ...
0059                           zerotol);