ufunction  [TP4_label,FP4_label,TP,TN,FP,FN] = MK(X,Y,p1,p,digest_index_temp,theta)
          p2 = 0;  pi0 = 0.9;  p_half = 0.12; p_none = 0.005;
          
      Xp = preprocess(X);
      delta_xcorr = Xp(:,1) + Xp(:,2);
      K_1 = TrainKernelMatrix('gaussian',1,Xp);
      K_02 = TrainKernelMatrix('gaussian',0.5,Xp);
      K_05 = TrainKernelMatrix('gaussian',0.2,Xp);


      True = find(Y==1);
      False = find(Y==-1);
      TrueNumber = length(True);
      FalseNumber = length(False);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%    stage   1   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%  centerlized process  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      K_1_False = K_1(Y==-1,Y==-1); 
      K_05_False = K_05(Y==-1,Y==-1); 
      K_02_False = K_02(Y==-1,Y==-1);
      K_1SumFalseTemp = sum(K_1_False);
      K_05SumFalseTemp = sum(K_05_False);
      K_02SumFalseTemp = sum(K_02_False);
      
      K_1SumFalse = sum(K_1SumFalseTemp);
      K_05SumFalse = sum(K_05SumFalseTemp);
      K_02SumFalse = sum(K_02SumFalseTemp);
      
      K_1_Mix = K_1(Y==1,Y==-1);
      K_05_Mix = K_05(Y==1,Y==-1);
      K_02_Mix = K_02(Y==1,Y==-1);
      
      K_1DistanceTrue = ones(TrueNumber,1) + K_1SumFalse/(FalseNumber*FalseNumber)*...
                         ones(TrueNumber,1) - 2*sum(K_1_Mix,2)/FalseNumber;
      K_05DistanceTrue = ones(TrueNumber,1) + K_05SumFalse/(FalseNumber*FalseNumber)*...
                             ones(TrueNumber,1) - 2*sum(K_05_Mix,2)/FalseNumber;
      K_02DistanceTrue = ones(TrueNumber,1) + K_02SumFalse/(FalseNumber*FalseNumber)*...
                              ones(TrueNumber,1) - 2*sum(K_02_Mix,2)/FalseNumber;
      DistanceTrue_temp = [K_1DistanceTrue,K_05DistanceTrue,K_02DistanceTrue];
      DistanceTrue = max(DistanceTrue_temp,[],2);
      % compte the distance between positive samples and center of negative samples
% DistanceN = ones(NegativeNumber,1) + SumNegative/(NegativeNumber*NegativeNumber)*...
%             ones(NegativeNumber,1) - 2*sum(Knegative,2)/NegativeNumber;
      DistanceTrueTemp = sort(DistanceTrue);
      DistanceThreshold = DistanceTrueTemp(round(p1*TrueNumber));        
      index1 = DistanceTrue > DistanceThreshold; 
      True1 = True(index1);
      c = [True1;False];   
