function vrepList = mergeVrepList(vrepList, zerotol)
% merge polytopes to decrease total number
%
% THIS FUNCTION IS NOT PUBLISHED

% Ducmentation draft for publishing:
%
% THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 
%   OF APPLICATIONS (decomposition case study; technical reports)
%
% Syntax: vrepList = mergeVrepList(vrepList)
%
% Reduce the set of simplices in vrepList by subsequent merging. The
%   algorithm that is used here is only implemented to generate some
%   figures in technical reports.
%
% See also: triangulateVrep, identifyPolytopeList, subtractPolytope

% 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
if ~exist('zerotol', 'var')
    zerotol = elkZerotol;
end

% The idea is the following:
%  1) cycle through all polytopes until the last is reached
%  2) merge this polytope with as much polytopes as possible:
%    2a) select a facet (or the next facet)
%    2b) find a neighbouring polytope that shares this facet
%    2c) if the convex hull contains vertices of other polytope, include
%        them in the current list
%    2d) if the volume of the convex hull of the current list equals the 
%        sum of the volumes of the polytopes in the list, merge polytope:
%       * substitute current polytope with merged polytope
%       * delete all merged polytopes in complete list
%    2e) update information and go to 2a
%    2f) if no merge was found, go to 2a anyway (without updating)
%    2g) if no facet is left, loop finishes
%  3) done

