
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  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  %%%%%%%%%%%%%%%%%%%%

numberExercise = ;  % Number of the simulation exercise
  paramVector1 = [];  paramVector2 = [];  % Parameters of the two models of the exercise

n1 = ;  % Size of the training sample of group 1
n2 = n1;  % Size of the training sample of group 2
m1 = ;  % Size of the testing sample of group 1
m2 = m1;  % Size of the testing sample of group 2

T = ;  % Length of the time series

contamination = false;  % Variable to indicate whether or not to introduce contaminating series
  contamType = 'A';  % Type of contamination
  n1c = 2; %ceil(n1*0.07);  % Number of contaminating series (in the training sample of group 1)
  alpha = 0.2;  % Proportion for the trim

B = ; % Number of runs

numbersTStoCLASSmethods = [1 1 1 1; 1 2 1 1];  % Numbers of the methods determined by each row
namesTStoCLASSmethods = {'DbC', 'DbC-alpha'};  % Names of the global methods

dyadicSplits = 0;  % Number of dyadic splits of the time series
numberBlocks = 2.^(0:dyadicSplits);  % Vector with the number of blocks

paramOptimization = false;  % To select the parameter(s) with a cross-validation loop. It is disabled if 'length(paramValues)=1'
  n1param = ;  % Training sample size of group 1
  n2param = n1param;  % Training sample size of group 2
  m1param = ;  % Testing sample size of group 1
  m2param = m1param;  % Testing sample size of group 2
  Bparam = ;  % Number of runs
  n1cParam = 1;  % Number of contaminating series in the training sample of group 1
methodSelection = false;
discretenessMargin = 0.25;


%%%% Setting some values... %%%%
if ~isequal(size(numbersTStoCLASSmethods,1),length(namesTStoCLASSmethods))
  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(numbersTStoCLASSmethods,1)
    namesTStoCLASSmethods{i1} = ['Method ' int2str(i1)];
  end
end
paramValues = numberBlocks;
if contamination && ~any(numbersTStoCLASSmethods(:,2) == 2)
  numbersTStoCLASSmethods = [numbersTStoCLASSmethods; 1 2 1 1 1 1];  % Method 2 is included
    namesTStoCLASSmethods = [namesTStoCLASSmethods, {'DbC-alpha'}];
  disp(' ')
  disp('Since there is contamination, method 2 (DbC-alpha) has been activated!')
  disp(' ')
