function growthRateMatrixNew = limitGrowthRate(growthRateMatrix, ...
    positionMatrix, partitionVector, dec, distanceLimit)
% limit growth rate so that hC-vector remains valid
%
% Syntax: growthRateMatrix = limitGrowthRate(growthRateMatrix, ...
%           positionMatrix, dec, distanceLimit, zerotol)
%
% During growth, hC vectors must be prevented that are outside the valid
%   cone. The first reason is that disappeared facets must retain touching
%   the crystal shape. Simulation results would be invalid when such a
%   facet reappears. The second reason is that measures cannot be
%   calculated for invalid vectors. However, this could be prevented by
%   mapping of the hC-vectors for the task of measure calculation.
%
% The input distanceLimit gives a maximal distance of a facet from its
%   limiting vertex before the growth rate of that facet is adopted. This
%   implies also that the growth rate is limited to zero growth at a size
%   that approximately corresponds to the distanceLimit. Hence,
%   disappearing crystals should be selected based on a pure-zero growth
%   rate if the growth rate was not pure-zero beforehand.
%
% See also: obtainDec, mapToValidCone, solvePopulationBalance

% The distance from a facet validity boundary A{V,i}*hC equals the distance
%   between the facet and the limiting vertex. Thus, selecting a margin
%   value for the facet distances of the cone equals a margin distance of
%   the facets from the limiting vertices.

%% input handling
if ~exist('distanceLimit', 'var') || isempty(distanceLimit)
    distanceLimit = 1e-6;
end
if isempty(growthRateMatrix)
    growthRateMatrixNew = growthRateMatrix;
    return
end

if ~dec.isComplete
    %% INCOMPLETE: no fixing at all
    growthRateMatrixNew = growthRateMatrix;
    return
    
elseif dec.isProper
    %% PROPER: handled in mapToValidityCone
    growthRateMatrixNew = mapToValidityCone(positionMatrix, dec, distanceLimit, ...
        growthRateMatrix);    

elseif ~dec.isProper && size(growthRateMatrix, 1) > size(positionMatrix, 1)
    %% IMPROPER: growth directly im embedding
    % use hE-vectors in proper embedding - we don't need to apply validity
    %   mapping to he-Vectors since the he-vectors are used only to
    %   determine which facets grow out additionally.
    heMatrix = dec.properData.mappingNewToProper * positionMatrix;
    % growth rates are already in proper embedding
    geMatrix = growthRateMatrix;
    % steps follow that are treated per partition:
    for thisPartIndex = unique(partitionVector)
        % filter affected vectors
        thisVectorFilter = partitionVector == thisPartIndex;
        if sum(thisVectorFilter) == 0,continue;end
        % apply validity mapping when growth violates the constraints given
        %   by the unified partition
        constraintIndexVector = dec.unifiedPartition{thisPartIndex...
            }.projection.invalidConstraintIndexVector;
        for iCon = 1:length(constraintIndexVector)
            thisConIndex = constraintIndexVector(iCon);
            applyFilter = thisVectorFilter;
            applyFilter(thisVectorFilter) = ...
                dec.properData.confinementMappingData.A(thisConIndex,:)*...
                geMatrix(:, thisVectorFilter) > 0;
            if sum(applyFilter) == 0,continue;end
            geMatrix(:, applyFilter) = ...
                dec.properData.confinementMappingData.projectionData(...
                dec.unifiedPartition{thisPartIndex}.projection.invalidGroupIndexVector(iCon)...
                ).projectionList{...               
                dec.unifiedPartition{thisPartIndex}.projection.projectionIndexVector(iCon)...
                } * geMatrix(:, applyFilter);
        end
        % otherwise, we ensure that the growth vector is consistent with
        %   the subspace by approximation:
        geMatrix(:, thisVectorFilter) = ...
            dec.unifiedPartition{thisPartIndex}.projection.orthogonalOntoSubspace*...
            geMatrix(:, thisVectorFilter);
        
        % we then apply "normal" growth rate mapping but apply a special
        %   filter for the constraints to be considered
        geMatrix(:, thisVectorFilter) = ...
            mapToValidityCone(heMatrix(:, thisVectorFilter), ...
              dec.properData, distanceLimit, ...
              geMatrix(:, thisVectorFilter), ...
              dec.unifiedPartition{thisPartIndex}.projection.applyConstraintFilter);
    end

    % transform to proper embedding according to known partition
    growthRateMatrixNew = mapFromProperEmbedding(geMatrix, dec, ...
        [], max(partitionVector,1));

else
    error('elk:pbeSolver:wrongInput', ['Please check the dimensions of ' ...
        'the growth rate function output and the position matrix']);
end
end
