%% Fig.4 in Paper: Mean Value of Normalized Power Allocation (Spherical Uncertainty Regions)
close all;
clear all;

% Parameters
varH = 2; % variance of channel matrix H
varE = 1; % variance of error matrix E
nT = 4; % number of transmit antennas
nR = 6; % number of receiver antennas
g = 0:0.1:0.9;
g = [g 0.95 0.99]; % parameter g

% Optimization problem with CVX (Eq.(23) in paper)
N = 300; % number of iterations
pDist = zeros(nT,12); % normalized power allocation for different radius
for a = 1:N
    H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
    E = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varE/2); % estimation error matrix
    estH = (H+E)*varH/(varH+varE); % estimated CSI matrix
    [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
    r = norm(estH(:))*g; % radius of spherical uncertainty region
    for b = 1:12
        cvx_begin
            variable t;
            variable delta(nR,nT) complex;     
            dual variables gamma{nT};
            minimize (t);
            norm(delta(:)) <= r(b);
            for c = 1:nT
                gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
            end
        cvx_end
        pDist(:,b) = pDist(:,b)+cell2mat(gamma);
    end
end
pDist = pDist/N; % mean value of normalized power allocation
save('data/powerDist', 'g', 'pDist'); % save data

% Figures
output_fig = 'figure/powerDist';
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 16);
hold on;
grid on;
plot(g, pDist(4,:), '-o', g, pDist(3,:), '--x', g, pDist(2,:), '-.*', g, pDist(1,:), ':^', 'linewidth', 2);
legend('Mean value of \gamma_1', 'Mean value of \gamma_2', 'Mean value of \gamma_3', 'Mean value of \gamma_4');
xlabel('Parameter g');
ylabel('Mean Value of \gamma_i');
eval(['print -depsc2 ' output_fig]);