%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Experiments with OFDMA Resource Allocation %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning('off','MATLAB:divideByZero')
CSI_Type = 0; % 0 - Perfect, 1 - Predicted
PLOT_CHAN = 0;
% close all
M = 2;                                  % No. of users
% Nfft = 2;                             % No. of subcarriers
% Nfftvec = [4,16,64,256];                  % No. of subcarriers
Nfftvec = 128;
Nu = 76;
CPfrac = 6/128;                          % Fraction of Nfft for CP
colorvec = ['m','b','r','g'];

Nchan = 10000;
BW = 1.25e6;                               % Bandwidth (Hz)
fc = 2.6e9;                             % Carrier frequency (Hz)
[tau,gnorm] = Pdp('vehicularA');        % Power Delay Profile
c = 3e8;                                % Speed of light (m/s)
Velvec = [30, 75, 120, 180, 250]*1000/3600;
v = Velvec(5);                          % Velocity (m/s)
fd = v/c*fc;                            % Doppler Frequency (Hz)
tcoh = 1/(2*fd);                        % Coherence time (s)
lambda = 3e8/fc;                        % Wavelength (m)
fs = 1.92e6;           % Sampling Frequency (Hz)


% Prediction params
Lpred = 4;
Dt = 1;

SNRvec = 5:5:15;

% SNR = 5;    % dB
% sigma = 10^(-SNR/20); % standard deviation
BER = 1e-3;
R = 0:2:6;
Eta = -log(5*BER)/1.6*(2.^R-1);
% Generate channel
wvec = linspace(0.1,0.9,11);
% wvec = 0;
Nfft = Nfftvec(1);
tsym = Nfft*(1+CPfrac)/fs;


  % [Hhat,sigma2ehat] = channelpred(M, 1, fs, tsym, Nu, tau, gnorm, fd, Lpred, Dt, CPfrac, sigma^2);

  C0 = zeros(M,length(wvec),length(SNRvec));
  P0 = zeros(M,length(wvec),length(SNRvec));
  G0 = zeros(1,length(wvec),length(SNRvec));
  I0 = G0;
  L0 = G0;

  C = zeros(M,length(wvec),length(SNRvec));
  P = zeros(M,length(wvec),length(SNRvec));
  G = zeros(1,length(wvec),length(SNRvec));
  I = G;
  L = G;

  C2 = zeros(M,length(wvec),length(SNRvec));
  G2 = zeros(1,length(wvec),length(SNRvec));
  I2 = G;
  L2 = G;
  C3 = zeros(M,length(wvec),length(SNRvec));
  G3 = zeros(1,length(wvec),length(SNRvec));
  I3 = G;
  L3 = G;
  C4 = zeros(M,length(wvec),length(SNRvec));
  G4 = zeros(1,length(wvec),length(SNRvec));
  I4 = G;
  L4 = G;
  C5 = zeros(M,length(wvec),length(SNRvec));
  G5 = zeros(1,length(wvec),length(SNRvec));
  I5 = G;
  L5 = G;
  C6 = zeros(M,length(wvec),length(SNRvec));
  G6 = zeros(1,length(wvec),length(SNRvec));
  I6 = G;
  L6 = G;
for s = 1 : length(SNRvec)
  SNR = SNRvec(s);    % dB
  sigma = 10^(-SNR/20); % standard deviation

  H = channel(M, 1000, fs, Nu, Nfft, tau, gnorm);
  Gamma = abs(H).^2/(sigma^2);
  for i = 1:length(wvec)
    w = [wvec(i),1-wvec(i)]
    %     [Gtmp0, Itmp0, Ltmp0, Ctmp0,
    %     Ptmp0]=ExpectedWeightedPerfectCSI2(w,[1,1]/(sigma^2),Nu,Nu);
    [Gtmp0, Itmp0, Ltmp0, Ctmp0, Ptmp0] = ExpectedDiscretePerfectCSI(Gamma,Gamma,w,Nu,R,Eta,0); 
    C0(:,i,s) = C0(:,i,s) + Ctmp0'/Nu;
    P0(:,i,s) = P0(:,i,s) + Ptmp0'/Nu;
%     G0(i) = G0(i) + Gtmp0;
%     I0(i) = I0(i) + Itmp0;
    L0(1,i,s) = L0(1,i,s) + Ltmp0;
  end;
  for c = 1 : Nchan
    c
    H = channel(M, 1, fs, Nu, Nfft, tau, gnorm);
    %   E = (randn(size(H))/sqrt(2)+sqrt(-1)*randn(size(H))/sqrt(2)).*repmat(sigma2ehat,1,M);
    Gamma = abs(H).^2/(sigma^2);
    %   Gammahat = abs(H-E).^2./(sigma^2);
    for i = 1:length(wvec)
      w = [wvec(i),1-wvec(i)];
