function [B,K,I] = hard_threshold(B,K,q,I)
% Hard-thresholding operator H_K(B,q) used in the compressive regularized
% (linear) discriminant analysis (CRDA). 
% 
% [B,K,I] = hard_threshold(B,K,q)
% [B,K,I] = hard_threshold(B,K,q,I)
% 
% Input
%   B       p x G coefficient matrix used in CRDA
%   K       integer in the interval [1,...,p] or empty set [] 
%           sparsity order used in the hardth
%           if emptyset is given then the function determines the sparsity 
%           order as K = sum_i=1^p I(l_i > mean(l)), where l_i are the row
%           lengths of B. The used length is determined by input q. 
%  q        scalar 0, inf or >= 1
%           If q = 0, then use sample variance as the row lenght of B 
%           if q = inf, then use L_Infinity-norm as the row length of B
%           If q >= 1, then use L_q-norm as the row length of B
%
% See also CRDACV and CRDA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Kmean = false; 
p = size(B,1);

if isempty(K)
    Kmean  = true;
end

if nargin < 4
    I = [];
end

if ~(( q >= 1 && q <= p) || q==0 || q == inf)
    error(sprintf('q  must value between 1 and %d or 0 or ''Inf''',p));
end

print_info = false;

switch q
    case 0
       len = var(B,[],2);   % Variance (rows as vectors)
    case inf
       len = max(abs(B), [],2); % Linf-norm (rows as vectors)
    otherwise
       len = sum(abs(B).^q, 2).^(1/q);	% Lq-norm (rows as vectors), q>=1
end
    
if print_info
    if q >= 1
       fprintf('CRDA: using L_{%.2f}-norm based hard thresholding\n',q);
    elseif q==0
       fprintf('CRDA: using sample variance based hard thresholding\n',q);
    else
       fprintf('CRDA: using sample variance based hard thresholding\n',q);
    end 
end

if Kmean
    K = sum(len>mean(len));
end
    
%% Hard threshold:
% keep the rows with  largest L_q-row norms and set others to zero  
if isempty(I)
    [~,I] = sort(len,'descend'); 
end
 
% Setting (p-K)-features of the coefficient matrix to zeros 
B(setdiff(1:p,I(1:K)),:) = 0;      
    
