%% Model-Data Comparison - stats

define_paths;
dirCompiledData = fullfile(dirTeamerHome, 'compiled_data');

%% LOAD DATA
model_run = 1; % [1 is v1, 2 is v2]

if ~exist('nRows','var')

% Full resolution with time shifting
if model_run == 1
    load('sim_estimates_model_v1_0min_sensor_avg_time_shifting.mat');
else
    load('sim_estimates_model_v2_0min_sensor_avg_time_shifting.mat');
end

load(fullfile(dirObservationalData, 'lib', 'AP2_metadata'));

% Correct testIDs 
SWIFT_60.testID = 0*SWIFT_60.X(:,1);
SWIFT_2m.testID = 0*SWIFT_2m.X(:,1);
for i = 1:length(meta.tTests)
    SWIFT_60.testID(SWIFT_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SWIFT_2m.testID(SWIFT_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

SK_30.testID = 0*SK_30.X(:,1);
SK_60.testID = 0*SK_60.X(:,1);
SK_2m.testID = 0*SK_2m.X(:,1);
SK_5m.testID = 0*SK_5m.X(:,1);
% For convenience, get test IDs for the averaged data too
for i = 1:length(meta.tTests)
    SK_30.testID(SK_30.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_30.TS(:,1) < meta.tTests(i)+minutes(25)) = i;   
    SK_60.testID(SK_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_2m.testID(SK_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_5m.testID(SK_5m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_5m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end


%   For time-shifting:
% N_[t_shift] = 31;
% N_modelAveraging = 1;
% N_dev = 12 [4*FS + 5*SK + 3*SWIFT]
% N_survey = 9; 

%   For model_averaging:
% N_[t_shift] = 1;
% N_modelAveraging = 4;
% N_dev = 12 [4*FS + 5*SK + 3*SWIFT]
% N_survey = 9; 

% Table
% {dev_type, dev_avg, model_avg, survey_id, time_shift, n_samp, iqr_ux, iqr_uy, iqr_uh} 

% FOR ALL DIFFERENCING: MODEL - OBSERVATION

nRows = 31*1*12*9 + 1*4*12*9;
n_samp = nan(nRows,1);
dev_type = n_samp; %[FS = 0; SWIFT = 1; SK = 2]
dev_avg = n_samp;
model_avg = n_samp; 
survey_id = n_samp;
time_shift = n_samp;
stat_UX = nan(nRows,7); 
stat_UY = stat_UX; 
stat_UH = stat_UX;

ind_full = 1;


dev = FS; 
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        stat_UX(ind_full,:) = quantile(dUX(:,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(:,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(:,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(:,:,j)), 'all');
        
        dev_type(ind_full) = 0;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;
        
        ind_full = ind_full + 1;
    end     
end

dev = FS_30; 
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        stat_UX(ind_full,:) = quantile(dUX(:,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(:,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(:,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(:,:,j)), 'all');
        
        dev_type(ind_full) = 0;
        dev_avg(ind_full) = 30;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;
        
        ind_full = ind_full + 1;
    end     
end

dev = FS_60; 
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        stat_UX(ind_full,:) = quantile(dUX(:,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(:,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(:,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(:,:,j)), 'all');
        
        dev_type(ind_full) = 0;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;
        
        ind_full = ind_full + 1;
    end     
end


dev = FS_2m; 
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        stat_UX(ind_full,:) = quantile(dUX(:,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(:,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(:,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(:,:,j)), 'all');
        
        dev_type(ind_full) = 0;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;
        
        ind_full = ind_full + 1;
    end     
end

dev = SWIFT;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey, :,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey, :,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey, :,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey, :,j)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end     

dev = SWIFT_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey, :,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey, :,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey, :,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey, :,j)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 


dev = SWIFT_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey, :,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey, :,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey, :,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey, :,j)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 


dev = SK;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:,j)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 

dev = SK_30;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:,j)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 30;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 

dev = SK_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:,j)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 

dev = SK_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:,j)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 


dev = SK_5m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
    for j = 1:length(t_shift)
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:,j), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:,j), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:,j), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:,j)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 300;
        model_avg(ind_full) = 0;
        time_shift(ind_full) = t_shift(j);
        survey_id(ind_full) = i;

        ind_full = ind_full + 1;
    end
end 

end

%% Model-averaging - 0 minute (full res, for comparison) - no time shifting
clear FS FS_30 FS_60 FS_2m 
clear SWIFT SWIFT_60 SWIFT_2m
clear SK SK_30 SK_60 SK_2m SK_5m

% Load data
if model_run == 1
    load('sim_estimates_model_v1_0min_sensor_avg.mat');  
elseif model_run == 2
    load('sim_estimates_model_v2_0min_sensor_avg.mat');  
end

model_avg_temp = 0;
time_shift_temp = 100; % Placeholder value to allow separate lookup for what should be redundant information. 

% Correct testIDs 
SWIFT_60.testID = 0*SWIFT_60.X(:,1);
SWIFT_2m.testID = 0*SWIFT_2m.X(:,1);
for i = 1:length(meta.tTests)
    SWIFT_60.testID(SWIFT_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SWIFT_2m.testID(SWIFT_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

SK_30.testID = 0*SK_30.X(:,1);
SK_60.testID = 0*SK_60.X(:,1);
SK_2m.testID = 0*SK_2m.X(:,1);
SK_5m.testID = 0*SK_5m.X(:,1);
% For convenience, get test IDs for the averaged data too
for i = 1:length(meta.tTests)
    SK_30.testID(SK_30.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_30.TS(:,1) < meta.tTests(i)+minutes(25)) = i;   
    SK_60.testID(SK_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_2m.testID(SK_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_5m.testID(SK_5m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_5m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

% Process data

dev = FS;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 0;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = time_shift_temp;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_30;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 30;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = time_shift_temp;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_60;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 60;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = time_shift_temp;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_2m;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 120;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = time_shift_temp;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = SWIFT;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_30;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 30;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_5m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 300;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = time_shift_temp;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

%% Model-averaging - 1 min - No time-shifting
clear FS FS_30 FS_60 FS_2m 
clear SWIFT SWIFT_60 SWIFT_2m
clear SK SK_30 SK_60 SK_2m SK_5m

% Load data
if model_run == 1
    load('sim_estimates_model_v1_1min_sensor_avg.mat');  
elseif model_run == 2
    load('sim_estimates_model_v2_1min_sensor_avg.mat');  
end

model_avg_temp = 60;

% Correct testIDs 
SWIFT_60.testID = 0*SWIFT_60.X(:,1);
SWIFT_2m.testID = 0*SWIFT_2m.X(:,1);
for i = 1:length(meta.tTests)
    SWIFT_60.testID(SWIFT_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SWIFT_2m.testID(SWIFT_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

SK_30.testID = 0*SK_30.X(:,1);
SK_60.testID = 0*SK_60.X(:,1);
SK_2m.testID = 0*SK_2m.X(:,1);
SK_5m.testID = 0*SK_5m.X(:,1);
% For convenience, get test IDs for the averaged data too
for i = 1:length(meta.tTests)
    SK_30.testID(SK_30.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_30.TS(:,1) < meta.tTests(i)+minutes(25)) = i;   
    SK_60.testID(SK_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_2m.testID(SK_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_5m.testID(SK_5m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_5m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

dev = FS;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 0;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_30;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 30;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_60;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 60;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_2m;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 120;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = SWIFT;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_30;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 30;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_5m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 300;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

%% Model-averaging - 2 minute - no time shifting
clear FS FS_30 FS_60 FS_2m 
clear SWIFT SWIFT_60 SWIFT_2m
clear SK SK_30 SK_60 SK_2m SK_5m

% Load data
if model_run == 1
    load('sim_estimates_model_v1_2min_sensor_avg.mat');  
elseif model_run == 2
    load('sim_estimates_model_v2_2min_sensor_avg.mat');  
end

model_avg_temp = 120;

% Correct testIDs 
SWIFT_60.testID = 0*SWIFT_60.X(:,1);
SWIFT_2m.testID = 0*SWIFT_2m.X(:,1);
for i = 1:length(meta.tTests)
    SWIFT_60.testID(SWIFT_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SWIFT_2m.testID(SWIFT_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

SK_30.testID = 0*SK_30.X(:,1);
SK_60.testID = 0*SK_60.X(:,1);
SK_2m.testID = 0*SK_2m.X(:,1);
SK_5m.testID = 0*SK_5m.X(:,1);
% For convenience, get test IDs for the averaged data too
for i = 1:length(meta.tTests)
    SK_30.testID(SK_30.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_30.TS(:,1) < meta.tTests(i)+minutes(25)) = i;   
    SK_60.testID(SK_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_2m.testID(SK_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_5m.testID(SK_5m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_5m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

% Process data

dev = FS;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 0;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_30;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 30;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_60;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 60;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_2m;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 120;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = SWIFT;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_30;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 30;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_5m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 300;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

%% Model-averaging - 5 minute - no time shifting
clear FS FS_30 FS_60 FS_2m 
clear SWIFT SWIFT_60 SWIFT_2m
clear SK SK_30 SK_60 SK_2m SK_5m

% Load data
if model_run == 1
    load('sim_estimates_model_v1_5min_sensor_avg.mat');  
elseif model_run == 2
    load('sim_estimates_model_v2_5min_sensor_avg.mat');  
end

model_avg_temp = 300;

% Correct testIDs 
SWIFT_60.testID = 0*SWIFT_60.X(:,1);
SWIFT_2m.testID = 0*SWIFT_2m.X(:,1);
for i = 1:length(meta.tTests)
    SWIFT_60.testID(SWIFT_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SWIFT_2m.testID(SWIFT_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SWIFT_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

SK_30.testID = 0*SK_30.X(:,1);
SK_60.testID = 0*SK_60.X(:,1);
SK_2m.testID = 0*SK_2m.X(:,1);
SK_5m.testID = 0*SK_5m.X(:,1);
% For convenience, get test IDs for the averaged data too
for i = 1:length(meta.tTests)
    SK_30.testID(SK_30.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_30.TS(:,1) < meta.tTests(i)+minutes(25)) = i;   
    SK_60.testID(SK_60.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_60.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_2m.testID(SK_2m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_2m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
    SK_5m.testID(SK_5m.TS(:,1) > meta.tTests(i)-minutes(5) &...
        SK_5m.TS(:,1) < meta.tTests(i)+minutes(25)) = i;
end

dev = FS;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 0;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_30;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 30;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_60;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 60;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = FS_2m;
for i = 1:9
    UX = dev(i).UX;
    UY = dev(i).UY;
    UH = (UX.^2 + UY.^2).^0.5;
    UXQ = dev(i).UXQ;
    UYQ = dev(i).UYQ;
    UHQ = (UXQ.^2 + UYQ.^2).^0.5;
    dUX = bsxfun(@minus, UXQ, UX);
    dUY = bsxfun(@minus, UYQ, UY);
    dUH = bsxfun(@minus, UHQ, UH);    
    metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
    stat_UX(ind_full,:) = quantile(dUX(:,:), metrics, 'all');
    stat_UY(ind_full,:) = quantile(dUY(:,:), metrics, 'all');
    stat_UH(ind_full,:) = quantile(dUH(:,:), metrics, 'all');
    n_samp(ind_full) = sum(~isnan(dUH(:,:)), 'all');
    dev_type(ind_full) = 0;
    dev_avg(ind_full) = 120;
    model_avg(ind_full) = model_avg_temp;
    time_shift(ind_full) = 0;
    survey_id(ind_full) = i;
    ind_full = ind_full + 1;
end 

dev = SWIFT;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SWIFT_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 1;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 0;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_30;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 30;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_60;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 60;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_2m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);
for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 120;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

dev = SK_5m;
UX = dev.UX;
UY = dev.UY;
UH = (UX.^2 + UY.^2).^0.5;
UXQ = dev.UXQ;
UYQ = dev.UYQ;
UHQ = (UXQ.^2 + UYQ.^2).^0.5;
dUX = bsxfun(@minus, UXQ, UX);
dUY = bsxfun(@minus, UYQ, UY);
dUH = bsxfun(@minus, UHQ, UH);

for i = 1:9
        metrics = [0, 0.1, 0.25, 0.5, 0.75, 0.9, 1];
        bSurvey = dev.testID(:,1) == i;
        stat_UX(ind_full,:) = quantile(dUX(bSurvey,:), metrics, 'all');
        stat_UY(ind_full,:) = quantile(dUY(bSurvey,:), metrics, 'all');
        stat_UH(ind_full,:) = quantile(dUH(bSurvey,:), metrics, 'all');
        n_samp(ind_full) = sum(~isnan(dUH(bSurvey,:)), 'all');
        dev_type(ind_full) = 2;
        dev_avg(ind_full) = 300;
        model_avg(ind_full) = model_avg_temp;
        time_shift(ind_full) = 0;
        survey_id(ind_full) = i;
        ind_full = ind_full + 1;
end 

%% SAVE CONSTRUCTED STATISTICS

% Load data
if model_run == 1
   save_name = 'assembled_comparison_statistics_v1';
elseif model_run == 2
   save_name = 'assembled_comparison_statistics_v2';
end


save(fullfile(dirCompiledData,save_name),...
    'stat_UX','stat_UY', 'stat_UH', 'n_samp',...
     'dev_type', 'dev_avg', 'model_avg', 'time_shift', 'survey_id',...
     't_shift','meta');





