%% Fig.8 in Paper: Probability of Service Provision vs. Probability of no QoS (TDD)
close all;
clear all;

% Parameters
varH = 2; % variance of channel matrix
noQoS = [0.00001 0.0001 0.001 0.01 0.05 0.1 0.2 0.4 0.6 0.8 1.0]; % probability of no QoS
qosP = 1 - noQoS; % probability of providing required QoS
estSNRdB = [9 7 6 5 3]; % estimation SNR (dB)
estSNR = 10.^(estSNRdB/10); % estimation SNR
varE = varH./estSNR; % variance of error matrix

% Algorithms
N = 1000; % number of iterations
pSev = N*ones(6,11); % number of service provision
for a = 1:N
    nT = 4; % number of transmit antennas
    nR = 4; % number of receive antennas
    H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
    E = randn(nR,nT) + i*randn(nR,nT); % estimation error matrix without specifying variance
    for b = 1:5
        estH = (H+E*sqrt(varE(b)/2))*varH/(varH+varE(b)); % estimated CSI matrix
        r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR(b))); % radius of spherical uncertainty region
        c = (norm(estH(:)) <= r);
        pSev(b,c) = pSev(b,c)-1; % actual CSI matrix H=0 belongs to the uncertainty region
    end
    
    nT = 2; % number of transmit antennas
    nR = 2; % number of receive antennas
    H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
    E = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varE(1)/2); % estimation error matrix
    estH = (H+E)*varH/(varH+varE(1)); % estimated CSI matrix
    r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR(1))); % radius of spherical uncertainty region
    b = (norm(estH(:)) <= r);
    pSev(6,b) = pSev(6,b)-1; % actual CSI matrix H=0 belongs to the uncertainty region
end
pSev = pSev/N; % probability of service provision
save('data/probServiceTDD.mat', 'noQoS', 'pSev'); % save data

% Figures
output_fig = 'figure/probServiceTDD';
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
hold on;
grid on;
plot(noQoS, pSev(1,:), '-o', noQoS, pSev(2,:), '-+', noQoS, pSev(3,:), '-*', 'linewidth', 2);
plot(noQoS, pSev(4,:), '-s', noQoS, pSev(5,:), '-v', noQoS, pSev(6,:), '--', 'linewidth', 2);
set(gca, 'XScale', 'log', 'Fontsize', 16);
legend('n_T=4, n_R=4, SNR_{est}=9dB', 'n_T=4, n_R=4, SNR_{est}=7dB', ...
    'n_T=4, n_R=4, SNR_{est}=6dB', 'n_T=4, n_R=4, SNR_{est}=5dB', ...
    'n_T=4, n_R=4, SNR_{est}=3dB', 'n_T=2, n_R=2, SNR_{est}=9dB', 'location', 'southeast');
xlabel('1 - P_{in} (Probability of no QoS)');
ylabel('Probability of Service Provision');
eval(['print -depsc2 ' output_fig]);