
function [G, DD, L, K]  = iakdaUpdate(K11, L11, K21, K22, T, opts)

% function [G, DD, L, K]  = iakdaUpdate(K11, L11, K21, K22, T, opts)
%
% Update step of incremental accelerated kernel discriminant analysis (IAKDA)
%
% Description:
%  Updating step of incremental and accelerated computation of expansion
%  coefficient matrix G and eigenvalue diagonal matrix DD for
%  dimensionality reduction
%
% IN: 
%   K11: n1 x n1 kernel matrix of training data at time t-1
%   K22: n2 x n2 kernel matrix of training data at time t
%   K21: n2 x n1 kernel matrix of training data Phi2' * Phi1
%   L11: Cholesky factor of K11
%   T: class labels
%   opts : additional options
%     regAlpha: regularization parameter for kernel matrix
%
% OUT:
%  G : N x D projection matrix
%  DD: eigenvalue diagonal matrix
%  L: Cholesky factor of the overall kernel data matrix K at time t
%
% Author: Nikolaos Gkalelis - CERTH-ITI
% Email: gkalelis@iti.gr
%
% First creation (Version 0.0), 1 Apr 2017.
%


fprintf('iakdaUpdate>>> Entering\n');

Clbl = unique(T); % class labels
if ~issorted(Clbl);
    fprintf('>>> Class labels are not sorted, sorting...\n');
    Clbl = sort(Clbl);
end

C = length(Clbl); % number of classes
n1 = size(K11,1);
n2 = size(K22,1);
N = n1 + n2; % number of observations

reguAlpha = opts.regAlpha; % regularization parameter for kernel matrices

% regularize K22
K22 = max(K22, K22');
K22 = K22 + eye(n2, n2) * reguAlpha;

% compute number of observations in each class and subclass
Ni = zeros(C,1); % number of observations per class
for i=1:C
    Ni(i) = sum(T == Clbl(i));
end

% construct elementary matrices
n_sqrt = zeros(C,1); % square root of cardinality vector
for i=1:C
    n_sqrt(i) = sqrt(Ni(i));
end

% construct core matrix
O = eye(C) -  n_sqrt* n_sqrt.' * (1/N); % (n_sqrt.' * n_sqrt) = N

N2 = zeros(C,C); % inverse of square root of cardinality matrix
for i=1:C
    N2(i,i) = 1/sqrt(Ni(i));
end

I1 = zeros(N,C); % I1 indicator matrix
for n=1:N
    for i=1:C % get class ID
        if T(n) == Clbl(i)
            I1(n,i) = 1;
            break;
        end
    end
end

opts_da.disp=0;
[KSI, DD] = eigs(full(O), rank(O), 'LM', opts_da); % usually rank(O) = D

TH = I1* N2* KSI;

L21 = K21/L11';
%KL22 = K22 - L21 * L21';
KL22 = K22 - (K21 / K11) * K21';

% regularize Gram matrix
KL22 = max(KL22, KL22'); %KL22 = (KL22 + KL22')/2;
KL22 = KL22 + eye(n2, n2) * reguAlpha/2;
%KL22 = KL22 + diag( 1 ./ diag(L21) ) * reguAlpha/2;

L22 = chol(KL22,'lower');

L = [L11, zeros(n1, n2); L21, L22];

K = [K11, K21'; K21, K22];

W = L'\(L\TH); % solve linear system

% normalize to unite matrix
tmpNorm = sqrt(sum((W'*K).*W',2));
G = W./repmat(tmpNorm',size(W,1),1);
DD = diag(DD);

fprintf('iakdaUpdate>>> Exiting\n');


