% cd ./cvx 
% cvx_setup
% 
% cd ..

clear all;
% clc;


%% System Parameter
Num_Cell = 7;
L_Macro = 1;
L_Pico = 3;
ISD = 0.8;%inter sector distance in km
Total_BS = L_Macro + L_Pico;
Num_TxAnt_Macro = 1;
Num_TxAnt_Pico = 1;
Num_Rx_Ant = 1;
K = 30; %number of users per cell
BW = 10e6;%bandwidth
P_max_Macro_dBm = 43 - dec2dB(BW); %dBm/Hz
P_max_Pico_dBm = 30 - dec2dB(BW); %dBm/Hz
P_max_total_dBm = dec2dB((dB2dec(P_max_Macro_dBm)*L_Macro + dB2dec(P_max_Pico_dBm)*L_Pico)*Num_Cell); %dBm/Hz
noise_dBm = -150; %dBm/Hz
snr_gap_dB = 9;
intf_2nd = zeros(Num_Rx_Ant,Num_Rx_Ant);    %2nd tier interference

% macro_const = [26.5075   68.3125  118.5453  181.3552];
% pico_const = [7.1892   10.6471   17.1229   29.4793];

macro_const = [4   8  16  24 32];
pico_const = [2   4   8   12 16];

backhaul_index = 2;
Macro_Bkhaul = macro_const(backhaul_index);%bps/Hz
Pico_Bkhaul = pico_const(backhaul_index);

per_BS_bkhaul_constraint = zeros(Num_Cell,Total_BS);
per_BS_power_constraint_matrix = zeros(Num_Cell,Total_BS);
for l = 1:Num_Cell
    for iBS = 1:Total_BS
        if iBS == 1 %Macro BS
            per_BS_power_constraint_matrix(l,iBS) = dB2dec(P_max_Macro_dBm);
            per_BS_bkhaul_constraint(l,iBS) = Macro_Bkhaul;
        else
            per_BS_power_constraint_matrix(l,iBS) = dB2dec(P_max_Pico_dBm);
            per_BS_bkhaul_constraint(l,iBS) = Pico_Bkhaul;
        end
    end
