clear all
num_class=79;                                      % 類別數
num_inclass=3;                                     % 每一類的樣本數
num_all=num_class*num_inclass;                     % 所有樣本數
m=100;                                             % 期望的行數 
n=60;                                             % 期望的列數
M=m*n;

min=1;
resultofratio=zeros(1,15);

%load K;                                  % 讀入訓練樣本的核函數(核函數為K)
%load K_new;                               % 讀入訓練樣本的核函數(核函數為K_new)
load K_rbf_minus;                        % 讀入訓練樣本的核函數(核函數為K_rbf_minus)

%load K_test;                             % 讀入測試樣本的核函數(核函數為K_test)
%load K_test_new;                          % 讀入測試樣本的核函數(核函數為K_test_new)
load K_test_rbf_minus;                   % 讀入測試樣本的核函數(核函數為K_test_rbf_minus)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PCA降維

K_one=ones(num_all,num_all);

% tempfait=K-(1/num_all*K*K_one)-(1/num_all*K_one*K)+(1/(num_all*num_all)*K_one*K*K_one);
%tempfait=K_new-(1/num_all*K_new*K_one)-(1/num_all*K_one*K_new)+(1/(num_all*num_all)*K_one*K_new*K_one);
tempfait=K_rbf_minus-(1/num_all*K_rbf_minus*K_one)-(1/num_all*K_one*K_rbf_minus)+(1/(num_all*num_all)*K_one*K_rbf_minus*K_one);

[V_tempfait,D_tempfait]=eig(tempfait);

Vtempst=zeros(num_all,(num_all-num_class));

for i=1:num_all-1
    k=i;
    for j=(i+1):num_all
        if(D_tempfait(k,k)<D_tempfait(j,j))           % 從大到小排
            k=j;
        end
        if(k~=i)
        
            temp_d=D_tempfait(i,i);
            temp_v=V_tempfait(:,i);
            
            D_tempfait(i,i)=D_tempfait(k,k);
            V_tempfait(:,i)=V_tempfait(:,k);
            
            D_tempfait(k,k)=temp_d;
            V_tempfait(:,k)=temp_v;
            
        end
    end
end
        

for i=1:(num_all-num_class)
    Vtempst(:,i)=V_tempfait(:,i);
end


% Y=K-(1/num_all*K_one*K);
%Y=K_new-(1/num_all*K_one*K_new);
Y=K_rbf_minus-(1/num_all*K_one*K_rbf_minus);
 
Z=Vtempst'*Y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sum_Z=zeros((num_all-num_class),1);

fai_b=zeros((num_all-num_class),num_class);
fai_w=zeros((num_all-num_class),num_all);

sum_inclass=zeros((num_all-num_class),num_class);
mean_inclass=zeros((num_all-num_class),num_class);


for i=1:num_all
   sum_Z=sum_Z+Z(:,i);                                                          % 所有訓練樣本的和 
end
mean_Z=sum_Z/num_all;                                                           % 計算所有樣本的平均值mean_Z

for i=1:num_class
    for j=1:num_inclass
        sum_inclass(:,i)=sum_inclass(:,i)+Z(:,(i-1)*num_inclass+j);
    end
    mean_inclass(:,i)=sum_inclass(:,i)/num_inclass;                             % 計算類內樣本的平均值mean_inclass
end


for i=1:num_class    
   fai_b(:,i)=mean_inclass(:,i)-mean_Z;                                          % 計算類間離散度sb的簡單算法                   
end

sb=fai_b*(fai_b');    

for i=1:num_class
    for j=1:num_inclass
       fai_w(:,(i-1)*num_inclass+j)=Z(:,(i-1)*num_inclass+j)-mean_inclass(:,i);   % 計算類內離散度sw的簡單算法
    end
end

sw=fai_w*(fai_w');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[V,D]=eig(inv(sw)*sb);       % 計算在sw的非零空間下，類內距離最小，類間距離最大的投影空間。 D降序排列

% dim_fisher=num_class-1;
dim_fisher_g=[10:10:150];

for index=1:15
    dim_fisher=dim_fisher_g(index);

project_fisher=zeros((num_all-num_class),dim_fisher);

for i=1:dim_fisher
    project_fisher(:,i)=V(:,i);
end

X_trainvector=project_fisher'*Z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K_test_one=ones(num_all,num_all);

% Y_test=K_test-(1/num_all*K_test_one*K_test);                          
%Y_test=K_test_new-(1/num_all*K_test_one*K_test_new);
Y_test=K_test_rbf_minus-(1/num_all*K_test_one*K_test_rbf_minus);
 
Z_test=Vtempst'*Y_test;

X_testvector=project_fisher'*Z_test;                              % 計算測試樣本的特征向量

for i=1:num_class
      
    for j=1:num_all
        temp=(X_testvector(:,i)-X_trainvector(:,j))'*(X_testvector(:,i)-X_trainvector(:,j));    % 計算測試樣本與訓練樣本的最近鄰距離
        distance(1,j)=sqrt(temp);
    end
    
    for j=1:num_all
        if  distance(1,j)<distance(1,min)
            min=j;                                                          % 找到與測試樣本距離最近的訓練樣本的編號
        end
    end
    
    remember(1,i)=min;                                                      % 將與相應的測試樣本距離最近的訓練樣本的編號存入remember一維行向量中
    
end

result_class=zeros(1,num_class);
num_recognition=0;

for i=1:num_class
    
    remainder=mod(remember(1,i),num_inclass);                               % 與相應的測試樣本距離最近的訓練樣本的編號除以每一類的樣本數的余數       
    quotient=fix(remember(1,i)/num_inclass);                                % 與相應的測試樣本距離最近的訓練樣本的編號除以每一類的樣本數的商       
    
    if(remainder==0)
        result_class(1,i)=quotient;                                         % 計算測試樣本的類別 （余數為零）
    else
        result_class(1,i)=quotient+1;                                       % 計算測試樣本的類別 （余數不為零）
    end
    
    if(result_class(1,i)==i)
        num_recognition=num_recognition+1;                                  % 計算識別正確的個數
    end
    
end

ratio_recognition(index)=100*num_recognition/num_class;                                % 計算識別率
end
ratio_recognition;
plot(10:10:150,ratio_recognition(1:15))
axis([0 150 0 100])
hold on