function [max_rate, B_opt, B_vector, rate_vector] = mud_vs_csi(T, Nt, SNR_dB, BFType, QuantizType, CQIType, SelectionMethod, MaxIter, BStart, BEnd)
% [rate_vector, B_vector, max_rate, B_opt] = mud_vs_csi(T, Nt, SNR_dB, BFType, QuantizType, CQIType, SelectionMethod, MaxIter, BStart, BEnd)
%
% Input Arguments:
%       T                         - Total feedback budget in bits
%       Nt                        - Number of base station antennas
%       SNR_dB             - The system SNR in dB
%       BFType               - The type of beamforming, can be 'ZF', 'PU2RC' or 'SUBF'
%       QuantizType       - The type of channel quantization used, can be 'None', 'RVQ', 'QUB', or 'Scalar' (ignored by PU2RC)
%                                     For perfect CSI, use 'None'
%       CQIType             - The type of channel quality indicator used, can be 'Norm' or 'SINR' (ignored for PU2RC and SUBF)
%       SelectionMethod - The type of user selection employed can be 'Greedy' or 'Norm' (ignored by PU2RC and SUBF)
%       MaxIter                - Optional. Number of Monte-carlo iterations (defaults to 100)
%       BStart/BEnd        - Optional. Range of values of B to simulate, omit for automatic calculation (ignored if QuantizType is 'None')
%
% Output Variables:
%       max_rate            - The maximum optimized sum rate (bps/hz)
%       B_opt                  - The feedback load per user that maximizes the sum rate (bits)
%       B_vector             - The values of feedback load per user corresponding to the rates in rate_vector (bits)
%       rate_vector          - A vector with sum rate values for each per-user feedback load in B_vector (bps/hz)
%
%
%   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: quantiz_channels.m, user_selection_ZF_greedy.m, compute_pu2rc_rate.m

%Check input arguments
if strcmpi(BFType, 'ZF')
    BFType = 0; %Zero-forcing
elseif strcmpi(BFType, 'PU2RC');
    BFType = 1; %RBF/PU2RC
else
    BFType = 2; %Single-user beamforming
end

if strcmpi(QuantizType, 'None')
    QuantizType = 0;    %Perfect CSI
elseif strcmpi(QuantizType, 'RVQ');
    QuantizType = 1; %Random vector quantization
elseif strcmpi(QuantizType, 'QUB');
    QuantizType = 2;    %Quantization upper bound
else
    QuantizType = 3;    %Scalar quantization
end

if strcmpi(CQIType, 'Norm')
    CQIType = 0;
else
    CQIType = 1;    %SINR feedback
end

if strcmpi(SelectionMethod, 'Greedy')
    SelectionMethod = 0;
else
    SelectionMethod = 1;    %Selection of users with largest channel norms
end
    
if nargin < 8
    MaxIter = 100; %Use a larger number for smoother curves
end
    
%Compute range of B values to simulate (automatically, if not specified)
if BFType == 0
    if nargin < 10
        K_vector = max(ceil(T/50), Nt):1:min(150, T);
    else
        K_vector = unique(round(T./(BStart:BEnd)));
    end
    B_vector = T./K_vector;
    [dummy, indices] = unique(floor(B_vector));
    K_vector = K_vector(indices);
    B_vector = B_vector(indices);
elseif BFType == 1
    if nargin < 10
        B_vector = [log2(Nt):1:7+log2(Nt)];
    else
        B_vector = [BStart:1:BEnd];
    end
    K_vector = ceil(T./B_vector);
else
    if nargin < 10
        K_vector = max(1, ceil(T/50)):1:min(200, T);
    else
        K_vector = unique(round(T./(BStart:BEnd)));
    end
    B_vector = T./K_vector;
    [dummy, indices] = unique(floor(B_vector));
    K_vector = K_vector(indices);
    B_vector = B_vector(indices);
end

P = 10^(SNR_dB/10);

rate_vector = zeros(1, length(B_vector));

%Display friendly status information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['Feedback budget = ' num2str(T) ' bits, SNR = ' num2str(SNR_dB) ' dB, BS Antennas = ' num2str(Nt,'%02d') ', Beamforming type = ' num2str(BFType)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for B_its = 1:length(B_vector);

    B = B_vector(B_its);
    K = ceil(K_vector(B_its));
    regions = 2^B;

    if (QuantizType ==3 && BFType ~= 1 && B<2*2*(Nt-1) )
        continue; %Too small a codebook for scalar quantization
    end

    rate_temp = 0;
    for its = 1:MaxIter

        H = 1/sqrt(2) * (randn(K, Nt) + i * randn(K, Nt)); %Downlink channels

        if BFType == 0  %ZF
            QuantH = quantiz_channels(H, regions, P, QuantizType, CQIType);
            rate_temp = rate_temp + user_selection_ZF_greedy(H, QuantH, P, SelectionMethod);
        elseif BFType == 1  %PU2RC
            NumBeams = Nt;
            NumOrthSets = ceil(2^(log2(regions)-log2(NumBeams)));
            BFCodebook = zeros(Nt, NumBeams, NumOrthSets);
            for OrthSet = 1:NumOrthSets
                BFCodebook(:, :, OrthSet) = orth(randn(Nt, NumBeams) + i*randn(Nt, NumBeams));
            end
            rate_temp = rate_temp + compute_pu2rc_rate(H, BFCodebook, P);
        else    %SUBF
            QuantH = quantiz_channels(H, regions, P, QuantizType, 0);
            temp_rate = 0;
            for k = 1:K
                SNR_eff = abs(H(k, :)*QuantH(k, :)')^2/norm(QuantH(k, :))^2;
                if log2(1 + P*SNR_eff) > temp_rate
                    temp_rate= log2(1 + P*SNR_eff);
                end
            end
            rate_temp = rate_temp + temp_rate;
        end

    end; %for MaxIter

    if QuantizType == 0 && BFType ~= 1
        rate_vector = rate_temp/MaxIter;
        B_vector = inf;
        break;  %Quit of perfect CSI simulation, no need to cycle through B_vector
    else
        rate_vector(B_its) = rate_temp/MaxIter;
        
        %Display friendly status information
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        disp(['   Feedback load per user = ' num2str(floor(B),'%03d') ' bits, Users = ' num2str(K, '%03d') ', Sum Rate = ' num2str(rate_vector(B_its)) ' bps/hz']);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    end

end

%Compute optimal B and throughput
[max_rate B_opt] = max(rate_vector);
B_opt = B_vector(B_opt)

%__EOF__