
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Title: StatisCLAS: A Package of Classification Methods   %%%%
%%%%  Author: David Casado de Lucas                            %%%%
%%%%  Version: 29th May 2013                                   %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%  Beginning Setting Parameters and Variables  %%%%%%%%%%%%%%%%%%%%

verbose = true;
dataFigures = true;
   numberElementsToPlot = ceil(m1param*0.05);
   modelsJointly = false;
   percentMargin = 5;
resultsFigures = true;
saveResults = false;

%%%%%%%%%%%%%%%%  Ending Setting Parameters and Variables  %%%%%%%%%%%%%%%%%%%%


clear all


%%%%%%%%%%%%%%%%  Beginning Statistical Parameters and Variables  %%%%%%%%%%%%%%%%%%%%

%%% REAL DATA

% Training samples
  allGroupFD1 = ;
  allGroupFD2 = ;

% New elements to be classified
  newDataFD = ;

% Is it necessary to smooth the curves? Take care to not destroying important details...


% %%% An example with pseudo-real data
%   numberExercise = 4;
%     paramVector1 = [4 1];  paramVector2 = [4.5 1]; 
%   T = 50;
%   % Training samples
%     N1= 48; N2 = 51;
%     allGroupFD1 = dataExercisesFD(numberExercise,1,N1,T,paramVector1);
%     allGroupFD2 = dataExercisesFD(numberExercise,2,N2,T,paramVector2);
%   % New data
%     M1 = 3; M2 = M1;
%     newDataFD = [dataExercisesFD(numberExercise,1,M1,T,paramVector1); dataExercisesFD(numberExercise,2,M2,T,paramVector2)];
%   clear numberExercise T t N N1 N2 paramVector1 paramVector2

%%%%%%%%%%%%%%%%  Ending Data Loading  %%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%  Beginning Statistical Parameters and Variables  %%%%%%%%%%%%%%%%%%%%

numbersFDtoCLASSmethods = [1 1 1 1; 1 2 1 1; 1 3 1 1];  % Numbers of the methods determined by each row
namesFDtoCLASSmethods = {'WI', 'WD', 'dKNN'};  % Names of the global methods
  diffOrders = [];  % Vector with the differentiations or derivatives
  diffsMode = 'derivatives';  % To choose between 'differentials' or 'derivatives'
  variaCutoff = sqrt(eps);  % MV: minimum between-group variability for a variable to be considered as nonconstant
  numberNeighbours = 3;  % MV: parameter of the KNN method
  
B = ;  % Number of runs
paramOptimization = false;  % To select the parameter(s) with a cross-validation loop. It is disabled if 'length(paramValues)=1'
  paramDataReuse = false;  % Whether to use all the training data (Ni-mi) firstly to select the parameter(s) and secondly as final training data
  n1 = -1; % Training sample size. It must be lower or equal than N1-m1 (strictly lower if 'depthBasedRobustness=true')
  n2 = -1; % Training sample size. It must be lower or equal than N2-m2 (strictly lower if 'depthBasedRobustness=true')
    % If 'paramDataReuse=true', n1 and n2 are assigned, respectively, the values N1-m1 and N2-m2 (possibly less the number of nondeep elements)
    % If 'paramDataReuse=false', n1=-1 and n2=-1, then n1 and n2 are set, respectively, 50% of N1-m1 and N2-m2 (possibly less the number of nondeep elements), respectively
  % For leave-one-out cross-validation
    m1param = ; % Testing sample size in the parameter cross-validation
    m2param = m1param; % Testing sample size in the parameter cross-validation
    leaveOneOutParam = false;
    Bparam = ;  % Number of cross-validation runs to select the best number of blocks. It is changed if 'leaveOneOutParam = true'
%   % For leave-one-out cross-validation
%     m1param = 1; % Testing sample size in the parameter cross-validation
%     m2param = 1; % Testing sample size in the parameter cross-validation
%     leaveOneOutParam = true;
%     Bparam = -1;  % Number of cross-validation runs to select the best parameter values. It is changed if 'leaveOneOutParam = true'
methodSelection = false;

depthBasedRobustness = false;  % To apply the approach of identifing the deepest elements and maintain only them in the training samples
  alpha = 0.2;  % Proportion for the trim


%%%% Setting some values... %%%%
M = size(newDataFD,1);
N1 = size(allGroupFD1,1);
N2 = size(allGroupFD2,1);
m1 = M;
m2 = M;
T = size(allGroupFD1,2);
t = 0:1/(T-1):1;  % Grid of equispaced points.
if ~isequal(size(numbersFDtoCLASSmethods,1),length(namesFDtoCLASSmethods))
  disp(' ')
  disp('The length of the vector with the names is different to the number of methods, generic names will be used')
  disp(' ')
  for i1 = 1:size(numbersFDtoCLASSmethods,1)
    namesFDtoCLASSmethods{i1} = ['Method ' int2str(i1)];
  end
end
paramValues = diffOrders;
% All the depth measurements
disp(' ')
disp('Measuring the depth of the functional data...')
disp(' ')
computationalTimesDepth = zeros(size(numbersFDtoCLASSmethods,1),length(paramValues));
if depthBasedRobustness
  deepestFpositions1 = zeros(N1-ceil(alpha*N1),size(numbersFDtoCLASSmethods,1),length(paramValues));
  deepestFpositions2 = zeros(N2-ceil(alpha*N2),size(numbersFDtoCLASSmethods,1),length(paramValues));
  for p = 1:length(paramValues)
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      if isequal(i1,1) || (i1>1 && ~any(numbersFDtoCLASSmethods(i1,1)==numbersFDtoCLASSmethods(1:(i1-1),1)))
        if isequal(p,1)  % The following lines are a trick to avoid a perturbation in measuring the computational time (see the helping file)
          [litter deepestFpositions1(:,i1,p)] = trimFunction(allGroupFD1, alpha);
          [litter deepestFpositions2(:,i1,p)] = trimFunction(allGroupFD2, alpha);
        end
        startClock = clock;
        % Identifying--but not removing--the deepest functions of each group
        [litter deepestFpositions1(:,i1,p)] = trimFunction(allGroupFD1, alpha);
        [litter deepestFpositions2(:,i1,p)] = trimFunction(allGroupFD2, alpha);
        computationalTimesDepth(i1,p) = computationalTimesDepth(i1,p) + etime(clock,startClock);
      else   % To use possible previous calculations, since i1-1, i1-2,... were also considered
        posts = find(numbersFDtoCLASSmethods(1:(i1-1),1)==numbersFDtoCLASSmethods(i1,1));
        computationalTimesDepth(i1,p) = computationalTimesDepth(posts(1),p);
        deepestFpositions1(:,i1,p) = deepestFpositions1(:,posts(1),p);
        deepestFpositions2(:,i1,p) = deepestFpositions2(:,posts(1),p);
      end
    end
  end
  for i1 = 1:size(numbersFDtoCLASSmethods,1)
    disp(' ')
    disp(['Mean computational times spent calculating the depth for each parameter values, method ' namesFDtoCLASSmethods{i1}])
    disp(['   '  num2str(computationalTimesDepth(i1,:))])
    disp(' ')
  end
