function outVar = computeEdgeSetVrep(vrep, varargin)
% identify the edges of a V-representation
%
% Syntax: vertexIndexMatrix = computeEdgeSetVrep(vrep)
%         outStruct = computeEdgeSetVrep(vrep, 'output', 'full')
%
% The matrix vertexIndexMatrix gives in each row the indices into vrep.V 
%   for the two vertices that define an edge. This algorithm is valid for 
%   arbitrary dimensions.
%
% computeEdgeSetVrep(.., 'output', out) lets you specify the type of
%   output that is provided. With out='normal', only the index matrix to
%   the vertices in vrep.V is returned. For out='full' a structure is
%   returned with the fields: 
%     vrep - the V-representation of the polytope
%     hrep - the H-representation of the polytope
%     nEdge - the total number of edges
%     vertexIndexMatrix - indices to the vertices in vrep.V for the edges
%     facetIndexMatrix - indices to the (n-1) facets that define the edge
%     incidenceMatrix - a matrix with a row for each extreme point and a
%       column for each facet that is assigned with 1 when the vertex
%       belongs to the facet and 0 otherwise.
%
% computeEdgeSetVrep(.., 'hrep', hrep) hands over the precalculated
%   H-representation of the polytope so that it is not calculated twice if
%   it is already known.
%
% computeEdgeSetVrep(.., 'zerotol', zerotol) allows to set the tolerance 
%   for underlying algorithms (determining flat polytopes; default in
%   elkZerotol).
%
% See also: obtainVrep, viewPolytope

% 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/>.

%% Input handling
[hrep, zerotol, output] = mapOption(varargin, ...
    'hrep', [], ...
    'zerotol', elkZerotol, ...
    'output', 'normal');

if isempty(hrep)
    hrep = convertVrepToHrep(vrep);
end

nVertex = size(vrep.V, 1);
nFacet = size(hrep.A, 1);
nDim = size(hrep.A, 2);

%% Incidence matrix
% generate an incidence matrix - rows for vertices and columns for facets -
%   a 1 is assigned when the vertex belongs to the corresponding facet - 0
%   is assigned otherwise.
incidenceMatrix = zeros(nVertex, nFacet, 'uint8');
for iVertex = 1:size(vrep.V, 1)
    for iFacet = 1:size(hrep.A, 1)
        incidenceMatrix(iVertex, iFacet) = abs(...
            hrep.A(iFacet, :)*vrep.V(iVertex,:)' - hrep.h(iFacet)) < zerotol;
    end
end

%% Evaluation of incidence
% search all vertex pairs for which the logical and-operation of the 
%   corresponding rows leaves exactly two ones. Less than two ones implies
%   that these two vertices do not have two facets in common (necessary
%   condition to form an edge). If more that two facets are adjacent for
%   both vertices, then it does no matter which two to take.
vertexIndexMatrix = [];
facetIndexMatrix = [];
for iVertexOne = 1:nVertex
    for iVertexTwo = (iVertexOne+1):nVertex
        thisFacetFilter = incidenceMatrix(iVertexOne, :) & ...
                          incidenceMatrix(iVertexTwo, :);
        nCommonFacet = sum(thisFacetFilter);
        if nCommonFacet >= (nDim-1)
            vertexIndexMatrix(end+1, :) = [iVertexOne, iVertexTwo];
            facetIndexMatrix(end+1, :) = find(thisFacetFilter, 2);
        end
%         if nCommonFacet > (nDim-1)
%             error('elk:polytope:numericalError', ...
%                 'You have two vertices that belong to the same (n-1) facets. This should not happen.');
%         end
    end
end

%% Output handling
switch lower(output)
    case 'normal'
        outVar = vertexIndexMatrix;
    case 'full'
        outVar.hrep = hrep;
        outVar.vrep = vrep;
        outVar.vertexIndexMatrix = vertexIndexMatrix;
        outVar.facetIndexMatrix = facetIndexMatrix;
        outVar.incidenceMatrix = incidenceMatrix;
        outVar.nEdge = size(vertexIndexMatrix, 1);
    otherwise
        error('elk:polytope:inputError', ...
            'Please use ''normal'' or ''full'' for the output parameter.');
end