function analyze_scenario(output_prefix, matrix_filename, captime_for_matrix)
% ANALYZE_SCENARIO: generates diagnostic plots to study an algorithm design
% scenario.
%
% Input: output_prefix: string with prefix for output plots (.eps & .fig)
%        matrix_filename: name of file with N lines for N instances, each
%        line <instance_name>, <performance of design1>, <performance of
%        design2>, ... In my experiments: design1 = default design (first
%        to be explored in IS-best)
%        captime_for_matrix: time after which runs for generating the
%        matrix were terminated -- matrix entries for terminated runs must
%        have values larger than this captime

rand('twister',1234);

data = csvread(matrix_filename,0,1);
data = data(randperm(size(data,1)),:);
K = 1000; % number of independent repetitions (sampled of IS-Best)
[P,M] = size(data);
N=P/2;

%=== Prepare data so mean is PAR with cutoff penalty 10.
cutoff_multiplier = 10;
for i=1:size(data,1)
    data(i, find(data(i,:)>captime_for_matrix-0.01)) = captime_for_matrix * cutoff_multiplier; 
end
data = max(data, 0.005);

%=== Compute train and test scores and plot train performance.
allscores = zeros(1,M);
for i=1:M
     allscores(i) = combine_performance(data(:,i));
end
[Y,ind] = sort(allscores);
def_score = allscores(1)
best_score = Y(1)

plotBasics = 1;
plot_traintest = 1;
plot_vary_N = 1;
plot_vary_captime = 1;
plot_vary_captime_and_N = 1;
plot_train_captime = 0; % not the best plot.

if plotBasics
    matrix = data'; % row: config, column: instance
    
    ind2 = ind;
    for i=1:P
         inst_hardness(i) = combine_performance(matrix(:,i));
    end
    [tmp,ind1] = sort(inst_hardness);

    figure;
    h=image(log10(matrix(ind2,ind1)), 'CDataMapping', 'scaled'); 
    colormap(gray); 
    hold on; 
