function [resultVrepList resultVolumeVector resultIndexVector] = subtractPolytope(...
    polyMinuend, polySubtrahend, varargin)
% calculates the set difference for a triangulated V-representation and any
%   H-representation
%
% THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 
%   OF APPLICATIONS (decomposition case study; technical reports)
%
% 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, isCone, minuendVolumeVector, scalingVolume, infoMinuend] = ...
    mapOption(varargin, ...
        'zerotol', elkZerotol, ...
        'isCone', 0, ...
        'volumeVector', [], ...
        'scalingVolume', 1, ...
        'infoMinuend', []);

% minuend (ensure triangulated vrep)
if isempty(infoMinuend)
    infoMinuend = identifyPolytope(polyMinuend);
end
if strcmpi(infoMinuend.rep, 'vrep')
    [minuendVrepList minuendVolumeVector] = triangulateVrep(...
        polyMinuend, zerotol, isCone);
    infoMinuend = identifyPolytope(minuendVrepList);
elseif strcmpi(infoMinuend.rep, 'hrep')
    polyMinuend = convertHrepToVrep(polyMinuend);
    [minuendVrepList minuendVolumeVector] = triangulateVrep(...
        polyMinuend, zerotol, isCone);
    infoMinuend = identifyPolytope(minuendVrepList);
elseif strcmpi(infoMinuend.rep, 'vrepList')
    if ~all(infoMinuend.m == infoMinuend.m(1))
        error('elk:polytope:wrongInput', ...
            'The given list of polytopes in V-representation is not a triangulation');
    end
    minuendVrepList = polyMinuend;
elseif strcmpi(infoMinuend.rep, 'hrepList')
    if ~all(infoMinuend.m == infoMinuend.m(1))
        error('elk:polytope:wrongInput', ...
            'The given list of polytopes in V-representation is not a triangulation');
    end
    for iSimplex = 1:infoMinuend.nPoly
        polyMinuend{iSimplex} = convertHrepToVrep(polyMinuend{iSimplex}, zerotol);
    end
    minuendVrepList = polyMinuend;
    infoMinuend = identifyPolytope(minuendVrepList);
else
    error('elk:polytope:wrongInput', ...
        'The input minuend should be vrep, hrep, vrepTri or hrepTri');
end

% subtrahend
infoSubtrahend = identifyPolytope(polySubtrahend);
if strcmpi(infoSubtrahend.rep, 'vrep')
    if ~isCone
        hrepSubtrahend = convertVrepToHrep(polySubtrahend, zerotol);
    else
        hrepSubtrahend = convertVrepToHrep(struct(...
            'V', [polySubtrahend.V; polySubtrahend.V(1,:)*0]), zerotol);
        closingFacetFilter = abs(hrepSubtrahend.h) > zerotol;
        hrepSubtrahend.A(closingFacetFilter, :) = [];
        hrepSubtrahend.h(closingFacetFilter) = [];
    end
    infoSubtrahend = identifyPolytope(hrepSubtrahend);
elseif strcmpi(infoSubtrahend.rep, 'hrep')
    hrepSubtrahend = polySubtrahend;
else
    error('elk:polytope:wrongInput', ...
        'The input subtrahend should be vrep or hrep');
end

%% Initial variables
nDim = infoMinuend.dim;
nMinuend = infoMinuend.nPoly;

resultVrepList = cell(0);
resultIndexVector = [];
resultVolumeVector = [];
remainingVrepList = cell(0);
remainingVolumeVector = [];
remainingIndexVector = [];


%% Algorithm overview

% This simple form of the set difference is based on the following term 
%   that represents the set difference of vrepMinuend by hrepSubtrahend:
%     Union over all {x | x in vrepMinuend * x <= 0 AND 
%                         hrepSubtrahend(i, :) * x > = 0}.
%   The inner terms of the union are all convex, contain only points of the
%   minuend but no point of the subtrahend (except boundaries). However,
%   they might have full-dimensional intersections. This problem can be
%   resolved.
% The outer algorithm works as follows:
%   1) From the simplices in in vrepMinuend and the last facet
%      hrepSubtrahend(end,:), calculate a list of simplices that are 
%      definitely contained in the difference and a remaining list. The
%      simplices in vrepResultList and vrepRemainingList have no full-
%      dimensional intersecntions and their union equals vrepMinuend.
%   2) Delete the last facet from hrepSubtrahend.
%   3) If hrepSubtrahend is empty, the algorithm terminates with
%      vrepResultList.
%   4) resultList = {resultList thisFunction(vrepRemainingList,
%                    hrepSubtrahend)}
%   4) Back to step 1 with vrepMinuendList=remainingList and
%      hrepSubtrahend.

