function [heMatrixOut targetUnifiedPartition] = projectVectorDec(heMatrix, dec, zerotol, ...
    useTrueMapping)
% approximize hE-vectors to the domain covered by an improper decomposition
%
% hcMatrix = projectHcVector(dec, heMatrix)
%
% see also: obtainDec, projectDec, mapToValidityCone

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

%% input handling
if ~exist('zerotol', 'var')
    zerotol = elkZerotol;
end
if ~exist('useTrueMapping', 'var')
    useTrueMapping = 0;
end

if useTrueMapping
    error('elk:development', 'not yet supported');
    %% true mapping
    heMatrixOut = decProj.properData.projectionMatrix * heMatrix;
    %! ToDo: The points in the heMatrix still need to be projected to the
    %    validity cone of the improper decomposition
    %/
elseif ~dec.isComplete
    %% incomplete decomposition data (only true possible)
    heMatrixOut = dec.properData.projectionMatrix * heMatrix;
else
    %% wrapped    
    % ensure input is valid
    heMatrix = mapToValidityCone(heMatrix, dec.properData, 0);
    
    % preset output matrix
    heMatrixOut = nan*heMatrix;

    % find closest point in hE-space that matches a unified partition
    minDistanceVector = inf(1, size(heMatrix, 2));
    thisDistanceVector = minDistanceVector;
    targetUnifiedPartition = nan(1, size(heMatrix, 2));
    for iPart = 1:length(dec.unifiedPartition)
        
        % orthogonal projection to selected subdomain
        orthogonalProjectionOnDomain = dec.unifiedPartition{iPart...
            }.projection.orthogonalOntoSubspace;
        heMatrixProjected = orthogonalProjectionOnDomain * heMatrix;
        % the resulting point must reside in the target domain to be
        %   accepted
        inDomainA = dec.unifiedPartition{iPart}.embeddingDomain.hrep.A;
        filterInDomain = all(inDomainA * heMatrixProjected < zerotol, 1);
        % distance to domain
        thisDistanceVector(filterInDomain) = sum((...
            heMatrixProjected(:,filterInDomain) - ...
            heMatrix(:, filterInDomain)).^2, 1);
        thisDistanceVector(~filterInDomain) = Inf;      
        
        % evaluate current distances
        thisFilter = thisDistanceVector < minDistanceVector;
        minDistanceVector(thisFilter) = thisDistanceVector(thisFilter);
        targetUnifiedPartition(thisFilter) = iPart;
        heMatrixOut(:, thisFilter) = heMatrixProjected(:,thisFilter);
    end
    
    %% second chance computations 
    % From the above, we have projected the heVectors directly onto the
    %   subspace in reach. We did considered the domain in reach in so far
    %   that the projected heVector must reside in that domain.
    % The following loop handles cases in which the projection onto the
    %   subspace in reach is not in the domain so that an additional
    %   projection onto the boundary of the domain is required. As an
    %   example think of an two-dimensional domain in 3D which has such a
    %   boundary.
    for iPart = 1:length(dec.unifiedPartition)
        % orthogonal projection to selected subdomain
        orthogonalProjectionOnDomain = dec.unifiedPartition{iPart...
            }.projection.orthogonalOntoSubspace;
        heMatrixProjected = orthogonalProjectionOnDomain * heMatrix;
        % distance to domain
        thisDistanceVector = sum((heMatrixProjected - heMatrix).^2, 1);
        
        % try to find closest vector in domain when there is a chance to be
        %   better than the current values
        filterSecondTry = thisDistanceVector < minDistanceVector;
        if ~any(filterSecondTry)
            continue
        end
        % find suitable x0
        seIndexVector = dec.unifiedPartition{iPart}.seIndexVector;
        x0 = sum(dec.se(seIndexVector, :), 1)';
        x0 = dec.unifiedPartition{iPart}.projection.validity*...
            dec.properData.mappingNewToProper*x0;
        inDomainA = dec.unifiedPartition{iPart}.embeddingDomain.hrep.A;
        % perform the projection
        heMatrixProjected(:, filterSecondTry) = ...
            closestProjectionOntoConfinementPartition(...
            heMatrix(:, filterSecondTry), ...
            x0, ...
            inDomainA, ...
            orthogonalProjectionOnDomain);
        if any(any(inDomainA * heMatrixProjected(:, filterSecondTry) > zerotol, 1))
            error('elk:decomposition:numericalProblem', ['The he-Vectors ' ...
                'after this step must reside in the given domain. This ' ....
                'should not happen.']);
        end
        
        % recompute distance vector.. ..this time only for the ones that
        %   ARE in the domain.. ..where only vectors can be excluded that
        %   already have a closest point assigned.
        thisDistanceVector(filterSecondTry) = sum((...
            heMatrixProjected(:,filterSecondTry) - ...
            heMatrix(:, filterSecondTry)).^2, 1);
        thisDistanceVector(~filterSecondTry) = Inf;
        
        % evaluate current distances
        thisFilter = thisDistanceVector < minDistanceVector;
        minDistanceVector(thisFilter) = thisDistanceVector(thisFilter);
        targetUnifiedPartition(thisFilter) = iPart;
        heMatrixOut(:, thisFilter) = heMatrixProjected(:,thisFilter);
    end
    
    % the above steps found the closest points:
    %    * by orthogonal projection onto a unified partition presented in
    %        the proper embedding (target domain)
    %    * the points actually REALLY reside in that target domain
    % BUT: it can happen that no such "closest" point exists at all because
    %   they must be projected ONTO some lowe dimensional subspace at once.
    %   (finding the closest point on a polytope surface)
    
    if any(isnan(targetUnifiedPartition))
        % If this fails, this could mean that the decomposition is not
        %   complete.
        error('elk:decomposition:numericalError', ['Could not find an ' ...
            'appropriate target partition for all input hE-vectors. ' ...
            'This should not happen.']);
    end    

    % handle all found target unified partitions subsequently
%     targetUnifiedPartition
%     hcMatrix = nan(dec.nC, size(heMatrix, 2));
%     for iPart = unique(targetUnifiedPartition)
%         thisFilter = targetUnifiedPartition == iPart;
%         if ~any(thisFilter),continue,end
%         % projection onto an nC-dimensional subspace of hE-vectors that are
%         %   described with the target hC-vectors
%         heMatrix(:, thisFilter) = ...
%             dec.unifiedPartition{iPart}.embedding.orthogonalProjection * ...
%             heMatrix(:, thisFilter);
%         
%         % ensure these vectors are valid: if it would not be valid after the
%         %   projection, it might not be confined at the end
% %         heMatrix(:, thisFilter) = mapToValidityCone(heMatrix(:, thisFilter), ...
% %             dec.properData, zerotol);
%         
%         % projection onto the domain of possibly invalid hE-vectors that
%         %   would be covered by mapNewToProper*hC
%         hcMatrix(:, thisFilter) = pinv(...
%             dec.unifiedPartition{iPart}.embedding.validityProjection*...
%             dec.properData.mappingNewToProper)*heMatrix(:, thisFilter);
%     end
%     
    % final projection into hC-space
%     hcMatrix = dec.properData.mappingProperToNew*heMatrix;
end

if nargout < 2
    clear targetUnifiedPartition
end
        