function [maxRs PP W W_v]=stats(RT1,T,T_v,P,grps)
%% selecting PC subspaces
% Input
%      RT: the resolutions of all the PCˇˇsubspaces
%      T: The scores of the training set after the rotation and  projection   
%      T_v: The scores of the test set after the rotation and  projection
%      P: the indices of all PC subspaces
%      grps: the structure of each class
% Output
%      maxRs: the Rs calculated after accumulation  
%      PP: the indices of the accepted PC subspaces
%      W: the scores of the accepted PCs in the training set
%      W_v: the scores of the accepted PCs in the test set

RTT = RT1; 
RT = []; PP = []; Ind = [];
m1 = find(RT1==max(RTT));
a0 = RT1(m1); T0 = T(:,m1); 
RT1(m1) = []; T(:,m1) = []; W = T0; PP = P(m1,:); P(m1,:) = [];
T_v0 = T_v(:,m1); T_v(:,m1) = []; W_v = T_v0; 
sumR = zeros(size(T(:,1)));

%% 
for i = 1:size(RT1,1)
i
d = size(T,2);
 for j = 1:d
   T2(:,j) = T0+T(:,j);
   sR = (T2(:,j)-min(T2(:,j)))./(max(T2(:,j))-min(T2(:,j)));
   [rt meanR1 meanR2 stdR1 stdR2] = dispsumRs(sR,grps);
   RT = [RT;rt];
 end
 n1 = find(RT>a0);
 if isempty(n1)==1, maxRs = a0; return
 else  ind = find(RT==min(RT(n1)));
       a0 = RT(ind); maxRs = a0;
       RT = []; W = [W T(:,ind)]; T(:,ind) = [];
       PP = [PP;P(ind,:)]; P(ind,:) = [];
       W_v =[W_v T_v(:,ind)]; T_v(:,ind) = [];
       T0 = T2(:,ind); T2 = [];
 end
end