% drop the positive samples with Distance<DistanceThreshold

      Y1 = Y(c);  
      m1 = length(c);
   
      
   MK = theta(1)*K_1+theta(2)*K_05+theta(3)*K_02;
   K1 = MK(c,c);
   model_precomputed1 = svmtrain(Y1, [(1:m1)', K1], '-t 4');
  % train the first classifier on sample S1
   [out1, accuracy_P1, dec_values_P1] = svmpredict(Y1, [(1:m1)', K1'], model_precomputed1);
  % accuracy on training set
   OutT1 = out1(1:length(True1));  
   OutF1 = out1(length(True1)+1:end);
   IndexTrueTemp = OutT1==1; IndexFalseTemp = OutF1==1;
   True1Number = length(True1);
   TP1_label = True1(IndexTrueTemp);
   TP1 = length(TP1_label);
   FP1_label = False(IndexFalseTemp);
   FP1 = length(find(OutF1==1));
   TN1 = length(find(OutF1==-1)); 
   FN1 = length(find(OutT1==-1));
   E_FDR1 = (pi0*length(True1)*FP1)/(FalseNumber*TP1);
   MS_FDR1 = 2*FP1/(FP1+TP1);                                   
   TPR1 = TP1/(TP1+FN1);
   % accuracy on true points 
   % according to prior knowledge, we expect it approach 10%
   SPC1 = TN1/(FP1+TN1);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%% stage 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%% denoise process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   TT = dec_values_P1(1:length(True1));
   aa = find (TT>p2);
   True2 = True1(aa); 
   cc = [True2;False];
   True2Number = length(True2);
   
 
   K2 = MK(cc,cc);
   Y2 = Y(cc);
   m2 = length(Y2);
     
   model_precomputed2 = svmtrain(Y2, [(1:m2)', K2], '-t 4');
   [out2, accuracy_P2, dec_values_P2] = svmpredict(Y2, [(1:m2)', K2], model_precomputed2);
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % predict the discriminate score of new sample ues the model model_prexomputed2 directly
   % 
   OutT2 = out2(1:length(True2));    OutF2 = out2(length(True2)+1:end);
   IndexTrueTemp = OutT2==1;  IndexFalseTemp = OutF2==1;
   TP2_label = True2(IndexTrueTemp);
   TP2 = length(TP2_label);
   FP2_label = False(IndexFalseTemp);
   FP2 = length(find(OutF2==1));
   TN2 = length(find(OutF2==-1));
   FN2 = length(find(OutT2==-1));
   E_FDR2 = (pi0*length(True2)*FP2)/(FalseNumber*TP2);
   MS_FDR2 = 2*FP2/(FP2+TP2);                                   
   TPR2 = TP2/(TP2+FN2);
   % accuracy on true points 
   % according to prior knowledge, we expect it approach 10%
   SPC2 = TN2/(FP2+TN2);
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%   Stage 3  %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%%%%%%%%%%%%%%%   Validation   %%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   OutValue_T2 = dec_values_P2(1:length(aa));
   OutValue_F = dec_values_P2(length(aa)+1:end);
   OutValue_F_Temp = sort(OutValue_F,'descend');
   TP3_label = cell(1,length(p)); %cell array storage the label of TP_3
   FP3_label = cell(1,length(p));
   TP3 = zeros(1,length(p)); FP3 = zeros(1,length(p));
   TN3 = zeros(1,length(p)); FN3 = zeros(1,length(p));
   TP3_full_label = cell(1,length(p)); FP3_full_label = cell(1,length(p));
   TP4_half_label = cell(1,length(p)); FP4_half_label = cell(1,length(p));
   TP3_none_label = cell(length(p),1); FP3_none_label = cell(length(p),1);
   TP4_none_label = cell(1,length(p)); FP4_none_label = cell(1,length(p));
   TP4 = zeros(1,length(p)); FP4 = zeros(1,length(p)); TN4 = zeros(1,length(p)); FN4 = zeros(1,length(p)); 
   TP3_half_label = cell(1,length(p)); FP3_half_label = cell(1,length(p));
   TP4_label = cell(1,length(p)); FP4_label = cell(1,length(p));
   for j = 1: length(p)
         p3 = p(j);
            
   t0 = OutValue_F_Temp(floor(length(False)*p3)+1);
   IndexTrueTemp = OutValue_T2>t0;   IndexFalseTemp = OutValue_F>t0;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % predict the new data according to the threshold t0
   True3Number = True2Number;
   
   
   TP3_label{1,j} = True2(IndexTrueTemp);
   FP3_label{1,j} = False(IndexFalseTemp);
   
   
   TP3(1,j) = length(TP3_label{1,j});
   FP3(1,j) = length(find(OutValue_F>t0));
   
   TN3 = length(False)*ones(1,length(p)) - FP3;
   FN3 = length(True2)*ones(1,length(p)) - TP3;
   
   %E_FDR3 = (pi0*length(True2)*FP3)/(FalseNumber*TP3);
   %MS_FDR3 = 2*FP3/(FP3+TP3);                                   
   %TPR3 = TP3/(TP3+FN3);
   % accuracy on true points 
   % according to prior knowledge, we expect it approach 10%
   %SPC3 = TN3/(FP3+TN3);
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%% add digest information %%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %    digest_index_temp
   
   TP3_full_label{1,j} = TP3_label{1,j}(digest_index_temp(TP3_label{1,j})==3);
   FP3_full_label{1,j} = FP3_label{1,j}(digest_index_temp(FP3_label{1,j})==3);
   
   
   TP3_half_label{1,j} = TP3_label{1,j}(digest_index_temp(TP3_label{1,j})==1|digest_index_temp(TP3_label{1,j})==2);
   FP3_half_label{1,j} = FP3_label{1,j}(digest_index_temp(FP3_label{1,j})==1|digest_index_temp(FP3_label{1,j})==2);
   
   
   TP3_none_label{1,j} = TP3_label{1,j}(digest_index_temp(TP3_label{1,j})==0);
   FP3_none_label{1,j} = FP3_label{1,j}(digest_index_temp(FP3_label{1,j})==0);
   
%    OutValueT2_half = OutValue_T2(digest_index_temp(TP3_label{1,j})==1|digest_index_temp(TP3_label{1,j})==2);
%    OutValueF_half = OutValue_F(digest_index_temp(FP3_label{1,j})==1|digest_index_temp(FP3_label{1,j})==2);
%    OutValue_half_temp = sort([OutValueT2_half; OutValueF_half],'descend');
   
   
%    if floor(length(OutValue_half_temp)*p_half) == 0  
%                    TP4_half_label{1,j} = [];
%                    FP4_half_label{1,j} = []; 
%    else
    
%    t0 = OutValue_half_temp(floor(length(OutValue_half_temp)*p_half));
%    IndexTrueHalf_Temp = OutValueT2_half>t0;   IndexFalseHalf_Temp = OutValueF_half>t0;
     delta_xcorr_half_temp = sort(delta_xcorr([TP3_half_label{1,j};FP3_half_label{1,j}]),'descend');
      if floor(length(delta_xcorr_half_temp)*p_half) == 0  
                    TP4_half_label{1,j} = [];
                    FP4_half_label{1,j} = []; 
      else
    
    t0 = delta_xcorr_half_temp(floor(length(delta_xcorr_half_temp)*p_half)+1);
    IndexTrueHalf_Temp = delta_xcorr(TP3_half_label{1,j})>t0;   
    IndexFalseHalf_Temp =  delta_xcorr(FP3_half_label{1,j})>t0;
       
   TP4_half_label{1,j} = TP3_half_label{1,j}(IndexTrueHalf_Temp);
   FP4_half_label{1,j} = FP3_half_label{1,j}(IndexFalseHalf_Temp);
       end
   

   
% %    OutValueT2_none = OutValue_T2(digest_index_temp(TP3_label{1,j})==0);
% %    OutValueF_none = OutValue_F(digest_index_temp(FP3_label{1,j})==0);
% %    OutValue_none_temp = sort([OutValueT2_none; OutValueF_none],'descend');
% %    
% %   
% %    if floor(length(OutValue_none_temp)*p_none) == 0  
% %       TP4_none_label{1,j} = []; 
% %       FP4_none_label{1,j} = [];
% %    else
% %    t0 = OutValue_none_temp(floor(length(OutValue_none_temp)*p_none));
% %    IndexTrueNone_Temp = OutValueT2_none>t0;   IndexFalseNone_Temp = OutValueF_none>t0;

 delta_xcorr_none_temp = sort(delta_xcorr([TP3_none_label{1,j};FP3_none_label{1,j}]),'descend');
      if floor(length(delta_xcorr_none_temp)*p_half) == 0  
                    TP4_none_label{1,j} = [];
                    FP4_none_label{1,j} = []; 
      else
    
    t0 = delta_xcorr_none_temp(floor(length(delta_xcorr_none_temp)*p_none)+1);
    IndexTrueNone_Temp = delta_xcorr(TP3_none_label{1,j})>t0;   
    IndexFalseNone_Temp =  delta_xcorr(FP3_none_label{1,j})>t0;
       
   TP4_none_label{1,j} = TP3_none_label{1,j}(IndexTrueNone_Temp);
   FP4_none_label{1,j} = FP3_none_label{1,j}(IndexFalseNone_Temp);
       end
   
   TP4_label{1,j} = [TP3_full_label{1,j}; TP4_half_label{1,j}; TP4_none_label{1,j}];
   FP4_label{1,j} = [FP3_full_label{1,j}; FP4_half_label{1,j}; FP4_none_label{1,j}];
   
   TP4(1,j) = length(TP4_label{1,j});
   FP4(1,j) = length(FP4_label{1,j});  
   TN4(1,j) = length(False) - FP4(1,j);
   FN4(1,j) = length(True2) - TP4(1,j);
   
   end
   
% %    TP = [TP1 TP2 TP3 TP4];
% %    TN = [TN1 TN2 TN3 TN4];
% %    FP = [FP1 FP2 FP3 FP4];
% %    FN = [FN1 FN2 FN3 FN4];
TP = TP4;
TN = TN4;
FP = FP4;
FN = FN4;


  