function [MGamma,newgamma,newGamma] = sHvectorize_indep_map(Gamma,gamma,tol_licols)
%function [MGamma,newgamma,newGamma] = sHvectorize_indep_map(Gamma,gamma,tol_licols)

%% vectorize (sHvec) the map Gamma(X)=gamma and sort out linearly dependent constraints
% Note: A,B are hermitian.
% Then <A,B>=trace(AB)=sum(A^T(:).*B(:))=Hvec(A)'*Hvec(B)
% Input: Gamma : cell array holding the data for the linear map
%        gamma : a vector holding the RHS of the constraints
% Output: MGamma : A matrix; vectorized map MGamma*Hvec(X) = b of Gamma(X)=gamma
%         newgamma : the RHS gamma after removing redundant constraints
%         newGamma : A cell array; Gamma after removing redundant constraints

verbose = false;

n2 = size(Gamma{1},1)^2;
MGamma = zeros(length(Gamma),n2); % each row is a constraint
for ii=1:length(gamma)
    MGamma(ii,:) = Hvec(Gamma{ii})';
end

if n2 ==1   % in this case, we may choose the first Gamma for the entire constraints 
    [~,indepid] = max(gamma);  
    newGamma = Gamma(indepid); % throw away the redundant constraints
    newgamma = gamma(indepid); % throw away the redundant constraints
    MGamma = MGamma(indepid,:);  % throw away the redundant constraints
    condMG = sqrt(condest(MGamma*MGamma'));
    if verbose
        fprintf('     in sHvectorize: licols removed: %i out of %i constraints; condest MGamma is %g \n',...
            length(gamma)-length(indepid),length(gamma),condMG)
    end
else
    [~,indepid] = licols(transpose(MGamma),tol_licols);  % pick the lin.indep constraints
    MGamma = MGamma(indepid,:);  % throw away the redundant constraints
    condMG = sqrt(condest(MGamma*MGamma'));
    if verbose
        fprintf('     in sHvectorize: licols removed: %i out of %i constraints; condest MGamma is %g \n',...
            length(gamma)-length(indepid),length(gamma),condMG)
    end
    newGamma = Gamma(indepid); % throw away the redundant constraints
    newgamma = gamma(indepid); % throw away the redundant constraints
    
end
end %of function
