clc
clear all
MC_run_1 = 1e4; %Monte Carlo Run for Fig 4,5
MC_run_2 = 2e4; %Monte Carlo Run for Fig 7
% Monte Carlo iterations will cost long time, especially for the CSI 
% scnario. Overall runtime may reach to a hour. You could play in a lower 
% iteration number. You may also skip to the graph section to see the 
% result with the 1e4 and 2e4 independent runs, which needs to uncomment 
% line 37, 38, 77
%% Fig 4 AWGN and Fig 5 Fading channel with full SI
% Waterfilling Solution

I_res_awgn(1,:) = FSI_SNR(1,MC_run_1,1);
I_res_awgn(2,:) = FSI_SNR(2,MC_run_1,1);
I_res_awgn(3,:) = FSI_SNR(4,MC_run_1,1);
I_res_ray(1,:) = FSI_SNR(1,MC_run_1,2);
I_res_ray(2,:) = FSI_SNR(2,MC_run_1,2);
I_res_ray(3,:) = FSI_SNR(4,MC_run_1,2);
save('Thrpt_FSI','I_res_awgn','I_res_ray');
fprintf('Finished: Fig 4 AWGN and Fig 5 Fading channel with full SI')

%% Fig 4 AWGN and Fig 5 Fading channel with Causal SI
% Bellman Equation
Thrpt_AWGN_CSI(1,:)=CSI_AWGN_SNR(1,MC_run_1);
Thrpt_AWGN_CSI(2,:)=CSI_AWGN_SNR(2,MC_run_1);
Thrpt_AWGN_CSI(3,:)=CSI_AWGN_SNR(4,MC_run_1);
Thrpt_Ray_CSI(1,:)=CSI_Ray_SNR(1,MC_run_1);
Thrpt_Ray_CSI(2,:)=CSI_Ray_SNR(2,MC_run_1);
Thrpt_Ray_CSI(3,:)=CSI_Ray_SNR(4,MC_run_1);

save('Thrpt_CSI','Thrpt_AWGN_CSI','Thrpt_Ray_CSI');
fprintf('Finished: Fig 4 AWGN and Fig 5 Fading channel with Causal SI')

%% GRAPH: Fig 4 and Fig 5
load('Thrpt_FSI','I_res_awgn','I_res_ray');
load('Thrpt_CSI','Thrpt_AWGN_CSI','Thrpt_Ray_CSI');
% load('Thrpt_FSI_1e4','I_res_awgn','I_res_ray');
% load('Thrpt_CSI_1e4','Thrpt_AWGN_CSI','Thrpt_Ray_CSI');

figure(1)
hold on
grid on
plot(0:5:20,I_res_awgn(1,:),'-or','LineWidth',2)
plot(0:5:20,I_res_awgn(2,:),'-og','LineWidth',2)
plot(0:5:20,I_res_awgn(3,:),'-ob','LineWidth',2)
plot(0:5:20,Thrpt_AWGN_CSI(1,:),'-xr','LineWidth',2)
plot(0:5:20,Thrpt_AWGN_CSI(2,:),'-xg','LineWidth',2)
plot(0:5:20,Thrpt_AWGN_CSI(3,:),'-xb','LineWidth',2)
legend('Full SI, K=1','Full SI, K=2','Full SI, K=4','Causal SI, K=1','Causal SI, K=2','Causal SI, K=4')
xlabel('SNR(dB)')
ylabel('Throughput pre slot')
title('AWGN channel (Fig. 4)')

figure(2)
hold on
grid on
plot(0:5:20,I_res_ray(1,:),'-or','LineWidth',2)
plot(0:5:20,I_res_ray(2,:),'-og','LineWidth',2)
plot(0:5:20,I_res_ray(3,:),'-ob','LineWidth',2)
plot(0:5:20,Thrpt_Ray_CSI(1,:),'-xr','LineWidth',2)
plot(0:5:20,Thrpt_Ray_CSI(2,:),'-xg','LineWidth',2)
plot(0:5:20,Thrpt_Ray_CSI(3,:),'-xb','LineWidth',2)
legend('Full SI, K=1','Full SI, K=2','Full SI, K=4','Causal SI, K=1','Causal SI, K=2','Causal SI, K=4')
xlabel('SNR(dB)')
ylabel('Throughput pre slot')
title('Fading channel (Fig. 5)')

%% Fig 6 Throughput comparision with optimal throughput and power-halving scheme
I_ph_Rayleigh = TS_ph(MC_run_2,2);
I_ph_AWGN = TS_ph(MC_run_2,1);
I_wf_AWGN = FSI_TS(MC_run_2,1);
I_wf_Rayleigh = FSI_TS(MC_run_2,2);
save('Comparision.mat','I_ph_Rayleigh','I_ph_AWGN','I_wf_AWGN','I_wf_Rayleigh');
fprintf('Finished: Fig 6 Throughput comparision with optimal throughput and power-halving scheme')
%% GRAPH: Fig 6
load('Comparision.mat','I_ph_Rayleigh','I_ph_AWGN','I_wf_AWGN','I_wf_Rayleigh');
%load('Comparision_2e4.mat','I_ph_Rayleigh','I_ph_AWGN','I_wf_AWGN','I_wf_Rayleigh');

figure(3)
hold on
grid on
plot(2.^(0:7),I_ph_Rayleigh,'-or','LineWidth',2)
plot(2.^(0:7),I_wf_Rayleigh,'-xb','LineWidth',2)
plot(2.^(0:7),I_ph_AWGN,'-.or','LineWidth',2)
plot(2.^(0:7),I_wf_AWGN,'-.xb','LineWidth',2)
legend('Fading Channel, full SI available','Fading Channel, power halving scheme','AWGN, full SI available','AWGN, power halving scheme')
xlabel('Time slot (K)')
ylabel('Throughput pre slot')
title('Throughput Comparision (Fig. 6)')
