function result = gatherValidityConeIneq(facetNormalMatrix, zerotol, verbose)
% generate the inequality constraints for valid polytopes in h-space
%
% THIS FUNCTION IS NO USER FUNCTION.

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

%% info
% This function generates inequality constraints in form of a
%   H-representation in the h-space that is associated with the provided
%   matrix of facet normals. For valid solutions the polytope is not empty
%   (at least contains a single point) and all support planes have a
%   nonempty intersection with the given polytope. The formal condition is
%   then: the polytope is valid for:
%     inequalityMatrix * h <= 0
% This function also returns in cone.responsibleFacetMatrix a matrix that
%   stores pairs of inequality indices (first column) and positive indices
%   to the facet that grows out at this boundary. Negative indices are
%   entered if the dissolution of that facet triggers an empty polytope.

%% Input handling
% some statistical numbers
nDim = size(facetNormalMatrix, 2);
nFacet = size(facetNormalMatrix, 1);

%% initialize matrices
% generate tuples to check
tupleMatrix = obtainOrderedSetMatrix(nFacet, nDim+1);
% using a sparse matrix was tested, but does not pay off.
inequalityMatrix = zeros(nDim, nFacet);
responsibleFacetVector = nan(nDim, 1);
thisInequalityIndex = 0;

%% start loop for all tuples
if verbose,fprintf('......');end
for iTuple = 1:size(tupleMatrix, 1)
    if verbose && mod(iTuple, 1000) == 1
        fprintf('\b\b\b\b\b\b%5.1f%%', iTuple / size(tupleMatrix, 1)*100);
    end
    
    thisTuple = tupleMatrix(iTuple, :);
     
    % find the n-tuple where the (n+1) element is in the positive OR 
    %   negative cone and the cone is full-dimensional. The last element is
    %   contained in the positive span, if:
    %     [a(1), .., a(n)] * lam = a(n+1)  and  lam >= 0
    %   and in the negative span for lam <= 0.
    coneType = 0; % (not found)
    
    for iFacet = thisTuple
        thisIndexVector = thisTuple(thisTuple ~= iFacet);
        [u,s,v] = svd(facetNormalMatrix(thisIndexVector, :));
        s = diag(s);
        isNotZeroVector = (s <= -zerotol) | (s >= zerotol);
        if all(isNotZeroVector)
            s(isNotZeroVector) = 1./s(isNotZeroVector);
            lambda = u * diag(s) * v' * facetNormalMatrix(iFacet, :)';
                           
            if all(lambda >= -zerotol)
                % tuplet combination found
                thisTuple = [thisIndexVector iFacet];
                coneType = 1; % (positive cone)
                break
            elseif all(lambda <= zerotol)
                % tuplet combination found
                thisTuple = [thisIndexVector iFacet];
                coneType = -1; % (negative cone)
                break
            end
        end
    end
    
    %! Info: There is an alternative implementation, that uses fewer lines
    %    of code and might be quicker to grasp, put it is not faster:
    