end
if isequal(size(numbersFDtoCLASSmethods,1),1)
  methodSelection = false;  % To disable the search of the best method
end
if methodSelection && ~paramOptimization
  paramOptimization = true;
  disp(' ')
  disp('The search of minimizing method has been disabled, since paramOptimization=false')
  disp(' ')
end
%%%%%%%%%%%%%%%%  Ending Statistical Parameters and Variables  %%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%  Beginning Main Computes  %%%%%%%%%%%%%%%%%%%%

% Hypermatrixes to store the misclassification error rates, the minimizing-power measures, the computational times and the coefficients
classificationOutputs = zeros(M*B,size(numbersFDtoCLASSmethods,1),length(paramValues)+1);  % N subhypermatrixes of height B
minimizingValuesClassificationOutputs = zeros(M*B,size(numbersFDtoCLASSmethods,1));  % N matrixes of height B
minimizingValuesMethodsClassificationOutputs = zeros(M*B,1);  % N vectors of height B
paramValuesMinpower = zeros(B,size(numbersFDtoCLASSmethods,1),length(paramValues));
methodsMinpower = zeros(B,size(numbersFDtoCLASSmethods,1));
computationalTimesRuns = zeros(B,size(numbersFDtoCLASSmethods,1),length(paramValues));
computationalTimesRunsWholeParamValues = zeros(B,size(numbersFDtoCLASSmethods,1));
coefficientsMVmethod = zeros(B,size(numbersFDtoCLASSmethods,1),length(paramValues));
computationalTimesRunsMVmethod = zeros(B,size(numbersFDtoCLASSmethods,1));
for i1 = 1:size(numbersFDtoCLASSmethods,1)
  for p = 1:length(paramValues)
    computationalTimesRuns(:,i1,p) = computationalTimesRuns(:,i1,p) + computationalTimesDepth(i1,p);
  end
  computationalTimesRunsWholeParamValues(:,i1) = computationalTimesRunsWholeParamValues(:,i1) + sum(computationalTimesDepth(i1,:));
  computationalTimesRunsMVmethod(:,i1) = computationalTimesRunsMVmethod(:,i1) + sum(computationalTimesDepth(i1,:));
end

% Runs
if verbose
  disp(' ')
  disp('START')
  disp(' ')
end
if verbose
  disp('Starting the iterations...')
  disp(' ')
