function rate = user_selection_ZF_greedy(H, QuantH, P, SelectionMethod)
% rate = user_selection_ZF_greedy(H, QuantH, P, SelectionMethod)
%
% Computes Zero-forcing rate for channels H, quantized channels QuantH and SNR P.
%
% H and QuantH must be of dimensions (Number of users) X (Number of base station antennas)
% SelectionMethod can be 0 (Greedy) or 1 (Norm)
%
%   Part of mud_vs_csi package.
%
%   Reference:  "Multi-User Diversity vs. Accurate Channel State Information in MIMO Downlink Channels"
%                       by Niranjay Ravindran and Nihar Jindal, Submitted to IEEE Transactions on Communications, 2009.
%

% Dependencies: compute_ZF_rate.m

Nt = size(H, 2);
K = size(H, 1);

CQI = abs(diag(QuantH*QuantH'));

if SelectionMethod == 0 %Greedy selection using QuantH
    % First select best single user
    [max_CQI selected_users] = max(CQI);
    rate = log2(1 + P * max_CQI);

    flag = 1;
    while (flag)
        temp_rate = zeros(K,1);
        for k = 1:K
            if (sum(k == selected_users) == 0)
                % Compute rate if k-th user added
                if (rank(QuantH([selected_users k], :))==numel(selected_users)+1)
                    temp_rate(k) = compute_ZF_rate(QuantH([selected_users k],:), P);
                end
            end
        end
        if (max(temp_rate) > rate)
            % If better than previous rate, add user
            [new_rate new_user] = max(temp_rate);
            selected_users = [selected_users new_user];
            rate = new_rate;
        else
            flag = 0;
        end;
        if (numel(selected_users) >= Nt)
            flag = 0;
        end
    end

else
    % Select users with largest CQI values
    [max_CQI, user_order]= sort(CQI, 'descend');
    selected_users = user_order(1);
    rate = log2(1+ P*max_CQI(1));
    for k = 2:Nt
        temp_rate = compute_ZF_rate(QuantH([selected_users user_order(k)],:), P);
        if (temp_rate > rate)
            rate = temp_rate;
            selected_users = [selected_users user_order(k)];
        else
            break;
        end
    end
end

num_users = numel(selected_users);

QuantH_ = QuantH(selected_users, :);
H_ = H(selected_users, :);

%Compute ZF beamformers
V = (QuantH_' * inv(QuantH_ * QuantH_'));
for k = 1:num_users
    V(:,k) = sqrt(P/num_users) * V(:,k)/norm(V(:,k)); %Normalize BF vectors with equal power allocation (waterfilling is optimal)
end;

%Compute sum rate using H for computed BF vectors
gain_matrix = abs(H_ * V).^2;
temp = 1+ sum(gain_matrix,2);
temp = diag(temp) * ones(num_users);
int_power = temp - gain_matrix;
SINR_matrix = gain_matrix ./ int_power;
rate = sum(log2(1 + diag(SINR_matrix)));

%__EOF__