function summandMatrix = decomposeMinkowskiHrep(hrep, varargin)
% decompose the given polytope into Minkowski summands
%
% Syntax: hrepList = decomposeMinkowskiHrep(hrep)
%
% THIS FUNCTION IS NO USER FUNCTION.

% ..event though, this function was originally designed to be. It might be
%   included in the polytope subsection in some later release. The main
%   missing feature is the generation of the summands itself.

% Draft documentation:
%
% Note that Minkowski decomposition is not unique. The result of this
%   algorithm gives indecomposable polytopes in hrepList so that it holds:
%       hrep = addMinkowskiPolytope(hrepList{1}, ..., hrepList{n}).
%   And because every hrepVector contains the original full matrix of facet
%   normals this is equivalent to:
%       hrep.h = hrepList{1}.h + .. + hrepList{n}.h
%
% The following options can be appended:
%   zerotol  - tolerance for singularities  and close to zero operations.
%              Default: 1e-13
%   verbose  - switch on (1) and off (0) some verbosity. Default: 0
%   matrix   - instead of a lost of H-representations a matrix with the 
%              facet distances is returned (1). One row for each summand, 
%              one column for each facet in hrep.A. Default: 0
%   result   - This algorithm is also the main component for the 
%              decomposition of the constrained h-space so that with the 
%              option set to 'indecomposable', all indecomposable summands
%              are returned. Default: 'summand'
%
% See also: addMinkowskiPolytope, addMinkowskiHrep

% 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/>.

%% sources
% The original sketch (and proof) of the algorithm can be found in Meyer,
%   1974: "Indecomposable Polytopes". It is modified and presented again in 
%   Mount and Silverman, 1991: "Combinatorial and Computational Aspects of 
%   Minkowski Decomposition". The algorithm here is adopted for constrained 
%   Minkowski decomposition. It also has an easy to implement approach to 
%   eliminate the degrees of freedom by translation. In Meyer 1974 the
%   Steiner-point was used.

%% Input
% some numbers
nDim = size(hrep.A, 2);
nFacet = size(hrep.A, 1);

% map options
[zerotol, verbose, outputResult, outputMatrix, constraintMatrix] = mapOption(varargin, ...
    'zerotol', elkZerotol, ...
    'verbose', 0, ...
    'result', 'summand',...
    'matrix', 0, ...
    'constraint', zeros(1, nFacet));

% check polytope vs. constraint
if ~isAllZero(constraintMatrix * hrep.h, zerotol)
    error('elk:decomposition:inputError', ...
        'The input polytope must fulfill the provided constraints');
end

%% Fix position & Scaling
% get additional constraint and move hrep
[positionConstraintMatrix, mapFacetDistanceToOffset] = ...
    fixPosition(hrep.A, constraintMatrix, zerotol);
positionFixOffset = mapFacetDistanceToOffset * hrep.h;
hrepPositionFixed = moveHrep(hrep, positionFixOffset);

% Scaling
% normal scaling cannot be applied according to position fixing. The only
%   thing that can be done is scaling the h-vector. We simply scale the
%   largest facet distance to one.
scaleFactor = 1/max(abs(hrepPositionFixed.h));
hrepScaled = hrepPositionFixed;
hrepScaled.h = hrepScaled.h * scaleFactor;

% check final polytope vs. input and position constraint
if ~isAllZero(constraintMatrix * hrepScaled.h, zerotol)
    error('elk:decomposition:wrongInput', ...
        'The input polytope must fulfill the given constraints');
end
if ~isAllZero(positionConstraintMatrix * hrepScaled.h, zerotol)
    error('elk:decomposition:numericalProblem', ...
        'The scaled polytope does not fulfill the fixed position constraint. This should not happen!');
end

if verbose,disp('Generating incidence matrix');end

%% Generate incidence matrix
% Facets, incident to a vertex are stored in an incidence matrix. Rows for 
% vertices, columns for facets.

% extreme points and initialization
vrepScaled = convertHrepToVrep(hrepScaled, zerotol);
nExtreme = size(vrepScaled.V, 1);
incidenceVerticesToFacets = zeros(nExtreme, nFacet, 'uint8');