%       [Gtmp, Itmp, Ltmp, Ctmp,P]=ExpectedWeightedPerfectCSI(Gamma,Gamma,w,Nu,L0(i));
%       [Gtmp2, Itmp2, Ltmp2, Ctmp2,P2]=ExpectedWeightedPerfectCSI(Gamma,Gamma,w,Nu,0);
      %       [Gtmp3, Itmp3, Ltmp3, Ctmp3, P3]=ExpectedWeightedImperfectCSI(Gammahat,Gamma,sigma^2,sigma2ehat,w,Nu);
      %       [Gtmp3, Itmp3, Ltmp3, Ctmp3, P3]=ExpectedWeightedPerfectCSI(Gammahat,Gamma,w,1);
      [Gtmp, Itmp, Ltmp, Ctmp,Ptmp] = ExpectedDiscretePerfectCSI(Gamma,Gamma,w,Nu,R,Eta,L0(1,i,s));
      [Gtmp2, Itmp2, Ltmp2, Ctmp2,P2] = ExpectedDiscretePerfectCSI(Gamma,Gamma,w,Nu,R,Eta,0);
     %       [Gtmp5, Itmp5, Ltmp5, Ctmp5, P5] = ExpectedDiscreteImperfectCSI(Gammahat,Gamma,sigma^2,sigma2ehat,w,Nu,R,Eta,BER);
      %       [Gtmp6, Itmp6, Ltmp6, Ctmp6, P6] = ExpectedDiscretePerfectCSI(Gammahat,Gamma,w,Nu,R,Eta);
      %       f = sum(w.*Ctmp);


      C(:,i,s) = C(:,i,s) + Ctmp'/Nu/Nchan;
      P(:,i,s) = P(:,i,s) + Ptmp'/Nu/Nchan;
      G(1,i,s) = G(1,i,s) + Gtmp/Nchan;
      I(1,i,s) = I(1,i,s) + Itmp/Nchan;
      L(1,i,s) = L(1,i,s) + Ltmp/Nchan;

      C2(:,i,s) = C2(:,i,s) + Ctmp2'/Nu/Nchan;
      G2(1,i,s) = G2(1,i,s) + Gtmp2/Nchan;
      I2(1,i,s) = I2(1,i,s) + Itmp2/Nchan;
      L2(1,i,s) = L2(1,i,s) + Ltmp2/Nchan;

      
      C3(:,i,s) = C2(:,i,s) + sum(Ctmp3)'/Nu/Nchan;
      
      %       f3 = sum(w.*Ctmp3);
      %       C3(:,i) = C3(:,i) + Ctmp3'/Nu/Nchan;
      %       G3(i) = G3(i) + Gtmp3/Nchan;
      %       I3(i) = I3(i) + Itmp3/Nchan;
      %       L3(i) = L3(i) + Ltmp3/Nchan;

      %       f4 = sum(w.*Ctmp4);
      %       C4(:,i) = C4(:,i) + Ctmp4'/Nchan;
      %       G4(i) = G4(i) + Gtmp4/Nchan;
      %       I4(i) = I4(i) + Itmp4/Nchan;
      %       L4(i) = L4(i) + Ltmp4/Nchan;
      %
      %       f5 = sum(w.*Ctmp5);
      %       C5(:,i) = C5(:,i) + Ctmp5'/Nchan;
      %       G5(i) = G5(i) + Gtmp5/Nchan;
      %       I5(i) = I5(i) + Itmp5/Nchan;
      %       L5(i) = L5(i) + Ltmp5/Nchan;
      %
      %       f6 = sum(w.*Ctmp5);
      %       C6(:,i) = C6(:,i) + Ctmp6'/Nchan;
      %       G6(i) = G6(i) + Gtmp6/Nchan;
      %       I6(i) = I6(i) + Itmp6/Nchan;
      %       L6(i) = L6(i) + Ltmp6/Nchan;
    end;
  end;
  figure(1)
  plot(C0(1,:,s),C0(2,:,s),'b-',C(1,:,s),C(2,:,s),'r-x',C2(1,:,s),C2(2,:,s),'go--')
  hold on
  grid on
  title('Capacity Region for 2 Users')
  xlabel('Capacity for user 1 (bps/Hz)')
  ylabel('Capacity for user 2 (bps/Hz)')
  legend('T-F WF: Theoretical','T-F WF: Simulation','F WF: Simulation')
end;

% save DiscRatePerfCSI
% figure(1)
% plot(wvec,abs(P(1,:)),'bo-')
% hold on
% plot(wvec,abs(P(2,:)),'rs-')
% plot(wvec,abs(sum(P)),'kx-')
% grid on
% legend('Power - User 1','Power - User 2','Sum Power')
% xlabel('w_1')
% ylabel('Power')
% figure(2)
% plot(C(1,:),C(2,:))
% grid on
% hold on
% title('Capacity Region for 2 Users')
% xlabel('Capacity for user 1 (bps/Hz)')
% ylabel('Capacity for user 2 (bps/Hz)')


% end;
% figure(1)
% plot(C(1,:),C(2,:),'-',C2(1,:),C2(2,:),'--')%,'--',C3(1,:),C3(2,:),'-.')
% % plot(C4(1,:),C4(2,:),'-',C5(1,:),C5(2,:),'--',C6(1,:),C6(2,:),'-.')
% grid on
% hold on
% title('Capacity Region for 2 Users')
% xlabel('Capacity for user 1 (bps/Hz)')
% ylabel('Capacity for user 2 (bps/Hz)')
% legend('P-CSI: WF','I-CSI: Opt','I-CSI: WF')
% figure(2)
% semilogy(wvec,G4)
% grid on
% hold on
% figure(3)
% plot(wvec,I)
% hold on
% grid on
% figure(4)
% plot(wvec,L)
% hold on
% grid on
% % end;
% figure(5)
% plot(wvec,G)
% grid on
% title('Duality Gap')
% xlabel('User 1 Weight')
% figure(6)
% plot(wvec,10*log10(I))
% title('No. of Iterations')

