%
% MULTIpLE benchmark evaluation.
%
% To evaluate MULTIpLE and the baselines w.r.t. Caltech-256, 5 related classes, run:
% >> NplusoneBenchmark.show
%
% This demo showcase takes approximately 2 hours to complete on Intel i7, 16G memory machine.
%
% To evaluate one of the following tasks:
% CODENAME                METHOD
% -------------------------------------------------
% nplusone                MULTIpLE
% priors_and_new          Source+1
% source_plus_one_hingel  Source+1 (hinge)
% no_transfer_all_data    Batch
% no_transfer             No transfer
% priors                  Source
% mktl                    MKTL
% multikt_ova             MultiKT-OVA
% pmtsvm_ova              PMT-SVM-OVA
% mtradaboost_ova         MultisourceTrAdaBoost-OVA
%
% run: >> NplusoneBenchmark.run('<name_of_the_task>')
%
%
% This code is part of the supplementary material to the CVPR 2013 paper
% "From N to N+1: Multiclass Transfer Incremental Learning",
% Ilja Kuzborskij, Francesco Orabona and Barbara Caputo
%
% Copyright (C) 2012-2013, Ilja Kuzborskij
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% Contact the author: ilja.kuzborskij [at] gmail.com
%

classdef NplusoneBenchmark < handle

    properties
        dataset = [];

        dataset_file = './data/n_5_related_14_features.mat';                    % dataset file
        support_dataset_file = '';                                            % unused in current version

        prior_location = 'tmp';                                               % source model location
        result_location = 'tmp';                                              % result location

        n_classes = [];
        new_class = [];
        task = '';

        prior_sample_slice = 1:30;                                            % source sample indices
        test_sample_slice = 31:60;                                            % test samples indices
        train_sample_slice = {61:60+5, 61:60+10, 61:60+15, 61:60+20};         % train sample progression indices
        splits=1:10;

        C_range = 10.^(-5:8);                                                 % C hyperparameter range
        gamma_range = {2.^(-5:8)};                                            % RBF gamma parameter ({} means combination)
        kernel_type = 'rbf';
        restrict_to_cue_ix = [];                                              % whether to restrict to feature (e.g. only SIFT) ( [] means all features)

        beta_iterations = 3000;                                              % # of \beta optimization iterations
        zero_beta = false;
        with_bias = true;

        show_conf_matrix = false;
    end

    methods

        function results = eval_priors(obj)
            disp('>>>>>>>>>>>>>>>>>>> PRIORS <<<<<<<<<<<<<<<<<<');

            results = {};

            for split=obj.splits
                priors = obj.get_priors(split);
                priors.parameters.new_class_label = obj.new_class; % required, so that evaluator knows new class

                [x, y] = obj.get_data(obj.get_classes, split, obj.test_sample_slice);
                rs = priors.evaluate(x, y);
                results{split, 1} = rs;

                fprintf('Split: %d; N+1 Accuracy: %.2f; N Accuracy: %.2f; +1 Accuracy: %.2f\n', ...
                        split, rs.test_accuracy, rs.n_accuracy, rs.npp_accuracy);

                if obj.show_conf_matrix
                    disp('Confusion matrix:');
                    disp(rs.conf_mat);
                end
            end
        end

        function rs = no_transfer(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "no transfer" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {obj.kernel_type};
            hyperr.gamma = obj.gamma_range;

            disp('tuning hyperparameters...');
            tuner = HyperSearch(MulticlassRLS(struct('optimize_l_tuning', true, 'with_bias', obj.with_bias)), hyperr);
            hyper = tuner.search(train_x, train_y)

            disp('training...');
            classifier = MulticlassRLS(struct('optimize_l_tuning', false, 'with_bias', obj.with_bias));
            classifier.train(hyper, train_x, train_y);

            classifier.parameters.new_class_label = obj.new_class; % required, so that evaluator knows new class
            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = no_transfer_all_data(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "Batch" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            [prior_train_x, prior_train_y] = obj.get_data(obj.n_classes, split_ix, obj.prior_sample_slice);

            train_x = arrayfun(@(i) [prior_train_x{i}; train_x{i}], 1:numel(prior_train_x), 'UniformOutput', 0);
            train_y = [prior_train_y; train_y];

            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {obj.kernel_type};
            hyperr.gamma = obj.gamma_range;

            disp('tuning hyperparameters...');
            tuner = HyperSearch(MulticlassRLS(struct('optimize_l_tuning', true, 'with_bias', obj.with_bias)), hyperr);
            hyper = tuner.search(train_x, train_y)

            disp('training...');
            classifier = MulticlassRLS(struct('optimize_l_tuning', false, 'with_bias', obj.with_bias));
            classifier.train(hyper, train_x, train_y);

            classifier.parameters.new_class_label = obj.new_class; % required, so that evaluator knows new class
            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = priors_and_new(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "Source+1" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {obj.kernel_type};
            hyperr.gamma = obj.gamma_range;

            parameters.new_class_label = obj.new_class;
            parameters.with_bias = obj.with_bias;
            parameters.priors = obj.get_priors(split_ix);
            parameters.optimize_l_tuning = true;

            tic
            disp('tuning hyperparameters...');
            classifier = SimpleNplusOne(parameters);
            tuner = HyperSearch(classifier, hyperr);
            hyper = tuner.search(train_x, train_y);
            toc
            disp(hyper)

            disp('training...');
            classifier = SimpleNplusOne(parameters);
            classifier.train(hyper, train_x, train_y);

            classifier.parameters.new_class_label = obj.new_class; % required, so that evaluator knows new class
            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = source_plus_one_hingel(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "Source+1 with hinge loss" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {obj.kernel_type};
            hyperr.gamma = obj.gamma_range;

            parameters.new_class_label = obj.new_class;
            parameters.with_bias = obj.with_bias;
            parameters.priors = obj.get_priors(split_ix);
            parameters.optimize_l_tuning = true;

            [prior_train_x, prior_train_y] = obj.get_data(obj.n_classes, split_ix, obj.prior_sample_slice);

            disp('tuning hyperparameters...');
            parameters.binary_classifier_class = 'Libsvm';
            cl = MulticlassOneVsRest(parameters);
            tuner = HyperSearch(cl, hyperr);
            hyper = tuner.search(prior_train_x, prior_train_y);

            disp('training source...');
            parameters.binary_classifier_class = 'Libsvm';
            source = MulticlassOneVsRest(parameters);
            source.train(hyper, prior_train_x, prior_train_y);

            parameters.source = source;
            classifier = SourcePlusOneHingeL(parameters);
            classifier.train(hyper, train_x, train_y);

            classifier.parameters.new_class_label = obj.new_class; % required, so that evaluator knows new class
            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = nplusone(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "MULTIpLE" given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            nplusone_parameters.iterations = obj.beta_iterations;
            nplusone_parameters.strongly_convex_objective = false;
            nplusone_parameters.new_class_label = obj.new_class;
            nplusone_parameters.priors = obj.get_priors(split_ix);

            disp('Reusing priors'' C');
            hyper.c_val = nplusone_parameters.priors.model.hyperparams.c_val;
            hyper.new_c_val = nplusone_parameters.priors.model.hyperparams.c_val;
            hyper.kernel_type = nplusone_parameters.priors.model.hyperparams.kernel_type;
            hyper.gamma = nplusone_parameters.priors.model.hyperparams.gamma;
            hyper.beta_l = 1;

            classifier = MULTIpLE(nplusone_parameters);

            disp('training...');
            classifier.train(hyper, train_x, train_y);
            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = mktl(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "MKTL" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            parameters.new_class_label = obj.new_class;
            parameters.with_bias = obj.with_bias;
            parameters.priors = obj.get_priors(split_ix);

            hyper = struct();
            hyper.c_val = parameters.priors.model.hyperparams.c_val;
            hyper.p_val = 1.5;
            hyper.gamma = obj.gamma_range{1};
            hyper.kernel_type = obj.kernel_type;

            disp('training...');
            classifier = MKTL(parameters);
            classifier.train(hyper, train_x, train_y);

            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = multikt_ova(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "MultiKT-OVA" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            parameters.new_class_label = obj.new_class;
            parameters.with_bias = obj.with_bias;

            % Predicting from priors once for efficiency reasons
            priors = obj.get_priors(split_ix);
            parameters.priors = priors;
            [junk, parameters.priors_margins_train] = priors.predict(train_x);
            [junk, parameters.priors_margins_test] = priors.predict(test_x);

            hyper = struct();
            hyper.c_val = parameters.priors.model.hyperparams.c_val;
            hyper.gamma = obj.gamma_range{1};
            hyper.kernel_type = obj.kernel_type;

            % Precomputing kernel matrices -- will use over all OVA instances
            parameters.Ktest = GenericClassifier.compute_mean_kernel_combinations_eff_(...
                hyper, train_x, test_x)';
            parameters.Ktrain = GenericClassifier.compute_mean_kernel_combinations_eff_(...
                hyper, train_x, train_x)';

            disp('training...');
            parameters.binary_classifier_class = 'MultiKT';
            classifier = MulticlassOneVsRest(parameters);
            classifier.train(hyper, train_x, train_y);

            rs = classifier.evaluate(test_x, test_y);
        end

        function rs = pmtsvm_ova(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "PMT-SVM-OVA" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
            parameters.new_class_label = obj.new_class;
            parameters.with_bias = obj.with_bias;
            parameters.cue_index = obj.restrict_to_cue_ix;

            % Getting prior data instead of priors (PmtSvm specific)
            [parameters.prior_train_x, parameters.prior_train_y] ...
                = obj.get_data(obj.n_classes, split_ix, obj.prior_sample_slice);

            % Tuning C for PMT-SVM: linear multiclass RLS on prior training data
            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {'linear'};
            classifier = MulticlassRLS(struct('optimize_l_tuning', true));
            tuner = HyperSearch(classifier, hyperr);
            disp('tuning hyperparameters for PMT-SVM...');
            tic; hyper = tuner.search({parameters.prior_train_x{parameters.cue_index}}, parameters.prior_train_y); toc
            disp(hyper)

            disp('training...');
            max_rs = [];
            max_pmt_gamma = 0;
            % Tuning gamma for PMT-SVM
            for pmt_gamma=[0,0.1,10]
                parameters.binary_classifier_class = 'PmtSvm';
                classifier = MulticlassOneVsRest(parameters);
                hyper.pmt_gamma = pmt_gamma;
                classifier.train(hyper, train_x, train_y);
                rs = classifier.evaluate(test_x, test_y);
                if isempty(max_rs)
                    max_rs = rs;
                    max_pmt_gamma = pmt_gamma; end
                if rs.test_accuracy > max_rs.test_accuracy
                    max_rs = rs;
                    max_pmt_gamma = pmt_gamma; end
            end

            max_pmt_gamma
            rs = max_rs;
        end

        function rs = mtradaboost_ova(obj, task_name, train_x, train_y, test_x, test_y, split_ix)
        % evaluates "MultiSourceTrAdaBoost" baseline given:
        % train_x -- training samples cell array, each cell contains row-sample matrix, which
        %            corresponds to a single feature (cue)
        % train_y -- train label vector
        % test_x -- testing samples (same format as in train_x)
        % test_y -- test label vector
        % split_ix -- current split index
        %
        % Returns result structure containing test accuracy, accuracy over N classes, over N+1 class,
        % confusion matrix, etc.
        %
            parameters.new_class_label = obj.new_class;
            parameters.with_bias = obj.with_bias;
            parameters.iterations = 10;

            % Getting prior data instead of priors (mtradaboost specific)
            [prior_train_x, prior_train_y] = obj.get_data(obj.n_classes, split_ix, obj.prior_sample_slice);

            % Stacking all features
            test_x_stack = [];
            train_x_stack = [];
            prior_train_x_stack = [];
            for i=1:numel(train_x)
                train_x_stack = [train_x_stack, train_x{i}];
                prior_train_x_stack = [prior_train_x_stack, prior_train_x{i}];
                test_x_stack = [test_x_stack, test_x{i}];
            end

            parameters.prior_train_x = {prior_train_x_stack};
            parameters.prior_train_y = prior_train_y;

            % Tuning C for mtradaboost: linear multiclass RLS on prior training data
            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {'linear'};
            classifier = MulticlassRLS(struct('optimize_l_tuning', true));
            tuner = HyperSearch(classifier, hyperr);
            disp('tuning hyperparameters for MultisourceTrAdaBoost...');
            tic; hyper = tuner.search(parameters.prior_train_x, parameters.prior_train_y); toc

            disp(hyper)

            disp('training...');
            classifier = MultisourceTrAdaBoost(parameters);
            data_obj.train_x = {train_x_stack};
            data_obj.train_y = train_y;
            data_obj.test_x = {test_x_stack};
            classifier.train_model(hyper, data_obj);

            rs = classifier.evaluate({test_x_stack}, test_y);
        end

        function results = evaluate(obj, task_name, task)
        % Evaluates a given task over obj.splits data splits.
        % task must be a MATLAB function handle of a prototype:
        % fn(task_name, train_x, train_y, test_x, test_y, split)
        % Look at no_transfer(...) as an example.
        %
        % Returns a cell matrix of result structures containing test accuracy,
        % accuracy over N classes, over N+1 class, confusion matrix, etc.
        %
            disp(sprintf('>>>>>>>>>>>>>>>>>>> %s <<<<<<<<<<<<<<<<<<', task_name));

            results = {};

            hyperr.c_val = obj.C_range;
            hyperr.kernel_type = {obj.kernel_type};
            hyperr.gamma = obj.gamma_range;

            for split=obj.splits % data splits
                [test_x, test_y] = obj.get_data(obj.get_classes, split, obj.test_sample_slice);

                for slice_i=1:numel(obj.train_sample_slice) % train progression
                    [train_x, train_y] = obj.get_data(obj.get_classes, split, obj.train_sample_slice{slice_i});

                    % evaluating the task
                    rs = task(task_name, train_x, train_y, test_x, test_y, split);
                    results{split, slice_i} = rs;

                    fprintf('Split: %d; Samples: %d; N+1 Accuracy: %.2f; N Accuracy: %.2f; +1 Accuracy: %.2f\n', ...
                            split, numel(obj.train_sample_slice{slice_i}), rs.test_accuracy, rs.n_accuracy, rs.npp_accuracy);

                    if obj.show_conf_matrix
                        disp('Confusion matrix:');
                        disp(rs.conf_mat);
                    end
                end
            end
        end

        function rs = evaluate_task(obj, task_name)
            if strfind(task_name, 'nplusone')
                rs = obj.evaluate(task_name, @obj.nplusone);
            elseif strfind(task_name, 'priors_and_new')
                rs = obj.evaluate(task_name, @obj.priors_and_new);
            elseif strfind(task_name, 'no_transfer_all_data')
                rs = obj.evaluate(task_name, @obj.no_transfer_all_data);
            elseif strfind(task_name, 'no_transfer')
                rs = obj.evaluate(task_name, @obj.no_transfer);
            elseif strfind(task_name, 'priors')
                rs = obj.eval_priors;
            elseif strfind(task_name, 'mktl')
                rs = obj.evaluate(task_name, @obj.mktl);
            elseif strfind(task_name, 'multikt_ova')
                rs = obj.evaluate(task_name, @obj.multikt_ova);
            elseif strfind(task_name, 'pmtsvm_ova')
                rs = obj.evaluate(task_name, @obj.pmtsvm_ova);
            elseif strfind(task_name, 'mtradaboost_ova')
                rs = obj.evaluate(task_name, @obj.mtradaboost_ova);
            elseif strfind(task_name, 'sourcepp_hingel')
                rs = obj.evaluate(task_name, @obj.source_plus_one_hingel);
            end
        end

        function [x_, y_, classes_rp] = get_data(obj, classes, split, slice, reduce_old_classes_to_num)
        %
        % Gets the data provided class indices, split number, slice to pick from randomly sampled data.
        %
        % Returns:
        % x_ -- features, cell array of row-sample matrices
        % y_ -- vector of labels
        % classes_rp -- random permutations of sample incides for all classes provided (cell array of vectors)
        %
            if isempty(obj.dataset)
                obj.dataset = load(obj.dataset_file);
            end

            % Random seed is fixed for reproducibility
            if strcmp(version('-release'), '2009a')
                randn('seed', 100+split);
            else
                rng(100+split);
            end
            x_ = {}; y_ = [];

            classes_rp = {};

            for c=classes'
                if c <= numel(obj.dataset.samples_per_class)
                    rp = randperm(obj.dataset.samples_per_class(c));
                    if rp
                        slice_ = slice;
                        if nargin == 5 && c ~= obj.new_class
                            slice_ = slice(1:reduce_old_classes_to_num);
                        end

                        ix = find(obj.dataset.y==c);
                        sample_ix = ix(rp(slice_));

                        y_ = [y_; obj.dataset.y(sample_ix)];

                        if isempty(x_)
                            x_ = cell(1, numel(obj.dataset.x)); end

                        for feature_i = 1:numel(obj.dataset.x)
                            x_{feature_i} = [x_{feature_i}; obj.dataset.x{feature_i}(sample_ix, :)];
                        end

                        classes_rp{end+1} = rp(slice_);
                    end
                end
            end

            if ~isempty(obj.restrict_to_cue_ix)
                if obj.restrict_to_cue_ix ~= -1
                    x_ = {x_{obj.restrict_to_cue_ix}};
                else
                    % Stacking all features
                    disp('Stacking features');
                    x_stack = [];
                    for i=1:numel(x_)
                        x_stack = [x_stack, x_{i}];
                    end
                    x_ = {x_stack};
                end
            end
        end

        function classes = get_classes(obj)
        % Get all classes IDs in the dataset.
        %
            if isempty(obj.dataset)
                obj.dataset = load(obj.dataset_file);
            end
            classes = unique(obj.dataset.y);
        end

        function [classifier, ds] = get_priors(obj, split, disposable)
        % Trains and stores OR loads source classifier provided data split.
        %
        % Returns a classifier, instance of class GenericClassifier.
        %
            if nargin < 3
                disposable = false; end % disposable = never load, don't save

            [junk, dataset_name] = fileparts(obj.dataset_file);
            [junk, support_dataset_name] = fileparts(obj.support_dataset_file);

            if support_dataset_name
                prior_file = fullfile(obj.prior_location, ...
                                      ['priors_', dataset_name, '_', support_dataset_name, '_balanced_split_', ...
                                    num2str(split), '_except_', num2str(obj.new_class), '.mat']);
                disp(sprintf('Using support data for priors from "%s"', support_dataset_name));
            else
                prior_file = fullfile(obj.prior_location, ...
                                      ['priors_', dataset_name, '_split_', ...
                                    num2str(split), '_except_', num2str(obj.new_class), '.mat']);
            end

            if ~exist(prior_file, 'file') || disposable
                disp('priors not found. training...');

                [train_x, train_y] = obj.get_data(obj.n_classes, split, obj.prior_sample_slice);

                %%% Supportive sources (not used in current evaluations)
                if obj.support_dataset_file
                     d = NplusoneBenchmark;
                     d.dataset_file = obj.support_dataset_file;
                     [supp_x, supp_y] = d.get_data(d.get_classes, split, obj.prior_sample_slice);
                     for i=1:numel(train_x)
                         train_x{i} = [train_x{i}; supp_x{i}];
                     end

                     if ~isempty(intersect(unique(train_y), unique(supp_y)))
                         warning('Priors'' and support priors'' classes intersect.');
                     end

                     train_y = [train_y; supp_y];
                end
                %%%

                hyperr.c_val = obj.C_range;
                hyperr.kernel_type = {obj.kernel_type};
                hyperr.gamma = obj.gamma_range;

                classifier = MulticlassRLS(struct('optimize_l_tuning', true));
                tuner = HyperSearch(classifier, hyperr);

                disp('tuning hyperparameters...');
                tic;
                hyper = tuner.search(train_x, train_y);
                toc

                disp(hyper)

                disp('training...');
                classifier = MulticlassRLS;
                classifier.train(hyper, train_x, train_y);

                if ~disposable
                    save(prior_file, 'classifier'); end
            else
                disp('loading priors...');
                classifier = load(prior_file);
                classifier = classifier.classifier;
                disp('loaded.');
            end
        end
    end

    methods (Static)
        function rs_set = run(task)
        % Evaluates specified task by name over obj.splits and
        % in leave-one-class out manner.
        %
        % Returns a cell array of result structure cell matrices
        % (see evaluate(...) for details).
        %
            e = NplusoneBenchmark;

            classes = e.get_classes;

            rs_set = {};

            for c=classes'
                e.n_classes = setdiff(classes, c);
                e.new_class = c;

                fprintf('Leaving out class: %d\n', c);
                rs = e.evaluate_task(task);
                rs_set{end+1} = rs;
            end
        end

        function show
            function acc = get_acc(rs)
                mx = cellfun(@(y) cellfun(@(x) x.test_accuracy, y), rs, 'UniformOutput', 0);
                split_train_class = reshape([mx{:}], size(mx{1}, 1), size(mx{1}, 2), numel(mx));
                acc = mean(mean(split_train_class, 3), 1);
            end

            tic
            figure(1); hold on;

            no_transfer_rs = NplusoneBenchmark.run('no_transfer');
            batch_rs = NplusoneBenchmark.run('no_transfer_all_data');
            source_one_rs = NplusoneBenchmark.run('priors_and_new');
            source_one_hinge_rs = NplusoneBenchmark.run('sourcepp_hingel');
            multikt_ova_rs = NplusoneBenchmark.run('multikt_ova');
            mktl_ova_rs = NplusoneBenchmark.run('mktl');
            mtradaboost_ova_rs = NplusoneBenchmark.run('mtradaboost_ova');
            nplusone_rs = NplusoneBenchmark.run('nplusone');
            toc

            h1 = NplusoneBenchmark.plot(5:5:20, get_acc(no_transfer_rs), 1);
            h2 = NplusoneBenchmark.plot(5:5:20, get_acc(batch_rs), 2);
            h3 = NplusoneBenchmark.plot(5:5:20, get_acc(source_one_rs), 3);
            h8 = NplusoneBenchmark.plot(5:5:20, get_acc(source_plus_one_hingel), 8);
            h4 = NplusoneBenchmark.plot(5:5:20, get_acc(multikt_ova_rs), 4);
            h5 = NplusoneBenchmark.plot(5:5:20, get_acc(mktl_ova_rs), 5);
            h6 = NplusoneBenchmark.plot(5:5:20, get_acc(mtradaboost_ova_rs), 6);
            h7 = NplusoneBenchmark.plot(5:5:20, get_acc(nplusone_rs), 7);

            legend([h1, h2, h3, h8, h4, h5, h6, h7], ...
                   'No transfer', 'Batch', 'Source+1', 'Source+1 (hinge)', 'MultiKT-OVA', 'MKTL', 'MultisourceTrAdaBoost-OVA', 'MULTIpLE');
            title('Caltech-256 5 related classes');
        end

        function h = plot(x, y, index)
            lines = {'-k', '-or', '-xm', '-sc', '-ob', '-xg', '-sm', '-kr', '-oy'};
            h = plot(x, y, lines{index}, 'LineWidth', 2 );
            set(gca,'YGrid','on');
            set(gca,'XGrid','on');
            set(gca,'XTick',x, 'FontSize', 16, 'fontweight','b')
            set(gca,'YTick',0.1:0.05:1, 'FontSize', 16, 'fontweight','b');
            xlabel('Number of Samples', 'FontSize', 16, 'fontweight','b');
            ylabel('N+1 accuracy', 'FontSize', 16, 'fontweight','b');
        end
    end
end
