% Setup simulation parameters.
clc; clear all;


simPar = setupSimParameters2();
simPar.v = 5;
simPar.vmax = 2*simPar.v;
N_sim       = 10;
N_scan      = 5;

SNR_dB_arr = flip(linspace(-10, 15, N_scan).'); 
% SNR_dB_arr = (linspace(-10, 15, N_scan).');

fprintf('SNR Range = ');
for idx = 1:N_scan
    fprintf('%.2f ', SNR_dB_arr(idx));
end
fprintf(' (dB)\n');
fprintf('Start simulation...\n');

NMSE_LS         = zeros(N_scan, 1);
NMSE_MMSE_ISO   = zeros(N_scan, 1);
NMSE_AR         = zeros(N_scan, 1);
NMSE_OMP        = zeros(N_scan, 1); 
NMSE_singleKernelGPR    = zeros(N_scan, 1);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 
NMSE_singleKernelGPRmu    = zeros(N_scan, 1);
NMSE_mixedKernelGPR     = zeros(N_scan, 1);
NMSE_GPRbes             = zeros(N_scan, 1);  

fprintf('+==========+===================================================================================================+\n');
fprintf('|          |                                              NMSE(dB)                                             |\n');
fprintf('| SNR(dB)  +---------------------------------------------------------------------------------------------------+\n');
fprintf('|          |  LS  | MMSE(ISO) | AR(1) | Prony |  OMP  |  GPR  | GPRmu  |  GPR(2) |  GPRbes |  Time Elapsed (s) |\n');
fprintf('+----------+------+-----------+-------+-------+-------+-------+--------+---------+---------+-------------------+\n');


for idx_scan = 1:N_scan

    simPar.SNR_dB  = SNR_dB_arr(idx_scan);
    simPar.sigma2_n = 1/db2pow(simPar.SNR_dB);
    
    innerNMSE_LS = zeros(N_sim, 1);
    innerNMSE_MMSE_ISO = zeros(N_sim, 1);
    innerNMSE_AR = zeros(N_sim, 1);
    innerNMSE_Prony = zeros(N_sim, 1);
    innerNMSE_OMP = zeros(N_sim, 1); 
    innerNMSE_singleKernelGPR = zeros(N_sim, 1);
    innerNMSE_singleKernelGPRmu = zeros(N_sim, 1);
    innerNMSE_mixedKernelGPR = zeros(N_sim, 1);
    innerNMSE_GPRbes            = zeros(N_sim, 1);
%     A = zeros(N_sim, simPar.N);
    tic; %
    rng(0); 
%     load('H_raytracing.mat')
    parfor idx_sim = 1:N_sim
%         rng(idx_sim+1);
        stream = RandStream('mt19937ar', 'Seed', idx_sim);
        RandStream.setGlobalStream(stream);
%         A(idx_sim,:) = (randn(stream,[simPar.N, 1]) + 1i*randn(stream,[simPar.N, 1]))/sqrt(2);
        user = [20;0;10];
        user = [20*rand()-10;0;10*rand()];
        user_theta = -pi/3+pi/3*rand(); %cdl not use this user
        user = 15*rand()*[sin(user_theta); 0; cos(user_theta)]; 

%         [h0, UE0] = get_channel0(simPar,user);
        [h1] = get_channel1(simPar,user);
        h = h1(:,simPar.T+1);
        h = h/norm(h)*sqrt(simPar.N);
%         h = simPar.h;
%         [h] = get_channel(simPar,user); 
        y = h + sqrt(simPar.sigma2_n)*(randn([simPar.N, 1]) + 1i*randn([simPar.N, 1]))/sqrt(2);
%         y0 = h0+sqrt(simPar.sigma2_n)*(randn([simPar.N, 1]) + 1i*randn([simPar.N, 1]))/sqrt(2);
        yg = zeros(simPar.N*simPar.T, 1);
        yg_prony = zeros(simPar.N, simPar.T);
        y_ = zeros(simPar.N, 1);
        
        for idx = 1:simPar.T
            yg((idx-1)*simPar.N+1:idx*simPar.N,1) = h1(:,idx)/norm(h1(:,idx))*sqrt(simPar.N)+sqrt(simPar.sigma2_n)*(randn([simPar.N, 1]) + 1i*randn([simPar.N, 1]))/sqrt(2);
%             yg((idx-1)*simPar.N+1:idx*simPar.N,1) = h1(:,idx)/norm(h1(:,idx))*simPar.N+sqrt(simPar.sigma2_n)*norm(h1(:,idx))/simPar.N*(randn([simPar.N, 1]) + 1i*randn([simPar.N, 1]))/sqrt(2);
            yg_prony(:,idx) = yg((idx-1)*simPar.N+1:idx*simPar.N,1);
            y_ = y_ + yg((idx-1)*simPar.N+1:idx*simPar.N,1);
        end
%         for idx = 1:simPar.T
%             hidx = squeeze(H(idx,:));
%             hidx = hidx.' *1e04;
%             yg((idx-1)*simPar.N+1:idx*simPar.N,1) = hidx+sqrt(simPar.sigma2_n)*(randn([simPar.N, 1]) + 1i*randn([simPar.N, 1]))/sqrt(2);
%             yg_prony(:,idx) = yg((idx-1)*simPar.N+1:idx*simPar.N,1);
%             y_ = y_ + yg((idx-1)*simPar.N+1:idx*simPar.N,1);
%         end
%         h = squeeze(H_U_pre(1,idx_sim,1,:,:));
%         h = h*7e05;
%         h = H(simPar.T+1,:)*1e04;
%         h = h.';
        y_  = y_ /simPar.T;
        
        h_hat_LS = predict_channel_v_20(yg_prony(:,2), simPar, "LS", user);
        innerNMSE_LS(idx_sim) = calc_NMSE(h_hat_LS, h);

        h_hat_MMSE_ISO = predict_channel_v_20(yg, simPar, "MMSE_ISO", user);
        innerNMSE_MMSE_ISO(idx_sim) = calc_NMSE(h_hat_MMSE_ISO, h);
        
%         h_hat_AR = predict_channel(yg((simPar.T-1)*simPar.N+1:simPar.T*simPar.N,1), simPar, "AR");
        h_hat_AR = predict_channel_v_20(yg, simPar, "AR", user);
        innerNMSE_AR(idx_sim) = calc_NMSE(h_hat_AR, h); 

        h_hat_Prony = predict_channel_v_20(yg_prony, simPar, "Prony", user);
        innerNMSE_Prony(idx_sim) = calc_NMSE(h_hat_Prony, h);

        h_hat_OMP = predict_channel_v_20(yg_prony(:,1), simPar, "OMP", user);
        innerNMSE_OMP(idx_sim) = calc_NMSE(h_hat_OMP, h);
        
        inner_simPar = simPar; 
        inner_simPar.channelEstimator.GPR.mixedKernel = false;
%         tic;
        h_hat_singleKernelGPR = predict_channel_v_20(yg, inner_simPar, "GPR", user);
        innerNMSE_singleKernelGPR(idx_sim) = calc_NMSE(h_hat_singleKernelGPR, h); 
%  timeElapsed = toc()
        h_hat_singleKernelGPRmu = predict_channel_v_20(yg, inner_simPar, "GPRmu", user);
        innerNMSE_singleKernelGPRmu(idx_sim) = calc_NMSE(h_hat_singleKernelGPRmu, h);
% tic;
        inner_simPar.channelEstimator.GPR.mixedKernel = true;
        h_hat_mixedKernelGPR = predict_channel_v_20(yg, inner_simPar, "GPR", user);
        innerNMSE_mixedKernelGPR(idx_sim) = calc_NMSE(h_hat_mixedKernelGPR, h);
%  timeElapsed = toc()        
%         h_hat_GPRbes                = estimate_channel(yg, simPar, "GPRbes"); 
%         innerNMSE_GPRbes(idx_sim)   = calc_NMSE(h_hat_GPRbes, h); 
        
    end

    timeElapsed = toc(); 

    NMSE_LS(idx_scan)       = mean(innerNMSE_LS);
    NMSE_MMSE_ISO(idx_scan) = mean(innerNMSE_MMSE_ISO);
    NMSE_AR(idx_scan)      = mean(innerNMSE_AR);
    NMSE_Prony(idx_scan)      = mean(innerNMSE_Prony);
    NMSE_OMP(idx_scan)      = mean(innerNMSE_OMP);
    NMSE_singleKernelGPR(idx_scan)      = mean(innerNMSE_singleKernelGPR);
    NMSE_singleKernelGPRmu(idx_scan)      = mean(innerNMSE_singleKernelGPRmu);
    NMSE_mixedKernelGPR(idx_scan)       = mean(innerNMSE_mixedKernelGPR); 
%     NMSE_GPRbes(idx_scan)           = mean(innerNMSE_GPRbes); 

    fprintf('|  % 6.2f  |% 6.2f|   % 6.2f  | % 6.2f| % 6.2f| % 6.2f| % 6.2f|  % 6.2f|  %6.2f |  %6.2f |      %7.2f      |\n', ...
        simPar.SNR_dB,...
        pow2db(NMSE_LS(idx_scan)), ...
        pow2db(NMSE_MMSE_ISO(idx_scan)), ...
        pow2db(NMSE_AR(idx_scan)), ...
        pow2db(NMSE_Prony(idx_scan)), ...
        pow2db(NMSE_OMP(idx_scan)), ...
        pow2db(NMSE_singleKernelGPR(idx_scan)), ...
        pow2db(NMSE_singleKernelGPRmu(idx_scan)), ...
        pow2db(NMSE_mixedKernelGPR(idx_scan)), ...
        pow2db(NMSE_GPRbes(idx_scan)),      ...
        timeElapsed);

end
fprintf('+==========+===================================================================================================+\n');
fprintf('Sim complete...\n');

%% Visualization.

set(0,'DefaultLineMarkerSize',  4);
set(0,'DefaultTextFontSize',    14);

set(0,'DefaultAxesFontSize',    12);
set(0,'DefaultLineLineWidth',   1.4);

set(0,'defaultfigurecolor',     'w');

figure('color', [1 1 1]);

% C = linspecer(5);
C = [   0.0000    0.0000    0.0000
        0.2000    0.4000    0.8494
        0.3718    0.7176    0.4000
        1.0000    0.5482    0.1000
        0.9650    0.1110    0.1330
        0.8650    0.1110    0.8330
        0.0000    0.0000    0.9000
    ]; 

plot(SNR_dB_arr, pow2db(NMSE_LS),       'Marker', 'diamond', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(1, :)); grid on; hold on;
plot(SNR_dB_arr, pow2db(NMSE_MMSE_ISO), 'Marker', 'square', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(3, :));
plot(SNR_dB_arr, pow2db(NMSE_AR),      'Marker', 'o', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(4, :));
plot(SNR_dB_arr, pow2db(NMSE_Prony),      'Marker', 'o', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(7, :));
% plot(SNR_dB_arr, pow2db(NMSE_OMP),      'Marker', '^', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(4, :)); 

% plot(SNR_dB_arr, pow2db(NMSE_GPRbes),   'Marker','o', 'MarkerSize',6, 'MarkerFaceColor', 'w', 'Color', C(7, :)); 
plot(SNR_dB_arr, pow2db(NMSE_singleKernelGPR),      'Marker', '>', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(5, :));
% plot(SNR_dB_arr, pow2db(NMSE_singleKernelGPRmu),      'Marker', '>', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(6, :));
plot(SNR_dB_arr, pow2db(NMSE_mixedKernelGPR),      'Marker', '>', 'MarkerSize', 6, 'MarkerFaceColor','w', 'color', C(6, :));

xlabel('SNR (dB)', 'Interpreter','latex');
ylabel('NMSE (dB)', 'Interpreter','latex');
legend({'No prediction LS', 'No learning EM', 'AR', "Prony", "EM kernel learning", "GEM kernel learning"});
set(gca,'FontName','Times New Roman');





