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

%%% REAL DATA

% Training samples
  allGroupTS1 = ;
  allGroupTS2 = ;

% New elements to be classified
  newDataTS = ;


% %%% An example with pseudo-real data
%   numberExercise = 2;
%     paramVector1 = [-0.1 0];  paramVector2 = [+0.1  0];
%   T = 64;
%   % Training samples
%     N1= 12; N2 = 13;
%     allGroupTS1 = dataExercisesTS(numberExercise,1,N1,T,paramVector1);
%     allGroupTS2 = dataExercisesTS(numberExercise,2,N2,T,paramVector2);
%   % New data (M elements of each group)
%     M = 3;
%     newDataTS = [dataExercisesTS(numberExercise,1,M,T,paramVector1); dataExercisesTS(numberExercise,2,M,T,paramVector2)];
%   clear numberExercise T N1 N2 M paramVector1 paramVector2


% Preprocess of the real data: standardization
for i1 = 1:size(newDataTS,1)
  newDataTS(i1,:) = (newDataTS(i1,:)-mean(newDataTS(i1,:),2))/std(newDataTS(i1,:));
end
for i1 = 1:size(allGroupTS1,1)
  allGroupTS1(i1,:) = (allGroupTS1(i1,:)-mean(allGroupTS1(i1,:),2))/std(allGroupTS1(i1,:));
end
for i1 = 1:size(allGroupTS2,1)
  allGroupTS2(i1,:) = (allGroupTS2(i1,:)-mean(allGroupTS2(i1,:),2))/std(allGroupTS2(i1,:));
end
 