% For step one, each simplex in vrepMinuendList is handled separately. For
%   each vertex of such a simplex, there are three possibilities according
%   to the facet in hrepSubtrahend(end, :):
%     ( 0)  - The vertex is on the facet
%     (-1)  - The vertex is on the opposite site of the facet and, hence,
%             contained in the difference
%     ( 1)  - The vertex is in the described halfspace and, hence, might be
%             contained in the difference, dependent on the other facets.
% For each simplex, there are three resulting cases:
%     (0) All vertices are case-(0) vertices. This should not happen, as this
%         implies that the simplex is not full-dimensional
%     (1) All vertices are case-(-1) or case-(0) vertices. The simplex is 
%         contained in the vrepResultList
%     (2) All vertices are case-(1) or case(0) vertices. The simplex might
%         be contained in the difference. They are stored in
%         vrepRemainingList.
%     (3) The vertices have mixed cases. This will be discussed in the
%         following.
% For case (3) of a simplex, we need to determine the extreme points of the
%   intersection between the plane and the simplex. These extreme points
%   can be described by:
%     xNew = lam*x_(-1) + (1-lam)*x_(1)
%   for any combination of simplex vertices x_(-1) of case-(-1) and x_(1)
%   of case-(1). These points must fulfill:
%     hrepSubtrahend.A(end, :) * xNew = hrepSubtrahend.h(end)
%   From here, the unknown lam can be determined:
%     a*(lam*x_(-1) + (1-lam)*x_(1)) = h
%     lam*(a*x_(-1) - a*x_(1)) = h - a*x_(1)
%     lam = (h - a*x_(1)) / (a*x_(-1) - a*x_(1)).
% All points xNew and the points x_(-1) are the V-representation of a
%   polytope that is definitely contained in the difference. All points
%   xNew and the points x_(1) are the V-representation of a polytope that
%   might be contained in the sum, dependent on the remaining facets. Both
%   polytopes are triangulated and the corresponding simplizes are added to
%   the appropriate lists.

% Along with the resultList, we store to which original simplex (from
%   minuendList) each simplex in resultList belongs. This makes it possible
%   to merge original simplices again that were splitted during the
%   algorithm.


% %% Determine parts of the difference
% % every combination of a hrep from hrepMinuend and a facet of
% %   hrepSubtrahend might create a new part. We initialize a cell array that
% %   contains the different hrep parts and a matrix that keeps the nonempty
% %   entries:
% hrepPartList = cell(nMinuend, nSubtrahendFacet);
% hrepPartNonemptyMatrix = zeros(nMinuend, nSubtrahendFacet, 'uint8');

