function rate = compute_pu2rc_rate(H, BFCodebook, P)
% rate = compute_pu2rc_rate(H, BFCodebook, P)
%
% Computes PU2RC rate for channels H and SNR P.
%
% H must be of dimension (Number of users) X (Number of base station antennas).
% BFCodebook is the codebook used and must be of dimension (Number of base station antennas) X (Number of transmit beams) X (Number of orthogonal sets),
% where the columns of BFCodebook(:, :, n) (the codebook for the n-th orthogonal set) consist of mutually
% orthogonal (Number of base station antennas) x 1 vectors with norm unity.
%
% If there is only 1 orthogonal set and BFCodebook is randomly generated, this scheme corresponds to RBF.
%
%   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.
%

num_users = size(H, 1);
NumBeams = size(BFCodebook, 2);
NumOrthSets = size(BFCodebook, 3);

Feedback_CQI = zeros(num_users, 1);
Feedback_OrthSet = NaN(num_users, 1);
Feedback_Beam = NaN(num_users, 1);

%Decide feedback information for each user
for k = 1:num_users
    channel_gain = norm(H(k, :));
    current_channel_dir = H(k, :)/channel_gain;
	%Find the "best beam" from among all orthogonal sets and feed back the corresponding SINR
    for OrthSet = 1:NumOrthSets
		[cos_sq, best_beam] = max(abs(current_channel_dir*BFCodebook(:, :, OrthSet)).^2);
        best_sinr = channel_gain^2*cos_sq/(NumBeams/P + channel_gain^2*sum(abs(current_channel_dir*BFCodebook(:, [1:best_beam-1 best_beam+1:NumBeams], OrthSet)).^2));
        if best_sinr > Feedback_CQI(k)
            Feedback_CQI(k) = best_sinr; %Feedback SINR
			%The following feedback information identifies an orthogonal set and a beam whithin that set
			%and requires log2(NumOrthSets)+log2(NumBeams)  bits for representation
            Feedback_Beam(k) = best_beam;
            Feedback_OrthSet(k) = OrthSet;
        end
    end
end

%Select users for downlink transmission and compute corresponding sum rate
rate = 0;
%Select "best" users for each orthogonal set, and decide on the set that has the highest sum rate
for OrthSet = 1:NumOrthSets
	%Begin computing sum rate for this orthogonal set
    sinr_values = NaN(NumBeams, 1);
    selected_user_indices = NaN(NumBeams, 1);
	%Select a user for each beam
    for m = 1:NumBeams
        users_on_this_beam = find(Feedback_OrthSet==OrthSet&Feedback_Beam==m);
        if ~isempty(users_on_this_beam)
            [best_sinr best_user_index] = max(Feedback_CQI(users_on_this_beam));
            sinr_values(m) = best_sinr;
            selected_user_indices(m) = users_on_this_beam(best_user_index);
        end
    end    

    if sum(isnan(sinr_values)) >= 1
		%If fewer users than than beams were selected, compute new SINR values and then compute the rate
        BFCodebook_Reduced = BFCodebook(:, :, OrthSet);
        for m = 1:NumBeams
            if isnan(sinr_values(m))
                BFCodebook_Reduced(:, m) = 0;
            end
        end
        rate_on_this_orth_set = 0;
        for m = 1:NumBeams
            if ~isnan(sinr_values(m))
                temp = abs(H(selected_user_indices(m), :)*BFCodebook_Reduced).^2;
                temp = (temp(m) / (sum(~isnan(sinr_values))/P + sum(temp([1:(m-1) (m+1):NumBeams]))));
                rate_on_this_orth_set = rate_on_this_orth_set + log2(1+temp);
            end
        end
    else
		%We already know the SINR values from the information that was fed back, so compute the rate
        rate_on_this_orth_set = sum(log2(1 + sinr_values));
    end 
    
    if rate_on_this_orth_set > rate
        rate = rate_on_this_orth_set;
        break;
    end
end

%__EOF__