x2zFunctions = {                            ...
    @prob_linear_RS_noise_x2z_test,         ...
    @prob_multiple_DP_x2z_test,             ...
    @prob_quadratic_RS_10term_x2z_test,     ...
    @prob_quadratic_RS_convex_x2z_test,     ...
    @prob_quadratic_RS_concave_x2z_test,    ...
    @prob_saddle_surface2_x2z_test,         ...
    @prob_highly_nonlinear_x2z_test,        ...
    @prob_highly_nonlinear2_x2z_test,       ...
    @prob_parallel_x2z_test,                ...
    @prob_multiple_DP_series_x2z_test,      ...
    @prob_multiple_DP_series2_x2z_test          };

randomVariables = {
    @prob_linear_RS_noise_stochast_test,        ...
    @prob_multiple_DP_stochast_test,            ...
    @prob_quadratic_RS_10term_stochast_test,    ...
    @prob_quadratic_RS_convex_stochast_test,    ...
    @prob_quadratic_RS_concave_stochast_test,   ...
    @prob_saddle_surface2_stochast_test,        ...
    @prob_highly_nonlinear_stochast_test,       ...
    @prob_highly_nonlinear2_stochast_test,      ...
    @prob_parallel_stochast_test,               ...
    @prob_multiple_DP_series_stochast_test,     ...
    @prob_multiple_DP_series2_stochast_test         };

exactResults = [    ...
    1.22e-2,        ...
    1.46e-7,        ...
    5.34e-3,        ...
    4.16e-3,        ...
    1.05e-1,        ...
    3.47e-2,        ...
    2.86e-3,        ...
    1.80e-4             ];

N = 100;
%results = cell(length(x2zFunctions), N);

clc;
warning off
for i = 6 %1:length(x2zFunctions)
    
    for j = 1:N
        
        if isempty(results{i,j})
            
            fprintf('%s, seed #%d\n', func2str(x2zFunctions{i}), j);
            
            try
                stochast = feval(randomVariables{i});

                results{i,j} = DS(...
                    'stochast',             stochast,           ...
                    'x2zFunction',          x2zFunctions{i},    ...
                    'minsamples',           0,                  ...
                    'DesignPointDetection', false,              ...
                    'ARS',                  false,              ...
                    'confidence',           .95,                ...
                    'seed',                 j                       );
            catch e
                disp(e.message)
            end
        
        end
    end
end

%save('GrootemanResults4.mat', 'x2zFunctions', 'randomVariables', 'exactResults', 'results');

%% print results

clc;

fprintf('%4s %-40s %10s %10s %10s %10s %10s  %6s %6s %6s %6s %6s\n', ...
    'ID', 'Name', 'Pf (exact)', 'Pf (mean)', 'Pf (std)', 'Error', ' ', ...
    'N', ' ', 'Dirs', ' ', 'Frac');

for i = 1:size(results, 1)
    result = [results{i,:}];
    
    if ~isempty(result)
        
        result = [result.Output];
        
        err      = [result.P_f] - exactResults(i);
        
        RMSE     = mean(abs(err));
        RMSE_rel = RMSE ./ exactResults(i) .* 100;
        nUn  = cellfun(@(x) size(x,1), {result.un});

        fprintf('#%03d %-40s %10.2e %10.2e %10.2e %10.2e %10.2f%% %6d %6d %6d %6d %6.2f\n', ...
            i, func2str(x2zFunctions{i}), exactResults(i), ...
            mean([result.P_f]), std([result.P_f]), ...
            RMSE, RMSE_rel, ...
            round(mean([result.Calc])), round(std([result.Calc])), ...
            round(mean(nUn)), round(std(nUn)), ...
            mean([result.Calc])/mean(nUn));
        
    end
end

%% plot results

close all;
for i = 4:5 % 1:size(results, 1)
    result = [results{i,:}];
    
    if ~isempty(result)
        
        result = [result.Output];
        
        err = abs(([result.P_f]-exactResults(i))./exactResults(i).*100);
        
        figure; hist(err, [00:10:100]); set(gca,'XLim',[00 100]);
        
        title(func2str(x2zFunctions{i}), 'Interpreter', 'none');
    end
end