% data must be loaded via SRread_script
close all
clc
%% Verify Measures
all_verified = 1;
for i_example = 1:length(example)
    try
    if sum(data.(example{i_example}).vermeas) ~= length(data.(example{i_example}).vermeas)
        disp([example{i_example} ' NOT VERIFIED! (vermeas)']);
        all_verified = 0;
    end
    for i = 1:length(data.(example{i_example}).meas_name)
        %[i_example i]
        if max(abs(data.(example{i_example}).meastest(i).relmin), ...
               abs(data.(example{i_example}).meastest(i).relmax)) > 1e-10
           disp([example{i_example} ', ' ...
                 data.(example{i_example}).meas_name{i} ...
                 ': RELATIVE ERRORS TOO HIGH! (meastest)']);
           disp([data.(example{i_example}).meastest(i).relmin ...
                 data.(example{i_example}).meastest(i).relmax]);             
            all_verified = 0;
        end
    end
    catch
        disp([example{i_example} ' field vermeas NOT FOUND!']);
    end
end
if all_verified
    disp('All examples verified (vermeas)');
end


%% Collect data
%res(1).meas = struct;
for i = 1:length(res)
res(i).meas.t_add = [];
res(i).meas.t_addSE = [];
res(i).meas.t_conv = [];
res(i).meas.t_get1 = [];
res(i).meas.t_get2 = [];
res(i).meas.t_get_dir = [];
res(i).meas.relerr = [];
end

nC_meas_count = zeros([nC_max, length(meas_filter)]);
nS_meas_count = zeros([nS_max, length(meas_filter)]);
%mean_t_get1 = nC_meas_count;
mean_t_get2 = nC_meas_count;
mean_t_add = nS_meas_count;
mean_t_total = sum(nC_meas_count, 2);
for i = 1:length(unique(ex_idx))
    i_ex = find(ex_idx==i);
%     res(i).meas.t_get1 = zeros(length(res(i).nC), length(meas_filter));
    for j = 1:length(i_ex)
        mean_t_total(res(i).nC(j)) = mean_t_total(res(i).nC(j)) ...
                                   + res(i).t_part(j) ...
                                   + res(i).t_verify(j);
        for k = 1:(length(meas_filter))
            for K = 1:(length(data.(example{i_ex(j)}).meas_name)+1)
                if strcmp(data.(example{i_ex(j)}).meas_name{k}, meas_filter{K})
                    break
                end
            end
            % K is appropriate measure in data for defined measure in
            % meas_filter(k)
        
            % res(i).nC(j) = size(data.(example{i_ex(j)}).SR_part.grouping, 2);
            res(i).meas.t_addSE(j,k) = data.(example{i_ex(j)}).time_meas1(K);
            
            res(i).meas.t_conv(j,k) = data.(example{i_ex(j)}).time_meas2(K);
            
            res(i).meas.t_add(j,k) = res(i).meas.t_conv(j,k) ...
                                   + res(i).meas.t_addSE(j,k);
            %[i,j,k,K]
            %data.(example{i_ex(j)}).meastest(K)
            res(i).meas.t_get1(j,k) = data.(example{i_ex(j)}).meastest(K).time_SR(1);
            res(i).meas.t_get2(j,k) = data.(example{i_ex(j)}).meastest(K).time_SR(2);
            res(i).meas.t_get_dir(j,k) = data.(example{i_ex(j)}).meastest(K).time_dir;
            
            res(i).meas.relerr(j,k) = max(abs(...
                data.(example{i_ex(j)}).meastest(K).relmin ...
                ), abs(...
                data.(example{i_ex(j)}).meastest(K).relmax));
            
            nC_meas_count(res(i).nC(j), k) = nC_meas_count(res(i).nC(j), k) + 1;
            nS_meas_count(res(i).nS(j), k) = nS_meas_count(res(i).nS(j), k) + 1;
            mean_t_get2(res(i).nC(j), k) = mean_t_get2(res(i).nC(j), k)...
                                         + res(i).meas.t_get2(j,k);
            mean_t_add(res(i).nS(j), k) = mean_t_add(res(i).nS(j), k)...
                                        + res(i).meas.t_add(j,k);
            mean_t_total(res(i).nC(j)) = mean_t_total(res(i).nC(j)) ...
                                        + res(i).meas.t_add(j,k) ...
                                        + data.(example{i_ex(j)}).time_meas3(K) ...
                                        + data.(example{i_ex(j)}).time_meas4(K);
        end
    end
end
sum(mean_t_total)/60/60
mean_t_get2 = mean_t_get2 ./ nC_meas_count;
filter = sum(nS_meas_count, 2) == 0;
nS_set = 1:nS_max;
nS_set(filter) = [];
mean_t_add(filter, :) = [];
nS_meas_count(filter, :) = [];
mean_t_add = mean_t_add ./ nS_meas_count;
mean_t_total = mean_t_total ./ sum(nC_meas_count, 2) * (1+1+4);

