
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  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(m1*0.05);
  modelsJointly = false;
  percentMargin = 5;
resultsFigures = true;
saveResults = false;

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


clear all


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

%%% REAL DATA

allGroupFD1 = ;
allGroupFD2 = ;

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


% %%% An example with pseudo-real data...
%   numberExercise = 2;
%     paramVector1 = [1 1];  paramVector2 = [ 1 1];
%   T = 224;
%   N1= 48; N2 = 51;
%   allGroupFD1 = dataExercisesFD(numberExercise,1,N1,T,paramVector1);
%   allGroupFD2 = dataExercisesFD(numberExercise,2,N2,T,paramVector2);
%   clear numberExercise T t N1 N2 paramVector1 paramVector2

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



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

% %%% For the leave-one-out cross-validation
% m1 = 1;
% m2 = 1;
% leaveOneOut = true;
%%% For the leave-m-out cross-validation
m1 = ; % It must be lower than N1
m2 = ; % It must be lower than N2
leaveOneOut = false;

B = ;  % Number of runs

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
  
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;
discretenessMargin = 0.25;

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... %%%%
N1 = size(allGroupFD1,1);
N2 = size(allGroupFD2,1);
if leaveOneOut
  B = N1+N2;
  disp(' ')
  disp('The value of B has been changed, since leave-one-out cross-validation is being run!')
  disp(' ')
end
if isequal(m1,1) && isequal(m2,1)
  if isequal(leaveOneOut,false)
    disp(' ')
    disp('Testing sample sizes are 1 and the elements will be chosen at random, not exhaustively!')
    disp(' ')
  end
