%% step1 固有値と直交基底の観察
clear all;
%% 訓練データと未知データの読み込み
%% 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); % 各クラスの直交基底 次元数×基底数×クラス数
figure(1), clf, axis square;
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に格納
  figure(1),subplot(2,5,j+1), plot(-value(1:100)); % 固有値をプロット
  s=sprintf('class %d', j);, title(s);
  xlabel('次元数');,ylabel('固有値');
  fprintf('class %d ... OK\n',j);
end
%% 各クラスの上位10本の固有ベクトルを表示
figure(2), clf, axis square;
IMG=zeros(16,16); %% 16 x 16 ピクセルの画像
count=1;
for j = 0 : 9
  for k = 1 : 10
    %% ベクトルの二次元化
    for i = 1 : 16, IMG(i,:)=W((i-1).*16+1:i.*16,k,j+1)';,end
    IMG=IMG-min(min(IMG));, IMG=IMG./max(max(IMG));
    %% 画像の表示
    figure(2),subplot(10,10,count),imshow(IMG);
    count=count+1;
  end
end
