function projectionData = generateOrthogonalProjectionDataForCone(A, ...
    extremeRayMatrix, zerotol)
% compute projection matrices and validity regions for orthogonal mapping
%
% Syntax: projectionData = generateOrthogonalProjectionDataForCone(A)
%
% This function generates a set of projection matrices P and validity 
%   regions {B*x <= 0} so that y=P*x is the closest point with (A*y <= 0)
%   for x. The output herefore is a data structure with fields:
%     .constraintMatrix - all possible rows of B
%     .projectionList - cell array with all projection matrices P
%     .responsibleMatrix - a matrix that assigns the rows of the constraint
%       matrix to the matrices B for each P. Each rows contains in the
%       first column a row-index of the constraintMatrix, in the second
%       column an index to a projection and in the third column +1 or -1
%       that might change the sign of the row in B (see example).
%     .targetDimensionVector - the target dimension of each projection
% 
% Alternatively, you can use the inputs ..(A, extremeRayMatrix, zerotol) to
%   provide the a prior calculated extreme rays of the cone A and set the
%   tolerance.
%
% See also: obtainDec, approximateDec, projectFromData

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

%% input handling

% zerotol
if ~exist('zerotol', 'var') || isempty(zerotol)
    zerotol = elkZerotol;
end

nDim = size(A,2);

% extreme ray matrix
if ~exist('extremeRayMatrix', 'var') || isempty(extremeRayMatrix)
    hrep.A = [A; mean(A, 1)];
    hrep.h = [zeros(size(A, 1), 1); 1];
    vrep = convertHrepToVrep(hrep);
    extremeRayMatrix = vrep.V;
    % eliminate origin
    extremeRayMatrix = eliminateRows(extremeRayMatrix, ...
        zeros(1, size(A, 2)), zerotol);
end

%% Catch empty cones
% The provided cone might be trivial so that only the origin itself
%   fulfills A*x <= 0. We notify this, if we try to calculate the
%   V-representation, which results in a single vertex (the origin) in that
%   case.
vrep = convertHrepToVrep(struct('A', A, 'h', zeros(size(A, 1), 1)));
if size(vrep.V, 1) == 1 && all(vrep.V < zerotol)
    projectionData.A = [1 zeros(nDim-1, 0)];
    projectionData.projectionList = {zeros(nDim, nDim), zeros(nDim, nDim)};
    projectionData.responsibleMatrix = [1 1 1; 1 2 -1];
    projectionData.targetDimensionVector = [0 0];
    return
end

%% Trivial cases
% especially simple are the cases where P is either the zero-matrix and
%   where P is the idendity matrix.

% In case of P being the identity matrix, the point x is already valid.
projectionData.A = A;
projectionData.responsibleMatrix = [(1:size(A,1))', ones(size(A, 1), 1), ...
    ones(size(A, 1), 1)];
projectionData.projectionList = {eye(nDim)};
projectionData.targetDimensionVector = nDim;

% In case of P being the zero-matrix, the closest point is the origin.
%   This projection is valid for any point in the positive cone of the
%   normal vectors in A.
vrep.V = [A; zeros(1, nDim)];
hrep = convertVrepToHrep(vrep);
% filter facets that do not hit 0
hrep.A(abs(hrep.h) > zerotol, :) = [];
% assign output (we accept double entries, they are unlikely)
projectionData = appendProjectionData(projectionData, hrep.A, ...
    zeros(nDim, nDim), zerotol);
projectionData.targetDimensionVector(end+1) = 0;

%% adjacency information
% which facet (row) is adjacent to which extreme ray (column)
adjacencyMatrix = abs(A * extremeRayMatrix') < zerotol;

% obtain all combinations of 0 up to nC selected elements from the
%   boundaries in A
permutationMatrix = concatenateOrderetSetMatrices(size(A,1), nDim, 1);
%! DEV: I've put concatenateOrderetSetMatrices as a seperate function under
%    private and adopted it to also return selecting nDim elements. This
%    should change things for here, even though it matches the expectation
%    of what concatenateOrderetSetMatrices should do.

% initialize
% constraintMatrix = zeros(0, nDim);
% responsibleMatrix = zeros(0, 3);
% projectionList = {};
% targetDimensionVector = [];

for iPerm = 1:size(permutationMatrix, 1)
    %% general case information
    % we use this permutation, now
    thisPerm = permutationMatrix(iPerm, :);
    % delete 'no facet selected' entries
    thisPerm(thisPerm == 0) = [];
    
    % basic information
    thisNullSpaceDim = length(thisPerm);
    thisTargetDim = nDim - length(thisPerm);
    thisNormalMatrix = A(thisPerm, :)';
    
    % skip invalid combinations: each set of facet normals must have 
    %   targetDim extreme rays in common. The intersection of nF facets
    %   leaves an (nDim - nF) dimensional subspace that is defined by (nDim
    %   - nF) points, the extreme points.
    filterMatrix = adjacencyMatrix(thisPerm, :);
    filterVector = sum(filterMatrix, 1) == thisNullSpaceDim;
    if ~isempty(thisPerm) && (sum(filterVector) < thisTargetDim)
        % the nF facets are not adjacent to (nDim-nF) extreme rays AND they
        % are not empty: skip this case
        continue
    end
    
    %% projection and condition
    % projection
    thisNullProjection = thisNormalMatrix * ...
        inv(thisNormalMatrix'*thisNormalMatrix)*thisNormalMatrix';
    thisProjection = eye(nDim) - thisNullProjection;
      
    % condition
    thisConstraintMatrix = [A*thisProjection; ...
                            -1*pinv(thisNormalMatrix)];
    % filter zero rows (happens for the facets that is projected onto)
    thisConstraintMatrix = eliminateRows(thisConstraintMatrix, ...
        zeros(1, nDim), zerotol);
    % normalize again
    thisConstraintMatrix = bsxfun(@times, thisConstraintMatrix, ...
        1./sqrt(sum(thisConstraintMatrix.^2, 2)));
    
    %% add information
    % add projection
    projectionData = appendProjectionData(projectionData, ...
        thisConstraintMatrix, thisProjection, zerotol);
    projectionData.targetDimensionVector(end+1) = thisTargetDim;
    
end