%%%%%%%%%%%%%%%%  Ending Data Loading  %%%%%%%%%%%%%%%%%%%%




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

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

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 number of blocks. 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(newDataTS,1);
N1 = size(allGroupTS1,1);
N2 = size(allGroupTS2,1);
m1 = M;
m2 = M;
T = size(allGroupTS1,2);
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;
% All functional transformations are calculated here
disp(' ')
disp('Calculating the functional data...')
disp(' ')
computationalTimesTStoFD = zeros(size(numbersTStoCLASSmethods,1),length(paramValues));
for p = 1:length(paramValues)
  infoTStoFDmethod = {paramValues(p)};
  for i1 = 1:size(numbersTStoCLASSmethods,1)
    if isequal(i1,1) || (i1>1 && ~any(numbersTStoCLASSmethods(i1,1)==numbersTStoCLASSmethods(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)
        eval(['allGroupFD1_' num2str(p) '_' num2str(i1) ' = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, allGroupTS1);'])
        eval(['allGroupFD2_' num2str(p) '_' num2str(i1) ' = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, allGroupTS2);'])
        eval(['newDataFD_' num2str(p) '_' num2str(i1) ' = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, newDataTS);'])
      end
      startClock = clock;
      eval(['allGroupFD1_' num2str(p) '_' num2str(i1) ' = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, allGroupTS1);'])
      eval(['allGroupFD2_' num2str(p) '_' num2str(i1) ' = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, allGroupTS2);'])
      eval(['newDataFD_' num2str(p) '_' num2str(i1) ' = fromTStoFD(numbersTStoCLASSmethods(i1,1), infoTStoFDmethod, newDataTS);'])
      computationalTimesTStoFD(i1,p) = computationalTimesTStoFD(i1,p) + etime(clock,startClock);
    else   % To use previous calculations, since i1-1, i1-2,... were also considered
      posts = find(numbersTStoCLASSmethods(1:(i1-1),1)==numbersTStoCLASSmethods(i1,1));
      computationalTimesTStoFD(i1,p) = computationalTimesTStoFD(posts(1),p);
      eval(['allGroupFD1_' num2str(p) '_' num2str(i1) ' = allGroupFD1_' num2str(p) '_' num2str(posts(1)) ';'])
      eval(['allGroupFD2_' num2str(p) '_' num2str(i1) ' = allGroupFD2_' num2str(p) '_' num2str(posts(1)) ';'])
      eval(['newDataFD_' num2str(p) '_' num2str(i1) ' = newDataFD_' num2str(p) '_' num2str(posts(1)) ';'])
    end
  end
end
% All the depth measurements
disp(' ')
disp('Measuring the depth of the functional data...')
disp(' ')
computationalTimesDepth = zeros(size(numbersTStoCLASSmethods,1),length(paramValues));
if depthBasedRobustness
  if any(numbersTStoCLASSmethods(:,2) == 2)  % If included, method 2 is removed
    namesTStoCLASSmethods(numbersTStoCLASSmethods(:,2) == 2) = [];
    numbersTStoCLASSmethods(numbersTStoCLASSmethods(:,2) == 2,:) = [];
    disp(' ')
    disp('Since the deepest elements have been identified, method 2 (DbC-alpha) has been disabled!')
    disp(' ')
  end
  deepestFpositions1 = zeros(N1-ceil(alpha*N1),size(numbersTStoCLASSmethods,1),length(paramValues));
  deepestFpositions2 = zeros(N2-ceil(alpha*N2),size(numbersTStoCLASSmethods,1),length(paramValues));
  for p = 1:length(paramValues)
    for i1 = 1:size(numbersTStoCLASSmethods,1)
      if isequal(p,1)  % The following lines are a trick to avoid a perturbation in measuring the computational time (see the helping file)
        eval(['[litter deepestFpositions1(:,i1,p)] = trimFunction(allGroupFD1_' num2str(p) '_' num2str(i1) ', alpha);']);
        eval(['[litter deepestFpositions2(:,i1,p)] = trimFunction(allGroupFD2_' num2str(p) '_' num2str(i1) ', alpha);']);
      end
      if isequal(i1,1) || (i1>1 && ~any(numbersTStoCLASSmethods(i1,1)==numbersTStoCLASSmethods(1:(i1-1),1)))
        startClock = clock;
        % Identifying--but not removing--the deepest functions of each group
        eval(['[litter deepestFpositions1(:,i1,p)] = trimFunction(allGroupFD1_' num2str(p) '_' num2str(i1) ', alpha);']);
        eval(['[litter deepestFpositions2(:,i1,p)] = trimFunction(allGroupFD2_' num2str(p) '_' num2str(i1) ', 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(numbersTStoCLASSmethods(1:(i1-1),1)==numbersTStoCLASSmethods(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(numbersTStoCLASSmethods,1)
    disp(' ')
    disp(['Mean computational times spent calculating the depth for each parameter values, method ' namesTStoCLASSmethods{i1}])
    disp(['   '  num2str(computationalTimesDepth(i1,:))])
    disp(' ')
  end
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  %%%%%%%%%%%%%%%%%%%%

% Hypermatrixes to store the misclassifications, the minimizing-power measures and the computational times
classificationOutputs = zeros(M*B,size(numbersTStoCLASSmethods,1),length(paramValues));  % N subhypermatrixes of height B
minimizingValuesClassificationOutputs = zeros(M*B,size(numbersTStoCLASSmethods,1));  % N matrixes of height B
minimizingValuesMethodsClassificationOutputs = zeros(M*B,1);  % N vectors of height B
paramValuesMinpower = zeros(B,size(numbersTStoCLASSmethods,1),length(paramValues));
methodsMinpower = zeros(B,size(numbersTStoCLASSmethods,1));
computationalTimesRuns = zeros(B,size(numbersTStoCLASSmethods,1),length(paramValues));
computationalTimesRunsWholeParamValues = zeros(B,size(numbersTStoCLASSmethods,1));
for i1 = 1:size(numbersTStoCLASSmethods,1)
  for p = 1:length(paramValues)
    computationalTimesRuns(:,i1,p) =...
        computationalTimesRuns(:,i1,p) + computationalTimesTStoFD(i1,p) + computationalTimesDepth(i1,p);
  end
  computationalTimesRunsWholeParamValues(:,i1) =...
      computationalTimesRunsWholeParamValues(:,i1) + sum(computationalTimesTStoFD(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(numbersTStoCLASSmethods,1),length(paramValues));
    errorRatesGroup2param = zeros(Bparam,size(numbersTStoCLASSmethods,1),length(paramValues));
    trustworthyRunsParam = true(Bparam,size(numbersTStoCLASSmethods,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);
      % Method calculations here with all possible values of the parameter(s)
      for p = 1:length(paramValues)
        infoTStoFDmethod = {paramValues(p)};
        % Calculating the functional data and applying the functional classification methods
        for i1 = 1:size(numbersTStoCLASSmethods,1)
          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));
            counter1 = length(indexes);
            positionsTestingGroup1paramTemp(indexes) = [];
            [litter, indexes] = setdiff(positionsTestingGroup2paramTemp,deepestFpositions2(:,i1,p));
            counter2 = length(indexes);
            positionsTestingGroup2paramTemp(indexes) = [];
            % Removing the untrustworthy training elements (maybe different for each method)
            [litter, indexes] = setdiff(positionsTrainingGroup1paramTemp,deepestFpositions1(:,i1,p));
            positionsTrainingGroup1paramTemp(indexes) = [];
            [litter, indexes] = setdiff(positionsTrainingGroup2paramTemp,deepestFpositions2(:,i1,p));
            positionsTrainingGroup2paramTemp(indexes) = [];
            % Checking if, after removing some elements, this run is reliable for this method
            if isempty([positionsTestingGroup1paramTemp; positionsTestingGroup2paramTemp]) ||...
                    isempty(positionsTrainingGroup1paramTemp)  || isempty(positionsTrainingGroup2paramTemp)
              trustworthyRunsParam(b2,i1,p) = false;  % To register this run
              disp(['      Untrustworthy run for [b b2 g i1] = ' num2str([b b2 g i1])])
              continue  % To skip this run
            end
            %%% From time series to functional data
            startClock = clock;
            eval(['groupFD1 = allGroupFD1_' num2str(p) '_' num2str(i1) '(positionsTrainingGroup1paramTemp,:);'])
            eval(['groupFD1test = allGroupFD1_' num2str(p) '_' num2str(i1) '(positionsTestingGroup1paramTemp,:);'])
            eval(['groupFD2 = allGroupFD2_' num2str(p) '_' num2str(i1) '(positionsTrainingGroup2paramTemp,:);'])
            eval(['groupFD2test = allGroupFD2_' num2str(p) '_' num2str(i1) '(positionsTestingGroup2paramTemp,:);'])
            timeTStoFD = etime(clock,startClock);
          else
            %%% From time series to functional data
            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;
              eval(['groupFD1 = allGroupFD1_' num2str(p) '_' num2str(i1) '(positionsTrainingGroup1paramTemp,:);'])
              eval(['groupFD1test = allGroupFD1_' num2str(p) '_' num2str(i1) '(positionsTestingGroup1paramTemp,:);'])
              eval(['groupFD2 = allGroupFD2_' num2str(p) '_' num2str(i1) '(positionsTrainingGroup2paramTemp,:);'])
              eval(['groupFD2test = allGroupFD2_' num2str(p) '_' num2str(i1) '(positionsTestingGroup2paramTemp,:);'])
              timeTStoFD = etime(clock,startClock);
            end
            counter1 = 0;
            counter2 = 0;
          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-counter1);
          errorRatesGroup2param(b2,i1,p) = errorG2temp/(m2param-counter2);
          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);    % Note that m1param = 0 if leaveOneOut = true
    % Mean error of the trustworthy runs
    meanErrorRatesParam = zeros(size(numbersTStoCLASSmethods,1),length(paramValues));
    for p = 1:length(paramValues)
      for i1 = 1:size(numbersTStoCLASSmethods,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(numbersTStoCLASSmethods,1),length(paramValues));
    methodMinError = zeros(1,size(numbersTStoCLASSmethods,1));  % To register, for each method, the minimum error rate
    minValue = zeros(1,size(numbersTStoCLASSmethods,1)); % To register, for each method, the smallest minimizing parameter value (as a reference to calculate some quantities later)
    for i1 = 1:size(numbersTStoCLASSmethods,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(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: one or several
    minMethod = posts(1);  % Smallest minimizing method (as a reference to calculate some quantities later)
    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
    % Selection of the final training groups
    positionsFinalTrainingGroup1 = positionsTrainingGroup1;
    positionsFinalTrainingGroup2 = positionsTrainingGroup2;
    methsXparams = true(size(numbersTStoCLASSmethods,1),length(paramValues));
  end
  if isequal(b,1) && dataFigures  % Figures: plotting the time series
    groupTS1 = allGroupTS1(positionsFinalTrainingGroup1,:);
    groupTS1test = newDataTS;
    groupTS2 = allGroupTS2(positionsFinalTrainingGroup2,:);
    groupTS2test = newDataTS; % = [];
    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))
    clear groupTS1 groupTS1test groupTS2 groupTS2test
  end
  n1 = length(positionsFinalTrainingGroup1); % New sample size 1
  n2 = length(positionsFinalTrainingGroup2); % New sample size 2
  % Hypermatrixes to store the multivariate variables
  groupMV1 = zeros(n1,size(numbersTStoCLASSmethods,1),length(paramValues));
  groupMV1test = zeros(m1,size(numbersTStoCLASSmethods,1),length(paramValues));
  groupMV2 = zeros(n2,size(numbersTStoCLASSmethods,1),length(paramValues));
  groupMV2test = zeros(m2,size(numbersTStoCLASSmethods,1),length(paramValues));
  % 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
    infoTStoFDmethod = {paramValues(p)};
    % Calculating the functional data and 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
      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 time series to functional data
      startClock = clock;
      eval(['groupFD1 = allGroupFD1_' num2str(p) '_' num2str(i1) '(positionsFinalTrainingGroup1Temp,:);'])
      eval(['groupFD1test = newDataFD_' num2str(p) '_' num2str(i1) ';'])
      eval(['groupFD2 = allGroupFD2_' num2str(p) '_' num2str(i1) '(positionsFinalTrainingGroup2Temp,:);'])
      eval(['groupFD2test = newDataFD_' num2str(p) '_' num2str(i1) ';'])
      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);
      switch numbersTStoCLASSmethods(i1,2)  % Necessary information for each method
        case 1  % Method DbC
          pseudoMisclass1 = outputFDtoCLASSmethod{1};
          pseudoMisclass2 = outputFDtoCLASSmethod{2};
        case 2  % Method DbC-alpha
          pseudoMisclass1 = outputFDtoCLASSmethod{1};
          pseudoMisclass2 = outputFDtoCLASSmethod{2};
      end
      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);
      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(indexesTrainingGroup1,i1,p), groupMV1test(:,i1,p),...
            groupMV2(indexesTrainingGroup2,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
  if paramOptimization
    % Registering the labels when the parameter value and the method are/were optimized
    for i1 = 1:size(numbersTStoCLASSmethods,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));
  end
  if isequal(b,1) && dataFigures  % Figures: results of this run
    % Results
    plotResults([],[],[],[],[],[],[],[],[],...
      paramValues, paramValuesMinpower(1,:,:), methodsMinpower(1,:), percentMargin, namesTStoCLASSmethods,...
      computationalTimesRuns(1,:,:), computationalTimesRunsWholeParamValues(1,:), paramOptimization, methodSelection,...
      [], [], true);
  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}])
      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))])
    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(numbersTStoCLASSmethods,1)
    disp(['METHOD: ' namesTStoCLASSmethods{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))])
  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, namesTStoCLASSmethods,...
    computationalTimesRuns, computationalTimesRunsWholeParamValues, paramOptimization, methodSelection,...
    [], [], true);
end

if saveResults
  save('resultsApplication','paramValues','methsXparams','paramValuesMinpower','methodsMinpower',...
      'percentMargin','namesTStoCLASSmethods','computationalTimesRuns','computationalTimesRunsWholeParamValues',...
      'paramOptimization','methodSelection')
  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 numberBlocks groupFD1 groupFD1test N1 groupFD2 numberElementsToPlot...
   N1param groupFD2test N2 groupMV1 N2param groupMV1test...
   outputFDtoCLASSmethod groupMV2 p groupMV2test paramDataReuse groupY1test groupY2test i1 indexes...
   infoFDreference infoFDtoCLASSmethod 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...
   errorG2temp errorRatesGroup1param verbose dataFigures resultsFigures
   for p = 1:length(paramValues)
     for i1 = 1:size(numbersTStoCLASSmethods,1)
       eval(['clear allGroupFD1_' num2str(p) '_' num2str(i1) ' allGroupFD2_' num2str(p) '_' num2str(i1)])
       eval(['clear newDataFD_' num2str(p) '_' num2str(i1)])
     end
   end
