function [edgeLengthFunctional vrepFunctional incidenceMatrix] = ...
    computeEdgeLengthFunctional(hrep, vrep, zerotol)
% compute function handles for edge lengths and vrep
%
% it is assumed that the edges in hrep exist and remain for inputs to the
%   function handles. edgeLengthFunctional(h) returns a vector of edge
%   lengths for P(H) and vrepFunctional(h) returns the corresponding vrep.
%
% THIS 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/>

if isempty(hrep.A)
    %% empty hrep is treated differently
    edgeLengthFunctional = @(h)([]);
    vrepFunctional = @(h)(struct('V', [0 0]));
    incidenceMatrix = zeros([0 1]);
else
    %% basic information
	% the procedure to identify vertex/edge pairs looks complicated but was
	%   necessary because of tolerance problems. The incidence matrix
	%   should only contain 2 entries in each row AND in each column.
    distanceMatrix = abs(bsxfun(@minus, hrep.A * vrep.V', hrep.h));
    % select first vertex
    tmpMinVec = min(distanceMatrix, [], 1);
    incidenceMatrix = bsxfun(@eq, distanceMatrix, tmpMinVec);
    % this distance will not be selected again
    distanceMatrix(incidenceMatrix) = 2*max(distanceMatrix(:));
    % try again
    tmpMinVec = min(distanceMatrix, [], 1);
    % but avoid selecting an additional vertex/edge when there are already
    %   2 selected
    tmpMinVec(sum(incidenceMatrix, 1) > 1) = 2*max(distanceMatrix(:));
    incidenceMatrix = incidenceMatrix | ...
        bsxfun(@eq, distanceMatrix, tmpMinVec);
    
%     incidenceMatrix = abs(bsxfun(@minus, hrep.A * vrep.V', hrep.h)) < zerotol;
%     % avoid multiple entries >> I can at least correct a single fault
%     if sum(sum(incidenceMatrix, 1) > 2) == 1
%         xIndex = find(sum(incidenceMatrix, 1) == 3);
%         yIndex = find(sum(incidenceMatrix, 2) == 3);
%         incidenceMatrix(yIndex, xIndex) = 0;
%     end
    if sum(sum(incidenceMatrix, 1) <= 2) == 1
        error('elk:boundary:internalError', ['Basically tolerance ' ...
            'problems. This should not happen']);
    end
    
    %% compute vrep functional
    coefficientMatrixX = zeros(size(incidenceMatrix, 2),...
                               size(incidenceMatrix, 1));
    coefficientMatrixY = coefficientMatrixX;
    for iVertex = 1:size(incidenceMatrix, 2)
        % from the incidence matrix, we easily know the adjacent edges:
        edgeIndexVector = find(incidenceMatrix(:, iVertex));
        % now, the vertex [x, y] is defined by:
        %     [a1, a2]' * [x, y]' = [h1, h2]'
        %   which becomes:
        %     [x y] = [h1, h2] * inv([a1, a2])
        %   and is rewritten by:
        %     PI = [p1, p2] = inv([a1, a2])
        %     x = p1' * [h1, h2]'
        %     y = p2' * [h1, h2]'
        PI = inv(hrep.A(edgeIndexVector, :)');
        coefficientMatrixX(iVertex, edgeIndexVector(1)) = PI(1, 1);
        coefficientMatrixX(iVertex, edgeIndexVector(2)) = PI(2, 1);
        coefficientMatrixY(iVertex, edgeIndexVector(1)) = PI(1, 2);
        coefficientMatrixY(iVertex, edgeIndexVector(2)) = PI(2, 2);
    end
    vrepFunctional = @(h)(struct('V', ...
        [coefficientMatrixX*h, coefficientMatrixY*h]));
    
    %% compute edgeLengthFunctional
    coefficientMatrix = zeros(size(hrep.A, 1));
    for iEdge = 1:size(hrep.A, 1)
        % identify bounding edges
        vertexIndexVector = find(incidenceMatrix(iEdge, :));
        thisFilter = ones(1, size(hrep.A, 1)) == 1;
        thisFilter(iEdge) = false;
        thisFilter = thisFilter & (...
            incidenceMatrix(:,vertexIndexVector(1))' | ...
            incidenceMatrix(:,vertexIndexVector(2))');
        edgeIndexVector = find(thisFilter);
        % assume, the whole polytope is moved such that the edge goes
        %   through the current support line:
        %      hM = h - A*(a0*h0)
        % now we can assign a distance of the two vertices from the origin,
        %   when they are projected to the support line (oblique
        %   projection):
        %      Px = a0Perp * inv(ax'*a0Perp) * ax'
        %         = (a0Perp*ax') / (ax'*a0Perp)
        %      d1 = a0Perp' * P1 * a1 * hM1
        %      d2 = a0Perp' * P2 * a2 * hM2
        %   where aPerp is perpendicular to (a). The edge length is then
        %   either (d1-d2) or (d2-d1). But is must be positive for lam=0.
        % This reduces to:
        %      d1 = (1/(a1'*a0Perp)) * (h1 - a1'*a0*h0)
        %      d2 = (1/(a2'*a0Perp)) * (h2 - a2'*a0*h0)
        %   and with an arbitrary sign assumed:
        %     d =  1/(a2'*a0Perp)*h2 ...
        %         -1/(a1'*a0Perp)*h1 ...
        %         +(a1'/(a1'*a0Perp) - a2'/(a2'*a0Perp))*a0*h0
        a0 = hrep.A(iEdge, :)';
        a1 = hrep.A(edgeIndexVector(1), :)';
        a2 = hrep.A(edgeIndexVector(2), :)';
        a0Perp = [a0(2); -1*a0(1)];
        coefficientMatrix(iEdge, iEdge) = (...
            a1'/(a1'*a0Perp) - a2'/(a2'*a0Perp))*a0;
        coefficientMatrix(iEdge, edgeIndexVector(1)) = -1/(a1'*a0Perp);
        coefficientMatrix(iEdge, edgeIndexVector(2)) = 1/(a2'*a0Perp);
        % correct sign of edge length
        if coefficientMatrix(iEdge, :) * hrep.h < 0
            coefficientMatrix(iEdge, :) = -1*coefficientMatrix(iEdge, :);
        end
    end
    
    edgeLengthFunctional = @(h)((coefficientMatrix * h)');
end


