global absolute;
absolute = false;    %adjusts plots to either absolute or relative errors

global relative_error_map;
relative_error_map =[0 0 0;
        0 0.75 0.75;
        0 0.5 0.75;
        0.75 0.75 0;
        0.75 0 0;
        0 0 0];

global error_diff_map;
error_diff_map =[0 0 0;
        0 0.75 0;
        0 0.5 0;
        0 0.25 0;
        0 0 0;
        0.25 0 0;
        0.5 0 0;
        0.75 0 0;
        0 0 0];
    
    
 [pattern_mtx, shift_mtx, comp_mtx, shift_diff_mtx, comp_diff_mtx]=read_data();
  
  ECDF_Prolog(pattern_mtx(12, :));
 plot_rel(pattern_mtx, '', '-');
 plot_rel(shift_mtx, 'o', "-.");
 plot_rel(comp_mtx, 'x', "--");
  
  
  
 diffs_Prolog(comp_diff_mtx(12, :));
 plot_diff(shift_diff_mtx, 'o', "-.");
 plot_diff(comp_diff_mtx, 'x', "--");
 
 [pavg, psd, p95, savg, ssd, s95, cavg, csd, c95]=read_scalars();
 plot_p95_prolog();
 plot_p95(pavg, psd, p95, [0, 0.75, 0], "Pattern Prediction");
 plot_p95(savg, ssd, s95, 'r', "Time Shift");
 plot_p95(cavg, csd, c95, 'b', "BGLP2020 Best");

 %from now, the real work begins...


function [pattern_row, shift_row, comp_row, shift_diff_row, comp_diff_row]=read_data_rows(ph)
    fname_parts = {pwd,'\output\bglp\',int2str(ph),'\all.csv'};
    fname = strjoin(fname_parts, '');
    data = readmatrix(fname);
   
    
    pattern_row = data(34:34, 2:7);   
    shift_row = data(35:35, 2:7);   
    comp_row = data(36:36, 2:7);   
    
    shift_diff_row = data(38:38, 2:10);
    comp_diff_row = data(39:39, 2:10);
end


function [pattern_mtx, shift_mtx, comp_mtx, shift_diff_mtx, comp_diff_mtx]=read_data()
    pattern_mtx = [];
    shift_mtx = [];
    comp_mtx = [];
    shift_diff_mtx = [];
    comp_diff_mtx = [];

    for ph=5:5:60
        
        [pattern_row, shift_row, comp_row, shift_diff_row, comp_diff_row] = read_data_rows(ph);

        pattern_mtx = [pattern_mtx; pattern_row];
        shift_mtx = [shift_mtx; shift_row];
        comp_mtx = [comp_mtx; comp_row];
        shift_diff_mtx = [shift_diff_mtx; shift_diff_row];      
        comp_diff_mtx = [comp_diff_mtx; comp_diff_row];
    end        
end


function ECDF_Prolog(err_pos)
    figure(1);
    view(2);
    
    title('Sliced ECDFs of the prediction signals at selected relative errors');
    xlabel('Cumulative probability [%]');
    ylabel('Prediction horizon [minutes]');
    
    hold on;
    
    ylim([5, 60]);
    xt = [35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 99];
    xticks(xt);
    xlim([35, 100]);
    

    grid on;
    rectangle('Position',[0, 5, 100, 55]);       

    clinical_clr = [0, 0.75, 0, 0.3];
    rectangle('Position',[95, 5, 100, 60], 'FaceColor', clinical_clr, 'EdgeColor', clinical_clr);
    clt = text(97.5, 30, '{To be clinically accurate at a given prediction horizon,\newline a prediction method must reach this area with 30% relative error.}');
    set(clt, 'Rotation',90);
    
    
    global relative_error_map;
    err_labels = ["None", "10%", "15%", "20%", "30%", "Max"];
    
    for i=2:length(err_pos)-1
        text(err_pos(i)*100+3, 59, err_labels(i), 'Color', relative_error_map(i, :));
    end
    
    %text(60, 62, 'Relative Error [%]');
    text(47, 63, '{\leftarrow Pattern prediction is worse                     Relative-error [%]                     Pattern prediction is better \rightarrow}');
    
    text(40, 27, {'{-- - Pattern Prediction}', '{o - Time Shift}', '{x - Best BGLP2020 Method}'}, 'EdgeColor', 'k', 'BackgroundColor', 'w');
end    

function plot_rel(mtx, ps, ls)
    
    global relative_error_map;
    
    
    y = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60];
    
    for i=2:length(mtx(1, :))-1
        x = mtx(:, i);
        x = x*100;
        xq = x;
        yq = makima(x,y,xq);
        
        if (ls =='none') 
            scatter(x, y, ps);
        else
            plot(x,y, ps,xq,yq, ls, 'MarkerEdgeColor', relative_error_map(i, :), 'Color', relative_error_map(i, :));
        end
        
        
    end

end