end
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
errorRatesGroup1 = zeros(B,size(numbersFDtoCLASSmethods,1),length(paramValues)+1);
errorRatesGroup2 = zeros(B,size(numbersFDtoCLASSmethods,1),length(paramValues)+1);
minimizingValuesErrorRates = zeros(B,size(numbersFDtoCLASSmethods,1));
minimizingValuesErrorRatesGroup1 = zeros(B,size(numbersFDtoCLASSmethods,1));
minimizingValuesErrorRatesGroup2 = zeros(B,size(numbersFDtoCLASSmethods,1));
paramValuesMinpower = zeros(B,size(numbersFDtoCLASSmethods,1),length(paramValues));
minimizingValuesMethodsErrorRates = zeros(B,1);
minimizingValuesMethodsErrorRatesGroup1 = zeros(B,1);
minimizingValuesMethodsErrorRatesGroup2 = zeros(B,1);
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
  % Selection of the final testing groups
  if leaveOneOut
    if b <= N1
      positionsFinalTestingGroup1 = b;
      m1 = 1;
      positionsFinalTestingGroup2 = [];
      m2 = 0;
    else
      positionsFinalTestingGroup1 = [];
      m1 = 0;
      positionsFinalTestingGroup2 = b-N1;
      m2 = 1;
    end
  else
    positionsFinalTestingGroup1 = randsample(N1,m1);
    positionsFinalTestingGroup2 = randsample(N2,m2);
  end
  positionsTrainingGroup1 = (1:N1)';
  positionsTrainingGroup2 = (1:N2)';
  % Rest of the data
  [litter, indexes] = setdiff(positionsTrainingGroup1,positionsFinalTestingGroup1);
  positionsTrainingGroup1 = positionsTrainingGroup1(indexes);
  [litter, indexes] = setdiff(positionsTrainingGroup2,positionsFinalTestingGroup2);
  positionsTrainingGroup2 = positionsTrainingGroup2(indexes);
  % 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((N1-m1)/2)  and  n2=ceil((N2-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(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
    for i1 = 1:size(numbersFDtoCLASSmethods,1)
      methodMinError(i1) = min(meanErrorRatesParam(i1,:));  % Minimum error rate
      methsXparams(i1,meanErrorRatesParam(i1,:)<=(methodMinError(i1)*(1+discretenessMargin))) = 1;  % Values with error rate in the interval
    end
    if methodSelection
      methsXparams(methodMinError>(min(methodMinError)*(1+discretenessMargin)),:) =...
          zeros(sum(methodMinError>(min(methodMinError)*(1+discretenessMargin))),length(paramValues));  % Disabling the other methods (for the results)
    end
    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 = allGroupFD1(positionsFinalTestingGroup1,:);
    groupFD2 = allGroupFD2(positionsFinalTrainingGroup2,:);
    groupFD2test = allGroupFD2(positionsFinalTestingGroup2,:);
    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
      startClock = clock;
      %%% From functional to multivariate data
      [groupMV1(indexesTrainingGroup1,i1,p) groupMV1test(:,i1,p) groupMV2(indexesTrainingGroup2,i1,p) groupMV2test(:,i1,p)] =...
          fromFDtoMV(numbersFDtoCLASSmethods(i1,1),infoFDtoMVmethod,...
                     allGroupFD1(positionsFinalTrainingGroup1Temp,:),allGroupFD1(positionsFinalTestingGroup1,:),...
                     allGroupFD2(positionsFinalTrainingGroup2Temp,:),allGroupFD2(positionsFinalTestingGroup2,:));
      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
      errorRatesGroup1(b,i1,p) = sum(outclass(1:m1));
      errorRatesGroup2(b,i1,p) = sum(1-outclass((m1+1):(m1+m2)));
      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
  % Calculations to select the parameter value(s) or the method(s)
  for i1 = 1:size(numbersFDtoCLASSmethods,1)
    posts = find(methsXparams(i1,:));  % Parameter values that have been considered
    if isempty(posts)  % This may happen only when 'methodSelection=true'
      continue
    end
    [minimizingValuesErrorRates(b,i1) minG] = min((errorRatesGroup1(b,i1,methsXparams(i1,:))+errorRatesGroup2(b,i1,methsXparams(i1,:)))/(m1+m2));
    posts = posts(minG);  % Minimizing parameter values
    minimizingValuesErrorRatesGroup1(b,i1) = errorRatesGroup1(b,i1,posts(1))/m1;
    minimizingValuesErrorRatesGroup2(b,i1) = errorRatesGroup2(b,i1,posts(1))/m2;
    if paramOptimization
      [litter,indexes] = setdiff(1:length(paramValues),posts);
      methsXparams(i1,indexes) = zeros(1,length(indexes));  % Disabling the other parameter values (for the results)
    end
    % paramValuesMinpower(b,i1,posts(1)) = 1;  % For the minimizing-power of the parameter values: only the smallest
    paramValuesMinpower(b,i1,posts) = 1/length(posts);  % For the minimizing-power of the parameter values: all the minimizing
  end
  posts = find(sum(methsXparams,2));  % Methods that have been considered
  [minimizingValuesMethodsErrorRates(b,1) minG] = min(minimizingValuesErrorRates(b,posts));  % Identifying the minimizing-in-mean method
  posts = posts(minG);  % Minimizing-in-mean methods
  minimizingValuesMethodsErrorRatesGroup1(b,1) = minimizingValuesErrorRatesGroup1(b,posts(1));
  minimizingValuesMethodsErrorRatesGroup2(b,1) = minimizingValuesErrorRatesGroup2(b,posts(1));
  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
  % Applying the multivariate method to the vector of discriminant variables
  if ~paramOptimization
    % 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
          infoMVtoCLASSmethod = {varsIncluded, stdvarsIncluded};
        case 2  % Method WD: classification based on weighted distance
          infoMVtoCLASSmethod = {varsIncluded, stdvarsIncluded, Sw};
        case 3  % Method: KNN (multivariate K-nearest neighbours)
          infoMVtoCLASSmethod = {numberNeighbours};
      end      
      [errorG1temp errorG2temp outputMVtoCLASSmethod] =...
          fromMVtoCLASS(numbersFDtoCLASSmethods(i1,2), infoMVtoCLASSmethod,...
            reshape(groupMV1(:,i1,:),n1,length(paramValues)), reshape(groupMV1test(:,i1,:),m1,length(paramValues)),...
            reshape(groupMV2(:,i1,:),n2,length(paramValues)), reshape(groupMV2test(:,i1,:),m2,length(paramValues)));
      errorRatesGroup1(b,i1,length(paramValues)+1) = errorG1temp;
      errorRatesGroup2(b,i1,length(paramValues)+1) = errorG2temp;
      switch numbersFDtoCLASSmethods(i1,2)  % Possible additional information of each method
        case 1  % Method WI: classification based on Fisher's method
          coefficientsMVmethod(b,i1,:) = outputMVtoCLASSmethod{1};  % Standardized version, for the figures
          coefficients = outputMVtoCLASSmethod{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;
        case 2  % Method WD: classification based on weighted distance
          coefficientsMVmethod(b,i1,:) = outputMVtoCLASSmethod{1};  % Standardized version, for the figures
          coefficients = outputMVtoCLASSmethod{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;
      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(errorRatesGroup1(1,:,:)/m1,errorRatesGroup2(1,:,:)/m2,...
        (errorRatesGroup1(1,:,:) + errorRatesGroup2(1,:,:))/(m1+m2),...
        minimizingValuesErrorRatesGroup1(1,:),minimizingValuesErrorRatesGroup2(1,:),...
        minimizingValuesErrorRates(1,:),minimizingValuesMethodsErrorRatesGroup1(1),...
        minimizingValuesMethodsErrorRatesGroup2(1),minimizingValuesMethodsErrorRates(1),...
        paramValues,paramValuesMinpower(1,:,:),methodsMinpower(1,:),...
        percentMargin,namesFDtoCLASSmethods,computationalTimesRuns(1,:,:),...
        computationalTimesRunsWholeParamValues(1,:),paramOptimization,methodSelection,...
        coefficientsMVmethod(1,:,:),computationalTimesRunsMVmethod(1,:));
  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}])
      disp(['      Error rates for group 1: ' '   ' num2str(errorRatesGroup1(b,i1,methsXparams(i1,:))/m1)])
      disp(['      Error rates for group 2: ' '   ' num2str(errorRatesGroup2(b,i1,methsXparams(i1,:))/m2)])
      disp(['      Overall error rates: ' '   ' num2str((errorRatesGroup1(b,i1,methsXparams(i1,:))+errorRatesGroup2(b,i1,methsXparams(i1,:)))/(m1+m2))])
      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(['        Error rates for group 1: ' '   ' num2str(errorRatesGroup1(b,i1,length(paramValues)+1)/m1)])
        disp(['        Error rates for group 2: ' '   ' num2str(errorRatesGroup2(b,i1,length(paramValues)+1)/m2)])
        disp(['        Overall error rates: ' '   ' num2str((errorRatesGroup1(b,i1,length(paramValues)+1)+errorRatesGroup2(b,i1,length(paramValues)+1))/(m1+m2))])
        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
  end
end

% Note that m1 = 0 if leaveOneOut = true
errorRates = (errorRatesGroup1 + errorRatesGroup2)/(m1+m2);
if leaveOneOut
  % Its original value is assigned to m1
  m1 = 1;
  % Taking the nonnull submatrices
  errorRatesGroup1 = errorRatesGroup1(1:N1,:);
  errorRatesGroup2 = errorRatesGroup2((N1+1):(N1+N2),:);
else
  errorRatesGroup1 = errorRatesGroup1/m1;
  errorRatesGroup2 = errorRatesGroup2/m2;
end

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

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




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

if verbose
  if methodSelection
    disp('MISCLASSIFICATION ERROR RATES (with both parameter optimization and method selection)')
      disp('-> Mean of the estimated error rates for group 1: ')
      disp(['   '  num2str(mean(minimizingValuesMethodsErrorRatesGroup1,1))])
      disp('  -> Standard deviation: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup1,1))])
      disp('  -> Standard error: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup1,1)/sqrt(B))])
      disp('-> Mean of the estimated error rates for group 2: ')
      disp(['   '  num2str(mean(minimizingValuesMethodsErrorRatesGroup2,1))])
      disp('  -> Standard deviation: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup2,1))])
      disp('  -> Standard error: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup2,1)/sqrt(B))])
      disp('-> Mean of the estimated overall error rates: ')
      disp(['   '  num2str(mean(minimizingValuesMethodsErrorRates,1))])
      disp('  -> Standard deviation: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRates,1))])
      disp('  -> Standard error: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRates,1)/sqrt(B))])
  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('  MISCLASSIFICATION ERROR RATES (with parameter optimization)')
          disp('  -> Mean of the estimated error rates for group 1: ')
          disp(['     '  num2str(mean(minimizingValuesErrorRatesGroup1(:,i1),1))])
          disp('    -> Standard deviation: ')
          disp(['       '  num2str(std(minimizingValuesErrorRatesGroup1(:,i1),1))])
          disp('    -> Standard error: ')
          disp(['       '  num2str(std(minimizingValuesErrorRatesGroup1(:,i1),1)/sqrt(B))])
          disp('  -> Mean of the estimated error rates for group 2: ')
          disp(['     '  num2str(mean(minimizingValuesErrorRatesGroup2(:,i1),1))])
          disp('    -> Standard deviation: ')
          disp(['       '  num2str(std(minimizingValuesErrorRatesGroup2(:,i1),1))])
          disp('    -> Standard error: ')
          disp(['       '  num2str(std(minimizingValuesErrorRatesGroup2(:,i1),1)/sqrt(B))])
          disp('  -> Mean of the estimated overall error rates: ')
          disp(['     '  num2str(mean(minimizingValuesErrorRates(:,i1),1))])
          disp('    -> Standard deviation: ')
          disp(['       '  num2str(std(minimizingValuesErrorRates(:,i1),1))])
          disp('    -> Standard error: ')
          disp(['       '  num2str(std(minimizingValuesErrorRates(:,i1),1)/sqrt(B))])
      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('  MISCLASSIFICATION ERROR RATES (without parameter optimization nor method selection)')
        disp('  -> Mean of the estimated error rates for group 1: ')
        disp(['     '  num2str(mean(errorRatesGroup1(:,i1,1:(end-1)),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(errorRatesGroup1(:,i1,1:(end-1)),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(errorRatesGroup1(:,i1,1:(end-1)),1)/sqrt(B))])
        disp('  -> Mean of the estimated error rates for group 2: ')
        disp(['     '  num2str(mean(errorRatesGroup2(:,i1,1:(end-1)),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(errorRatesGroup2(:,i1,1:(end-1)),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(errorRatesGroup2(:,i1,1:(end-1)),1)/sqrt(B))])
        disp('  -> Mean of the estimated overall error rates: ')
        disp(['     '  num2str(mean(errorRates(:,i1,1:(end-1)),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(errorRates(:,i1,1:(end-1)),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(errorRates(:,i1,1:(end-1)),1)/sqrt(B))])
      disp('  MISCLASSIFICATION ERROR RATES (as if parameter optimization had been applied)')
        disp('  -> Mean of the estimated error rates for group 1: ')
        disp(['     '  num2str(mean(minimizingValuesErrorRatesGroup1(:,i1),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(minimizingValuesErrorRatesGroup1(:,i1),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(minimizingValuesErrorRatesGroup1(:,i1),1)/sqrt(B))])
        disp('  -> Mean of the estimated error rates for group 2: ')
        disp(['     '  num2str(mean(minimizingValuesErrorRatesGroup2(:,i1),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(minimizingValuesErrorRatesGroup2(:,i1),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(minimizingValuesErrorRatesGroup2(:,i1),1)/sqrt(B))])
        disp('  -> Mean of the estimated overall error rates: ')
        disp(['     '  num2str(mean(minimizingValuesErrorRates(:,i1),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(minimizingValuesErrorRates(:,i1),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(minimizingValuesErrorRates(:,i1),1)/sqrt(B))])
    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 of the estimated error rate for group 1: ')
        disp(['     '  num2str(mean(errorRatesGroup1(:,i1,end),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(errorRatesGroup1(:,i1,end),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(errorRatesGroup1(:,i1,end),1)/sqrt(B))])
        disp('  -> Mean of the estimated error rate for group 2: ')
        disp(['     '  num2str(mean(errorRatesGroup2(:,i1,end),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(errorRatesGroup2(:,i1,end),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(errorRatesGroup2(:,i1,end),1)/sqrt(B))])
        disp('  -> Mean of the estimated overall error rates: ')
        disp(['     '  num2str(mean(errorRates(:,i1,end),1))])
        disp('    -> Standard deviation: ')
        disp(['       '  num2str(std(errorRates(:,i1,end),1))])
        disp('    -> Standard error: ')
        disp(['       '  num2str(std(errorRates(:,i1,end),1)/sqrt(B))])
        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))])
  else
    disp('MISCLASSIFICATION ERROR RATES (as if both parameter optimization and method selection had been applied)')
      disp('-> Estimation of the mean of the estimated error rates for group 1: ')
      disp(['   '  num2str(mean(minimizingValuesMethodsErrorRatesGroup1,1))])
      disp('  -> Standard deviation: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup1,1))])
      disp('  -> Standard error: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup1,1)/sqrt(B))])
      disp('-> Estimation of the mean of the estimated error rates for group 2: ')
      disp(['   '  num2str(mean(minimizingValuesMethodsErrorRatesGroup2,1))])
      disp('  -> Standard deviation: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup2,1))])
      disp('  -> Standard error: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRatesGroup2,1)/sqrt(B))])
      disp('-> Estimation of the mean of the estimated global error rates: ')
      disp(['   '  num2str(mean(minimizingValuesMethodsErrorRates,1))])
      disp('  -> Standard deviation: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRates,1))])
      disp('  -> Standard error: ')
      disp(['     '  num2str(std(minimizingValuesMethodsErrorRates,1)/sqrt(B))])
  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(errorRatesGroup1,errorRatesGroup2,errorRates,...
      minimizingValuesErrorRatesGroup1,minimizingValuesErrorRatesGroup2,minimizingValuesErrorRates,...
      minimizingValuesMethodsErrorRatesGroup1,minimizingValuesMethodsErrorRatesGroup2,minimizingValuesMethodsErrorRates,...
      paramValues,paramValuesMinpower,methodsMinpower,percentMargin,namesFDtoCLASSmethods,...
      computationalTimesRuns,computationalTimesRunsWholeParamValues,...
      paramOptimization,methodSelection,coefficientsMVmethod,computationalTimesRunsMVmethod);
end

if saveResults
  save('resultsRealData','errorRatesGroup1','errorRatesGroup2','errorRates',...
      'minimizingValuesErrorRatesGroup1','minimizingValuesErrorRatesGroup2','minimizingValuesErrorRates',...
      'minimizingValuesMethodsErrorRatesGroup1','minimizingValuesMethodsErrorRatesGroup2','minimizingValuesMethodsErrorRates',...
      'paramValues','methsXparams','paramValuesMinpower','methodsMinpower','percentMargin','namesFDtoCLASSmethods',...
      'computationalTimesRuns','computationalTimesRunsWholeParamValues',...
      'paramOptimization','methodSelection','coefficientsMVmethod','computationalTimesRunsMVmethod')
  if verbose
    disp('The results have been saved')
    disp(' ')
  end
end
% load('resultsRealData')

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


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

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