function resultHrepList = subtractPolytopeHrep(...
    hrepMinuendList, hrepSubtrahend, zerotol)
% calculates the set difference for a list of hrep's versus a single hrep
%
% THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 
%   OF APPLICATIONS (decomposition case study; technical reports)
%

%! ToDo: redo documentation
% Syntax: vrepList = subtractPolytope(polyMinuend, polySubtrahend)
%
% The input minuend can be a V-representation, H-representation or a cell
%   array of these representations. While in the latter case it is expected
%   that this list contains only simplices. The subtrahend should be a
%   H-representation of V-representation only.
%
% subtractPolytope(.., zerotol) lets you specify the zerotol that is used 
%   to identify values that are close to zero by zero. Default: elkZerotol.
%
% subtractPolytope(.., zerotol, isCone)  is a specific logical parameter that
%   identifies the input as a convex cone that contains the origin as one
%   point. All other points must be at a hyperplane. This parameter is
%   passed to triangulateVrep.
%
% See also: triangulateVrep, identifyPolytope

% 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
% zerotol
if ~exist('zerotol', 'var')
    zerotol = elkZerotol;
end
% sort the constraints of the subtrahend
[hrepSubtrahend.A idxVector] = sortrows(hrepSubtrahend.A);
hrepSubtrahend.h = hrepSubtrahend.h(idxVector);

% The principle is that every part of the minuend is treated seperately and
%   then we go through the disjoint convex compartements for the complement
%   of the subtrahend.
% some numbers:
nSubtrahendFacet = length(hrepSubtrahend.h);
nMinuend = length(hrepMinuendList);
nResult = 0;
% common options for Chebychev ball computation:
% problem.solver = 'linprog';
% problem.options = optimset('Display', 'off', 'largeScale', 'off', ...
%                            'TolFun', 0.01*zerotol, 'simplex', 'off');
% problem.f = [-1 zeros(1, size(hrepSubtrahend.A, 2))];
% problem.lb = -1;
% problem.ub = 1;
% problem.x0 = zeros(1 + size(hrepSubtrahend.A, 2), 1);
problem.obj = [-1 zeros(1, size(hrepSubtrahend.A, 2))];
problem.lin = [];


resultHrepList = cell(1, nSubtrahendFacet*nMinuend);

for iSubtrahendFacet = 1:length(hrepSubtrahend.h)
    for iMinuend = 1:length(hrepMinuendList)
        % construct hrep
        thisHrep.A = [hrepMinuendList{iMinuend}.A; ...
                   -1*hrepSubtrahend.A(iSubtrahendFacet, :); ...
                      hrepSubtrahend.A(1:(iSubtrahendFacet-1), :)];
        thisHrep.h = [hrepMinuendList{iMinuend}.h; ...
                   -1*hrepSubtrahend.h(iSubtrahendFacet, :); ...
                      hrepSubtrahend.h(1:(iSubtrahendFacet-1), :)];
%         thisHrep = reduceHrep(thisHrep);
        % find chebychev radius
        
%         problem.Aineq = [ones(size(thisHrep.A, 1), 1) thisHrep.A];
%         problem.bineq = thisHrep.h;
%         [x, ~, exitFlag] = linprog(problem);
%         if (exitFlag == 1) && (x(1) < zerotol)
%             % empty polytope, skip
%             continue
%         elseif exitFlag ~= 1
%             % no solution
%             error('elk:decomposition:numericalError', ['Chebychev ball ' ...
%                 'radius computation failed. This should not happen.']);
%         end
        problem.A = [ones(size(thisHrep.A, 1), 1) thisHrep.A];
        problem.B = thisHrep.h;
        OPT =  callCddLib('solve_lp', problem);
        if OPT.how == 1 && OPT.xopt(1) < 0.5*min(zerotol, zerotolCddlib)
            continue
        elseif OPT.how ~= 1
            error('elk:decomposition:numericalError', ['Chebychev ball ' ...
                'radius computation failed. This should not happen.']);
        end
        
        nResult = nResult + 1;
        resultHrepList{nResult} = thisHrep;
    end
end

% cleanup hrep
resultHrepList((nResult+1):end) = [];