%
% MULTIpLE (extends MulticlassRLS)
%
% MULTIpLE algorithm.
%
% Usage example:
%
% source = obj.get_priors(split_ix);
%
% % settings
% parameters.iterations = obj.beta_iterations;
% parameters.new_class_label = obj.new_class;
% parameters.priors = source;

% % hyperparameters
% hyper.c_val = parameters.source.model.hyperparams.c_val;
% hyper.new_c_val = parameters.source.model.hyperparams.c_val;
% hyper.kernel_type = parameters.source.model.hyperparams.kernel_type;
% hyper.gamma = parameters.priors.source.hyperparams.gamma;
% hyper.beta_l = 1;

% % Instantiating, training & testing, here:
% % 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
% classifier = MULTIpLE(parameters);
% classifier.train(hyper, train_x, train_y);
% rs = classifier.evaluate(test_x, test_y);
%
%
% 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 MULTIpLE < MulticlassRLS
    properties
        invGd = [];
        loss_track = [];
        new_class_pos = [];
        old_class_pos = [];
        non_supp_pos = [];
    end

    methods
        function obj = MULTIpLE(parameters)
        % Constructs MULTIpLE
        %
        % Required 'parameters' structure values:
        %
        % parameters.priors                              -- 'priors' classifier
        % parameters.iterations                          -- number of iterations
        % parameters.new_class_label                     -- +1 class label (numeric)
        % (optional) parameters.zero_B                   -- whether to use beta value = 0 if < 0 heuristic
        % (optional) parameters.use_priors_for_priors    -- use don't use new prior models from prior classes
        %
            parameters.with_bias = true;
            parameters.strongly_convex_objective = false;

            obj@MulticlassRLS(parameters);
        end

        function num = get_num_priors(obj)
            num = numel(obj.parameters.priors.model.class_map);
        end

        function prior_margins = predict_from_priors(obj, samples)
            [junk, prior_margins] = obj.parameters.priors.predict(samples);
        end

        function precompute_loo_model(obj)
            Gd = diag(obj.model.Ginv);

            params = obj.model.Ginv * ...
                     [obj.model.encoded_train_y_orig; ...
                      zeros(1, size(obj.model.encoded_train_y_orig, 2))]; % alpha'
            obj.model.alpha_loo = params(1:end-1, :);
            obj.model.bias_loo = params(end, :);

            y_from_priors = obj.predict_from_priors(obj.model.SVs);
            params = obj.model.Ginv * ...
                     [y_from_priors; zeros(1, size(y_from_priors, 2))]; % alpha''
            obj.model.alpha_prior = params(1:end-1, :);
            obj.model.bias_prior = params(end, :);

            obj.invGd = diag(1./Gd(1:end-1));
            obj.model.A = (obj.invGd * obj.model.alpha_prior)'; % precomputed for efficiency
        end

        function optimize_beta(obj)
            n = size(obj.model.train_y, 1); % number of samples
            m = numel(obj.model.class_map); % number of classes
            k = obj.get_num_priors;         % number of priors

            beta = zeros(k, 1);
            obj.loss_track = [];

            % LOO estimates for new priors (doesn't change during optimization)
            loo_old = obj.model.encoded_train_y_orig(:, obj.old_class_pos) - ...
                      obj.invGd * (obj.model.alpha_loo(:, obj.old_class_pos) - obj.model.alpha_prior);

            % Matrix for all LOO estimates
            combined_loo = zeros(n, m);
            combined_loo(:, obj.old_class_pos) = loo_old;

            sum_beta = zeros(k, 1);
            sum_step = 0;

            for iter=1:obj.parameters.iterations

                % New class LOO estimates given beta
                loo_new = obj.model.encoded_train_y_orig(:, obj.new_class_pos) - ...
                          obj.invGd * (obj.model.alpha_loo(:, obj.new_class_pos) - obj.model.alpha_prior * beta);

                % Forming (sample X model) LOO pred. matrix
                combined_loo(:, obj.new_class_pos) = loo_new;

                beta_subg = zeros(k, 1);
                sum_losses = 0;

                sum_losses2 = 0;

                for i=1:n % iterating over samples
                    yi = obj.model.train_y(i) == obj.model.class_map;

                    %%%%%%%
                    % Computing Crammer-Singer loss given sample
                    losses = 1 + combined_loo(i, ~yi) - combined_loo(i, yi);
                    max_r2 = max(losses);

                    if max_r2 > 0
                        sum_losses2 = sum_losses2 + max_r2;
                    end
                    %%%%%%%

                    %%%%%%%
                    cost = obj.model.sample_weights(i);
                    %%%%%%%

                    yi = find(yi);

                    if yi ~= obj.new_class_pos
                        pw_loss = 1 + combined_loo(i, obj.new_class_pos) - combined_loo(i, yi);
                        pw_loss = cost * pw_loss;

                        if pw_loss > 0
                            sum_losses = sum_losses + pw_loss;

                            beta_subg = beta_subg + cost * obj.model.A(:, i);
                        end
                    else
                        pw_loss = 1 + combined_loo(i, setdiff(1:m, obj.new_class_pos)) - combined_loo(i, obj.new_class_pos);
                        pw_loss = cost * pw_loss;
                        [max_r, r] = max(pw_loss);

                        if max_r > 0
                            sum_losses = sum_losses + max_r;

                            beta_subg = beta_subg - cost * obj.model.A(:, i);
                        end
                    end
                end

                obj.loss_track(end+1) = sum_losses;

                if obj.parameters.strongly_convex_objective
                    eta = sqrt(iter) - sqrt(iter-1);
                    beta = beta * (1 - 1/iter) - 1/(iter * eta) * beta_subg;
                else
                    beta = beta - beta_subg * (1 / (n*sqrt(iter))); % updating beta
                end

                % Non-negative heuristic step
                if isfield(obj.parameters, 'zero_B') && obj.parameters.zero_B
                    beta = (beta > 0) .* beta;
                end

                if isfield(obj.parameters, 'l1_projection') && obj.parameters.l1_projection
                    beta = ProjectOntoL1Ball(beta, 1/obj.model.hyperparams.beta_l); % l1 projection
                else
                    beta = min(1, 1/ ( norm(beta) * obj.model.hyperparams.beta_l )) * beta; % l2 projection
                end

                if obj.parameters.strongly_convex_objective
                    sum_beta = sum_beta + beta;
                    beta = (1/iter) * sum_beta;
                end
            end

            obj.model.beta = beta;
        end

        function train_model(obj, hyperparameters, data_obj)

            % Querying priors
            prior_margins = obj.predict_from_priors(data_obj.train_x);

            % Detecting which labels are from "supportive priors" (non-evaluated)
            supp_priors = setdiff(unique(obj.parameters.priors.model.class_map), unique(data_obj.train_y));

            % Coding labels (OVA)
            [train_y, class_map] = obj.coding_ova(data_obj.train_y, supp_priors);
            obj.model.encoded_train_y_orig = train_y;
            obj.model.train_y = data_obj.train_y;

            % New (+1), old (priors) and supportive prior class positions
            obj.new_class_pos = find(class_map==obj.parameters.new_class_label);
            obj.old_class_pos = setdiff(1:numel(class_map), obj.new_class_pos);
            obj.non_supp_pos = arrayfun(@(x) ismember(x, unique(data_obj.train_y)), class_map);

            % Non-uniform class ratio handling
            weights = zeros(size(data_obj.train_y));
            for l=1:numel(class_map)
                num_samples = numel(data_obj.train_y);
                weights(data_obj.train_y==class_map(l)) = num_samples / ...
                    (numel(class_map) * sum(data_obj.train_y==class_map(l)));
            end
            obj.model.sample_weights = weights;

            K = obj.compute_mean_kernel_combinations(hyperparameters, data_obj.train_x);

            % Separate C tuning (if required)
            if isfield(hyperparameters, 'new_c_val')
                new_class_I = diag(obj.parameters.new_class_label == data_obj.train_y);
                old_class_I = diag(obj.parameters.new_class_label ~= data_obj.train_y);
                L = (1/hyperparameters.c_val) * ( diag(1./weights) .* old_class_I ) ...
                    + (1/hyperparameters.new_c_val) * ( diag(1./weights) .* new_class_I );
            else
                L = (1/hyperparameters.c_val) * diag(1./weights);
            end

            obj.model.SVs = data_obj.train_x;
            obj.model.hyperparams = hyperparameters;
            obj.model.class_map = class_map;

            % Inverting K + lamb*I
            if obj.with_bias
                G = [K + L                 ones(size(K, 1), 1);
                     ones(1, size(K, 2))   0                  ];
            else
                G = K + L;
            end

            try
                obj.model.Ginv = inv(G);
            catch
                obj.model.Ginv = pinv(G);
                warning('Gram matrix was badly scaled. Used pinv instead of inv.');
            end

            % Optimizing beta
            obj.precompute_loo_model;
            obj.optimize_beta;

            % Solving RLS
            train_y(:, obj.old_class_pos) = train_y(:, obj.old_class_pos) - prior_margins;
            train_y(:, obj.new_class_pos) = train_y(:, obj.new_class_pos) - prior_margins * obj.model.beta;

            if obj.with_bias
                params = obj.model.Ginv * [train_y; ...
                                    zeros(1, size(train_y, 2))];
                obj.model.alpha = params(1:end-1, :);
                obj.model.bias = params(end, :);
            else
                obj.model.alpha = obj.model.Ginv * double(train_y);
            end
        end

        function [pred, margins] = predict(obj, samples)
            K = obj.compute_mean_kernel_combinations(obj.model.hyperparams, samples, obj.model.SVs);
            margins = bsxfun(@plus, K * obj.model.alpha, obj.model.bias);
            prior_margins = obj.predict_from_priors(samples);

            margins(:, obj.old_class_pos) = margins(:, obj.old_class_pos) + prior_margins;
            margins(:, obj.new_class_pos) = margins(:, obj.new_class_pos) + (prior_margins * obj.model.beta);

            margins = margins(:, obj.non_supp_pos);
            obj.model.predict_class_map = obj.model.class_map(obj.non_supp_pos);

            [junk, class_idx] = max(margins, [], 2);
            pred = obj.model.predict_class_map(class_idx);
        end

        function clear(obj)
            obj.model = [];
        end

    end
end
