function extremePointMatrix = ...
    computeExtremePointsOfCone(equationMatrix, inequalityMatrix, zerotol)
% 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.

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

%% Singular value decomposition
% Given are homogenous equations and inequalities:
%   eq * h = 0 and ineq * h = 0
% With a singular value decomposition, we obtain:
%   eq = U*S*V' with U'*U=E and V'*V=E.
% If we transform the h-vector into a new space by using:
%   hDash = V'*h
% the inequality conditions can be solved in this lower-dimensional space:
%   ineqDash = ineq*V (from ineq*h = 0 over ineq*V*V'*h = 0)
% where only the last r coordinates are free that have eigenvalues of 0:
%   V = V(:, r).
% When the inequalities are solved for this subspace that give summands
%   in hDash. These are finally transformed back to h-values by:
%     h = V*hDash

projectToReduced = computeProjectionFromConstraint(equationMatrix, zerotol);
reducedDimension = size(projectToReduced, 2);


% transfer inequalities
reducedInequalityMatrix = inequalityMatrix * projectToReduced;

% remove zero rows and trim
reducedInequalityMatrix = trimConstraintMatrix(reducedInequalityMatrix, ...
    zerotol, 1);
% remove duplicate rows
reducedInequalityMatrix = unique(reducedInequalityMatrix, 'rows');

% find extreme points of polytope
extremePointMatrixReduced = callCddLib('extreme', struct(...
    'A', reducedInequalityMatrix, ...
    'B', zeros(size(reducedInequalityMatrix, 1), 1)...
    ));
if ~isempty(extremePointMatrixReduced.V)
    error('elk:decomposition:numericalProblem', ...
         ['Conversion of cone should only return rays, no extreme points. ', ...
          'This should not happen.']);
end
% these are the rays:
extremePointMatrixReduced = extremePointMatrixReduced.R;

% normalize rays to points on plane
planeNormal = mean(extremePointMatrixReduced, 1);
for iSummand = 1:size(extremePointMatrixReduced, 1)
    extremePointMatrixReduced(iSummand, :) = ...
        extremePointMatrixReduced(iSummand, :) / ...
        (planeNormal * extremePointMatrixReduced(iSummand, :)');
end

extremePointMatrix = (projectToReduced * extremePointMatrixReduced')';
% projectionMatrix = projectToReduced';
