function [User,sumrate, wgt_sumrate,Per_BS_Bkhaul, Cells, BS] = WSR_Solver_WMMSE(Cells,User,Chn,Sim,intf_2nd,Power_Mode,BS)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This function solves WSR problem using WMMSE approach with
%%% 1. SumPower constraint in P_max_total_dBm (Power_Mode = 1) or
%%%     per BS power constraint (Power_Mode = 2)
%%% 2. Given set of weights for WSR stored in User.wgt_wsr
%%% 3. Given scheduled users stored in Cells.Scheduled_User
%%% 4. Given clustering scheme for each scheduled user
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% initialize TX beamformers
% for l = 1:Num_Cell
%     for kk = 1:length(Cells(l).Scheduled_User)
%         k = Cells(l).Scheduled_User(kk); %scheduled user in cell l
%         Num_ServingAnt = User(l,k).ServingAnt;
%         temp_tx = randn(Num_ServingAnt,1) + 1i*randn(Num_ServingAnt,1);
%         User(l,k).beam_tx = sqrt(per_BS_power_constraint(2)/10)*temp_tx/norm(temp_tx);
%     end
% end
Total_BS = Sim.Num_BS;
Num_Cell = Sim.Num_Cell;
noise_dB = Sim.noise_dBm;
snr_gap_dB = Sim.snr_gap_dB;
Num_TxAnt_Macro = Sim.Num_Ant_Macro;
Num_TxAnt_Pico = Sim.Num_Ant_Pico;
Num_Rx_Ant = Sim.Num_Rx_Ant;
L_Macro = Sim.Num_Macro;
L_Pico = Sim.Num_Pico;

previous_rate = 1;
current_rate = 2;
max_iteration = 50;
iteration = 1;
max_diff = 1;
bkhaul_flag = 1;
%%
fprintf('\n $$$ Solve WSR via WMMSE $$$ \n')
while abs(current_rate - previous_rate)/current_rate > 0.01 || iteration <=Sim.MinWMMSEIter
    % stop criterion: WSR converges or max_iteration reaches
    
    previous_rate = current_rate;
    fprintf('\n$$$$$$$$$ WMMSE iteration= %d   $$$$$$$$$$ \n', iteration);
    
    fprintf('$$$$$$$$$  Update Rx BF and Weight   $$$$$$$$$$ \n');
    [User] = RXBF_Wgt(Num_Cell,Cells,Num_Rx_Ant,User,Chn,noise_dB,snr_gap_dB,intf_2nd,L_Macro,L_Pico,Num_TxAnt_Pico,Sim);   
    
    fprintf('$$$$$$$$$  Compute CVX parameters   $$$$$$$$$$ \n');
    [User, Total_Ant, Kz, b,Per_BS_Coop,Per_BS_Coop_Bkhaul,Per_BS_ServAnt, Per_BS_CompressionAnt] = ...
        Compute_Kz_UL(Num_Cell,User,Cells,Num_TxAnt_Macro,Num_TxAnt_Pico,L_Macro,L_Pico,Chn,Sim);

    fprintf('$$$$$$$$$  Compute Tx BF using CVX   $$$$$$$$$$ \n');
    if Power_Mode == 3
        [User, max_diff] = CVX_QCQP_TxBF(User,Cells,BS,Num_Cell,Power_Mode,Chn,...
    Total_BS,Total_Ant, Kz, b,Per_BS_Coop,Per_BS_Coop_Bkhaul,Per_BS_ServAnt,Per_BS_CompressionAnt,Sim);
    end
    
    fprintf('$$$$ Tx BF diff = %f   $$$$ \n',max_diff);
    
    %check power constraint
    [~,per_BS_power,User,max_power_diff,BS] = Sum_Power(Num_Cell,Cells,User,L_Macro,L_Pico,Num_TxAnt_Macro,Num_TxAnt_Pico,BS);
    
    
    %fprintf('\n');
    fprintf('$$$ max power diff = %f $$$\n',max_power_diff);
    % sumpower = Sum_Power(Num_Cell,Cells,User,L_Macro,L_Pico,Num_TxAnt_Macro,Num_TxAnt_Pico);
    
    [User, sumrate, wgt_sumrate, Per_BS_Bkhaul, curr_user_num, Cells, max_rate_diff] = ...
        rate_DL(Num_Cell,Num_Rx_Ant,User,Cells,Chn,noise_dB,intf_2nd,L_Macro,L_Pico,Num_TxAnt_Pico,Sim);
    
    fprintf('$$$ max rate diff = %f $$$\n',max_rate_diff);
%     bkhaul_gap = (Per_BS_Bkhaul - per_BS_bkhaul_constraint)./per_BS_bkhaul_constraint;
%     max_bkhaul_gap = max(bkhaul_gap);
    
    max_power_gap = -inf;
    max_bkhaul_gap = -inf;
    for l=1:Num_Cell
        for b=1:Total_BS
            power_gap = (per_BS_power(l,b) - BS(l,b).total_power_budget)/BS(l,b).total_power_budget;
            bkhaul_gap = (Per_BS_Bkhaul(l,b) - BS(l,b).total_bkhaul_budget)/BS(l,b).total_bkhaul_budget;
            if power_gap > max_power_gap
                max_power_gap = power_gap;
            end
            
            if bkhaul_gap > max_bkhaul_gap
                max_bkhaul_gap = bkhaul_gap;
            end
        end
    end
    fprintf('$$$ max power gap = %f $$$ \n',max_power_gap);    
    fprintf('$$$ max bkhaul gap = %f $$$ \n',max_bkhaul_gap);   
    
    if max(max_bkhaul_gap) < 1
        bkhaul_flag = 0;
    else
        bkhaul_flag = 1;
    end
    
    current_rate = wgt_sumrate;
    
    fprintf('$$$ Current Survived User = %d  $$$ \n',curr_user_num);
    fprintf('$$$$$$$$$  Updated Wgted Sum Rate= %f   $$$$$$$$$$ \n', current_rate);
       
    iteration = iteration + 1;
    
    if iteration > max_iteration
        break;
    end
    
%     if mod(iteration,10) == 0
%         iteration
%     end
   
end