for iExtreme = 1:nExtreme
    incidenceVerticesToFacets(iExtreme, :) = isZero( ...
	    hrepScaled.A * vrepScaled.V(iExtreme, :)' - hrepScaled.h, ...
	    zerotol);
end

%% Equations and Inequalities
% The following loop goes through all extreme point / facet 
%   combinations. While it follows three states:
%    (1) finding the representative facets
%    (2) adding equality constraints
%    (3) adding inequality constraints
% Both matrices are stored in the form that each row is a equation
%   (or inequality) and each column corresponds to a facet.

if verbose,disp('Extracting equations and inequalities');end

% initialize equality/inequality matrices
nEquation = sum(sum(incidenceVerticesToFacets, 2) - nDim);
nInequality = sum(sum(~incidenceVerticesToFacets, 2));
equationMatrix = zeros(nEquation, nFacet);
inequalityMatrix = zeros(nInequality, nFacet);
thisEquationIndex = 0;
thisInequalityIndex = 0;

% loop over extreme points
for iExtreme = 1:nExtreme

	thisAdjacentFacetVector = find(incidenceVerticesToFacets(iExtreme, :));
	thisNonadjacentFacetVector = find(~incidenceVerticesToFacets(iExtreme, :));
	
	% find representative facets
	thisRank = 0;
	repFacetMatrix = zeros(nDim, nDim);
	repFacetIndexVector = zeros(1, nDim);
	for iFacet = thisAdjacentFacetVector
		repFacetMatrix(thisRank+1, :) = ...
			hrepScaled.A(iFacet, :);
		if rank(repFacetMatrix) > thisRank
			thisRank = thisRank + 1;
			repFacetIndexVector(thisRank) = iFacet;
			if thisRank == nDim, break;end
		end
	end
	
	% error if representatives not found
	if thisRank < nDim
		error('elk:decomposition:numericalProblem', ...
			'Failed to find representative facets. This should not happen!');
	end
	
	% the inversion of the representative facet normals
	% this is used together with the fourth facet normal (a4). For all 
	%   four one can find Li that satisy (L4 = 1):
	%     A' * [L1, L2, L3]' = -L4*a4
	%   by using:
    %     [L1, L2, L3] = -L4 * a4' * inv(A)
	invertedRepresentativeNormalMatrix = inv(repFacetMatrix);
	
	% eliminating representatives from adjacency vector
	thisAdjacentFacetVector = eliminateRows(thisAdjacentFacetVector', ...
        repFacetIndexVector', eps)';
	
	% adding equations
	% Knowing L1, L2, L3 and L4 we gather the equation:
    %   [L1 L2 L3 L4] * [h1; h2; h3; h4] = 0
    for iFacet = thisAdjacentFacetVector
        thisEquationIndex = thisEquationIndex + 1;
        equationMatrix(thisEquationIndex, iFacet) = 1;
        lambdaVector = -1 * hrepScaled.A(iFacet, :) * ...
            invertedRepresentativeNormalMatrix;
        equationMatrix(thisEquationIndex, ...
            repFacetIndexVector) = lambdaVector;
    end

	% adding inequalities
    % Knowing L1, L2, L3 and L4 we gather the equation:
    %   -1 * [L1 L2 L3 L4] * [h1; h2; h3; h4] <= 0
    for iFacet = thisNonadjacentFacetVector
        thisInequalityIndex = thisInequalityIndex + 1;
        inequalityMatrix(thisInequalityIndex, iFacet) = -1;
        lambdaVector = 1 * hrepScaled.A(iFacet, :) * ...
            invertedRepresentativeNormalMatrix;
        inequalityMatrix(thisInequalityIndex, ...
            repFacetIndexVector) = lambdaVector;
    end
	
end

if thisEquationIndex < nEquation
	error('elk:decomposition:numericalProblem', ...
		'Failed to add all equations. This should not happen!');
end
if thisInequalityIndex < nInequality
	error('elk:decomposition:numericalProblem', ...
		'Failed to add all inequalities. This should not happen!');
end

if verbose,disp(['Generated ' num2str(nEquation) ' equations and ' ...
        num2str(nInequality) ' inequalities']);end

if verbose,disp('Solving equations and inequalities');end

%% Solve equations and inequalities
summandMatrix = computeExtremePointsOfCone(...
    [equationMatrix; constraintMatrix; positionConstraintMatrix], ...
    inequalityMatrix, zerotol);
nSummand = size(summandMatrix, 1);

%% Output handling

% generate summands or not
switch lower(outputResult)
    case 'indecomposable'
        %! ToDo: there should be some scaling. But this does not seem to be
        %        trivial with fixed position and in hrep.
        %  Update: Scaling is now possible in a reasonable manner, using
        %    the mean width (implemented in scaleToReal)
        %/
    case 'summand'
        error('this option is currently switched off');
        % get scalar factors
        % For any number of summands we search a solution for:
        %   summandMatrix' * lam = h with lam >= 0.
        % We simply use an optimizer with the function to be minimized:
        %   (1:nSummands) * lam
%         cddlibOutput = callCddLib('solve_lp', struct(...
%             'A', [summandMatrix'; eye(nSummand)], ...
%             'B', [hrepScaled.h; zeros(nSummand, 1)], ...
%             'lin', 1:nFacet, ...
%             'obj', nSummand:-1:1 ...
%             ));
%         scalarMultiples = cddlibOutput.xopt;

%! ToDo: This part does currently not work. I could need to triangulate the
%        cone to find an appropriate simplex partition for which the
%        condition lam > 0 is automatically fulfilled.
        
%         summandMatrixReduced = (projectionMatrix * summandMatrix')';
%         facetDistanceVectorReduced = (projectionMatrix * hrepScaled.h)';
        
%         triangulationIndexMatrix = computeTriangulationOfCone(summandMatrixReduced);
%         for iSimplex = 1:size(triangulationIndexMatrix, 1)
%             thisRayMatrix = summandMatrixReduced(triangulationIndexMatrix(iSimplex, :), :);
%             thisMultipleVector = inv(thisRayMatrix') * facetDistanceVectorReduced';
%             if any(thisMultipleVector < 0)
%                 scalarMultipleVector = [];
%                 continue
%             else
%                 scalarMultipleVector = zeros(size(summandMatrixReduced, 1), 1);
%                 scalarMultipleVector(triangulationIndexMatrix(iSimplex, :)) = ...
%                     thisMultipleVector;
%                 break;
%             end
%         end
%         
%         if isempty(scalarMultipleVector)
%             error('elk:decomposition:numericalProblem', ...
%                 'Could not determine scalar coefficients for summands. This should not happen.');
%         end
%         
%         for iSummand = 1:nSummand
%             summandMatrix(iSummand, :) = summandMatrix(iSummand, :) * ...
%                                          scalarMultipleVector(iSummand);
%         end
%         
%         % throw away almost zero rows
%         summandMatrix(...
%             prod(...
%             double(isZero(summandMatrix, zerotol))...
%             , 2) > 0, :) = [];
%         nSummand = size(summandMatrix, 1);
    otherwise
        error('elk:decomposition:inputIssue', ...
            'The option ''result'' is either ''summand'' or ''indecomposable''');
end

% revert scaling
summandMatrix = summandMatrix / scaleFactor;

% revert position fix
% this should only be done for the first summand and only for result type
%   'summand'
if strcmpi(outputResult, 'summand')
    summandMatrix(1,:) = (summandMatrix(1,:)' - hrep.A * positionFixOffset)';
end

switch outputMatrix
    case 1
        % nothing to do.
        
    case 0
        % convert to a cell array of h-representations
        hrepTemplate.A = hrep.A;
        hrepList = cell(nSummand, 1);
        for iSummand = 1:nSummand
            hrepList{iSummand} = hrepTemplate;
            hrepList{iSummand}.h = summandMatrix(iSummand, :)';
        end
        summandMatrix = hrepList;
        
    otherwise
        error('elk:decomposition:wrongInput', ...
            'The option ''matrix'' is either 1 or 0');
end

