function [PercentCorrect, MuMeasure, FeatureMatrix, ElapsedTime] = ...
    test_satellite(m,Nruns,PlotFlag,InitCondition,Method)

%Input arguments:
%---------------
%
%m - (1 x 1) the size of feature space (chosen by the user)
%
%Nruns - (1 x 1) the number of optimization runs
%  
%PlotFlag - (1 x 1) if PlotFlag = 1 create plot, otherwise don't
%  
%InitCondition - (string): 'random' random initial feature extraction matrix
%                             'lda' linear discriminant analysis matrix
%                             'che' the matrix proposed by Loog & Duin in
%                             IEEE TPAMI, 26, 732-739, 2004.   
%  
%Method - (string): 'tr' for the trust-region method
%                   'cg' for the conjugate gradient method

to = clock;

if PlotFlag == 1  
  close all
end

load ./Satimage.mat

Train = Data(1:4435,1:36);
Valid = Data(4436:6435,1:36);
TrainMemship = Data(1:4435,37);
ValidMemship = Data(4436:6435,37);

[T Mu] = ...
    ida_feature_extraction_matrix(m,Train,TrainMemship, ...
                                  num2str(Method),10^(-8)*[1 1], ...
                                  2000,num2str(InitCondition),Nruns);

MuMeasure = Mu;
FeatureMatrix = T;

%feature extraction done; classification ensues

FeatureTrain = (Train * T');

if PlotFlag == 1
  
  C = length(unique(TrainMemship));

  Kolor = colormap;
  Ind = 1:floor(64/C):64;

  for i = 1:C
    ind = find(TrainMemship == (i-1));
  
    if m == 1
      ph = plot(FeatureTrain(ind), zeros(length(ind),1),'ko');
      hold on
    elseif m == 2
      ph = plot(FeatureTrain(ind,1), FeatureTrain(ind,2),'ko');
      hold on
    elseif m >= 3
      ph = plot3(FeatureTrain(ind,1),FeatureTrain(ind,2), ...
                 FeatureTrain(ind,3),'ko');
      hold on
    end
    set(ph,'Color',Kolor(Ind(i),:),'MarkerSize',4, ...
           'MarkerFaceColor', Kolor(Ind(i),:))
  end
  
  if m > 1
    set(gca,'DataAspectRatio',[1 1 1])
  end
  
end

FeatureValid = (Valid * T');
Class_l = classify(FeatureValid,FeatureTrain,TrainMemship,'linear', ...
                 'empirical');
Class_q = classify(FeatureValid,FeatureTrain,TrainMemship,'quadratic', ...
                 'empirical');
PercentCorrect_l = ...
    100 * sum((Class_l - ValidMemship) == 0)/length(ValidMemship);

PercentCorrect_q = ...
    100 * sum((Class_q - ValidMemship) == 0)/length(ValidMemship);

fprintf('Percent correct (linear): %1.3f\n',PercentCorrect_l)
fprintf('Percent correct (quadratic): %1.3f\n',PercentCorrect_q)

et = etime(clock,to);
disp(['Elapsed time: ' num2str(et)])
ElapsedTime = et;

PercentCorrect = [PercentCorrect_l PercentCorrect_q];