%% Treat last facet of hrepSubtrahend
thisNormal = hrepSubtrahend.A(end,:);
thisDistance = hrepSubtrahend.h(end);
% start loop for minuends
% disp(['subtractPolytope: nMinuend = ' num2str(nMinuend)])
for iMinuend = 1:nMinuend
    thisSimplex = minuendVrepList{iMinuend};
    
    %% Classify vertices according to last facet
    distanceVector = thisNormal * thisSimplex.V' - thisDistance;
    filterMinus = distanceVector < (-1 * zerotol);
    filterPlus = distanceVector > zerotol;
    filterZero = ~filterMinus & ~filterPlus;
    
    %% Case (0) - Invalid
    if all(filterZero)
        error('elk:decomposition:numericalProblem', ...
              'A flat polytope (partition) was generated. This should not happen.');
    end
    
    %% Case (1) - Simplex fully contained
    if all(filterPlus | filterZero)
        resultVrepList{end+1} = thisSimplex;
        resultIndexVector(end+1) = iMinuend;
        if ~isempty(minuendVolumeVector)
            resultVolumeVector(end+1) = minuendVolumeVector(iMinuend);
        else
            resultVolumeVector(end+1) = computeVolumeVrep(thisSimplex, zerotol);
        end
    end
    
    %% Case (2) - Simplex fully a "maybe"
    if all(filterMinus | filterZero)
        remainingVrepList{end+1} = thisSimplex;
        remainingIndexVector(end+1) = iMinuend;
        if ~isempty(minuendVolumeVector)
            remainingVolumeVector(end+1) = minuendVolumeVector(iMinuend);
        else
            remainingVolumeVector(end+1) = computeVolumeVrep(thisSimplex, zerotol);
        end
    end
    
    %% Case (3) - Partially contained
    if any(filterPlus) && any(filterMinus)
        % prepare mid-point matrix
        midPointMatrix = nan([sum(filterMinus)*sum(filterPlus), nDim]);
        iMidPoint = 0;
        
        % determine mid-points
        for iCaseMinus = find(filterMinus)
        for iCasePlus = find(filterPlus)
            iMidPoint = iMidPoint + 1;
            
            % lam = (h - a*x_(1)) / (a*x_(-1) - a*x_(1))
            lam = (thisDistance - thisNormal * thisSimplex.V(iCasePlus, :)') / ...
                  (thisNormal * thisSimplex.V(iCaseMinus, :)' - ...
                   thisNormal * thisSimplex.V(iCasePlus, :)');
               
            % xNew =  lam*x_(-1) + (1-lam)*x_(1)
            midPointMatrix(iMidPoint, :) = ...
                lam * thisSimplex.V(iCaseMinus, :) + ...
                (1-lam) * thisSimplex.V(iCasePlus, :);
        end
        end
        
        % generate result simplex/simplices
        [plusVrepList plusVolumeVector] = triangulateVrep(struct('V', ...
            [thisSimplex.V(filterPlus, :); ...
             thisSimplex.V(filterZero, :); ...
             midPointMatrix...
             ]), zerotol, isCone); 
         % Note: DEGENERATE SPLITTING
         % When a boundary cuts the simplex only slighty because of
         %   numerical issues, then the simplex is split into a part that
         %   is roughly itself and a part that is flat/degenerate. Hence,
         %   the following error is not reasonable.
%         if isempty(thisPartPlus)
%             error('elk:decomposition:numericalProblem', ...
%                   'Polytope in set difference not splitted properly. This should not happen.');
%         end
        % But we still can ignore these "flat" polytopes:
        plusVrepList(plusVolumeVector < zerotol*scalingVolume) = [];
        plusVolumeVector(plusVolumeVector < zerotol*scalingVolume) = [];
        %/
        nPlus = length(plusVrepList);
        [resultVrepList{end + (1:nPlus)}] = plusVrepList{:};
        resultIndexVector(end + (1:nPlus)) = iMinuend;
        resultVolumeVector(end + (1:nPlus)) = plusVolumeVector(:);
        
        % generate remaining simplex/simplices
        [minusVrepList minusVolumeVector] = triangulateVrep(struct('V', ...
            [thisSimplex.V(filterMinus, :); ...
             thisSimplex.V(filterZero, :); ...
             midPointMatrix...
             ]), zerotol, isCone);
        
        % Commented out, see note "DEGENERATE SPLITTING", above.
%         if isempty(thisPartMinus)
%             error('elk:decomposition:numericalProblem', ...
%                   'Polytope in set difference not splitted properly. This should not happen.');
%         end
        minusVrepList(minusVolumeVector < zerotol*scalingVolume) = [];
        minusVolumeVector(minusVolumeVector < zerotol*scalingVolume) = [];
        %/
        nMinus = length(minusVrepList);
        [remainingVrepList{end + (1:nMinus)}] = minusVrepList{:};
        remainingIndexVector(end + (1:nMinus)) = iMinuend;
        remainingVolumeVector(end + (1:nMinus)) = minusVolumeVector(:);
    end
end

%% Recursively handle facet (end-1) from hrepSubtrahend
if length(hrepSubtrahend.h) > 1 && ~isempty(remainingVrepList)
    % delete last facet (already treated)
    hrepSubtrahend.A(end, :) = [];
    hrepSubtrahend.h(end) = [];
    
    % only the remaining part needs a decision whether they are contained 
    %   or not. This is a recursive call of this function.
    [vrepNextList nextVolumeVector nextIndexVector] = ...
        subtractPolytope(remainingVrepList, hrepSubtrahend, ...
        'zerotol', zerotol, 'isCone', isCone, ...
        'volumeVector', remainingVolumeVector, ...
        'scalingVolume', scalingVolume);
    
    % assign result
    [resultVrepList{end + (1:length(vrepNextList))}] = vrepNextList{:};
    resultIndexVector(end + (1:length(vrepNextList))) = ...
        remainingIndexVector(nextIndexVector);
    resultVolumeVector(end + (1:length(vrepNextList))) = nextVolumeVector(:);
end

% %% recover uneccesary splitting of simplices
% % the following is only performed, when the volume vector of the
% %   minuendList is provided

%! DEV: reprogram the following such that the lists/vectors are reduced in
% only one step. This step seems to take far too long. Same for the lists
% above.
%/
if ~isempty(minuendVolumeVector)
filterRemove =  false(1, length(minuendVrepList));
for iMinuend = 1:length(minuendVrepList)
    thisFilter = (resultIndexVector == iMinuend);
   
    if sum(thisFilter) > 1 && ...
            abs(sum(resultVolumeVector(thisFilter)) - ...
            minuendVolumeVector(iMinuend)) < zerotol*scalingVolume
        
        % get real indices
        indexVector = find(thisFilter);
        
        % we substitute the first element with the non-splitted original
        resultVrepList{indexVector(1)} = minuendVrepList{iMinuend};
        resultVolumeVector(indexVector(1)) = minuendVolumeVector(iMinuend);
        % .. index to original remains unchanged        
        
        % remove all other simplices
        filterRemove(indexVector(2:end)) = 1;
    end
end
resultVrepList(filterRemove) = [];
resultVolumeVector(filterRemove) = [];
resultIndexVector(filterRemove) = [];
end