function heMatrixProjected = estimateHeMatrix(heMatrix, ...
    projectionMatrix, sampleProjectionData)
% compute the "closest" hE-vector in reach of the given hC-representation
%
% Distinguishing confined/valid projections and the PCA approach is ensured
%   via the construction of sampleProjetionData and/or the evaluation
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2014 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/>.


%% comments (based on confinement projection)
% We assume appropriate confinement partitions exist so that the sample
%   points are approximated by shapes for which the same validity
%   projection is valid. The target confinement partition is concluded 
%   directly from the sample vector. 
% If this assumption is not proper, we cannot simply ignore these points 
%   when the error measure is calculated. Furthermore, we have to avoid 
%   jumps in the objective function.
% To resolve the first issue, a ratio is defined that describes, how
%   many columns of heMatrix contribute to the error:
%      ratio = sum(notProperVector(hcMatrix)) / nCol.
%   Invalid vectors in the hcMatrix are substituted by their original
%   values so that these cannot contribute to the error anymore.
% To resolve the second issue, we have to switch continuously from the
%   proper state to the improper state. This is, fortunately relatively
%   simple when the distance to the constreints is evaluated accordingly.
%   A parameter "distanceLimit" controls this behavior, analogously to the 
%   growth rate validity mapping. However, this time we use an approximate 
%   of the outer distance so that the fade-out does start when the vectors 
%   already have become invalid. (This outer distance approximate is NOT the
%   distance from the validity cone).
% to additionally avoid that NO vectors are proper, a penalty scheme is
%   added at the very end. It fixes optimizations that do converge towards
%   "no vectors", but require a good knowledge or tuning to set the slope
%   according to the expected error values. This penalty scheme could also
%   very simply be used to optimize for "ensure that the existent facets
%   match for original and approximation for as much as posisble vectors".
    
% mapping matrices
mappingFullToReduced = orth(projectionMatrix)';
mappingReducedToFull = pinv(mappingFullToReduced);

% function handle to convert distances to "smoothing-filter" values
% smoothFilter = @(distanceVector)(1*(distanceVector <= 0) + ...
%     (distanceVector <= par.distanceLimit).*(distanceVector > 0).*...
%             (par.distanceLimit - distanceVector)/par.distanceLimit);

% we cycle through the groups of projections to generate the
%   approximated vectors
heMatrixProjected = nan*heMatrix;
% fadeOutScaling = zeros(1, size(heMatrix, 2));
for iProjection = 1:length(sampleProjectionData.list)
    heFilter = sampleProjectionData.assignedProjectionVector == iProjection;
    
    % generate the closest vectors in the subspace given by the current
    %   projection, assuming a confinement region where the
    %   pre-calculated validity projection is valid
    % when the validityProjection is the unit matrix, we should use the
    %   initial projectionMatrix instead of the mapping matrices for
    %   accuracy reasons
    thisValidityProjection = sampleProjectionData.list(iProjection).validityProjection;
    if all(diag(thisValidityProjection) == 1) && ...
       all(all(abs(thisValidityProjection - diag(diag(thisValidityProjection))) <= eps))
        heMatrixInSubspace = projectionMatrix * heMatrix(:, heFilter);
    else
        heMatrixInSubspace = mappingReducedToFull*pinv(...
            thisValidityProjection * mappingReducedToFull) * ...
            heMatrix(:, heFilter);
    end
%     projectionApplicableDistance = max(...
%         sampleProjectionData.list(iProjection).projectionApplicableConstraint * ...
%         heMatrixInSubspace, [], 1);
    % the projection is smoothened out while we assume this projection
    %   is valid for a few MORE points than exact.
%     if isempty(projectionApplicableDistance)
%         projectionApplicableDistance = ones(1, sum(vectorFilter));
%     else
%         projectionApplicableDistance =  smoothFilter(projectionApplicableDistance);
%     end
    
    % generate the corresponding valid hC-vectors
    heMatrixValid = sampleProjectionData.list(iProjection).validityProjection * ...
        heMatrixInSubspace;
%     facetsValidDistance = max(...
%         sampleProjectionData.list(iProjection).facetPresentConstraint * ...
%         heMatrixValid , [], 1);
%     if isempty(facetsValidDistance)
%         facetsValidDistance = ones(1, sum(vectorFilter));
%     else
%         facetsValidDistance = smoothFilter(facetsValidDistance);
%     end
    
    % join applicability
%     vectorSmoothFilter = min(projectionApplicableDistance, ...
%         facetsValidDistance);
    
    % assign values
%     heMatrixProjected(:, vectorFilter) = ...
%         bsxfun(@times, heMatrixValid, vectorSmoothFilter) + ...
%         bsxfun(@times, heMatrix(:, vectorFilter), (1-vectorSmoothFilter));
    heMatrixProjected(:, heFilter) = heMatrixValid;
%     fadeOutScaling(vectorFilter) = vectorSmoothFilter;
end

if isnan(heMatrixProjected)
    error('elk:approximation:internalError', ['Estimated projection of ' ...
        'hE-matrix contains nan. This should not happen.']);
end

% fadeOutScaling = ones(1, size(heMatrix, 2));