%% Examine assembled comparison data

define_paths;
model_run = 1;

%% Load data
if model_run == 1
    load('assembled_comparison_statistics_v1.mat');
elseif model_run == 2
    load('assembled_comparison_statistics_v2.mat');
end
% Description:
%   stat_UX - 7-num description of Ux difference distribution
%   stat_UY - 7-num description of Uy difference distribution
%   stat_UH - 7-num description of Uh difference distribution

%   Rows are all the assembled combinations of devices, device averaging,
%   model averaging, and time-shifting. Use logical indexing to call subset. 

%   dev_type: [0 =
%   dev_avg: [0, 30, 60, 120, 300]
%       FS: 0 = full; 30; 60; 120
%       SWIFT: 0 = full = 30 sec; 60; 120
%       SK: 0 = full = 15 sec; 30; 60; 120; 300 = all from location
%   model_avg: [0, 60, 120, 300]
%       0 is full, but really 30 seconds interpolated
%   time_shift: full array, unique values are in t_shift = [-15:15]

% meta - survey information

% For the report, we used UH statistics. To plot specificcaly for UX or UY,
% change the following line:
stat = stat_UH;

%% Assess Optimal Time Shift

opt_shift = nan(3,9);
best_med = opt_shift;
for j = 1:9
    for k = 1:3
        bRow = dev_type == k-1 & survey_id == j & dev_avg == 0 & model_avg == 0;
        [best_med(k,j), iShift] = min(abs(stat(bRow,4)));
        temp_shift = time_shift(bRow);
        opt_shift(k,j) = temp_shift(iShift);
    end
end

opt_shift_iqr = nan(3,9);
best_iqr =  nan(3,9);
opt_shift_idr = nan(3,9);
best_idr =  nan(3,9);
for j = 1:9
    for k = 1:3
        bRow = dev_type == k-1 & survey_id == j & dev_avg == 0 & model_avg == 0;
        temp_shift = time_shift(bRow);
        temp_iqr = stat(bRow,5)-stat(bRow,3);
        [best_iqr(k,j), iShift] = min(temp_iqr);
        opt_shift_iqr(k,j) = temp_shift(iShift);
        
        temp_idr = stat(bRow,6)-stat(bRow,2);
        [best_idr(k,j), iShift] = min(temp_idr);
        opt_shift_idr(k,j) = temp_shift(iShift);
    end
end

%
cla;
for i = 1:3
    bSub = dev_type == i-1 & t_shift(1);
    med = stat(bSub,3);
    iqr = stat(bSub,5)-stat(bSub,3);
    plot(med, iqr, '.')
end


%% Supplementary Figues - Optimal Time Shifting
f2 = figure(2); clf(f2); f2.Color = 'w';
ax = gca; hold(ax, 'on');
cm = lines(7);
for i = 1:3
    plot(1:9, opt_shift(i,:), '.:', 'MarkerSize', 14, 'Color', cm(i,:));
    %plot(1:9, opt_shift_iqr(i,:), '+:', 'MarkerSize', 8, 'Color', cm(i,:));
    %plot(1:9, opt_shift_idr(i,:), '^:', 'MarkerSize', 8, 'Color', cm(i,:));
end
    
legend(ax, {'floats', 'SWIFTs', 'station-keeping'})


%% Report Figure 20 - Plot Time Shifting
f1 = figure(1); clf(f1); f1.Color = 'w';
ax1 = buildAxes(f1, 3, 9); set(ax1, 'NextPlot','add');
f1.Units = 'inches';
f1.Position = [8.0000   1   9.0    5.5];

for j = 1:9
    for k = 1:3
    bRow = dev_type == k-1 & survey_id == j & dev_avg == 0 & model_avg == 0; 
    ax = ax1(k,j);
        for i = 1:length(t_shift)
            cBox = 0.5*[1,1,1];
            if opt_shift(k,j) == t_shift(i)
                cBox = 0.5*[1,0,0];
            end
            % Plot interdecile range
            % FOR HORIZONTAL BOX PLOTS
            bSub = bRow & time_shift == t_shift(i);
            temp = stat(bSub, :);
            n_temp = mode(n_samp(bSub));
            % Plot interdecile
            x = temp(1, [2,6]);
            y = t_shift(i)*[1, 1];
            plot(ax, x, y, '-', 'Color', cBox, 'LineWidth', 1);

            % Plot interquartile
            x = temp(1, [3,5]);
            plot(ax, x, y, '-', 'Color', cBox, 'LineWidth', 4);

            % Plot median
            x = temp(1, 4) * [1, 1];
            y = y + 0.3*[-1, 1];
            plot(ax, x, y, '-', 'Color', [1,1,1], 'LineWidth', 1);


        end
    % Number of samples
    if n_temp < 1000
        n_str = sprintf('%0.1fk',n_temp/1000);
    else
        n_str = sprintf('%0.0fk',round(n_temp/1000));
    end
    text(ax, 1, -14, n_str, 'HorizontalAlignment', 'right',...
        'FontSize', 8, 'FontWeight', 'normal');
    
    grid(ax, 'on');
    end
end

set(ax1, 'XLim', [-1, 1], 'YLim',[-16, 16]);

for i = 1:numel(ax1)
    ax1(i).XAxis.FontSize = 10;
    ax1(i).YAxis.FontSize = 10;
end

ylabel(ax1(2,1), 'time shift [minutes]',...
    'FontSize',12);
xlabel(ax1(3,5), 'horizontal velocity difference (model - obs.) [m/s]',...
    'FontSize',12);

for i = 1:9
    title(ax1(1, i), meta.testNames{i},...
    'FontSize',12);
end

dev_names = {'Float', 'SWIFT', 'Station-Keeping'}; 
for k = 1:3
    tx(k) = text(ax1(k,9), 1.3, 0, dev_names{k});
    set(tx(k), 'HorizontalAlignment','center','Rotation',90,'FontWeight','bold',...
    'FontSize',12);
end

%{
annotation(f1,'textbox',...
    [0.2146 0.9594 0.5170 0.038],...
    'String','Impact of time shifting on error with full resolution model and device data [no averaging]',...
    'FitBoxToText','off',...
    'EdgeColor','none');

an = annotation(f1,'textbox',...
    'Position', [0.7+(model_run-1)*0.15 0.95 0.2 0.0500],...
    'String',sprintf('Model run v%i',model_run),...
    'FontSize',12,...
    'FontWeight', 'bold',...
    'FitBoxToText','off',...
    'EdgeColor','none');
%}

set(ax1(:,2:end), 'YTickLabel', '');
set(ax1(1:2,:), 'XTickLabel', '');

% Adjust axes size
x_marg = 0.05; 
x_sep = 0.02;
w = 0.08;

y_marg = 0.08; 
y_sep = 0.02;
h = 0.28;


for j = 1:9
    for k = 1:3
        x_temp = x_marg + (j-1)*(w+x_sep);
        y_temp = y_marg + (3-k)*(h+y_sep);
        ax1(k,j).Position = [x_temp, y_temp, w, h];
    end
end

if model_run == 1
    save_name = 'time_shifting_full_res_v1';
elseif model_run == 2
    save_name = 'time_shifting_full_res_v2';
end

print(f1, fullfile(dirFig, save_name), '-djpeg', '-r600');

%% Report Figure 19 - Plot Averaging Impact
f1 = figure(1); clf(f1); f1.Color = 'w';
ax1 = buildAxes(f1, 3, 9); set(ax1, 'NextPlot','add');
f1.Units = 'inches';
f1.Position = [8.0000   1   9.0    5.5];

c_model = [20, 20, 20;...
           44,17,87;...
           17,23,87;...
           17,87,74]/256 + 0.2*[1,1,1];

model_opt = [0, 60, 120, 300];
for j = 1:9   % Each surveys
    for k = 1:3   % Each device
        if k == 1    
            % Floats
            dev_avg_opt = [0, 30, 60, 120];      
        elseif k == 2
            % SWIFTS
            dev_avg_opt = [0, 60, 120];
        elseif k == 3
            % Station-keeping
            dev_avg_opt = [0, 30, 60, 120, 300];
        end
        
        n_data_rows = length(model_opt)*(length(dev_avg_opt)+1);
        ax = ax1(k,j);
        
        for m = 1:4  % Each model averaging
            % Model averaging   
            for n = 1:length(dev_avg_opt)  % Each device averaging
            bRow = dev_type == k-1 & survey_id == j & dev_avg == dev_avg_opt(n) & model_avg == model_opt(m) & time_shift == 0;
            if m == 1
            bRow = dev_type == k-1 & survey_id == j & dev_avg == dev_avg_opt(n) & model_avg == model_opt(m) & time_shift == 0;
            end
            temp = stat(bRow, :);
            %n_temp = mode(n_samp(bRow));
            cBox = c_model(m,:);
            
            % Plot interdecile
            x = temp(1, [2,6]);
            y = (n_data_rows-((m-1)*(length(dev_avg_opt)+1)+n))*[1, 1];
            plot(ax, x, y, '-', 'Color', cBox, 'LineWidth', 1);

            % Plot interquartile
            x = temp(1, [3,5]);
            plot(ax, x, y, '-', 'Color', cBox, 'LineWidth', 4);

            % Plot median
            x = temp(1, 4) * [1, 1];
            y = y + 0.3*[-1, 1];
            plot(ax, x, y, '-', 'Color', [1,1,1], 'LineWidth', 1);

            end
        end
    % Number of samples
    %if n_temp < 1000
    %    n_str = sprintf('%0.1fk',n_temp/1000);
    %else
    %    n_str = sprintf('%0.0fk',round(n_temp/1000));
    %end
    %text(ax, 1, -14, n_str, 'HorizontalAlignment', 'right',...
    %    'FontSize', 8, 'FontWeight', 'normal');
    
    grid(ax, 'on');
    ax.YTickLabel = '';
    end
end

%
set(ax1, 'XLim', [-0.5, 0.5]);
set(ax1(1,:), 'YLim', [0, 20]);
set(ax1(2,:), 'YLim', [0, 17]);
set(ax1(3,:), 'YLim', [0, 24]);

%set(ax1(1:3,9), 'XLim', [-1, 0]);
set(ax1(1:3,5:9), 'XLim', [-1, 1]);


for i = 1:numel(ax1)
    ax1(i).XAxis.FontSize = 8;
    ax1(i).YAxis.FontSize = 8;
end

set(ax1(1:2,:), 'XTickLabel', '');

for i = 1:9
    title(ax1(1, i), meta.testNames{i},...
    'FontSize',12);
end

dev_names = {'Float', 'SWIFT', 'Station-Keeping'}; 
for k = 1:3
    ax = ax1(k,9);
    tx(k) = text(ax, ax.XLim(2)+0.15*range(ax.XLim), ax.YLim(1)+0.5*range(ax.YLim), dev_names{k});
    set(tx(k), 'HorizontalAlignment','center','Rotation',90,'FontWeight','bold',...
    'FontSize',12);
end


xlabel(ax1(3,5), 'horizontal velocity difference (model - obs.) [m/s]',...
    'FontSize',12);

% Legend plots 


for m = 1:4
    cBox = c_model(m,:);
    p_leg(m) = plot(ax1(1), NaN, NaN, '-', 'Color', cBox, 'LineWidth', 8);
end

leg_names = {'30 s - full res.', '60 s', '120 s', '300 s'}; 

lg = legend(ax1(1), p_leg, leg_names,...
    'Orientation', 'Horizontal',...
    'Position', [0, 0.95, 1, 0.05],...
    'Box', 'off',...
    'FontSize', 10);

annotation(f1,'textbox',...
    [0, 0.95 0.35 0.05],...
    'String','Model Averaging Window',...
    'FontSize',10,...
    'FitBoxToText','off',...
    'EdgeColor','none',...
    'VerticalAlignment', 'middle',...
    'HorizontalAlignment', 'center',...
    'FontWeight', 'bold'); 


%{
annotation(f1,'textbox',...
    [0.8993 0.46454 0.099664 0.14],...
    'String',{'Model Averaging','black = none = full res.','','purple = 60 s','','blue = 120 s','','green = 300 s'},...
    'FontSize',8,...
    'FitBoxToText','off',...
    'EdgeColor','none');


% Create textbox
annotation(f1,'textbox',...
    [0.094 0.8094 0.0752 0.1400],...
    'String',{'1 s = full res.','30 s','60 s','120 s'},...
    'FontSize',8,...
    'FitBoxToText','off',...
    'EdgeColor','none');

annotation(f1,'textbox',...
    [0.0995 0.509 0.07529 0.140],...
    'String',{'30 s = full res.','60 s','120 s'},...
    'FontSize',8,...
    'FitBoxToText','off',...
    'EdgeColor','none');

an = annotation(f1,'textbox',...
    [0.0978 0.2121 0.0752 0.1400],...
    'String',{'15 s = full res.','30 s','60 s','120 s','300 s'},...
    'FontSize',8,...
    'FitBoxToText','off',...
    'EdgeColor','none');

%}    
    
%{
an = annotation(f1,'textbox',...
    [0.2+(model_run-1)*0.2 0.95 0.2 0.0500],...
    'String',sprintf('Model run v%i',model_run),...
    'FontSize',12,...
    'FontWeight', 'bold',...
    'FitBoxToText','off',...
    'EdgeColor','none');
%}

% Adjust axes size
x_marg = 0.03; 
x_sep = 0.025;
w = 0.08;

y_marg = 0.08; 
y_sep = 0.02;
h_norm = 0.31;
h(1) = range(ax1(1).YLim);
h(2) = range(ax1(2,1).YLim);
h(3) = range(ax1(3,1).YLim);
h = h_norm.*h./max(h);

for j = 1:9
    for k = 1:3
        x_temp = x_marg + (j-1)*(w+x_sep);
        y_temp = y_marg + sum(h((k+1):end)) + y_sep*(3-k);
        ax1(k,j).Position = [x_temp, y_temp, w, h(k)];
    end
end


if model_run == 1
    save_name = 'compiled_impact_averaging_v1';
elseif model_run == 2
    save_name = 'compiled_impact_averaging_v2';
end

print(f1, fullfile(dirFig, save_name), '-djpeg', '-r600');