end
for b = 1:B
  if verbose
    disp(['   ITERATION: ' int2str(b)]) 
  end
  positionsTrainingGroup1 = (1:N1)';
  positionsTrainingGroup2 = (1:N2)';
  % Optimization of the parameter(s)
  if paramOptimization
    if verbose
      disp('    Optimization loop...')
    end
    if paramDataReuse
      % Selection of the final training groups
      positionsFinalTrainingGroup1 = positionsTrainingGroup1;
      positionsFinalTrainingGroup2 = positionsTrainingGroup2;
    else
      if all([n1 n2]==-1)
        n1Temp = ceil(length(positionsTrainingGroup1)/2);
        n2Temp = ceil(length(positionsTrainingGroup2)/2);
        disp('      The values  n1=ceil((M1-m1)/2)  and  n2=ceil((M2-m2))/2  have temporarily been set!')
      else
        n1Temp = n1;
        n2Temp = n2;
      end
      % Selection of the final training groups
      positionsFinalTrainingGroup1 = positionsTrainingGroup1(randsample(length(positionsTrainingGroup1),n1Temp));
      positionsFinalTrainingGroup2 = positionsTrainingGroup2(randsample(length(positionsTrainingGroup2),n2Temp));
      % Rest of the data
      [litter, indexes] = setdiff(positionsTrainingGroup1,positionsFinalTrainingGroup1);
      positionsTrainingGroup1 = positionsTrainingGroup1(indexes);
      [litter, indexes] = setdiff(positionsTrainingGroup2,positionsFinalTrainingGroup2);
      positionsTrainingGroup2 = positionsTrainingGroup2(indexes);
    end
    N1param = length(positionsTrainingGroup1);
    N2param = length(positionsTrainingGroup2);
    if leaveOneOutParam && ~isequal(Bparam,N1param+N2param)
      Bparam = N1param+N2param;
      if verbose && isequal(b,1)
        disp('      The value of Bparam may been changed, since leave-one-out cross-validation is being run!')
      end
    end
    % Hypermatrixes to store the misclassification error rates and the reliability of the runs
    errorRatesGroup1param = zeros(Bparam,size(numbersFDtoCLASSmethods,1),length(paramValues));
    errorRatesGroup2param = zeros(Bparam,size(numbersFDtoCLASSmethods,1),length(paramValues));
    trustworthyRunsParam = true(Bparam,size(numbersFDtoCLASSmethods,1),length(paramValues));  % To register the reliable runs
    for b2 = 1:Bparam
      % Selection of the final testing groups
      if leaveOneOutParam
        if b2 <= N1param
          positionsTestingGroup1param = positionsTrainingGroup1(b2);
          m1param = 1;
          positionsTestingGroup2param = [];
          m2param = 0;
        else
          positionsTestingGroup1param = [];
          m1param = 0;
          positionsTestingGroup2param = positionsTrainingGroup2(b2-N1param);
          m2param = 1;
        end
      else
        positionsTestingGroup1param = positionsTrainingGroup1(randsample(N1param,m1param));
        positionsTestingGroup2param = positionsTrainingGroup2(randsample(N2param,m2param));
      end
      % Rest of the data
      [litter, indexes] = setdiff(positionsTrainingGroup1,positionsTestingGroup1param);
      positionsTrainingGroup1param = positionsTrainingGroup1(indexes);
      [litter, indexes] = setdiff(positionsTrainingGroup2,positionsTestingGroup2param);
      positionsTrainingGroup2param = positionsTrainingGroup2(indexes);
      n1param = length(positionsTrainingGroup1param); % New sample size 1
      n2param = length(positionsTrainingGroup2param); % New sample size 2
      % Hypermatrixes to store the multivariate variables
      groupMV1 = zeros(n1param,size(numbersFDtoCLASSmethods,1),length(paramValues));
      groupMV2 = zeros(n2param,size(numbersFDtoCLASSmethods,1),length(paramValues));
      groupMV1test = zeros(m1param,size(numbersFDtoCLASSmethods,1),length(paramValues));      
      groupMV2test = zeros(m2param,size(numbersFDtoCLASSmethods,1),length(paramValues));
      % Method calculations here with all possible values of the parameter(s)
      for p = 1:length(paramValues)
        for i1 = 1:size(numbersFDtoCLASSmethods,1)
          switch numbersFDtoCLASSmethods(i1,4)
            case 1  % For the mean
              infoFDreference = {};
            case 2  % For the alpha-trimmed mean
              alpha = 0.2;  % Parameter of the trim
              infoFDreference = {alpha};
          end
          switch numbersFDtoCLASSmethods(i1,1)
            case 1  % Necessary information for our discriminant variables
              infoFDtoMVmethod = {paramValues(p),diffsMode,numbersFDtoCLASSmethods(i1,3),numbersFDtoCLASSmethods(i1,4),infoFDreference};
          end
          positionsTestingGroup1paramTemp = positionsTestingGroup1param;
          positionsTestingGroup2paramTemp = positionsTestingGroup2param;
          positionsTrainingGroup1paramTemp = positionsTrainingGroup1param;
          positionsTrainingGroup2paramTemp = positionsTrainingGroup2param;
          if depthBasedRobustness
            % Removing the untrustworthy testing elements (maybe different for each method)
            [litter, indexes] = setdiff(positionsTestingGroup1paramTemp,deepestFpositions1(:,i1,p));
            positionsTestingGroup1paramTemp(indexes) = [];
            [litter, indexesTestingGroup1param] = intersect(positionsTestingGroup1param,positionsTestingGroup1paramTemp);
            [litter, indexes] = setdiff(positionsTestingGroup2paramTemp,deepestFpositions2(:,i1,p));
            positionsTestingGroup2paramTemp(indexes) = [];
            [litter, indexesTestingGroup2param] = intersect(positionsTestingGroup2param,positionsTestingGroup2paramTemp);
            % Removing the untrustworthy training elements (maybe different for each method)
            [litter, indexes] = setdiff(positionsTrainingGroup1paramTemp,deepestFpositions1(:,i1,p));
            positionsTrainingGroup1paramTemp(indexes) = [];
            [litter, indexesTrainingGroup1param] = intersect(positionsTrainingGroup1param,positionsTrainingGroup1paramTemp);
            [litter, indexes] = setdiff(positionsTrainingGroup2paramTemp,deepestFpositions2(:,i1,p));
            positionsTrainingGroup2paramTemp(indexes) = [];
            [litter, indexesTrainingGroup2param] = intersect(positionsTrainingGroup2param,positionsTrainingGroup2paramTemp);
            % Checking if, after removing some elements, this run is reliable for this method
            if isempty([positionsTestingGroup1paramTemp; positionsTestingGroup2paramTemp]) ||...
                    isempty(positionsTrainingGroup1paramTemp)  || isempty(positionsTrainingGroup2paramTemp)
              disp(['      Untrustworthy run for [b b2 g i1] = ' num2str([b b2 g i1])])
              continue  % To skip this run
            end
            startClock = clock;
            %%% From functional to multivariate data
            % [groupMV1(:,i1,p) groupMV1test(:,i1,p) groupMV2(:,i1,p) groupMV2test(:,i1,p)] =...
            [groupMV1(indexesTrainingGroup1param,i1,p) groupMV1test(indexesTestingGroup1param,i1,p)...
                groupMV2(indexesTrainingGroup2param,i1,p) groupMV2test(indexesTestingGroup2param,i1,p)] =...
                fromFDtoMV(numbersFDtoCLASSmethods(i1,1),infoFDtoMVmethod,...
                  allGroupFD1(positionsTrainingGroup1paramTemp,:),allGroupFD1(positionsTestingGroup1paramTemp,:),...
                  allGroupFD2(positionsTrainingGroup2paramTemp,:),allGroupFD2(positionsTestingGroup2paramTemp,:));
            timeFDtoMV = etime(clock,startClock);
          else
            if isequal(i1,1) || (i1>1 && ~isequal(numbersFDtoCLASSmethods(i1,1),numbersFDtoCLASSmethods(i1-1,1)))
              startClock = clock;
              %%% From functional to multivariate data
              [groupMV1(:,i1,p) groupMV1test(:,i1,p) groupMV2(:,i1,p) groupMV2test(:,i1,p)] =...
                fromFDtoMV(numbersFDtoCLASSmethods(i1,1),infoFDtoMVmethod,...
                  allGroupFD1(positionsTrainingGroup1paramTemp,:),allGroupFD1(positionsTestingGroup1paramTemp,:),...
                  allGroupFD2(positionsTrainingGroup2paramTemp,:),allGroupFD2(positionsTestingGroup2paramTemp,:));
              timeFDtoMV = etime(clock,startClock);
            else  % Previous calculations can be used, since i1-1, i1-2,... were also considered
              groupMV1(:,i1,p) = groupMV1(:,i1-1,p);
              groupMV2(:,i1,p) = groupMV2(:,i1-1,p);
              groupMV1test(:,i1,p) = groupMV1test(:,i1-1,p);
              groupMV2test(:,i1,p) = groupMV2test(:,i1-1,p);
            end
          end
          computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + timeFDtoMV;
          computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + timeFDtoMV;
          startClock = clock;
          %%% The classifications (by using the multivariate data) 
          % Univariate methods DFMp: classification based on D^p(f)
          % Option 1: Minimizing the distance to the training centroids (in this case it's not necessary to look at the sign of the variable)
          % Option 2: Looking at the sign of the variable
          outclass = 1-([groupMV1test(:,i1,p); groupMV2test(:,i1,p)] < 0);
          equalDist = ([groupMV1test(:,i1,p); groupMV2test(:,i1,p)] == 0);
          if any(equalDist)  % The elements in the frontier are allocated at random
            disp(['      For method ' namesFDtoCLASSmethods{i1} ', number of the ' int2str(sum(m1param+m2param)) ' testing elements such that x' int2str(p) ' = 0: ' int2str(sum(equalDist))])
            outclass(equalDist) = randi(2,sum(equalDist),1)-1;
          end
          errorRatesGroup1param(b2,i1,p) = mean(outclass(1:length(positionsTestingGroup1paramTemp)),1);
          errorRatesGroup2param(b2,i1,p) =...
              mean(1-outclass((length(positionsTestingGroup1paramTemp)+1):(length(positionsTestingGroup1paramTemp)+length(positionsTestingGroup2paramTemp))),1);
          computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + etime(clock,startClock);
          computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + etime(clock,startClock);
        end
      end
    end
    errorRatesParam = (m1param*errorRatesGroup1param + m2param*errorRatesGroup2param)/(m1param+m2param);    % Note that m1param = 0 if leaveOneOut = true
    % Mean error of the trustworthy runs
    meanErrorRatesParam = zeros(size(numbersFDtoCLASSmethods,1),length(paramValues));
    for p = 1:length(paramValues)
      for i1 = 1:size(numbersFDtoCLASSmethods,1)
        meanErrorRatesParam(i1,p) = mean(errorRatesParam(trustworthyRunsParam(:,i1,p),i1,p),1);
      end
    end
    % Calculations to select the parameter value(s) or the method(s)
    methsXparams = false(size(numbersFDtoCLASSmethods,1),length(paramValues));
    methodMinError = zeros(1,size(numbersFDtoCLASSmethods,1));  % To register, for each method, the minimum error rate
    minValue = zeros(1,size(numbersFDtoCLASSmethods,1)); % To register, for each method, the smallest minimizing parameter value (as a reference to calculate some quantities later)
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      methodMinError(i1) = min(meanErrorRatesParam(i1,:));  % Minimum error rate
      posts = find(meanErrorRatesParam(i1,:)==methodMinError(i1));  % Values for which the minimum is reached
      methsXparams(i1,posts) = 1;  % To enable the minimizing parameter values
      paramValuesMinpower(b,i1,posts) = 1/length(posts);  % To register the minimizing-power of the parameter values
      minValue(i1) = posts(1);  % The smallest minimizing parameter value (as a reference to calculate some quantities later)
    end
    posts = find(methodMinError==min(methodMinError));  % Methods for which the minimum-in-mean is reached    
    if methodSelection
      [litter,indexes] = setdiff(1:size(numbersFDtoCLASSmethods,1),posts);
      methsXparams(indexes,:) = zeros(length(indexes),length(paramValues));  % Disabling the other methods (for the results)
    end
    methodsMinpower(b,posts) = 1/length(posts);  % Minimizing-power of the methods: one or several
    minMethod = posts(1);  % Smallest minimizing method (as a reference to calculate some quantities later)
    if verbose
      for i1 = 1:size(numbersFDtoCLASSmethods,1)
        disp(['      Information for method ' namesFDtoCLASSmethods{i1}])
        disp(['        Mean overall error rates for the parameter values: ' '   ' num2str(meanErrorRatesParam(i1,:))])
      end
      disp(['      Minimum error rate for the methods: ' '   ' num2str(methodMinError)])
    end
    if verbose
      disp('    ...end of the optimization loop')
    end
  else
    % Selection of the final training groups
    positionsFinalTrainingGroup1 = positionsTrainingGroup1;
    positionsFinalTrainingGroup2 = positionsTrainingGroup2;
    methsXparams = true(size(numbersFDtoCLASSmethods,1),length(paramValues));
  end
  if isequal(b,1) && dataFigures  % Crude data
    groupFD1 = allGroupFD1(positionsFinalTrainingGroup1,:);
    groupFD1test = newDataFD;
    groupFD2 = allGroupFD2(positionsFinalTrainingGroup2,:);
    groupFD2test = newDataFD; % = [];
    figure  % Figures: plotting the functional data
    plotFunctions(t,groupFD1test(1:numberElementsToPlot,:),groupFD2test(1:numberElementsToPlot,:),...
        'Some testing curves and the mean of the training samples',percentMargin,modelsJointly,mean(groupFD1,1),mean(groupFD2,1))
  end
  n1 = length(positionsFinalTrainingGroup1); % New sample size 1
  n2 = length(positionsFinalTrainingGroup2); % New sample size 2
  % Hypermatrixes to store the discriminant variables and the discriminant function
  groupMV1 = zeros(n1,size(numbersFDtoCLASSmethods,1),length(paramValues));
  groupMV2 = zeros(n2,size(numbersFDtoCLASSmethods,1),length(paramValues));
  groupMV1test = zeros(m1,size(numbersFDtoCLASSmethods,1),length(paramValues));      
  groupMV2test = zeros(m2,size(numbersFDtoCLASSmethods,1),length(paramValues));
  groupY1test = zeros(m1,size(numbersFDtoCLASSmethods,1));
  groupY2test = zeros(m2,size(numbersFDtoCLASSmethods,1));
  % Method calculations here with the selected value(s) of the parameter(s)
  for p = 1:length(paramValues)
    if all(methsXparams(:,p)==0)
      continue  % Not to consider this value
    end
    % Calculating the functional data and applying the functional classification methods    
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      if all(methsXparams(i1,:)==0) || isequal(methsXparams(i1,p),0)
        continue  % Not to consider this method or this value
      end
      switch numbersFDtoCLASSmethods(i1,4)
        case 1  % For the mean
          infoFDreference = {};
        case 2  % For the alpha-trimmed mean
          alpha = 0.2;  % Parameter of the trim
          infoFDreference = {alpha};
      end
      switch numbersFDtoCLASSmethods(i1,1)
        case 1  % Necessary information for our discriminant variables
          infoFDtoMVmethod = {paramValues(p),diffsMode,numbersFDtoCLASSmethods(i1,3),numbersFDtoCLASSmethods(i1,4),infoFDreference};
      end
      positionsFinalTrainingGroup1Temp = positionsFinalTrainingGroup1;
      positionsFinalTrainingGroup2Temp = positionsFinalTrainingGroup2;
      if depthBasedRobustness
        [litter, indexes] = setdiff(positionsFinalTrainingGroup1Temp,deepestFpositions1(:,i1,p));
        positionsFinalTrainingGroup1Temp(indexes) = [];
        [litter, indexesTrainingGroup1] = intersect(positionsFinalTrainingGroup1,positionsFinalTrainingGroup1Temp);
        [litter, indexes] = setdiff(positionsFinalTrainingGroup2Temp,deepestFpositions2(:,i1,p));
        positionsFinalTrainingGroup2Temp(indexes) = [];
        [litter, indexesTrainingGroup2] = intersect(positionsFinalTrainingGroup2,positionsFinalTrainingGroup2Temp);
      else
        indexesTrainingGroup1 = 1:length(positionsFinalTrainingGroup1);
        indexesTrainingGroup2 = 1:length(positionsFinalTrainingGroup2);
      end
      %%% From functional to multivariate data
      startClock = clock;
      [groupMV1(indexesTrainingGroup1,i1,p) groupMV1test(:,i1,p) groupMV2(indexesTrainingGroup2,i1,p) groupMV2test(:,i1,p)] =...
          fromFDtoMV(numbersFDtoCLASSmethods(i1,1),infoFDtoMVmethod,...
                     allGroupFD1(positionsFinalTrainingGroup1Temp,:),newDataFD,...
                     allGroupFD2(positionsFinalTrainingGroup2Temp,:),newDataFD);
      timeFDtoMV = etime(clock,startClock);
      computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + timeFDtoMV;
      computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + timeFDtoMV;
      computationalTimesRunsMVmethod(b,i1) = computationalTimesRunsMVmethod(b,i1) + timeFDtoMV;
      startClock = clock;
      %%% The classifications (by using the multivariate data)  
      % Univariate methods DFMp: classification based on D^p(f)
      % Option 1: Minimizing the distance to the training centroids (in this case it's not necessary to look at the sign of the variable)
      % Option 2: Looking at the sign of the variable
      outclass = 1-([groupMV1test(:,i1,p); groupMV2test(:,i1,p)] < 0);
      equalDist = ([groupMV1test(:,i1,p); groupMV2test(:,i1,p)] == 0);
      if any(equalDist)  % The elements in the frontier are allocated at random
        disp(['For method ' namesFDtoCLASSmethods{i1} ', number of the ' int2str(sum(m1+m2)) ' testing elements such that x' int2str(p) ' = 0: ' int2str(sum(equalDist))])
        % The elements such that groupMVtest(:,i1,p) == 0 are modified so as to assign them at random...
        outclass(equalDist) = randi(2,sum(equalDist),1)-1;
      end
      pseudoMisclass1 = outclass(1:m1);
      pseudoMisclass2 = 1-outclass((m1+1):(m1+m2));
     if ~isequal(sum(pseudoMisclass1 ~= pseudoMisclass2),M)  % This must hold
       disp(['Different classification, when (p,i1)=(' num2str(p) ',' num2str(i1) ') for new data ' num2str(find(pseudoMisclass1 == pseudoMisclass2)) ])
     end
      classificationOutputs(((M*(b-1)+1):(M*b)),i1,p) = pseudoMisclass1;  % Labels: 0 for group 1 and 1 for group 2
      computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + etime(clock,startClock);
      computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + etime(clock,startClock);
      if isequal(b,1) && dataFigures && ~paramOptimization
        if ~isequal(paramValues(p),0)
          switch diffsMode
            case 'derivatives'
              groupFD1Temp = ((T-1)^paramValues(p))*diff(groupFD1,paramValues(p),2);
              groupFD1testTemp = ((T-1)^paramValues(p))*diff(groupFD1test,paramValues(p),2);
              groupFD2Temp = ((T-1)^paramValues(p))*diff(groupFD2,paramValues(p),2);
              groupFD2testTemp = ((T-1)^paramValues(p))*diff(groupFD2test,paramValues(p),2);
            case 'differentials'
              groupFD1Temp = diff(groupFD1,paramValues(p),2);
              groupFD1testTemp = diff(groupFD1test,paramValues(p),2);
              groupFD2Temp = diff(groupFD2,paramValues(p),2);
              groupFD2testTemp = diff(groupFD2test,paramValues(p),2);
          end
          figure  % Functional data
          xMin = 0;
          xMax = 1;
          % xRange = xMax - xMin;
          yMin = min(min([groupFD1Temp; groupFD1testTemp; groupFD2Temp; groupFD2testTemp]));
          yMax = max(max([groupFD1Temp; groupFD1testTemp; groupFD2Temp; groupFD2testTemp]));
          yRange = yMax - yMin;
          if isequal(yRange,0)
            yRange = 0.01;
          end
          subplot(2,1,1)
            refGroupFD1 = referenceElement(numbersFDtoCLASSmethods(i1,4),infoFDreference,groupFD1Temp);
            refGroupFD2 = referenceElement(numbersFDtoCLASSmethods(i1,4),infoFDreference,groupFD2Temp);
            plotFunctions(t(1:(end-paramValues(p))),groupFD1Temp(1:numberElementsToPlot,:),groupFD2Temp(1:numberElementsToPlot,:),...
              ['Some functions of the training and testing samples, method ' namesFDtoCLASSmethods{i1}],percentMargin,true,...
              refGroupFD1,refGroupFD2)
              xlabel('Training elements and the mean')
              ylabel(['Order = ' num2str(paramValues(p))])
              set(gca,'XLim',[xMin xMax])
              set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
          subplot(2,1,2)
            refGroupFD1 = referenceElement(numbersFDtoCLASSmethods(i1,4),infoFDreference,groupFD1testTemp);
            refGroupFD2 = referenceElement(numbersFDtoCLASSmethods(i1,4),infoFDreference,groupFD2testTemp);
            plotFunctions(t(1:(end-paramValues(p))),groupFD1testTemp(1:numberElementsToPlot,:),groupFD2testTemp(1:numberElementsToPlot,:),...
              ' ',percentMargin,true,refGroupFD1,refGroupFD2)
              xlabel('Testing elements and the mean')
              ylabel(['Order = ' num2str(paramValues(p))])
              set(gca,'XLim',[xMin xMax])
              set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        end  
        figure  % Discriminant variable
        xMin = 0;
        yMin = min(min([groupMV1(:,i1,p); groupMV1test(:,i1,p); groupMV2(:,i1,p); groupMV2test(:,i1,p)]));
        yMax = max(max([groupMV1(:,i1,p); groupMV1test(:,i1,p); groupMV2(:,i1,p); groupMV2test(:,i1,p)]));
        yRange = yMax - yMin;
        if isequal(yRange,0)
          yRange = 0.01;
        end
        subplot(3,1,1)
          plotDiscrimVariables(groupMV1(:,i1,p),groupMV1(:,i1,p),groupMV2(:,i1,p),groupMV2(:,i1,p),...
              ['Discriminant variable, method ' namesFDtoCLASSmethods{i1}],percentMargin)
          xlabel('Training elements')
          ylabel(['Order = ' num2str(paramValues(p))])
          xMax = max([n1 n2])+1;
          % xRange = xMax - xMin;
          set(gca,'XLim',[xMin xMax])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        subplot(3,1,2)
          plotDiscrimVariables(groupMV1test(:,i1,p),groupMV1test(:,i1,p),...
              groupMV2test(:,i1,p),groupMV2test(:,i1,p),' ',percentMargin)
          xlabel('Testing elements')
          ylabel(['Order = ' num2str(paramValues(p))])
          xMax = max([m1 m2])+1;
          % xRange = xMax - xMin;
          set(gca,'XLim',[xMin xMax])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        subplot(3,1,3)
          hist([groupMV1test(:,i1,p); groupMV2test(:,i1,p)])
          xlabel(['Histogram of the variable (testing elements), method ' namesFDtoCLASSmethods{i1}])
          ylabel(['Order = ' num2str(paramValues(p))])
      end
    end
  end
  clear groupFD1 groupFD1test groupFD2 groupFD2test...
      groupFD1Temp groupFD1testTemp groupFD2Temp groupFD2testTemp
  if paramOptimization
    % Registering the labels when the parameter value and the method are/were optimized
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      minimizingValuesClassificationOutputs(((M*(b-1)+1):(M*b)),i1) = classificationOutputs(((M*(b-1)+1):(M*b)),i1,minValue(i1));
    end
    minimizingValuesMethodsClassificationOutputs(((M*(b-1)+1):(M*b)),1) = classificationOutputs(((M*(b-1)+1):(M*b)),minMethod,minValue(minMethod));
  else
    % Applying the multivariate method to the vector of discriminant variables
    % To avoid a classification based on numerical errors.
    groupMV1 = groupMV1.*(abs(groupMV1) > eps^(1/3));
    groupMV2 = groupMV2.*(abs(groupMV2) > eps^(1/3));
    groupMV1test = groupMV1test.*(abs(groupMV1test) > eps^(1/3));
    groupMV2test = groupMV2test.*(abs(groupMV2test) > eps^(1/3));
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      startClock = clock;
      % Preprocessing of the discriminant variables
      Sw = ((n1-1)*cov(reshape(groupMV1(:,i1,:),n1,length(paramValues))) +...
          (n2-1)*cov(reshape(groupMV2(:,i1,:),n2,length(paramValues))))/(n1+n2-2);  % The pooled covariance matrix
      varsIncluded = find(diag(Sw) > variaCutoff);  % Checking the variances
      varsExcluded = find(diag(Sw) <= variaCutoff);
      if verbose
        disp(['    Information for method ' namesFDtoCLASSmethods{i1}])
        vars = [];
        for i2 = 1:length(varsExcluded)
          vars = [vars ' ' ['x' int2str(varsExcluded(i2))]];
        end
        if isempty(vars)
          disp('      Variables with null variability in the whole data set: None')        
        else
          disp(['      Variables with null variability in the whole data set: ' vars])
        end
        vars = [];
        for i2 = 1:length(paramValues)
          vars = [vars ' ' ['x' int2str(i2)]];
        end
        disp(['      Discriminant variables available: ' vars])
        vars = [];
        for i2 = 1:length(varsIncluded)
          vars = [vars ' ' ['x' int2str(varsIncluded(i2))]];
        end
        disp(['      Discriminant variables considered: ' vars])
      end
      % Considering the included variables
      Sw = Sw(varsIncluded,varsIncluded);
      stdvarsIncluded = std([reshape(groupMV1(:,i1,varsIncluded),n1,length(varsIncluded));...
          reshape(groupMV2(:,i1,varsIncluded),n2,length(varsIncluded))])';
      % Multivariate method
      switch numbersFDtoCLASSmethods(i1,2)  % Necessary information for each method
        case 1  % Method WI: classification based on Fisher's method
          infoMVtoCLASmethod = {varsIncluded, stdvarsIncluded};
        case 2  % Method WD: classification based on weighted distance
          infoMVtoCLASmethod = {varsIncluded, stdvarsIncluded, Sw};
        case 3  % Method: KNN (multivariate K-nearest neighbours)
          infoMVtoCLASmethod = {numberNeighbours};
      end
      [errorG1temp errorG2temp outputMVtoCLASmethod] =...
          fromMVtoCLASS(numbersFDtoCLASSmethods(i1,2), infoMVtoCLASmethod,...
            reshape(groupMV1(:,i1,:),n1,length(paramValues)), reshape(groupMV1test(:,i1,:),m1,length(paramValues)),...
            reshape(groupMV2(:,i1,:),n2,length(paramValues)), reshape(groupMV2test(:,i1,:),m2,length(paramValues)));
      switch numbersFDtoCLASSmethods(i1,2)  % Possible additional information of each method
        case 1  % Method WI: classification based on Fisher's method
          coefficientsMVmethod(b,i1,:) = outputMVtoCLASmethod{1};  % Standardized version, for the figures
          coefficients = outputMVtoCLASmethod{2};  % Nonstandardized version, for the projection
          % Projecting the testing samples
          groupY1test(:,i1) = reshape(groupMV1test(:,i1,:),m1,length(paramValues))*coefficients;
          groupY2test(:,i1) = reshape(groupMV2test(:,i1,:),m2,length(paramValues))*coefficients;
          % Projecting the centroids of the training samples
          groupY1c = mean(reshape(groupMV1(:,i1,:),n1,length(paramValues)),1)*coefficients;
          groupY2c = mean(reshape(groupMV2(:,i1,:),n2,length(paramValues)),1)*coefficients;
          pseudoMisclass1 = outputMVtoCLASmethod{3};
          pseudoMisclass2 = outputMVtoCLASmethod{4};
        case 2  % Method WD: classification based on weighted distance
          coefficientsMVmethod(b,i1,:) = outputMVtoCLASmethod{1};  % Standardized version, for the figures
          coefficients = outputMVtoCLASmethod{2};  % Nonstandardized version, for the projection
          % Projecting the testing samples
          groupY1test(:,i1) = reshape(groupMV1test(:,i1,:),m1,length(paramValues))*coefficients;
          groupY2test(:,i1) = reshape(groupMV2test(:,i1,:),m2,length(paramValues))*coefficients;
          % Projecting the centroids of the training samples
          groupY1c = mean(reshape(groupMV1(:,i1,:),n1,length(paramValues)),1)*coefficients;
          groupY2c = mean(reshape(groupMV2(:,i1,:),n2,length(paramValues)),1)*coefficients;
          pseudoMisclass1 = outputMVtoCLASmethod{3};
          pseudoMisclass2 = outputMVtoCLASmethod{4};
        case 3  % Method KNN: k-nearest neighbours
          pseudoMisclass1 = outputMVtoCLASmethod{1};
          pseudoMisclass2 = outputMVtoCLASmethod{2};
      end
      if ~isequal(sum(pseudoMisclass1 ~= pseudoMisclass2),M)  % This must hold
        disp(['Different MV classification, when i1=' num2str(i1) ', for new data ' num2str(find(pseudoMisclass1 == pseudoMisclass2)) ])
      end
      computationalTimesRunsMVmethod(b,i1) = computationalTimesRunsMVmethod(b,i1) + etime(clock,startClock);
    end
  end
  if isequal(b,1) && dataFigures  % Figures: results of this run
    if ~paramOptimization 
      for i1 = 1:size(numbersFDtoCLASSmethods,1)
        switch length(paramValues)
          case 2
            figure
            plotDiscrimVariables(reshape(groupMV1(:,i1,1:2),n1,2),reshape(groupMV1test(:,i1,1:2),m1,2),...
                reshape(groupMV2(:,i1,1:2),n2,2),reshape(groupMV2test(:,i1,1:2),m2,2),...
                ['Two-dimensional representation, testing samples and centroids of the samples (testing in thinner line), method ' namesFDtoCLASSmethods{i1}],...
                percentMargin,coefficientsMVmethod(b,i1,:))
          case 3
            figure
            plotDiscrimVariables(reshape(groupMV1(:,i1,1:3),n1,3),reshape(groupMV1test(:,i1,1:3),m1,3),...
                reshape(groupMV2(:,i1,1:3),n2,3),reshape(groupMV2test(:,i1,1:3),m2,3),...
                ['Three-dimensional representation, testing samples and centroids of the samples (testing in thinner line), method ' namesFDtoCLASSmethods{i1}],...
                percentMargin,coefficientsMVmethod(b,i1,:))
        end
        % A matrix of scatter plots
        figure
        gplotmatrix([reshape(groupMV1test(:,i1,:),m1,length(paramValues)); reshape(groupMV2test(:,i1,:),m2,length(paramValues))],...
            [],[zeros(m1,1); ones(m2,1)],'br','os',3)
        title(['Array of scatter plots for the discriminant variables of the testing samples, method ' namesFDtoCLASSmethods{i1}]);
        if any([1 2] == numbersFDtoCLASSmethods(i1,2))
          plotDiscrimFunctions(groupY1test(:,i1),groupY2test(:,i1),groupY1c,groupY2c,...
              ['Discriminant function and centroids proyection (testing in dotted), method ' namesFDtoCLASSmethods{i1}],percentMargin)
        end
      end
    end
    % Results
    plotResults([],[],[],[],[],[],[],[],[],...
        paramValues, paramValuesMinpower(1,:,:), methodsMinpower(1,:), percentMargin, namesFDtoCLASSmethods,...
        computationalTimesRuns(1,:,:), computationalTimesRunsWholeParamValues(1,:), paramOptimization, methodSelection,...
        coefficientsMVmethod(1,:,:), computationalTimesRunsMVmethod(1,:,:), true);
  end
  if verbose
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      if all(methsXparams(i1,:)==0)  % Not to consider this method
        continue
      end
      disp(['    Information for method ' namesFDtoCLASSmethods{i1}])
      for p = 1:length(paramValues)
        if isequal(methsXparams(i1,p),0)  % Not to consider this value for this method
          continue
        end
        disp(['      Classification labels for the new elements, for p = ' num2str(p) ',   ' num2str(1+classificationOutputs(((M*(b-1)+1):(M*b)),i1,p)')])
      end
      if paramOptimization
        disp(['      Classification labels for the new elements, when the parameter value is optimized: ' '   ' num2str(1+minimizingValuesClassificationOutputs(((M*(b-1)+1):(M*b)),i1)')])
      end
      disp(['      Time spent for each value of the parameter values: ' num2str(computationalTimesRuns(b,i1,:))])
      disp(['      Joint time spent for all values of the parameter values: ' num2str(computationalTimesRunsWholeParamValues(b,i1))])
      if ~paramOptimization
        disp('      Multivariate method')
        disp(['        Classification labels for the new elements: ' '   ' num2str(1+classificationOutputs(((M*(b-1)+1):(M*b)),i1,length(paramValues)+1)')])
        disp(['        Coefficients of the parameter values: ' num2str(coefficientsMVmethod(b,i1,methsXparams(i1,:)))])
        disp(['        Time spent with the vector of variables: ' num2str(computationalTimesRunsMVmethod(b,i1))])
      end
    end
    if paramOptimization
      disp(['      Classification labels for the new elements, when both the parameter value and the method are optimized: ' '   ' num2str(1+minimizingValuesMethodsClassificationOutputs(((M*(b-1)+1):(M*b)),1)')])
    end
  end
end

if verbose
  disp(' ')
  disp('... ending the iterations') 
  disp(' ')
end

%%%%%%%%%%%%%%%%  Ending Main Computes  %%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%  Beginning Results Processing  %%%%%%%%%%%%%%%%%%%%

if verbose
  if methodSelection
    disp('CLASSIFICATION LABELS (with both parameter optimization and method selection)')
      disp('-> Labels for the new data: ')
      for b = 1:B
        disp(['  Run b = ' num2str(b) ':  '  num2str(1+minimizingValuesMethodsClassificationOutputs(((M*(b-1)+1):(M*b)),1)')])
      end
      if B>1  % Majority vote information
        matrixOutput = zeros(B,M);
        for b = 1:B
          matrixOutput(b,:) = 1+minimizingValuesMethodsClassificationOutputs(((M*(b-1)+1):(M*b)),1)';
        end
        vectorOutput = zeros(1,M);
        for m = 1:M
          if sum(matrixOutput(:,m)==1)>sum(matrixOutput(:,m)==2)
            vectorOutput(1,m) = 1;
          elseif sum(matrixOutput(:,m)==1)<sum(matrixOutput(:,m)==2)
            vectorOutput(1,m) = 2;
          end
        end
        disp(['  Majority vote (0 for ties), ' num2str(vectorOutput)])
      end
  end
  for i1 = 1:size(numbersFDtoCLASSmethods,1)
    disp(['METHOD: ' namesFDtoCLASSmethods{i1}])
    if paramOptimization    % Errors have been calculated for possibly different values of the parameter
      if ~methodSelection
        disp('  CLASSIFICATION LABELS (with parameter optimization)')
          disp('  -> Labels for the new data: ')
          for b = 1:B
            disp(['    Run b = ' num2str(b) ':  '  num2str(1+minimizingValuesClassificationOutputs(((M*(b-1)+1):(M*b)),i1)')])
          end
          if B>1  % Majority vote information
            matrixOutput = zeros(B,M);
            for b = 1:B
              matrixOutput(b,:) = 1+minimizingValuesClassificationOutputs(((M*(b-1)+1):(M*b)),i1)';
            end
            vectorOutput = zeros(1,M);
            for m = 1:M
              if sum(matrixOutput(:,m)==1)>sum(matrixOutput(:,m)==2)
                vectorOutput(1,m) = 1;
              elseif sum(matrixOutput(:,m)==1)<sum(matrixOutput(:,m)==2)
                vectorOutput(1,m) = 2;
              end
            end
            disp(['    Majority vote (0 for ties), ' num2str(vectorOutput)])
          end
      end
      disp('  MINIMIZING-POWER OF THE PARAMETER VALUES')
        disp('  -> Mean: ')
        disp(['     '  num2str(mean(paramValuesMinpower(:,i1,:),1))])
    else  % Errors for all values of the parameter have been calculated
      disp('  CLASSIFICATION LABELS (without parameter optimization nor method selection)')
        disp('  -> Labels for the new data: ')
        for p = 1:length(paramValues)
          disp(['    Parameter value i-th = ' num2str(paramValues(p))])
          for b = 1:B
            if isequal(methsXparams(i1,p),0)  % Not to consider this value for this method
              continue
            end
            disp(['      Run b = ' num2str(b) ':  '  num2str(1+classificationOutputs(((M*(b-1)+1):(M*b)),i1,p)')])
          end
          if B>1  % Majority vote information
            matrixOutput = zeros(B,M);
            for b = 1:B
              matrixOutput(b,:) = 1+classificationOutputs(((M*(b-1)+1):(M*b)),i1,p)';
            end
            vectorOutput = zeros(1,M);
            for m = 1:M
              if sum(matrixOutput(:,m)==1)>sum(matrixOutput(:,m)==2)
                vectorOutput(1,m) = 1;
              elseif sum(matrixOutput(:,m)==1)<sum(matrixOutput(:,m)==2)
                vectorOutput(1,m) = 2;
              end
            end
            disp(['      Majority vote (0 for ties), ' num2str(vectorOutput)])
          end
        end
    end
    disp('  COMPUTATIONAL TIMES')
      disp('  -> Mean time spent for each value of the parameter values: ')
      disp(['     '  num2str(mean(computationalTimesRuns(:,i1,:),1))])
      disp('  -> Mean joint time spent for all values of the parameter values: ')
      disp(['     '  num2str(mean(computationalTimesRunsWholeParamValues(:,i1),1))])
    if ~paramOptimization
      disp('  MULTIVARIATE METHOD')
        disp('  -> Mean coefficients of the parameter values: ')
        disp(['     '  num2str(mean(coefficientsMVmethod(:,i1,:),1))])
        disp('  -> Mean joint time spent for all values of the parameter values: ')
        disp(['     '  num2str(mean(computationalTimesRunsMVmethod(:,i1),1))])
    end
  end
  if methodSelection
    disp('MINIMIZING-POWER OF THE METHODS')
      disp('-> Mean: ')
      disp(['   '  num2str(mean(methodsMinpower,1))])
  end
end

% Boxplot with the results: error rates and the coefficients
if ~(dataFigures && isequal(B,1)) && resultsFigures
  if verbose
    disp(' ')
    disp('Result figures...')
    disp(' ')
  end
  plotResults([],[],[],[],[],[],[],[],[],...
      paramValues, paramValuesMinpower, methodsMinpower, percentMargin, namesFDtoCLASSmethods,...
    computationalTimesRuns, computationalTimesRunsWholeParamValues, paramOptimization, methodSelection,...
    coefficientsMVmethod, computationalTimesRunsMVmethod, true);
end

if saveResults
  save('resultsApplication','paramValues','methsXparams','paramValuesMinpower','methodsMinpower',...
      'percentMargin','namesFDtoCLASSmethods','computationalTimesRuns','computationalTimesRunsWholeParamValues',...
      'paramOptimization','methodSelection','coefficientsMVmethod','computationalTimesRunsMVmethod')
  if verbose
    disp('The results have been saved')
    disp(' ')
  end
end
% load('resultsApplication')

%%%%%%%%%%%%%%%%  Ending Results Processing  %%%%%%%%%%%%%%%%%%%%


if verbose
  disp(' ')
  disp('STOP    (And remember: `Everything and nothing is possimpible´, Barney Stinson) ')
  disp(' ')
  disp(' ')
end

clear errorRatesGroup2param Bparam errorRatesParam groupFD1 groupFD1test N1 groupFD2 numberElementsToPlot...
    N1param groupFD2test numberNeighbours N2 groupMV1 N2param groupMV1test...
    outputFDtoCLASmethod groupMV2 p groupMV2test paramDataReuse groupY1test groupY2test i1 indexes...
    infoFDreference infoFDtoCLASmethod infoTStoFDmethod positionsFinalTrainingGroup1 leaveOneOutParam...
    positionsFinalTrainingGroup1Temp litter positionsFinalTrainingGroup2 m1 positionsFinalTrainingGroup2Temp...
    m1param positionsTestingGroup1param m2 positionsTestingGroup1paramTemp m2param positionsTestingGroup2param...
    meanErrorRatesParam positionsTestingGroup2paramTemp b positionsTrainingGroup1 b2...
    positionsTrainingGroup1param methodMinError positionsTrainingGroup1paramTemp minMethod...
    positionsTrainingGroup2 positionsTrainingGroup2param minValue positionsTrainingGroup2paramTemp...
    posts pseudoMisclass1 modelsJointly pseudoMisclass2 counter1 n1 counter2 n1Temp n2 startClock...
    depthBasedRobustness n2Temp timeTStoFD dyadicSplits trustworthyRunsParam errorG1temp variaCutoff...
    errorG2temp errorRatesGroup1param verbose dataFigures resultsFigures
