function collectiveStruct = computeCollectiveShapeStruct(hrep, zerotol)
% compute information to quickly operate on lambda-shapes
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 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/>.

%% initial structure
collectiveStruct(1).hrep = hrep;
collectiveStruct(1).lambda = 0;

while ~isempty(collectiveStruct(end).hrep.A)
    %% ensure nonredundant representation
    collectiveStruct(end).hrep = reduceHrep(collectiveStruct(end).hrep);
    
    %% check emptyness
    % compute vrep with large infinity box
    collectiveStruct(end).vrep = ...
        convertHrepToVrep(collectiveStruct(end).hrep, [], 1e5);
    
    % polytope is not sufficiently closed if a vertex exceeds the following
    %   threshold; it is also considered final when its volume does not
    %   contain a pixel anymore
    if any(abs(collectiveStruct(end).vrep.V(:)) > 9e4) || ...
       (computeVolumeVrep(collectiveStruct(end).vrep) < 1e-3)
        % if not closed, there is no hrep and vrep equals the origin
        collectiveStruct(end).hrep = struct('A', [], 'h', []);
        collectiveStruct(end).vrep.V = [0 0];
    end
    
    %% compute vrep and edgeLength functionals
    [collectiveStruct(end).edgeLengthFunctional, ...
     collectiveStruct(end).vrepFunctional, ...
     collectiveStruct(end).incidenceMatrix] = ...
     computeEdgeLengthFunctional(collectiveStruct(end).hrep, ...
                                 collectiveStruct(end).vrep, zerotol);
    
    %% find next lambda
    % now we investigate the edge lengths of P(h-lam) which are calculated
    %   by:
    %     e = M*(h-lam)
    %   while we demand e>=0 and actually search for the largest lambda
    %   that fulfills this inequality:
    %     <maximize lam> with (M*1)*lam<=M*h.
    %   We can directly calculate the candidate lambda values for each edge
    %   and subequently have to take the smallest lambda
    nEdge = size(collectiveStruct(end).hrep.A, 1);
    leftVector = collectiveStruct(end).edgeLengthFunctional(ones(nEdge, 1));
    rightVector = collectiveStruct(end).edgeLengthFunctional(...
        collectiveStruct(end).hrep.h);
    lambdaVector = rightVector ./ leftVector;
    [thisLambda thisRemIndex] = min(lambdaVector);
    
    %% assign next hrep
    if  ~isempty(collectiveStruct(end).hrep.A)
        collectiveStruct(end+1).hrep = collectiveStruct(end).hrep;
        collectiveStruct(end).hrep.A(thisRemIndex, :) = [];
        collectiveStruct(end).hrep.h(thisRemIndex) = [];
        collectiveStruct(end).hrep.h = collectiveStruct(end).hrep.h - thisLambda;
        collectiveStruct(end).lambda = collectiveStruct(end-1).lambda + ...
            thisLambda;
    end
    
end
