function data = generateSampleProjectionData(heMatrix, dec, par)
% determine active facet validity boundaries, corresponding validity
%   projection and applicable region for projection
%
% THIS IS NO USER FUNCTION

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

% The input hcMatrix should not contain zero vectors. The approximation
%   will never contribute to any error (always perfectly approximated.
% The input does also not contain invalid vectors. But vectors at the
%   boundary of the validity cone for which the information in this
%   function is actually generated.

%! DEV: This code lines are somewhat badly written, it is sometimes not
%       clear what projection index is meant between the different lists
%       and purposes (last for loop).
%/

% For the initialization of the projection data, we already assume that
%   some vectors will be valid. But even if this is not the case and the
%   index 1 is not used in the end, this is OK.
nE = size(heMatrix, 1);
nSample = size(heMatrix, 2);

data.list.validityProjection = eye(nE);
data.list.facetPresentConstraint = zeros(0, nE);
data.list.projectionApplicableConstraint = zeros(0, nE);
data.assignedProjectionVector = ones(1, nSample);

% shortcut for PCA objective, incomplete representations and valid
%   projections

%! DEV: removed options errorType=PCA and ~dec.isComplete, for these cases
%       we have to apply the 'true' projection somewhere else. The
%       arguments were:
% 1) everey vector stays projected as it is, no fancy things to do
% 2) For an incomplete representation, we can even have invalid vectors
%   during measure calculation as we have only one unified partition.
%   The aim of extended embeddings (incomplete representations) IS to
%   use no shape-dependent stuff like boundaries.
%/
if strcmpi(par.projectionType, 'true')
    return
elseif strcmpi(par.projectionType, 'valid')
    % the resulting vectors must be valid
    data.list.facetPresentConstraint = dec.validityCone.A;
    return
    %! DEV: When is that required??
elseif strcmpi(par.projectionType, 'confined')
    % the gibberish below will be executed. ..nothing to do
else
    error('elk:approximation:internalError', ['The projection type is ' ...
          'unknown. This should not happen.'])
end

% Cycling through all facet groups might save some time, even though the
%   handling of new projections seems to be more complicated.
for iGroup = 1:nE
    nProjection = max(data.assignedProjectionVector);
    
    % the constraint for this facet group being present
    facetPresentConstraint = dec.confinementMappingData.A(...
        dec.confinementMappingData.responsibleVector == iGroup, :);  
    % find vectors with corresponding disappeared facets while they are
    %   slightly more quickly assumed to be at this boundary
    heFilter = any(facetPresentConstraint * heMatrix ...
        > -1*par.zerotol, 1);
    heIndexVector = find(heFilter);
    
    % if the current facet group is present for all vectors, we make a
    %   short circuit, here
    if all(~heFilter)
        for iProjection = 1:nProjection
            data.list(iProjection).facetPresentConstraint = [...
                data.list(iProjection).facetPresentConstraint; ...
                facetPresentConstraint];
        end
        continue;
    end
    
    % obtain projection matrices and constraints for the applicability of
    %   the projection - we need to apply some larger zerotol than for
    %   filtering the initial vectors because this constraint is also
    %   contained there
    [dump projectionIndexVector constraintList] = projectFromData(...
        heMatrix(:, heFilter), ...
        dec.confinementMappingData.projectionData(iGroup), ...
        2*par.zerotol, 1);
    groupProjectionVector = unique(projectionIndexVector);
    nProjectionNew = length(groupProjectionVector);
    
    % cycle through the original projections and see, if they need to be
    %   split up
    for iProjectionOriginal = 1:nProjection
        % copy original projection and constraints
        listOriginal = data.list(iProjectionOriginal);
        % affected vectors
        thisVectorFilter = heFilter & (...
            data.assignedProjectionVector == iProjectionOriginal);
        
        if any(~thisVectorFilter)
            % no projection is concatenated for this projection since
            %   all assigned vectors have the current facet group being
            %   present
            
            % ..projection remains
            % ..necessarily present facets are expanded
            data.list(iProjectionOriginal).facetPresentConstraint = [...
                listOriginal.facetPresentConstraint; ...
                facetPresentConstraint];
            % ..no contraint on the applicability of the projection
            % ..the assigned vectors change, but this is done when the new
            %   projections are created
            
            % new projection entries must be created for each new
            % projection
            thisProjectionIndexVector = nProjection + (1:nProjectionNew);
        elseif all(thisVectorFilter)
            % this entry can be used to contain the very first projection
            thisProjectionIndexVector = [iProjectionOriginal ...
                (nProjection + 1:(nProjectionNew-1))];
        end
        
        % expand projection list, but make sure to do this last entry to
        %   first entry as the first entry overwrites the original data.
        for iProjectionNew = nProjectionNew:-1:1
            thisProjectionListIndex = thisProjectionIndexVector(iProjectionNew);
            thisGroupProjectionIndex = groupProjectionVector(iProjectionNew);
            
            % ..projection is extended
            data.list(thisProjectionListIndex).validityProjection = ...
                listOriginal.validityProjection * ...
                dec.confinementMappingData.projectionData(iGroup...
                ).projectionList{thisGroupProjectionIndex};
            % ..necessarily present facets remain (this group is not
            %   present)
            data.list(thisProjectionListIndex).facetPresentConstraint = ...
                listOriginal.facetPresentConstraint;
            % ..constraints for the applicability of the projection are
            %   expanded
            data.list(thisProjectionListIndex).projectionApplicableConstraint = ...
              [listOriginal.projectionApplicableConstraint; ...
               constraintList{find(...
                 projectionIndexVector==thisGroupProjectionIndex, ...
                 1, 'first')}...
               ];
            % ..vectors need to be assigned, here
            data.assignedProjectionVector(heIndexVector(...
                projectionIndexVector == thisGroupProjectionIndex)) = ...
                thisProjectionListIndex;
            
        end
    end
end