%% step2 テストパターンを直交展開
r=30; % 部分空間の次元数+1，最大100まで
x=2;  % テストサンプルの番号 (1---2007)

if exist('W')==0 % すでに計算済みか？
  %% 訓練データと未知データの読み込み
  %% D：次元数×訓練データ数の行列, trai_label  ：訓練データのラベル
  %% Q：次元数×テストデータ数の行列, test_label：未知データのラベル
  D=load('./USPS/trai_data.txt');, trai_label=load('./USPS/trai_label.txt');
  Q=load('./USPS/test_data.txt');, test_label=load('./USPS/test_label.txt');
  %% Dim：次元数，xx_num：データ数
  [Dim,trai_num]=size(D);,[Dim,test_num]=size(Q);
  %% 各パターンのノルムを1に正規化
  for i = 1 : trai_num, D(:,i)=D(:,i)./norm(D(:,i));, end
  for i = 1 : test_num, Q(:,i)=Q(:,i)./norm(Q(:,i));, end
  %% 各クラスで部分空間を求める
  W=zeros(Dim,100,10); % 各クラスの直交基底 次元数×基底数×クラス数
  for j = 0 : 9
    X=D(:,find(trai_label==j)); % ラベルがjの訓練パターンをXに格納
    C=X*X';                     % 自己相関行列を計算
    [eig_vec, eig_val]=eig(C);  % 固有ベクトル，固有値を計算
    [value index]=sort(-diag(eig_val)); % 固有値を降順に並べ替え
    W(:,:,j+1)=eig_vec(:,index(1:100)); % 固有値の大きい上位100本に対応する固有ベクトルをWに格納
    fprintf('class %d ... OK\n',j);
  end
end
%% 直交基底による線型近似
figure(1), clf, hold on;
IMG=zeros(16,16);
% 未知パターンの表示
for i = 1 : 16, IMG(i,:)=Q((i-1).*16+1:i.*16,x)';,end
IMG=IMG-min(min(IMG));, IMG=IMG./max(max(IMG));
figure(1),subplot(2,10,5),imshow(IMG);
title('test sample');
% 線型近似パターンの表示
for j = 0 : 9
  c=W(:,1:r,j+1)'*Q(:,x); % 係数を計算
  QA=W(:,1:r,j+1)*c;       % 線型近似パターンの作成
  for i = 1 : 16, IMG(i,:)=QA((i-1).*16+1:i.*16,1)';,end
  IMG=IMG-min(min(IMG));, IMG=IMG./max(max(IMG));
  figure(1),subplot(2,10,11+j),imshow(IMG);,
  s=sprintf('class %d',j);, title(s);
end
fprintf('次元数 r=%d\n',r);