%% Creating some data (only determining p)
detp.nPU = [];
detp.tp = [];
for i = 1:length(unique(ex_idx))
    for j = 1:length(res(i).nPU)
        detp.nPU(end+1) = res(i).nPU(j);
        detp.tp(end+1,:) = [res(i).meas.t_get1(j,1)-res(i).meas.t_get2(j,1), ...
                            res(i).meas.t_get1(j,2)-res(i).meas.t_get2(j,2), ...
                            res(i).meas.t_get1(j,3)-res(i).meas.t_get2(j,3)];
    end
end
disp('linear fit for t_determine_p');
p = polyfit([detp.nPU';detp.nPU';detp.nPU'], [detp.tp(:,1);detp.tp(:,2);detp.tp(:,3)], 1);
p = p*1e3

%% Plot getMeas (with determining p)
figure,hold on,grid off,box on
plot([1:150], [1:150]*p(1)+p(2), 'k-', 'linewidth', 2, 'color', [1 1 1]*0.5)

plot(res(1).nPU, 1e3*(res(1).meas.t_get1(:,1)-res(1).meas.t_get2(:,1)), 'kx', 'MarkerSize', 5)
plot(res(1).nPU, 1e3*(res(1).meas.t_get1(:,2)-res(1).meas.t_get2(:,2)), 'kx', 'MarkerSize', 4)
plot(res(1).nPU, 1e3*(res(1).meas.t_get1(:,3)-res(1).meas.t_get2(:,3)), 'kx', 'MarkerSize', 4)
for i = 2:length(unique(ex_idx))
    plot(res(i).nPU, 1e3*(res(i).meas.t_get1(:,1)-res(i).meas.t_get2(:,1)), 'kx', 'MarkerSize', 5)
    plot(res(i).nPU, 1e3*(res(i).meas.t_get1(:,2)-res(i).meas.t_get2(:,2)), 'kx', 'MarkerSize', 4)
    plot(res(i).nPU, 1e3*(res(i).meas.t_get1(:,3)-res(i).meas.t_get2(:,3)), 'kx', 'MarkerSize', 4)
end
xlabel('nPU');ylabel('tpunknown')
xlim([1 2.5e2]);ylim([2e-5,1e-3]*1e3);
set(gca, 'yscale', 'log', 'xscale', 'log');
standardplot('t_determine_p')



%% Plot getMeas (without determining p)
figure,hold on,grid off,box on
plot(res(1).nC, res(1).meas.t_get2(:,1), 'kx', 'MarkerSize', 5)
plot(res(1).nC, res(1).meas.t_get2(:,2), 'kd', 'MarkerSize', 4)
plot(res(1).nC, res(1).meas.t_get2(:,3), 'k*', 'MarkerSize', 4)
for i = 2:length(unique(ex_idx))
    plot(res(i).nC, res(i).meas.t_get2(:,1), 'kx', 'MarkerSize', 5)
    plot(res(i).nC, res(i).meas.t_get2(:,2), 'kd', 'MarkerSize', 4)
    plot(res(i).nC, res(i).meas.t_get2(:,3), 'k*', 'MarkerSize', 4)
end
% plot(1:nC_max, mean_t_get2(:, 1), 'k-', 'linewidth', 2)
% plot(1:nC_max, mean_t_get2(:, 2), 'k-.', 'linewidth', 1)
% plot(1:nC_max, mean_t_get2(:, 3), 'k-', 'linewidth', 1)
xlabel('nC');ylabel('tpknown');
xlim([1 10]);%ylim([1e-4,1e-3]);
set(gca, 'yscale', 'log', 'xscale', 'log');
standardplot('t_p_known')

%% Plot getMeas (direct)
figure,hold on,grid off,box on
%p = p*1e-3;
%plot([1:150], [1:150]*p(1)+p(2), 'k-', 'linewidth', 2, 'color', [1 1 1]*0.5);
plot(res(1).nPU, (res(1).meas.t_get1(:,3)-res(1).meas.t_get2(:,3)), 'kx', 'MarkerSize', 4)
%plot(res(1).nPU, res(1).meas.t_get_dir(:,1), 'kx')
%plot(res(1).nPU, res(1).meas.t_get_dir(:,2), 'kd')
plot(res(1).nPU, res(1).meas.t_get_dir(:,3), 'kd', 'MarkerSize', 4)
xlim([1 2.5e2]);
legend('new method', 'direct (MPT-toolbox: extreme, convhulln)', 'Location', 'SouthEast')
for i = 2:length(unique(ex_idx))
%    plot(res(i).nPU, res(i).meas.t_get_dir(:,1), 'kx')
%    plot(res(i).nPU, res(i).meas.t_get_dir(:,2), 'kd')
    plot(res(i).nPU, (res(i).meas.t_get1(:,3)-res(i).meas.t_get2(:,3)), 'kx', 'MarkerSize', 4)
    plot(res(i).nPU, res(i).meas.t_get_dir(:,3), 'kd', 'MarkerSize', 4)
end
set(gca, 'yscale', 'log', 'xscale', 'log');

%% Fitting for createMes (add)
%f = fittype('a*x.^3');
% p1 = fit(nS_set', mean_t_add(:, 1), 'a*x.^3')
% p2 = fit(nS_set', mean_t_add(:, 2), 'a*x.^2')
% p3 = fit(nS_set', mean_t_add(:, 3), 'a*x.^1')
p1.a = 0.007774;
p2.a = 0.03;
p3.a = 0.1697;


%% Plot createMeas (add)
figure,hold on,grid off,box on
plot(nS_set, p1.a*nS_set.^3, 'k-', 'linewidth', 2, 'color', [1 1 1]*0.5)
plot(nS_set, p2.a*nS_set.^2, 'k--', 'linewidth', 1, 'color', [1 1 1]*0.5)
plot(nS_set, p3.a*nS_set.^1, 'k-', 'linewidth', 1, 'color', [1 1 1]*0.5)
%plot(nS_set, mean_t_add(:, 1), 'k-', 'linewidth', 2, 'color', [1 1 1]*0.5)
%plot(nS_set, mean_t_add(:, 2), 'k--', 'linewidth', 1, 'color', [1 1 1]*0.5)
%plot(nS_set, mean_t_add(:, 3), 'k-', 'linewidth', 1, 'color', [1 1 1]*0.5)
plot(res(1).nS, res(1).meas.t_add(:,1), 'kx')
plot(res(1).nS, res(1).meas.t_add(:,2), 'kd')
plot(res(1).nS, res(1).meas.t_add(:,3), 'k*')
for i = 2:length(unique(ex_idx))
    plot(res(i).nS, res(i).meas.t_add(:,1), 'kx')
    plot(res(i).nS, res(i).meas.t_add(:,2), 'kd')
    plot(res(i).nS, res(i).meas.t_add(:,3), 'k*')
end
xlabel('nS');ylabel('tadd');ylim([10^-2, 10^4]);
set(gca, 'yscale', 'log', 'xscale', 'log');
standardplot('t_add_measure', [1 1.5])

%% Plot createMeas (conv)
figure,hold on,grid off,box on
plot3(res(1).nC, res(1).nPU, res(1).meas.t_conv(:,1), 'kx')
plot3(res(1).nC, res(1).nPU, res(1).meas.t_conv(:,2), 'kd')
plot3(res(1).nC, res(1).nPU, res(1).meas.t_conv(:,3), 'k*')
legend(meas_filter)
for i = 2:length(unique(ex_idx))
    plot3(res(i).nC, res(i).nPU, res(i).meas.t_conv(:,1), 'kx')
    plot3(res(i).nC, res(i).nPU, res(i).meas.t_conv(:,2), 'kd')
    plot3(res(i).nC, res(i).nPU, res(i).meas.t_conv(:,3), 'k*')
end
xlabel('nC');ylabel('tcon')
set(gca, 'zscale', 'log');
% standardplot('t_p_known')

% %% Plot getMeas - mean
% figure,hold on
% plot(1:nC_max, mean_t_get(:,1), 'k-x')
% plot(1:nC_max, mean_t_get(:,2), 'k-d')
% plot(1:nC_max, mean_t_get(:,3), 'k-*')
% legend(meas_filter)


% %% Plot addMeas times
% figure,hold on
% plot(res(1).nC, res(1).meas.t_add(:,1), 'kx')
% plot(res(1).nC, res(1).meas.t_add(:,2), 'kd')
% plot(res(1).nC, res(1).meas.t_add(:,3), 'k*')
% legend(meas_filter)
% for i = 2:length(unique(ex_idx))
%     plot(res(i).nC, res(i).meas.t_add(:,1), 'kx')
%     plot(res(i).nC, res(i).meas.t_add(:,2), 'kd')
%     plot(res(i).nC, res(i).meas.t_add(:,3), 'k*')
% end

% %% Plot addMeas - mean
% figure,hold on
% plot(1:nC_max, mean_t_add(:,1), 'k-x')
% plot(1:nC_max, mean_t_add(:,2), 'k-d')
% plot(1:nC_max, mean_t_add(:,3), 'k-*')
% legend(meas_filter)

% %% Mean Total Times
% figure
% plot(1:nC_max, mean_t_total, 'linewidth', 3)