% counters and loop control
iPoly = 1;
iFacet = nan;
filterNeighbour = 0;
updateThisPoly = 1;
while iPoly < length(vrepList)
    disp(['mergeVrepList: nPart = ' num2str(length(vrepList))]);
    % catch change of current polytope
    if updateThisPoly
        % polytope was merged, or the next polytope in list was chosen
        vrepMain = vrepList{iPoly};
        hrepMain = convertVrepToHrep(vrepMain);
        updateThisPoly = 0;
        iFacet = 0;
        filterNeighbour = 0;
        continue;
    end
    
    % catch change of facet, when neighbourList is empty
    if sum(filterNeighbour) <= 0
        % no merging happened and no neighbours on facet remain
        iFacet = iFacet + 1;
        % if all facets are searched, chose next polytope
        if iFacet > length(hrepMain.h)
            iPoly = iPoly + 1;
            updateThisPoly = 1;
            continue
        end
        
        % determine new vrepNeighbourList
        % a neighbour must:
        %  * must have n extreme points that can be expressed as a linear
        %    combination of the points of the current facet
        facetVertexMatrix = vrepMain.V(...
            isZero(hrepMain.A(iFacet, :) * vrepMain.V' - hrepMain.h(iFacet), zerotol), ...
            :);
        % singular value decomposition for inner-loop calculations
        [u,s,v] = svd(facetVertexMatrix, 0);
%             uSum = sum(u, 1);
        % loop init
        filterNeighbour = ones(1, length(vrepList)) == 1;
        filterNeighbour(iPoly) = 0;
        for iNeighbour = find(filterNeighbour)
            thisVrepNeighbour = vrepList{iNeighbour};
            
            thisFilterVector = isZero(...
                hrepMain.A(iFacet, :) * thisVrepNeighbour.V' - hrepMain.h(iFacet), ...
                zerotol);
                        
            % quick check, that at least n points of the neighbour must
            %   reside on the supporting hyperplane of the facet
            if sum(thisFilterVector) < size(vrepMain.V, 2)
                filterNeighbour(iNeighbour) = 0;
%                 disp('first drop out')
                continue
            end
%             disp('accepted')
%             hrepMain.A(iFacet, :) * thisVrepNeighbour.V' - hrepMain.h(iFacet)
            
            % for the vertices of the neighbour that are on the facet of
            %   the current polytope, it must hold:
            %     lam * vFacet = v; lam > 0
            %   where lam is [lam1, .., lamN], vFacet are all facets of the
            %   current polytope and v (column vector) is a single vertex
            %   that already resides on the facet.
            % accorindg to singular value decomposition, we might have:
            %     vFacet = U*S*(V')
            %   for which we introduce:
            %     gam = lam*U OR lam = gam*(U')
            %   to obtain:
            %     gam * S = v*V; gam*(U') > 0
            %   or finally:
            %     v*V*(S^-1)*(U') > 1
            nVertexOnFacet = 0;
            % for the vertices that are already on the supporting
            %   hyperplane, count the ones that are on the facet
            for iVertex = find(thisFilterVector)
                thisVertex = thisVrepNeighbour.V(iVertex, :);
                if all(thisVertex*v*(diag(1./diag(s)))*(u') > -zerotol)
                    nVertexOnFacet = nVertexOnFacet + 1;
                end
            end
            
            if nVertexOnFacet >= size(vrepMain.V, 2)
%                 disp('taken')
                filterNeighbour(iNeighbour) = 1;
            else
%                 disp('second drop out')
                filterNeighbour(iNeighbour) = 0;
            end
            
        end
        
%         displayNeighbourFilter = filterNeighbour
        
        % enter next part while reentering the loop
        continue;
    end
    
    % check merging last entry in neighbour list
    if sum(filterNeighbour) > 0 && iFacet <= length(hrepMain.h)
        % skip that one for now:
%         filterNeighbour = 0;
%         continue
        
        iNeighbour = find(filterNeighbour, 1, 'last');
        
        % add single polytope
        vrepMerged.V = [vrepMain.V; vrepList{iNeighbour}.V];
        hrepMerged = convertVrepToHrep(vrepMerged);
        
        % check for other polytopes
        filterPolyAlso = 1:length(vrepList) ~= iNeighbour;
        filterPolyAlso(iPoly) = 0;
        for iPolyAlso = find(filterPolyAlso)
            if any(prod(double(...
                 bsxfun(@minus, ...
                   hrepMerged.A * vrepList{iPolyAlso}.V', ...
                   hrepMerged.h) < -zerotol*10), 1)) || ...
               all(prod(double(...
                 bsxfun(@minus, ...
                   hrepMerged.A * vrepList{iPolyAlso}.V', ...
                   hrepMerged.h) < zerotol*10), 1))
%                 disp('added polytope')
                filterPolyAlso(iPolyAlso) = 1;
                vrepMerged.V = [vrepMerged.V; vrepList{iPolyAlso}.V];
            else
                filterPolyAlso(iPolyAlso) = 0;
            end
        end
        
        % the very final check for the volume
        totalVolume = 0;
        for iPart = [iPoly, iNeighbour, find(filterPolyAlso)]
            totalVolume = totalVolume + computeVolumeVrep(vrepList{iPart});
        end
        
        if abs(totalVolume - computeVolumeVrep(vrepMerged)) < zerotol
            % accept merging
%             listLengthBefore = length(vrepList)
%             mergedPolys = 2 + sum(filterPolyAlso)
            vrepList{iPoly} = reduceVrep(vrepMerged);
            filterDelete = filterPolyAlso;
            filterDelete(iNeighbour) = 1;
            vrepList(filterDelete) = [];
%             listLengthAfter = length(vrepList)
            updateThisPoly = 1;
%             disp('polytopes merged');
        else
            if false
            totalVolume - computeVolumeVrep(vrepMerged)
            figure(66), hold on
            viewPolytope(vrepList{iPoly}, 'color', [1 0 0])
            viewPolytope(vrepList{iNeighbour}, 'color', [0 1 0])
            for iTmp = find(filterPolyAlso)
%                 pause
                viewPolytope(vrepList{iTmp})
            end
%             figure(67)
%             viewPolytope(vrepMerged)
            pause
            close(66)
%             close(67)
            end
            % neighbour did not succeed, delete it and try next one
            filterNeighbour(iNeighbour) = 0;
%             disp('polytopes rejected');
        end
    else
        error('elk:polytope:numericalError', ...
              'Neighbour list is empty, this should not happen at this stage');
    end
    
end