%     % switch some warnings to errors to detect lower rank during mldivide
%     warningStateOne = warning('error', 'MATLAB:nearlySingularMatrix');
%     warningStateTwo = warning('error', 'MATLAB:singularMatrix');
%     
%     for iFacet = thisTuple
%         thisIndexVector = thisTuple(thisTuple ~= iFacet);
%         % use mldivide directly and catch low rank
%         try
%               lambda = facetNormalMatrix(thisIndexVector, :)' \ ...
%                 facetNormalMatrix(iFacet, :)';
%         catch
%               continue
%         end
%         %! Assignement as above!
%     end
%     % reset warning states
%     warning(warningStateOne);
%     warning(warningStateTwo);
    %/
    
    % result: we have evaluated the tuple to form a negative cone 
    %   (coneType=-1), positve cone (coneType=1) or none of these
    %   (coneType=0). The array thisIndexVector also have the indces to the
    %   cone-forming facet normals in the first three elements (that
    %   defines the vertex). It contains the forth facet as the last
    %   element.
    
    % continue with next tuple (this tuple beeing invalid) or continue with
    %   checking locality
    if coneType == 0
        %disp([num2str(iTuple) ': no full dimensional positive/negative span'])
        %thisNormalMatrix
        continue
    else
        % this tuple now contains the same (n+1) normal vectors that we
        %   started with, but the first n generate the positive/negative 
        %   cone where the (n+1)st element is an element of.
        thisNormalMatrix = facetNormalMatrix(thisTuple, :);
    end
    
    % the disabled code represents improvements by Karl (1995, "local 
    %   conditions") that were found to give no significant improvement.
%     % check if the n-tuple is local
%     % Another facet normal is part of the positive cone, whenever we have:
%     %   [a(1), .., a(n)] * lam = a(i)  and  lam >= 0
%     % which can be rewritten to:
%     %   inv([a(1), .., a(n)]) * a(i) >= 0
%     % the same facet normal is part of the negative cone, for:
%     %   -inv([a(1), .., a(n)]) * a(i) >= 0
%     notInTuple = 1:nFacet;
%     notInTuple(thisTuple) = [];
    if coneType == 1
        testMatrix = inv(thisNormalMatrix(1:nDim, :)');
    else
        testMatrix = inv(thisNormalMatrix(1:nDim, :)');
    end
%     exitFlag = 0;
%     for iFacet = notInTuple
%         lambda = testMatrix * facetNormalMatrix(iFacet, :)';
%         % as this is a "we save equations" condition, less samples should
%         %   be valid to fulfill this condition.
%         if all(lambda >= zerotol)
%             % there is another facet in the positive/negative span. Hence,
%             %   the tuple is not local.
%             exitFlag = 1;
%             break
%         end
%     end
%     if exitFlag
%         %disp([num2str(iTuple) ': not local positive span'])
%         %thisNormalMatrix
%         continue;
%     else
%          %disp([num2str(iTuple) ': add constraint'])
%     end
    
    % now, add the condition
    % from Karl (1995) we have:
    %   a(n+1)' * inv([a(1)'; ...; a(n)']) * [h(1); ..; h(n)] >= h(n+1)
    % which is equal to:
    %   (inv([a(1), ..., a(n)]) * a(n+1))' * [h(1); ..; h(n)] >= h(n+1)
    % where the inverse of the second equation is the testMatrix from above
    % for positive cones. 
    % if coneType == -1
    %     testMatrix = -1*testMatrix;
    % end
    % Further more, for negative cones the inequality sign changes so that
    %   we have for negative cones:
    %     (inv([a(1), ..., a(n)]) * a(n+1))' * [h(1); ..; h(n)] - h(n+1) <= 0
    %
    % The common H-representation of the cone uses <= instead of >= so that
    %   for the positive cone, we obtain:
    %   -1*(inv([a(1), ..., a(n)]) * a(n+1))' * [h(1); ..; h(n)] + h(n+1) <= 0
    
    thisInequalityIndex = thisInequalityIndex + 1;
     
    % the following equals inv([a(1), ..., a(n)]) * a(n+1)
    coefficientVector = (testMatrix * thisNormalMatrix(end, :)')';
    
    if coneType == 1
        inequalityMatrix(thisInequalityIndex, thisTuple(end)) = 1;
        inequalityMatrix(thisInequalityIndex, thisTuple(1:(end-1))) = ...
           -1 * coefficientVector;
        responsibleFacetVector(thisInequalityIndex) = thisTuple(end);
    elseif coneType == -1 
        inequalityMatrix(thisInequalityIndex, thisTuple(end)) = -1;
        inequalityMatrix(thisInequalityIndex, thisTuple(1:(end-1))) = ...
           coefficientVector;
        responsibleFacetVector(thisInequalityIndex) = -1*thisTuple(end);
    else
        error('elk:decomposition:numericalError', ...
            'cone is neither positive nor negative. This should not happen');
    end
    
    % reserve more space in sparce matrix, if required
    if thisInequalityIndex >= size(inequalityMatrix, 1)
        % allocate more space
        estimatedRowSize = thisInequalityIndex / iTuple * size(tupleMatrix, 1);
        % round and take some more
        estimatedRowSize = round(estimatedRowSize) + nDim;
        % enlarge matrix
        inequalityMatrix(estimatedRowSize, :) = 0;
        responsibleFacetVector((end+1):estimatedRowSize) = nan;
    end
    
end
if verbose, fprintf('\b\b\b\b\b\b');end

%% remove the reservated rows
inequalityMatrix((thisInequalityIndex+1):end, :) = [];
responsibleFacetVector((thisInequalityIndex+1):end, :) = [];

result.A = inequalityMatrix;
result.responsibleFacetVector = responsibleFacetVector;