end
if isequal(size(numbersTStoCLASSmethods,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  %%%%%%%%%%%%%%%%%%%%

if verbose
  disp(' ')
  disp('START')
end
if verbose
  disp(' ')
  disp('Generating the data...')
  disp(' ')
end
try
  % Generating all the time series
  allGroupTS1 = dataExercisesTS(numberExercise,1,n1*B,T,paramVector1);
  allGroupTS1test = dataExercisesTS(numberExercise,1,m1*B,T,paramVector1);
  allGroupTS2 = dataExercisesTS(numberExercise,2,n2*B,T,paramVector2);
  allGroupTS2test = dataExercisesTS(numberExercise,2,m2*B,T,paramVector2);
  if contamination
    % Generating all the contaminants
    allContamGroupTS1 = dataExercisesTS(numberExercise,1,n1c*B,T,paramVector1,contamType);
  end
  memoryProblem = false;
catch
  % If there is any problem with the memory of the computer, the functions are generated later
  memoryProblem = true;
end

if contamination
  n1c2 = ceil(n1c*0.1);
end

% Hypermatrixes to store the misclassification error rates, the minimizing-power measures and the computational times
errorRatesGroup1 = zeros(B,size(numbersTStoCLASSmethods,1),length(paramValues));
errorRatesGroup2 = zeros(B,size(numbersTStoCLASSmethods,1),length(paramValues));
minimizingValuesErrorRates = zeros(B,size(numbersTStoCLASSmethods,1));
minimizingValuesErrorRatesGroup1 = zeros(B,size(numbersTStoCLASSmethods,1));
minimizingValuesErrorRatesGroup2 = zeros(B,size(numbersTStoCLASSmethods,1));
paramValuesMinpower = zeros(B,size(numbersTStoCLASSmethods,1),length(paramValues));
minimizingValuesMethodsErrorRates = zeros(B,1);
minimizingValuesMethodsErrorRatesGroup1 = zeros(B,1);
minimizingValuesMethodsErrorRatesGroup2 = zeros(B,1);
methodsMinpower = zeros(B,size(numbersTStoCLASSmethods,1));
computationalTimesRuns = zeros(B,size(numbersTStoCLASSmethods,1),length(paramValues));
computationalTimesRunsWholeParamValues = zeros(B,size(numbersTStoCLASSmethods,1));

% Runs
if verbose
  disp('Starting the iterations...')
  disp(' ')
end
for b = 1:B
  if verbose
    disp(['  ITERATION: ' int2str(b)]) 
  end
  if paramOptimization
    if verbose
      disp('    Optimization loop...')
    end
    % Hypermatrixes to store the misclassification error rates
    errorRatesGroup1param = zeros(Bparam,size(numbersTStoCLASSmethods,1),length(paramValues));
    errorRatesGroup2param = zeros(Bparam,size(numbersTStoCLASSmethods,1),length(paramValues));
    for b2 = 1:Bparam
      % Generation of the data      
      groupTS1 = dataExercisesTS(numberExercise,1,n1param,T,paramVector1);
      groupTS1test = dataExercisesTS(numberExercise,1,m1param,T,paramVector1);
      groupTS2 = dataExercisesTS(numberExercise,2,n2param,T,paramVector2);
      groupTS2test = dataExercisesTS(numberExercise,2,m2param,T,paramVector2);
      if contamination
        groupTS1(1:n1cParam,:) = dataExercisesTS(numberExercise,1,n1cParam,T,paramVector1,contamType);
      end
      for p = 1:length(paramValues)
        infoTStoFDmethod = {paramValues(p)};
        % Calculating the functional data and applying the functional classification methods
        for i1 = 1:size(numbersTStoCLASSmethods,1)
          if isequal(i1,1) ||...  % Not to use the previous calculations for i1-1
                  (i1>1 && ~isequal(numbersTStoCLASSmethods(i1,1),numbersTStoCLASSmethods(i1-1,1)))
            startClock = clock;
            %%% From time series to functional data
            groupFD1 = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS1);
            groupFD1test = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS1test);
            groupFD2 = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS2);
            groupFD2test = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS2test);
            timeTStoFD = etime(clock,startClock);
          end
          computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + timeTStoFD;
          computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + timeTStoFD;              
          startClock = clock;
          switch numbersTStoCLASSmethods(i1,2)  % Necessary information for each method
            case 1  % Method DbC
              infoFDtoCLASSmethod = {numbersTStoCLASSmethods(i1,3),numbersTStoCLASSmethods(i1,4),{}};
            case 2  % Method DbC-alpha
              infoFDtoCLASSmethod = {numbersTStoCLASSmethods(i1,3),numbersTStoCLASSmethods(i1,4),alpha,true,alpha};
          end
          [errorG1temp errorG2temp outputFDtoCLASSmethod] =...
              fromFDtoCLASS(numbersTStoCLASSmethods(i1,2),infoFDtoCLASSmethod,groupFD1,groupFD1test,groupFD2,groupFD2test);
          errorRatesGroup1param(b2,i1,p) = errorG1temp/m1param;
          errorRatesGroup2param(b2,i1,p) = errorG2temp/m2param;
          computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + etime(clock,startClock);
          computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + etime(clock,startClock);
          equalDist = outputFDtoCLASSmethod{3};
          if any(equalDist)  % Elements in the frontier
            disp(['      For method ' namesTStoCLASSmethods{i1} ', number of the ' int2str(sum(m1param+m2param)) ' testing elements such that x' int2str(p) ' = 0: ' int2str(sum(equalDist))])
          end
        end
      end
    end
    errorRatesParam = (m1param*errorRatesGroup1param + m2param*errorRatesGroup2param)/(m1param+m2param);
    % Mean error of the runs
    meanErrorRatesParam = zeros(size(numbersTStoCLASSmethods,1),length(paramValues));
    for p = 1:length(paramValues)
      for i1 = 1:size(numbersTStoCLASSmethods,1)
        meanErrorRatesParam(i1,p) = mean(errorRatesParam(:,i1,p),1);
      end
    end
    % Calculations to select the parameter value(s) or the method(s)
    methsXparams = false(size(numbersTStoCLASSmethods,1),length(paramValues));
    methodMinError = zeros(1,size(numbersTStoCLASSmethods,1));  % To register, for each method, the minimum error rate
    for i1 = 1:size(numbersTStoCLASSmethods,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(numbersTStoCLASSmethods,1)
        disp(['      Information for method ' namesTStoCLASSmethods{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
    methsXparams = true(size(numbersTStoCLASSmethods,1),length(paramValues));
  end
  if memoryProblem
    % Generating the time series for this run
    groupTS1 = dataExercisesTS(numberExercise,1,n1,T,paramVector1);
    groupTS1test = dataExercisesTS(numberExercise,1,m1,T,paramVector1);
    groupTS2 = dataExercisesTS(numberExercise,2,n2,T,paramVector2);
    groupTS2test = dataExercisesTS(numberExercise,2,m2,T,paramVector2);
    if contamination
      % Generating the contaminant for this run
      contamGroupTS1 = dataExercisesTS(numberExercise,1,n1c,T,paramVector1,contamType);
    end
  else
    % Taking the time series for this run
    groupTS1 = allGroupTS1(((b-1)*n1+1):(b*n1),:);
    groupTS1test = allGroupTS1test(((b-1)*m1+1):(b*m1),:);
    groupTS2 = allGroupTS2(((b-1)*n2+1):(b*n2),:);
    groupTS2test = allGroupTS2test(((b-1)*m2+1):(b*m2),:);
    if contamination
      % Taking the contaminant for this run
      contamGroupTS1 = allContamGroupTS1(((b-1)*n1c+1):(b*n1c),:);
    end
  end
  if contamination
    % Allocating part of the contaminants in the first rows (for the figures)
    groupTS1(1:n1c2,:) = contamGroupTS1(1:n1c2,:);
    if n1c > 1
      % Substituting the last series of the training sample of group 1 by contaminants
      groupTS1((end-n1c+n1c2+1):end,:) = contamGroupTS1((n1c2+1):end,:);
    end
  end
  if isequal(b,1) && dataFigures  % Figures: plotting the time series
    figure
    plotTimeSeries(1:T,groupTS1test(1:numberElementsToPlot,:),groupTS2test(1:numberElementsToPlot,:),...
        'Some testing series and the mean of the training samples',...
        percentMargin,modelsJointly,mean(groupTS1,1),mean(groupTS2,1))
  end
  % Hypermatrixes to store the multivariate variables
  groupMV1 = zeros(n1,size(numbersTStoCLASSmethods,1),length(paramValues));
  groupMV2 = zeros(n2,size(numbersTStoCLASSmethods,1),length(paramValues));
  groupMV1test = zeros(m1,size(numbersTStoCLASSmethods,1),length(paramValues));
  groupMV2test = zeros(m2,size(numbersTStoCLASSmethods,1),length(paramValues));
  % Main calculations
  for p = 1:length(paramValues)
    if all(methsXparams(:,p)==0)
      continue  % Not to consider this value
    end
    infoTStoFDmethod = {paramValues(p)};
    % Applying the functional classification methods
    for i1 = 1:size(numbersTStoCLASSmethods,1)
      if all(methsXparams(i1,:)==0) || isequal(methsXparams(i1,p),0)
        continue  % Not to consider this method or this value
      end
      switch numbersTStoCLASSmethods(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
      startClock = clock;
      %%% From time series to functional data
      groupFD1 = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS1);
      groupFD1test = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS1test);
      groupFD2 = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS2);
      groupFD2test = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, groupTS2test);
      timeTStoFD = etime(clock,startClock);
      computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + timeTStoFD;
      computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + timeTStoFD;              
      startClock = clock;
      switch numbersTStoCLASSmethods(i1,2)  % Necessary information for each method
        case 1  % Method DbC
          infoFDtoCLASSmethod = {numbersTStoCLASSmethods(i1,3),numbersTStoCLASSmethods(i1,4),{}};
        case 2  % Method DbC-alpha
          infoFDtoCLASSmethod = {numbersTStoCLASSmethods(i1,3),numbersTStoCLASSmethods(i1,4),alpha,true,alpha};
      end
      [errorG1temp errorG2temp outputFDtoCLASSmethod] =...
          fromFDtoCLASS(numbersTStoCLASSmethods(i1,2),infoFDtoCLASSmethod,groupFD1,groupFD1test,groupFD2,groupFD2test);
      errorRatesGroup1(b,i1,p) = errorG1temp;
      errorRatesGroup2(b,i1,p) = errorG2temp;
      computationalTimesRuns(b,i1,p) = computationalTimesRuns(b,i1,p) + etime(clock,startClock);
      computationalTimesRunsWholeParamValues(b,i1) = computationalTimesRunsWholeParamValues(b,i1) + etime(clock,startClock);
      equalDist = outputFDtoCLASSmethod{3};
      if any(equalDist)  % Elements in the frontier
        disp(['For method ' namesTStoCLASSmethods{i1} ', number of the ' int2str(sum(m1+m2)) ' testing elements such that x' int2str(p) ' = 0: ' int2str(sum(equalDist))])
      end
      % Calculating the vector of discriminant variables
      if ~paramOptimization
        infoFDtoMVmethod = {0,'',numbersTStoCLASSmethods(i1,3),numbersTStoCLASSmethods(i1,4),infoFDreference};
        [groupMV1(:,i1,p), groupMV1test(:,i1,p), groupMV2(:,i1,p), groupMV2test(:,i1,p)] =...
            fromFDtoMV(1, infoFDtoMVmethod, groupFD1, groupFD1test, groupFD2, groupFD2test);
      end
      % Figures: the first training curves and the mean of both training samples, and the univariate variable
      % measuring the difference between the distances from some training curves to the mean curves
      if isequal(b,1) && dataFigures && ~paramOptimization
        figure  % Functional data
        xMin = 1;
        xMax = size(groupFD1,2);
        % xRange = xMax - xMin;
        yMin = min(min([groupFD1; groupFD1test; groupFD2; groupFD2test]));
        yMax = max(max([groupFD1; groupFD1test; groupFD2; groupFD2test]));
        yRange = yMax - yMin;
        if isequal(yRange,0)
          yRange = 0.01;
        end
        subplot(2,1,1)
          refGroupFD1 = referenceElement(numbersTStoCLASSmethods(i1,4),infoFDreference,groupFD1);
          refGroupFD2 = referenceElement(numbersTStoCLASSmethods(i1,4),infoFDreference,groupFD2);
          plotFunctions(1:size(groupFD1,2),groupFD1(1:numberElementsToPlot,:),groupFD2(1:numberElementsToPlot,:),...
            ['Some functions of the training and testing samples, method ' namesTStoCLASSmethods{i1}],percentMargin,true,...
            refGroupFD1,refGroupFD2)
          xlabel('Training elements and the reference function')
          ylabel(['For G = ' 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(numbersTStoCLASSmethods(i1,4),infoFDreference,groupFD1test);
          refGroupFD2 = referenceElement(numbersTStoCLASSmethods(i1,4),infoFDreference,groupFD2test);
          plotFunctions(1:size(groupFD1test,2),groupFD1test(1:numberElementsToPlot,:),groupFD2test(1:numberElementsToPlot,:),...
            ' ',percentMargin,true,refGroupFD1,refGroupFD2)
          xlabel('Testing elements and the reference function')
          ylabel(['For G = ' num2str(paramValues(p))])
          set(gca,'XLim',[xMin xMax])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        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 ' namesTStoCLASSmethods{i1}],percentMargin)
          xlabel('Training elements')
          ylabel(['For G = ' 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(['For G = ' 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 ' namesTStoCLASSmethods{i1}])
          ylabel(['For G = ' num2str(paramValues(p))])
      end
    end
  end
  % Calculations to select the parameter value(s) or the method(s)
  for i1 = 1:size(numbersTStoCLASSmethods,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(numbersTStoCLASSmethods,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
  if isequal(b,1) && dataFigures  % Figures: results of this run
    % 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,namesTStoCLASSmethods,computationalTimesRuns(1,:,:),...
        computationalTimesRunsWholeParamValues(1,:),paramOptimization,methodSelection,[],[]);
  end
  if verbose
    for i1 = 1:size(numbersTStoCLASSmethods,1)
      if all(methsXparams(i1,:)==0)  % Not to consider this method
        continue
      end
      disp(['    Information for method ' namesTStoCLASSmethods{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))])
    end
  end
end

errorRates = (errorRatesGroup1 + errorRatesGroup2)/(m1 + m2);
errorRatesGroup1 = errorRatesGroup1/m1;
errorRatesGroup2 = errorRatesGroup2/m2;

if verbose
  disp(' ') 
  disp('... ending 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(numbersTStoCLASSmethods,1)
    disp(['METHOD: ' namesTStoCLASSmethods{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))])
  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,namesTStoCLASSmethods,...
      computationalTimesRuns,computationalTimesRunsWholeParamValues,...
      paramOptimization,methodSelection);
end

if saveResults
  save(['resultsSimulationExercise' int2str(numberExercise)], 'errorRatesGroup1','errorRatesGroup2',...
      'errorRates','minimizingValuesErrorRatesGroup1','minimizingValuesErrorRatesGroup2','minimizingValuesErrorRates',...
      'minimizingValuesMethodsErrorRatesGroup1','minimizingValuesMethodsErrorRatesGroup2','minimizingValuesMethodsErrorRates',...
      'paramValues','methsXparams','paramValuesMinpower','methodsMinpower','percentMargin',...
      'namesTStoCLASSmethods','computationalTimesRuns','computationalTimesRunsWholeParamValues','paramOptimization',...
      'methodSelection')
  if verbose
    disp('The results have been saved')
    disp(' ')
  end
end
% load(['resultsSimulationExercise' int2str(numberExercise)])

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


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

clear i2 Bparam infoFDreference p numberBlocks Sw allGroupTS1test minG paramVector1 paramVector2...
    m1 allGroupTS2test m1param alpha m2 groupY2test coefficients m2param...
    memoryProblem percentMargin posts postsMin refGroupFD1 contamType...
    refGroupFD2 contamination dyadicSplits stdvarsIncluded errorG1temp timeTStoFD errorG2temp modelsJointly...
    n1 vars n1c varsExcluded n1cParam varsIncluded groupMV1 n1param groupMV1test n2 xMax groupMV2...
    n2param xMin groupMV2test yMax groupY1c numberElementsToPlot yMin groupY1test numberExercise yRange...
    groupY2c verbose dataFigures resultsFigures methodMinError