function diffs_Prolog(diff_pos)
    global error_diff_map;

    figure(2);
    view(2);
    title('Sliced ECDFs of the prediction signals at selected relative-errors differences against the pattern prediction');

    xlabel('Cumulative probability [%]');
    ylabel('Prediction horizon [minutes]');

    hold on;
    
    ylim([5, 60]);
    xt = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100];
    xticks(xt);
 
    grid on;
    rectangle('Position',[0, 5, 100, 55]);
    
    
    diff_labels = ["-Max", "-30%", "-20%", "-10%", "0%", "10%", "20%", "30%", "Max"];
    
    for i=2:length(diff_pos)-1
        text(diff_pos(i)*100-3, 59, diff_labels(i), 'Color', error_diff_map(i, :));
    end
    
    text(18, 63, '{\leftarrow Pattern prediction is better                     Relative-Error Difference [%]                     Pattern prediction is worse \rightarrow}');
    
    text(65, 20, {'{o - Time Shift}', '{x - Best BGLP2020 Method}'}, 'EdgeColor', 'k', 'BackgroundColor', 'w');
end


function plot_diff(mtx, ps, ls)
    
    global error_diff_map;
    
    
    y = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60];
    
    for i=2:length(mtx(1, :))-1
        x = mtx(:, i);
        
        f_idx = 2;
        l_idx = length(x);
        
        first = x(1);
        while (first == x(f_idx)) && (f_idx<=l_idx)
            f_idx = f_idx+1;
        end
        f_idx = f_idx-1;
        
       
        x = x(f_idx:l_idx);
                       
        x = x*100;
        xq = x;
        yq = makima(x,y(f_idx:l_idx),xq);
        
        if (ls =='none') 
            scatter(x, y(f_idx:l_idx), ps);
        else
            plot(x,y(f_idx:l_idx), ps,xq,yq, ls, 'MarkerEdgeColor', error_diff_map(i, :), 'Color', error_diff_map(i, :));
        end
        
        
    end

end



function [pavg, psd, p95, savg, ssd, s95, cavg, csd, c95]=read_data_scalars(ph)
    fname_parts = {pwd,'\output\bglp\',int2str(ph),'\all.csv'};
    fname = strjoin(fname_parts, '');
    data = readmatrix(fname);
   
    
    savg = data(1:1, 2:2);   
    ssd = data(2:2, 2:2);   
    s95 = data(9:9, 2:2);   
    
    pavg = data(1:1, 3:3);   
    psd = data(2:2, 3:3);   
    p95 = data(9:9, 3:3);   
    
    cavg = data(1:1, 4:4);   
    csd = data(2:2, 4:4);   
    c95 = data(9:9, 4:4);   
end


function [pavg, psd, p95, savg, ssd, s95, cavg, csd, c95]=read_scalars()
    pavg= [];
    psd= [];
    p95= [];
    savg= [];
    ssd= [];
    s95= [];
    cavg= [];
    csd= [];
    c95= [];
    

    for ph=5:5:60
        
        [r_pavg, r_psd, r_p95, r_savg, r_ssd, r_s95, r_cavg, r_csd, r_c95] = read_data_scalars(ph);

        pavg = [pavg; r_pavg];
        psd = [psd; r_psd];
        p95 = [p95; r_p95];
        
        savg = [savg; r_savg];
        ssd = [ssd; r_ssd];
        s95 = [s95; r_s95];
        
        cavg = [cavg; r_cavg];        
        csd = [csd; r_csd];        
        c95 = [c95; r_c95];        
    end        
end


function plot_p95_prolog
    figure(3);
    view(2);
    hold on;
    grid on;
   
    global absolute;
    
    if absolute 
        title('Absolute error markers');
        xlabel('Absolute Error [mmol/L]');
        ylabel('Prediction horizon [minutes]');
        xlim([0, 5]);
    else
        xlabel('Relative Error [%]');
        title('Relative error markers');
        
        xt = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60];
    
        xlim([0, 60]);
        xticks(xt);    
    end
    
     

    
    
    
    
    ylim([5, 60]);
%   text(10, 50, {
%        '{\color{green}-- - Pattern Prediction}', '{\color{red}-- - Time Shift}', '{\color{blue}-- - Best BGLP2020 Method}'}, 'EdgeColor', 'k', 'BackgroundColor', 'w'); 
    
  lg = legend();
  lg_rect = [0.7, 0.2, 0.2, 0.3];
  set(lg, 'Position', lg_rect)
end

function plot_p95(avg, sd, p95, clr, name) 
  y = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60];
  
  global absolute;
  
  if absolute 
      sd = avg+sd;
  else
      avg = avg*100;
      sd = avg + sd*100;
      p95 = p95*100;
  end
  
  
  xq = avg;
  yq = makima(avg,y,xq);
  plot(avg,y, "-",xq,yq, "-", 'Color', clr, 'DisplayName', append(name,' avg'));
  
  
  xq = sd;
  yq = makima(sd,y,xq);
  plot(sd,y, "-.", xq,yq, "-.", 'Color', clr, 'DisplayName', append(name,' avg+std.dev.'));
  
  
  xq = p95;
  yq = makima(p95,y,xq);
  plot(p95,y, "--",xq,yq, "--", 'Color', clr, 'DisplayName', append(name,' 95th percentile'));
end