% kernel mutual subspace method (KMSM.m) by Seiji Hotta
% 19, June, 2010

clear all;
close all;

ri=5; % dimensionality of input subspace
rd=6; % dimensionality of dictionaly subspace
scale=1.0; %scale of kernel parameter

nclass=30; % total number of classes
load('doubutsu.mat'); % load training and test data

ntrai=size(trai,2);
ntest=size(test,2);

% mean distance for kernel parameter
meandist=zeros(nclass,1);
%% kernel computing
for jj = 1 : nclass
  X=trai(jj).X;
  [dim ndata]=size(X);
  K=zeros(ndata);
  for ii = 1 : ndata
    for ll = ii : ndata
      K(ii,ll)=norm(X(:,ii)-X(:,ll)).^2;
      K(ll,ii)=K(ii,ll);
    end
  end
  trai(jj).K=K;
  meandist(jj)=1./(mean(K(:)));
end
% for Gaussian kernel parameter
sig=mean(meandist).*scale;
% Gramian matrix for training data
for jj = 1 : nclass, trai(jj).K=exp(-sig.*trai(jj).K);,end

for jj = 1 : ntest
  X=test(jj).X;
  [dim ndata]=size(X);
  K=zeros(ndata);
  for ii = 1 : ndata
    for ll = ii : ndata
      K(ii,ll)=norm(X(:,ii)-X(:,ll)).^2;
      K(ll,ii)=K(ii,ll);
    end
  end
  test(jj).K=K;
  % Gramian matrix for test data
  test(jj).K=exp(-sig.*test(jj).K);
end
%% eigenvalue decomposition for each class
for jj = 1 : nclass
  K=trai(jj).K;
  [eig_vec eig_val]=eig(K);
  [value index]=sort(diag(eig_val),'descend');
  mat_rank=rank(K);
  for l = 1 : mat_rank
    eig_vec(:,index(l))=eig_vec(:,index(l))./sqrt(value(l));
  end
  trai(jj).V=eig_vec(:,index(1:mat_rank));
end
%% eigenvalue decomposition for each test sample
for jj = 1 : ntest
  K=test(jj).K;
  [eig_vec eig_val]=eig(K);
  [value index]=sort(diag(eig_val),'descend');
  mat_rank=rank(K);
  for l = 1 : mat_rank
    eig_vec(:,index(l))=eig_vec(:,index(l))./sqrt(value(l));
  end
  test(jj).U=eig_vec(:,index(1:mat_rank));
end
%% kernel MSM
sim=zeros(nclass,1);
% confusion matrix
CONF=zeros(nclass);
for query = 1 : ntest
  U=test(query).U(:,1:ri);
  Q=test(query).X;
  nQ=size(Q,2);
  for jj = 1 : nclass
    V=trai(jj).V(:,1:rd);
    X=trai(jj).X;
    nX=size(X,2);
    % (test images x class images) kernel matrix
    K=zeros(nX,nQ);
    for ii = 1 : nX
      for ll = 1 : nQ
        K(ii,ll)=norm(X(:,ii)-Q(:,ll)).^2;
      end
    end
    K=exp(-sig.*K);
    S=U'*K'*V*V'*K*U;
    [eig_vec eig_val]=eig(S);
    sim(jj)=max(diag(eig_val));
  end
  [value index]=max(sim);
  CONF(test(query).label,index)=CONF(test(query).label,index)+1;
end
accuracy=(sum(diag(CONF))./ntest).*100;
fprintf(2,'ri=%d rd=%d accuracy=%3.1f %%\n',ri,rd,accuracy);
