%
% MulticlassRLS (extends GenericClassifier)
%
% Multiclass kernelized regularized least squares.
% Classes are coded using OVA coding.
%
% 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 MulticlassRLS < GenericClassifier
    properties
        Q=[];
        L=[];
        Ktest=[];
    end

    methods
        function obj = MulticlassRLS(parameters)
        % Constructs MulticlassRLS
        %
        % (optional, default=false) parameters.with_bias -- wether to model bias term (b)
        %
            if nargin < 1
                parameters = struct;
            end

            obj@GenericClassifier(parameters);

            warning('error', 'MATLAB:nearlySingularMatrix');
        end

        function flag = with_bias(obj)
            flag = isfield(obj.parameters, 'with_bias') && obj.parameters.with_bias;
        end

        function flag = with_weighted_classes(obj)
            flag = isfield(obj.parameters, 'with_weighted_classes') && obj.parameters.with_weighted_classes;
        end

        function flag = with_optimize_l_tuning(obj)
            flag = isfield(obj.parameters, 'optimize_l_tuning') && obj.parameters.optimize_l_tuning;
        end

        function train_model(obj, hyperparameters, data_obj)
        % hyperparameters.kernel_type
        % hyperparameters.gamma
        % hyperparameters.c_val
        % data_obj.train_x - cell array with row-sample cues (kernels averaged)
        % data_obj.train_y - labels

            [train_y, class_map] = obj.coding_ova(data_obj.train_y);

            if obj.with_weighted_classes
                weights = zeros(size(data_obj.train_y));
                for l=1:numel(class_map)
                    weights(data_obj.train_y==class_map(l)) = numel(data_obj.train_y) / ...
                        (numel(class_map) * sum(data_obj.train_y==class_map(l)));
                end
            else
                weights = ones(size(data_obj.train_y));
            end

            Lh = (1/hyperparameters.c_val) * diag(1./weights);

            try
                if ~obj.with_optimize_l_tuning
                    K = GenericClassifier.compute_mean_kernel_combinations_eff_(hyperparameters, data_obj.train_x, data_obj.train_x);

                    if ~obj.with_bias
                        obj.model.alpha = (K + Lh) \ train_y;
                    else
                        G = [K + Lh                 ones(size(K, 1), 1);
                             ones(1, size(K, 2))   0                  ];
                        params = G \ [train_y; zeros(1, size(train_y, 2))];
                        obj.model.alpha = params(1:end-1, :);
                        obj.model.bias = params(end, :);
                    end
                else
                    if isempty(obj.Q)
                        %disp('Computing kernel and diagonalizing...');
                        %tic;
                        K = GenericClassifier.compute_mean_kernel_combinations_eff_(hyperparameters, data_obj.train_x, data_obj.train_x);
                        [obj.Q,obj.L] = eig(K);
                        %toc
                    end

                    %disp(sprintf('Computing parameters for C: %f', hyperparameters.c_val));
                    %tic;

                    Kinv = (obj.Q * inv(obj.L + Lh) * obj.Q');
                    if ~obj.with_bias
                        obj.model.alpha = Kinv * train_y;
                    else
                        % Using matrix block inversion
                        z=inv(-sum(sum(Kinv)));
                        A = Kinv + Kinv * ones(size(Kinv, 1), 1) * z * ones(1, size(Kinv, 2)) * Kinv;
                        B = - Kinv * ones(size(Kinv, 1), 1) * z;
                        C = -z * ones(1, size(Kinv, 2)) * Kinv;
                        G = [A B; C z];
                        params = G * [train_y; zeros(1, size(train_y, 2))];
                        obj.model.alpha = params(1:end-1, :);
                        obj.model.bias = params(end, :);
                    end

                    %toc
                end
            catch e
                error('MulticlassRLS:terminated', 'MulticlassRLS:terminated');
            end

            obj.model.SVs = data_obj.train_x;
            obj.model.train_y = data_obj.train_y;
            obj.model.encoded_train_y = train_y;
            obj.model.hyperparams = hyperparameters;
            obj.model.class_map = class_map;
        end

        function [pred, margins] = predict(obj, test_x)
            if (obj.with_optimize_l_tuning && isempty(obj.Ktest)) || ~obj.with_optimize_l_tuning
                obj.Ktest = GenericClassifier.compute_mean_kernel_combinations_eff_(obj.model.hyperparams, test_x, obj.model.SVs);
            end

            if ~obj.with_bias
                margins = obj.Ktest * obj.model.alpha;
            else
                margins = bsxfun(@plus, obj.Ktest * obj.model.alpha, obj.model.bias);
            end

            [junk, class_idx] = max(margins, [], 2);
            pred = obj.model.class_map(class_idx); % Decoding to class label values
        end

        function clear(obj)
            obj.Q = [];
            obj.L = [];
            obj.Ktest = [];
        end
    end

    methods (Static)
        function test
            rng(100);

            [train_x, train_y, class_stat] = SyntheticData.gen_synth_data(5, 20, [1 1 1 1 1 1 1], 3);
            [test_x, test_y] = SyntheticData.gen_synth_data_(class_stat, 100, [1 1 1 1 1 1 1]);

            tic;
            cl = MulticlassRLS(struct('optimize_l_tuning', true, 'with_bias', true));
            hyperr.c_val = 10.^(-20:20);
            hyperr.gamma = {2.^(-10:10)};
            hyperr.kernel_type = {'rbf'};
            tuner = HyperSearch(cl, hyperr);
            hyper = tuner.search({train_x, train_x}, train_y)
            toc

            cl = MulticlassRLS(struct('optimize_l_tuning', false, 'with_bias', true));
            cl.train(hyper, {train_x, train_x}, train_y);
            cl.evaluate({test_x, test_x}, test_y)
        end
    end

end