%    title(output_prefix); 
    ylabel('design (sorted by PAR)', 'FontSize', 26); 
    xlabel('instance (sorted by hardness)', 'FontSize', 26);
    hold on
    h=colorbar;
    mini = min(min(log10(data)));
    maxi = max(max(log10(data)));
    
    set(gca, 'FontSize', 16);
    set(h, 'YLim', [mini,maxi])
    
    filename_prefix = strcat([output_prefix, '-wholematrix.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));


    %=== Now, plotting results on single instances (columns of matrix).
    figure
    handles = [];
    legendstrings = {};

    q50 = ind1(ceil(length(ind1)/2));
    q25 = ind1(ceil(length(ind1)/4));
    q75 = ind1(ceil(length(ind1)/4*3));
    hardest = ind1(length(ind1));
    easiest = ind1(1);

    idxs = [easiest, q25, q50, q75, hardest];
    fonts = {'k-.', 'b:', 'r-', 'c--', 'kx:'};
    widths = [2,2,3,2,1];
    
%    legendsentries = {'easiest instance', 'q_{0.25} instance', 'q_{0.50} instance', 'q_{0.75} instance', 'hardest instance'};
    legendsentries = {'easiest', 'q_{0.25}', 'q_{0.50}', 'q_{0.75}', 'hardest'};
    
    for j=1:length(idxs)
        i = idxs(j);
        numSolved = length(find(matrix(:,i) < captime_for_matrix));
        if numSolved > 0
            sorted = sort(matrix(:,i));
            h=semilogx([sorted(1), sorted(1:numSolved)', captime_for_matrix],([0, 1:numSolved,numSolved])./M*100, fonts{j});
            set(h, 'LineWidth', widths(j));
            hold on
            handles(end+1) = h;
            legendstrings{end+1} = legendsentries{j};
        end
    end

    xlabel('CPU time[s]','FontSize',26);
    ylabel('% designs successful', 'FontSize', 26);
    set(gca, 'FontSize', 16);
      
    set(gca,'XTick',[0.01, 0.1, 1, 10, 100, 1000])
    axis([0.001, captime_for_matrix*2, 0, 100]);
%    title(output_prefix);

    legend(handles, legendstrings, 'Location', 'NorthWest');    

    filename_prefix = strcat([output_prefix, '-single_insts_times.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));


    %=== Now, plotting results averaged over configs (horizontal marginal of matrix).
    figure

    h=semilogx(allscores(ind),(1:M)./M*100, 'k-');
    hold on
    xlabel('PAR score','FontSize',26);
    ylabel('% designs', 'FontSize', 26);

    set(h, 'LineWidth', 2);
    hold on
%    axis([0.1, 100, 0, 100]);
    axis([0.03, 2*captime_for_matrix*cutoff_multiplier, 0, 100]);
%    title(output_prefix);
    hDef = semilogx(allscores(1), find(ind==1)./M*100, 'rx', 'MarkerSize', 15);
    set(hDef, 'LineWidth', 2);
    handles = [hDef];
    legendstrings = {'Default'};
    set(gca, 'FontSize', 16);
    
    legend(handles, legendstrings, 'Location', 'SouthEast');
    
    filename_prefix = strcat([output_prefix, '-mean_across_insts_times.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));
 
    
    
    %=== Now, plotting results on single configs (rows of matrix).
    figure
    handles = [];
    legendstrings = {};

    default_idx = 1;
    q50 = ind2(ceil(length(ind2)/2));
    q25 = ind2(ceil(length(ind2)/4));
    q75 = ind2(ceil(length(ind2)/4*3));
    worst = ind2(length(ind2));
    best = ind2(1);

    idxs = [default_idx, best, q25, q50, q75, worst];
    fonts = {'k-.', 'b:', 'r-', 'c--', 'kx:', 'b-'};
    widths = [2,2,3,2,1,3];
    legendsentries = {'default_.', 'best', 'q_{0.25}', 'q_{0.50}', 'q_{0.75}', 'worst'};
    
    % need to have the _. since 'default' seems to be reserved in Matlab
    % legends!
    
    for j=1:length(idxs)
        i = idxs(j);
        numSolved = length(find(matrix(i,:) < captime_for_matrix));
        if numSolved > 0
            sorted = sort(matrix(i,:));
            h=semilogx([sorted(1), sorted(1:numSolved), captime_for_matrix],([0, 1:numSolved,numSolved])./P*100, fonts{j});
            set(h, 'LineWidth', widths(j));
            hold on
            handles(end+1) = h;
            legendstrings{end+1} = legendsentries{j};
        end
    end
    
    xlabel('CPU time[s]','FontSize',26);
    ylabel('% instances solved', 'FontSize', 26);
    set(gca, 'FontSize', 16);
    
    set(gca,'XTick',[0.01, 0.1, 1, 10, 100, 1000])
      
    axis([0.001, captime_for_matrix*2, 0, 100]);
%    title(output_prefix);

    legend(handles, legendstrings, 'Location', 'NorthWest');    
    filename_prefix = strcat([output_prefix, '-single_configs_times.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));
end

for i=1:M
     trainscores(i) = combine_performance(data(1:N,i));
     testscores(i) = combine_performance(data(N+1:2*N,i));
end


if plot_traintest
    numTrain = N;
    numTest = N;
    numTrainToUse = min(10, N);
    [change_times, change_incs, change_testincs, defaultTrain, defaultTest] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, numTrainToUse);

    %=== FIGURE FOR TRAIN PERFORMANCE
    figure
    loglog([1e-10,1e10], [defaultTrain,defaultTrain], 'r--', 'LineWidth', 2);

    hold on
    L = quantile(change_incs,0.1,1);
    U = quantile(change_incs,0.9,1);
    MED = quantile(change_incs,0.5,1);
    h = confplot(change_times, MED, MED-L, U-MED, 'k', 'LineWidth', 2);
%    h2 = loglog([1e-10,1e10], [trainscores(1),trainscores(1)],'r--');
    h2 = loglog([1e-10,1e10], [defaultTrain,defaultTrain],'r--', 'LineWidth', 2);

    mini = min(min(change_times))
    maxi = max(max(change_times))
    %axis([mini/2,maxi*2, min(min(change_incs))/2,max(max(change_incs))*2]);

    %    legend([h(1),h(2),h2],{'Best configuration so far (median)', 'Best configuration so far (variation)', 'Default configuration'}, 'Location', 'SouthWest')
    legend([h(1),h2],{strcat(['Incumbent, N=', num2str(numTrainToUse)]), 'Default'}, 'Location', 'SouthWest')
    ylabel('PAR score (train)','FontSize',26)
    xlabel('CPU time used for tuner [s]','FontSize',26)
%    title(output_prefix)
%    axis([10,500000,0.1,5*cutoff_multiplier])
    axis([10,captime_for_matrix*M*100,0.03,captime_for_matrix*cutoff_multiplier])
    set(gca,'XTick',[10,100,1000,10000,100000,1000000,10000000])
    set(gca,'YTick',[0.1,1,10,100,1000])
    set(gca, 'FontSize', 16);
    
    filename_prefix = strcat([output_prefix, '-trainperf.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));

    
    %=== ANOTHER FIGURE FOR TEST PERFORMANCE
    figure

    mean_train = mean(change_incs,1);
    mean_test = mean(change_testincs,1);
    h4 = loglog(change_times, mean_train, 'r--', 'LineWidth', 3);
    hold on
    h = loglog(change_times, mean_test, 'k', 'LineWidth', 2);    

    legend([h(1), h4],{'Test', 'Training, N=10'}, 'Location', 'SouthWest');
    ylabel('PAR score','FontSize',26)
    xlabel('CPU time used for tuner [s]','FontSize',26);
%    title(output_prefix)
    %axis([10,500000,0.1,5*cutoff_multiplier])
    axis([10,captime_for_matrix*M*100,0.005,captime_for_matrix*cutoff_multiplier])
    set(gca,'XTick',[10,100,1000,10000,100000,1000000,10000000])
    set(gca,'YTick',[0.01,0.1,1,10,100,1000])
    set(gca, 'FontSize', 16);

    filename_prefix = strcat([output_prefix, '-testperf.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));
end


if plot_vary_N
    %=== FIGURE WITH TEST PERFORMANCE FOR VARYING TRAINING N.
    figure
    numTrain = N;
    numTest = N;

    [change_times1, change_incs1, change_testincs1, defaultTrain1, defaultTest1] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, 1);
    [change_times10, change_incs10, change_testincs10, defaultTrain10, defaultTest10] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, 10);
    [change_times100, change_incs100, change_testincs100, defaultTrain100, defaultTest100] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, min(100, N));
    
    h4 = loglog([1e-10,1e10], [defaultTest1,defaultTest1],'k-','LineWidth', 1);
    hold on
    h1 = loglog(change_times1, mean(change_testincs1,1), 'b-', 'LineWidth', 2);
    h2 = loglog(change_times10, mean(change_testincs10,1), 'g-.', 'LineWidth', 2);
    h3 = loglog(change_times100, mean(change_testincs100,1), 'r--', 'LineWidth', 3);

%     h1 = loglog(change_times1, quantile(change_testincs1,0.5,1), 'k:');
%     h2 = loglog(change_times10, quantile(change_testincs10,0.5,1), 'r-');
%     h3 = loglog(change_times100, quantile(change_testincs100,0.5,1), 'g.-');
    
    defaultTest1
    defaultTest10
    defaultTest100
    
    nStr = num2str(min(100, N));
    legend([h4, h1, h2, h3],{'Default', 'Incumbent, N=1', 'Incumbent, N=10', strcat(['Incumbent, N=',nStr])}, 'Location', 'SouthWest');
%    ylabel('Mean runtime of config.','FontSize',26)
    ylabel('PAR score (test)','FontSize',26)
    xlabel('CPU time used for tuner [s]','FontSize',26)
%    title(output_prefix)
    mini = min(min(change_times1));
    maxi = max(max(change_times100));
    %axis([0.01,500000,0.1,5*cutoff_multiplier])
    axis([0.01,captime_for_matrix*M*100,0.03,captime_for_matrix*cutoff_multiplier])
    set(gca,'XTick',[0.01,0.1,1,10,100,1000,10000,100000,1000000,10000000])
    set(gca,'YTick',[0.1,1,10,100,1000,10000])
    set(gca, 'FontSize', 16);
    
    filename_prefix = strcat([output_prefix, '-varyN.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));
   
%    axis([mini/2,maxi*2, min(min(change_testincs100))/2,max(max(change_testincs1))*2]);
    
%     %=== OVERCONFIDENCE/OVERTUNING PLOT.
%     figure
%     loglog(change_times10, quantile(change_incs10,0.5,1), 'b:');
%     hold on
%     L = quantile(change_testincs10,0.1,1);
%     U = quantile(change_testincs10,0.9,1);
%     MED = quantile(change_testincs10,0.5,1);
% 
%     h = confplot(change_times10, MED, MED-L, U-MED, 'k');
%     h2 = loglog(change_times10, quantile(change_incs10,0.5,1), 'b:');
% 
% %     L = quantile(change_incs10,0.1,1);
% %     U = quantile(change_incs10,0.9,1);
% %     MED = quantile(change_incs10,0.5,1);
% % 
% %     h = confplot(change_times10, MED, MED-L, U-MED, 'r');
% 
%     legend([h(1), h2],{'Incumbent(test)', 'Incumbent(train)'}, 'Location', 'SouthWest');
%     ylabel('PAR score');
%     xlabel('CPU time used for tuner [s]');
%     title(output_prefix)
%     mini = min(min(change_times1));
%     maxi = max(max(change_times1));
% % 
% %     axis([mini/2,maxi*2, min(min(change_incs1))/2,max(max(change_testincs1))*2]);
%     clear change_times1;
%     clear change_incs1;
%     clear change_testincs1;
% 
%     clear change_times10;
%     clear change_incs10;
%     clear change_testincs10;
%     
%     clear change_times100;
%     clear change_incs100;
%     clear change_testincs100;
end


if plot_vary_captime
    %=== FIGURE WITH TEST PERFORMANCE FOR VARYING training captime_for_matrix.
    figure
    numTrain = N;
    numTest = N;

    [change_times1, change_incs1, change_testincs1, defaultTrain1, defaultTest1] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, min(100, N), captime_for_matrix/100);
    [change_times10, change_incs10, change_testincs10, defaultTrain10, defaultTest10] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, min(100, N), captime_for_matrix/10);
    [change_times100, change_incs100, change_testincs100, defaultTrain100, defaultTest100] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, min(100, N), captime_for_matrix);

    h4 = loglog([1e-10,1e10], [defaultTest1,defaultTest1],'k-','LineWidth', 1);
    hold on
    h1 = loglog(change_times1, mean(change_testincs1,1), 'b-', 'LineWidth', 2);
    h2 = loglog(change_times10, mean(change_testincs10,1), 'g-.', 'LineWidth', 2);
    h3 = loglog(change_times100, mean(change_testincs100,1), 'r--', 'LineWidth', 3);
    
    defaultTest1
    defaultTest10
    defaultTest100
    
    legend([h4, h1, h2, h3],{'Default', strcat(['Incumbent, \kappa=', num2str(captime_for_matrix/100), 's']), strcat(['Incumbent, \kappa=', num2str(captime_for_matrix/10), 's']), strcat(['Incumbent, \kappa=', num2str(captime_for_matrix), 's'])}, 'Location', 'SouthWest');
    ylabel('PAR score (test)','FontSize',26)
    xlabel('CPU time used for tuner [s]','FontSize',26)
%    title(output_prefix)
    mini = min(min(change_times1));
    maxi = max(max(change_times100));
    %axis([0.01,500000,0.1,5*cutoff_multiplier])
    axis([0.01,captime_for_matrix*M*100,0.03,captime_for_matrix*cutoff_multiplier])
    set(gca,'XTick',[0.01,0.1,1,10,100,1000,10000,100000,1000000,10000000])
    set(gca,'YTick',[0.1,1,10,100,1000,10000])
    set(gca, 'FontSize', 16);
    %axis([mini/2,maxi*2, min(min(change_testincs100))/2,max(max(change_testincs1))*2]);

    filename_prefix = strcat([output_prefix, '-varyC.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));
end


if plot_vary_captime_and_N
    figure

    Ns = [];
    n = N;%P;
    while 1
        Ns = [Ns, n];
        if n == 1
            break;
        end
        n = ceil(n/2);
    end

    Cs = [];
    kappa = captime_for_matrix;
    while kappa > max([captime_for_matrix / 400, 0.025])
        Cs = [Cs, kappa];
        kappa = kappa/2;
    end

    performance_mean = zeros(length(Ns),length(Cs));
    ratio_correct_comparisons = zeros(length(Ns),length(Cs));
    for n_i=1:length(Ns)
        n = Ns(n_i)
        for c_i=1:length(Cs)
            c = Cs(c_i)

            %=== FIGURE WITH TEST PERFORMANCE FOR VARYING training captime and N.
            numTrain = N;
            numTest = N;

            [mean_perf, stddev_perf, ratio_correct_mean, ratio_correct_Std] = get_perf(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, n, c)
            performance_mean(n_i,c_i) = mean_perf;
            ratio_correct_comparisons(n_i, c_i) = ratio_correct_mean;
        end
    end
    h=surf(Cs, Ns, performance_mean); xlabel('Captime, \kappa','FontSize',26); ylabel('#instances, N   ','FontSize',26); zlabel('PAR score (test)','FontSize',26)
    set(get(h,'Parent'),'XScale','log'); 
    set(get(h,'Parent'),'YScale','log'); 
    
    axis tight
    set(gca, 'FontSize', 16);
    set(gca,'YTick',[1;10;100;1000])
    set(gca,'YTickLabel',{'1';'10';'100';'1000'})
    if captime_for_matrix ==5000
        set(gca,'XTick',[0.05;0.5;5;50;500;5000])
        set(gca,'XTickLabel',{'0.05';'0.5';'5';'50';'500';'5000'})
    elseif captime_for_matrix == 300
        set(gca,'XTick',[0.3;3;30;300])
        set(gca,'XTickLabel',{'0.3';'3';'30';'300'})
    else
        set(gca,'XTick',[0.05;0.5;5])
        set(gca,'XTickLabel',{'0.05';'0.5';'5'})
    end
    view(-37.5+180,30);
    
    filename_prefix = strcat([output_prefix, '-varyNandC.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));
    
    %=== Plot ratio of correct comparisons.
    figure
    h=surf(Cs, Ns, 100*(1-ratio_correct_comparisons)); xlabel('\kappa','FontSize',26); ylabel('N   ','FontSize',26); zlabel('% wrong comparisons','FontSize',26)
    set(get(h,'Parent'),'XScale','log'); 
    set(get(h,'Parent'),'YScale','log'); 
    
    axis tight
    set(gca, 'FontSize', 16);
    set(gca,'YTick',[1;10;100;1000])
    set(gca,'YTickLabel',{'1';'10';'100';'1000'})
    if captime_for_matrix ==5000
        set(gca,'XTick',[0.05;0.5;5;50;500;5000])
        set(gca,'XTickLabel',{'0.05';'0.5';'5';'50';'500';'5000'})
    elseif captime_for_matrix == 300
        set(gca,'XTick',[0.3;3;30;300])
        set(gca,'XTickLabel',{'0.3';'3';'30';'300'})
    else
        set(gca,'XTick',[0.05;0.5;5])
        set(gca,'XTickLabel',{'0.05';'0.5';'5'})
    end
    view(-37.5+180,30);
    
    filename_prefix = strcat([output_prefix, '-varyNandC-ratio.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));    
end


if plot_train_captime
    %=== This plot is not very good since different orderings imply
    % different overall runtimes, so many runs stop early and are just
    % counted as the final performance afterwards.
    % What would be better would be a plot quantifying the speedup of
    % every run (e.g. distribution of speedups, depending on the ordering
    % in which configurations are visited.
    
    numTrain = N;
    numTest = N;
    numTrainToUse = 1000;

    [change_times, change_times_ac, change_incs, change_incs_ac, change_testincs, change_testincs_ac, defaultTrain, defaultTest] = get_train_and_test_ac(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, numTrainToUse);

    %=== FIGURE FOR TRAIN PERFORMANCE
    figure
    loglog([1e-10,1e10], [defaultTrain,defaultTrain], 'r--', 'LineWidth', 2);
    hold on
    L = quantile(change_incs,0.1,1);
    U = quantile(change_incs,0.9,1);
    MED = quantile(change_incs,0.5,1);
    h = loglog(change_times, MED, 'k--', 'LineWidth', 2);
%    h = confplot(change_times, MED, MED-L, U-MED, 'k');
%    h2 = loglog([1e-10,1e10], [trainscores(1),trainscores(1)],'r--');
%    h2 = loglog([1e-10,1e10], [defaultTrain,defaultTrain],'r--');

    MED_ac = quantile(change_incs_ac,0.5,1);
    h2 = loglog(change_times_ac, MED_ac, 'b-', 'LineWidth', 2);
    

    mini = min(min(change_times))
    maxi = max(max(change_times))
    %axis([mini/2,maxi*2, min(min(change_incs))/2,max(max(change_incs))*2]);

    %    legend([h(1),h(2),h2],{'Best configuration so far (median)', 'Best configuration so far (variation)', 'Default configuration'}, 'Location', 'SouthWest')
    legend([h,h2],{'No capping', 'Adaptive Capping'}, 'Location', 'NorthEast')
    ylabel('PAR score (test)','FontSize',26)
    xlabel('CPU time used for tuner [s]','FontSize',26)
%    title(output_prefix)
    %axis([10,500000,0.1,5*cutoff_multiplier])
    axis([10,captime_for_matrix*M*100,0.05,captime_for_matrix*cutoff_multiplier])
    set(gca,'XTick',[10,100,1000,10000,100000,1000000,10000000])
    set(gca,'YTick',[0.1,1,10,100,1000,10000])
    set(gca, 'FontSize', 16);

    filename_prefix = strcat([output_prefix, '-adaptiveCap.']);
    saveas(gcf, strcat(filename_prefix,'fig'));
    set(gcf, 'PaperPositionMode', 'auto');
    print('-depsc2', strcat(filename_prefix,'eps'));    
end


function [change_times, change_incs, change_testincs, defaultTrain, defaultTest] = get_train_and_test(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, numTrainToUse, traincaptime)
%%%%%%%% RANDOMIZE OVER ORDERING OF CONFIGURATIONS, returning training and
%%%%%%%% test qualities over time for each random ordering.
if nargin == 7
    traincaptime = captime_for_matrix;
end
assert(numTrain >= numTrainToUse);

[P,M] = size(data);
data = max(data,0.001);

times = zeros(K,M);
incumbents = zeros(K,M);
testincumbents = zeros(K,M);

%=== Artificially censor the data with the specified training captime.
trainData = data;
for i=1:size(trainData,1)
    trainData(i, find(trainData(i,:)>traincaptime-0.01)) = cutoff_multiplier*traincaptime;
end

%=== Do K repetitions, keeping track of each change_point of the incumbent.
for i=1:K
    %=== Use different test set, and different training subset in each repetition.
    instperm = randperm(P);
    trainind = instperm(1:numTrainToUse);
    testind = instperm(numTrain+1:end);
    
    %=== Compute scores for training and test phase.
%     trainscores = zeros(1,M);
%     testscores = zeros(1,M);

    %=== Faster in this vectorized form.
    trainscores = combine_performance(trainData(trainind,:));
    testscores  = combine_performance(data(testind,:));
    
    for m=1:M
%         trainscores(m) = combine_performance(trainData(trainind,m));
%         testscores(m)  = combine_performance(data(testind,m));

         tmp = min(trainData(trainind,m), traincaptime);
         traintimes(m) = sum(tmp);
    end

    defaultTrain(i) = trainscores(1);
    defaultTest(i) = testscores(1);
    %=== For initialization with the default.
    randind = [1, 1+randperm(M-1)];

    %=== For random initialization.
    %randind = randperm(M);
    
    time = traintimes(randind(1));
    incumbent = trainscores(randind(1));
    testincumbent = testscores(randind(1));

    times(i,1) = time;
    incumbents(i,1) = incumbent;
    testincumbents(i,1) = testincumbent;
    for j=2:M
        time = time + traintimes(randind(j));
        %=== Break ties randomly (for each tie, change with 0.5 probability).
        if trainscores(randind(j)) <= incumbent + 1e-5
            if trainscores(randind(j)) < incumbent
                incumbent = trainscores(randind(j));
                testincumbent = testscores(randind(j));
            elseif rand < 0.5
                incumbent = trainscores(randind(j));
                testincumbent = testscores(randind(j));
            end
        end
        incumbents(i,j) = incumbent;
        testincumbents(i,j) = testincumbent;
        times(i,j) = time;
    end
end
defaultTrain = mean(defaultTrain);
defaultTest = mean(defaultTest);

%startTime = min(times(:,1));
startTime = max(times(:,1));
%quantile(times(:,1),0.99,1);

%=== Create an array of qualities for the change points of all repetitions.
change_times = sort(unique(times(:)))';

%=== Don't start until the 90% quantile is defined.
change_times = change_times(find(change_times >= startTime));

%=== If there are too many change_times, subsample 1000 of them.
if length(change_times) > 1000
    sub = 0:log(length(change_times))/1000:log(length(change_times));
    sub = [1, floor(exp(sub)), length(change_times)];
    sub = unique(sub);
    change_times = change_times(sub);
end
P = length(change_times);
change_incs = zeros(K,P);
change_testincs = zeros(K,P);


change_times = [change_times(1), change_times];
change_incs = [ones(K,1)*defaultTrain, change_incs];
change_testincs = [ones(K,1)*defaultTest, change_testincs];

for i=1:K
%    inc = captime_for_matrix*10;
%    testinc = captime_for_matrix*10;
    inc = defaultTrain;
    testinc = defaultTest;
    p=2;
    for j=1:M
        while p<=P && times(i,j) > change_times(p)
            change_incs(i,p) = inc;
            change_testincs(i,p) = testinc;
            p=p+1;
        end
        inc = incumbents(i,j);
        testinc = testincumbents(i,j);
    end
    while p<=P
        change_incs(i,p) = inc;
        change_testincs(i,p) = testinc;
        p=p+1;
    end
end


function a=combine_performance(scores)
% Other performance metrics, such as median etc could be implemented here.
a=mean(scores,1);


function [mean_perf, stddev_perf, ratio_correct_mean, ratio_correct_std] = get_perf(K, data, captime_for_matrix, cutoff_multiplier, numTrain, numTest, numTrainToUse, traincaptime)
%%%%%%%% RANDOMIZE OVER ORDERING OF CONFIGURATIONS, returning training and
%%%%%%%% test qualities at end.

if nargin == 7
    traincaptime = captime_for_matrix;
end
[P,M] = size(data);
data = max(data,0.001);

trainData = data;
for i=1:size(trainData,1)
    trainData(i, find(trainData(i,:)>traincaptime-1e-5)) = cutoff_multiplier*traincaptime;
end

train_perfs = zeros(K,1);
test_perfs = zeros(K,1);
ratio_correct_means = zeros(K,1);

%=== Do K repetitions.
for i=1:K
    instperm = randperm(P);
    trainind = instperm(1:numTrainToUse);
    testind = instperm(numTrain+1:end);
    
%     %=== Compute scores for training and test phase.
%     trainscores = zeros(1,M);
%     testscores = zeros(1,M);

    %=== Faster in this vectorized form.
    trainscores = combine_performance(trainData(trainind,:));
    testscores  = combine_performance(data(testind,:));
    
    train_perf = min(trainscores);
    %Break ties randomly.
    minInd = find(trainscores <= train_perf + 1e-5);
    minIndRand = randperm(length(minInd));
    ind = minInd(minIndRand(1));
    
    test_perf = testscores(ind);
    train_perfs(i) = train_perf;
    test_perfs(i) = test_perf;
    
    
    %=== Get percentage of correct comparions.

    %=== Shuffle scores first, so we break ties at random.
    [tmp, test_sort_idx] = mysort(testscores);

    %=== Sort trainscores by test_scores
    trainscores = trainscores(test_sort_idx);

    %=== Again, shuffle to avoid systematic bias (sort is order-preserving)
    [tmp, train_sort_idx] = mysort(trainscores); % train_sort_idx(i) = j <=> the ith smallest entry in trainscores is at index j

    
    a(train_sort_idx) = 1:M;
    num_correct_comparisons = 0;
    for m=1:M
        train_idx = find(a==m);
        num_correct_comparisons = num_correct_comparisons + length(find(a(train_idx+1:end) > m));
    end
    num_correct_comparisons = num_correct_comparisons*2;

%     Straight-forward, slow implementation.    
%     num_correct_comparisons = 0;
%     for m=1:M
%         
%         train_idx = find(train_sort_idx==m);
%         test_idx = find(test_sort_idx==m);
%         num_correct_comparisons = num_correct_comparisons + length( intersect(train_sort_idx(1:train_idx-1), test_sort_idx(1:test_idx-1)) );
%         num_correct_comparisons = num_correct_comparisons + length( intersect(train_sort_idx(train_idx+1:end), test_sort_idx(test_idx+1:end)) );
%     end
    ratio_correct_means(i) = num_correct_comparisons/(M*(M-1));
end
mean_perf = mean(test_perfs);
stddev_perf = std(test_perfs);
ratio_correct_mean = mean(ratio_correct_means);
ratio_correct_std = std(ratio_correct_means);



function [sorted, sort_idx] = mysort(array)
%=== Just like sort, but NOT preserving order -- instead random tie breaking
perm = randperm(length(array));
arrayperm = array(perm);
[sorted, sort_idx_perm] = sort(arrayperm);
sort_idx(:,1) = perm(sort_idx_perm);



function H = confplot(varargin)
%CONFPLOT Linear plot with continuous confidence/error boundaries.
%
%   CONFPLOT(X,Y,L,U) plots the graph of vector X vs. vector Y with
%   'continuous' confidence/error boundaries specified by the vectors
%   L and U.  L and U contain the lower and upper error ranges for each
%   point in Y. The vectors X,Y,L and U must all be the same length.  
%
%   CONFPLOT(X,Y,E) or CONFPLOT(Y,E) plots Y with error bars [Y-E Y+E].
%   CONFPLOT(...,'LineSpec') uses the color and linestyle specified by
%   the string 'LineSpec'.  See PLOT for possibilities.
%
%   H = CONFPLOT(...) returns a vector of line handles.
%
%   For example,
%      x = 1:0.1:10;
%      y = sin(x);
%      e = std(y)*ones(size(x));
%      confplot(x,y,e)
%   draws symmetric continuous confidence/error boundaries of unit standard deviation.
%
%   See also ERRORBAR, SEMILOGX, SEMILOGY, LOGLOG, PLOTYY, GRID, CLF, CLC, TITLE,
%   XLABEL, YLABEL, AXIS, AXES, HOLD, COLORDEF, LEGEND, SUBPLOT, STEM.
%
%     © 2002 - Michele Giugliano, PhD (http://www.giugliano.info) (Bern, Monday Nov 4th, 2002 - 19:02)
%    (bug-reports to michele@giugliano.info)
%   $Revision: 1.1 $  $Date: 2007/12/07 22:02:37 $
%                        

if (nargin<2)
 disp('ERROR: not enough input arguments!');
 return;
end % if

x = [];  y = [];  z1 = [];  z2 = [];  spec = '';

switch nargin
 case 2
  y  = varargin{1};
  z1 = y + varargin{2};
  z2 = y - varargin{2};
  x  = 1:length(y);
 case 3
  x  = varargin{1};
  y  = varargin{2};
  z1 = y + varargin{3};
  z2 = y - varargin{3};
 case 4
  x  = varargin{1};
  y  = varargin{2};
  z2 = y - varargin{3};
  z1 = y + varargin{4};
end % switch

if (nargin >= 5)
 x  = varargin{1};
 y  = varargin{2};
 z1 = y + varargin{4};
 z2 = y - varargin{3};
 spec = 'ok';
end % 


p = plot(x,y,x,z1,x,z2);    YLIM = get(gca,'YLim');    delete(p);
a1 = area(x,z1,min(YLIM)); 
hold on;
%set(a1,'LineStyle','none');     set(a1,'FaceColor',[0.9 0.9 0.9]);
set(a1,'LineStyle','none');     set(a1,'FaceColor',[0.7 0.7 0.7]);
%fh: set(a1,'LineStyle','none');     set(a1,'FaceColor',[0.8 0.8 1]);
a2 = area(x,z2,min(YLIM)); 
set(a2,'LineStyle','none');     set(a2,'FaceColor',[1 1 1]);
if (~isempty(spec)),     
 spec = sprintf('p = plot(x,y,varargin{5}');
 for i=6:nargin,  spec = sprintf('%s,varargin{%d}',spec,i); end % for
 spec = sprintf('%s);',spec);
 eval(spec);
else                     p = plot(x,y); 
end;
%hold off;

%%set(gca,'Layer','top','XGrid','on','YGrid','on');               
%fh: commented: set(gca,'Layer','top');               
set(gca,'Layer','bottom');

H = [p, a1, a2];

%if (nargout>1) varargout{1} = H; end;
%if (nargout>1) varargout{1} = H; end;
