clear all;, close all;

r=13;      % dimensionality of each subspace
imgnum=1;  % test sample number for displaying Uc

nclass=10; % number of classes

% loading data-file
load('./usps_resampled/usps.mat');
[d,ndata]=size(trai);

% normalization
for ii = 1 : ndata
  trai(:,ii)=trai(:,ii)./norm(trai(:,ii));
end

figure(1),clf
%% forming subspaces
for c = 1 : nclass
  X=trai(:,find(trai_label==c-1));
  D=X*X';
  maxr=rank(D);
  [eig_vec, eig_val]=eig(D);
  [value index]=sort(diag(eig_val),'descend');
  C(c).U=eig_vec(:,index(1:maxr));
  C(c).eigval=value(1:maxr);
  %% displaying U of each class
  for ii = 1 : 10
    IMG=reshape(C(c).U(:,ii),[16 16]);
    IMG=IMG-min(IMG(:));, IMG=IMG./max(IMG(:));
    figure(1),subplot(10,10,(c-1).*10+ii),imshow(IMG);
  end
  fprintf(2,'class %d ... OK\n',c-1);
end

figure(2),clf
%% displaying Uc
Q=test(:,imgnum)./norm(test(:,imgnum));
IMG=reshape(Q,[16 16]);
IMG=IMG-min(IMG(:));, IMG=IMG./max(IMG(:));
figure(2),subplot(2,10,5),imshow(IMG);
title('test sample');
for c = 1 : nclass
  a=C(c).U(:,1:r)'*Q;
  IMG=reshape(C(c).U(:,1:r)*a,[16 16]);
  IMG=IMG-min(IMG(:));, IMG=IMG./max(IMG(:));
  figure(2),subplot(2,10,10+c),imshow(IMG);
  s=sprintf('class %d',c-1);, title(s);
end

%% CLAFIC
S=zeros(nclass,1);
CONF=zeros(nclass);
tic
for ii = 1 : ndata
  test(:,ii)=test(:,ii)./norm(test(:,ii));
  for c = 1 : nclass
    S(c)=norm(C(c).U(:,1:r)'*test(:,ii));
  end
  [value index]=max(S);
  CONF(index,test_label(ii)+1)=CONF(index,test_label(ii)+1)+1;
end
finish=toc;
accuracy=(sum(diag(CONF))./ndata).*100;
fprintf(2,'accuracy=%3.2f\n',accuracy);
fprintf(2,'classification time per sample [s]: %f\n',finish./ndata);