end
per_BS_power_constraint = reshape(per_BS_power_constraint_matrix',[],1);
%per BS power constraint in dBm/Hz

%% Simulation parameter
Power_Mode = 3; 
%1:sum power constraint ; 
%2: per BS power constraint using subgraident method to find optimal dual
%variable
%21: per BS power constraint using gradient descent method to find optimal dual
%variable
%3: two total power constraint, one for Macro BS's, the other for Pico BS's
load Chn.mat;
Cells = struct([]);
User = struct([]);
%SumBackhaulConstraint = 80;

%%
%initialization
for iCell = 1:Num_Cell
    for k = 1:K
        User(iCell,k).long_avg_rate = 0.1;
        User(iCell,k).avg_rate = 0;
    end
    Cells(iCell).UnScheduled_User = 1:K;
end

M = 8; %maximum number of serving BSs per User

[User Asso_User] = Neighbour(Num_Cell, K, path_loss_dB, per_BS_power_constraint_matrix,User,M,noise_dBm);

out_iter = 1;
previous_utility = 1;
current_utility = 2;
Log_Data = struct([]);

Num_Scheduled_User = 4;

%while abs((current_utility - previous_utility)/current_utility) > 1e-2 || out_iter <=70
    % each time schedule 10 users, 5 iterations is needed for all users
    % being served. we wish every user being served at least 10 times
    
    previous_utility = current_utility;
    
    tic;
    
    fprintf('$$$$$$ User Scheduling iter = %d $$$$$$$ \n', out_iter);
    %User scheduler
%     for iCell = 1:Num_Cell
%         [Cells(iCell).Scheduled_User Cells(iCell).UnScheduled_User]= Round_Robin(Cells(iCell).UnScheduled_User,Num_Scheduled_User);
%        
% %         Log_Data(out_iter).Cells(iCell).Scheduled_User = Cells(iCell).Scheduled_User;
% %         Log_Data(out_iter).Cells(iCell).UnScheduled_User = Cells(iCell).UnScheduled_User; % log user scheduling
%         if isempty(Cells(iCell).UnScheduled_User)
%             Cells(iCell).UnScheduled_User = 1:K;
%         end
%   %      Cells(iCell).initial_lambda = zeros(1,Total_BS);
%     end
    
Cells(1).Scheduled_User = [8 11 12 19];
Cells(2).Scheduled_User = [3 7 21 30];
Cells(3).Scheduled_User = [6 12 19 28];
Cells(4).Scheduled_User = [3 6 13 27];
Cells(5).Scheduled_User = [13 16 19 21];
Cells(6).Scheduled_User = [7 9 13 24];
Cells(7).Scheduled_User = [2 4 24 27];

Cells(1).UnScheduled_User = [1     2     3     4     5     6     7     9    10    13    14    15    16    17    18    20    21    22    23    24    25    26    27    28    29    30];
Cells(2).UnScheduled_User = [1     2     4     5     6     8     9    10    11    12    13    14    15    16    17    18    19    20    22    23    24    25    26    27    28    29];
Cells(3).UnScheduled_User = [1     2     3     4     5     7     8     9    10    11    13    14    15    16    17    18    20    21    22    23    24    25    26    27    29    30];
Cells(4).UnScheduled_User = [1     2     4     5     7     8     9    10    11    12    14    15    16    17    18    19    20    21    22    23    24    25    26    28    29    30];
Cells(5).UnScheduled_User = [1     2     3     4     5     6     7     8     9    10    11    12    14    15    17    18    20    22    23    24    25    26    27    28    29    30];
Cells(6).UnScheduled_User = [1     2     3     4     5     6     8    10    11    12    14    15    16    17    18    19    20    21    22    23    25    26    27    28    29    30];
Cells(7).UnScheduled_User = [1     3     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19    20    21    22    23    25    26    28    29    30];
    
    %initial Clustering 
    %every scheduled user connects with all the candidate BS's
    Num_Serving_User = zeros(Num_Cell,Total_BS); %num of users each BS supports at the begining
    for iCell = 1:Num_Cell
        for ik = 1:length(Cells(iCell).Scheduled_User)
            k = Cells(iCell).Scheduled_User(ik);
            User(iCell,k).ServingCluster = User(iCell,k).NeighbourCluster; 
            User(iCell,k).ServingAnt = 0;
            for iBS = 1:length(User(iCell,k).ServingCluster)
                
                curr_cell = User(iCell,k).Neighbour(iBS).Cell_BS(1);
                curr_BS = User(iCell,k).Neighbour(iBS).Cell_BS(2);
                User(iCell,k).bkhaul_wgt(curr_cell,curr_BS) = 0.01;
                Num_Serving_User(curr_cell,curr_BS) = Num_Serving_User(curr_cell,curr_BS) + 1;
                if curr_BS == 1 %Macro BS
                    User(iCell,k).ServingAnt = User(iCell,k).ServingAnt + Num_TxAnt_Macro;
                else
                    User(iCell,k).ServingAnt = User(iCell,k).ServingAnt + Num_TxAnt_Pico;
                end
            end
        end
    end
    
    %Initialize wgts
    max_wgt = 0;
    for iCell = 1:Num_Cell
        for ik = 1:length(Cells(iCell).Scheduled_User)
            k = Cells(iCell).Scheduled_User(ik);
            User(iCell,k).wgt_wsr = 1/User(iCell,k).long_avg_rate;
            User(iCell,k).inst_rate = 0.01;
            User(iCell,k).prev_inst_rate = 0.01;
            if User(iCell,k).wgt_wsr > max_wgt
                max_wgt = User(iCell,k).wgt_wsr;
            end
        end
    end
    %normalize user wgts by max_wgt
    for iCell = 1:Num_Cell
        for ik = 1:length(Cells(iCell).Scheduled_User)
            k = Cells(iCell).Scheduled_User(ik);
            User(iCell,k).wgt_wsr = User(iCell,k).wgt_wsr/max_wgt;
            
        end
    end
    
    % initialize TX beamformers and sparse Matrix
    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);
            
            ant_head = 1;
            temp_tx2 = [];
            for m = 1:length(User(l,k).ServingCluster)
                
                curr_cell = ceil(User(l,k).ServingCluster(m)/Total_BS);
                curr_BS = User(l,k).ServingCluster(m) - (curr_cell-1)*Total_BS;
                if curr_BS == 1
                    Num_Tx_Ant = Num_TxAnt_Macro;
                else
                    Num_Tx_Ant = Num_TxAnt_Pico;
                end
                ant_tail = ant_head + Num_Tx_Ant - 1;
                norm_tx = temp_tx(ant_head:ant_tail)/norm(temp_tx(ant_head:ant_tail));
                temp_tx2 = [temp_tx2;sqrt(per_BS_power_constraint_matrix(curr_cell,curr_BS)/Num_Serving_User(curr_cell,curr_BS))*norm_tx];
                ant_head = ant_tail+1;
                User(l,k).BS_Power(m) = per_BS_power_constraint_matrix(curr_cell,curr_BS)/Num_Serving_User(curr_cell,curr_BS);
            end

            User(l,k).beam_tx = temp_tx2;
            if length(User(l,k).beam_tx) ~= User(l,k).ServingAnt
                fprintf('$$$ Tx Initialization Not Right $$$ \n');
            end
            User(l,k).SparseMatrix = eye(User(l,k).ServingAnt);
        end
    end
    
    %check per BS power constraint
    
 %   [sumpower, per_BS_power, User] = Sum_Power(Num_Cell,Cells,User,L_Macro,L_Pico,Num_TxAnt_Macro,Num_TxAnt_Pico);
    
    [User,sumrate, wgt_sumrate, Per_BS_Bkhaul, Cells] = WSR_Solver_WMMSE(Num_Cell,Cells,User,P_max_total_dBm,Chn,noise_dBm,snr_gap_dB,intf_2nd,...
        L_Macro,L_Pico,Num_TxAnt_Macro,Num_TxAnt_Pico,Num_Rx_Ant,Power_Mode,per_BS_power_constraint,per_BS_power_constraint_matrix,per_BS_bkhaul_constraint); 

    Log_Data(out_iter).Per_BS_Bkhaul = Per_BS_Bkhaul;
    
    [Log_Utility,User,Log_Data] = Compute_Utility(Num_Cell,K,User,Cells,out_iter,Log_Data);

    fprintf('$$$$ Current Log Utility = %f $$$$ \n\n\n\n',Log_Utility);
    
    current_utility = Log_Utility;
    out_iter = out_iter + 1;
    
%    save Log_Data_proposed1.mat Log_Data User;
    
    toc;
    
%end

