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

% PURPOSE: This function generates boxplots from matrices with information
%          like the error rates, minimizing-power, computational times.

% INPUT VARIABLES:
%    errorRatesGroup1: Hypermatrix with the misclassification error rates of the testing sample of group 1.
%    errorRatesGroup2: Hypermatrix with the misclassification error rates of the testing sample of group 2.
%    errorRates: Hypermatrix with the misclassification error rates of all the testing data.
%    minimizingValuesErrorRatesGroup1: For group 1, matrix with the (minimum) errors when the parameter has been optimized
%    minimizingValuesErrorRatesGroup2: For group 2, matrix with the (minimum) errors when the parameter has been optimized
%    minimizingValuesErrorRates: Matrix with the (minimum) global errors when the parameter has been optimized
%    minimizingValuesMethodsErrorRatesGroup1: For group 1, matrix with the (minimum) errors when both the method and the parameter have been optimized
%    minimizingValuesMethodsErrorRatesGroup2: For group 2, matrix with the (minimum) errors when both the method and the parameter have been optimized
%    minimizingValuesMethodsErrorRates: Matrix with the (minimum) global errors when both the method and the parameter have been optimized
%    paramValues: Vector with the values of the parameter.
%    paramValuesMinpower: Hypermatrix with the minimizing-power measure of the parameter values.
%    methodsMinpower: Hypermatrix with the minimizing-power of the methods.
%    percentMargin: Optional variable to control the distance from the more
%                   external data to the margin of the axes.
%    namesDATAtoCLASSmethods: Cell array with the names of the global methods, to be used in the figures.
%    computationalTimesRuns: Hypermatrix with the computational times for each value of the parameter.
%    computationalTimesRunsWholeParamValues: Matrix with the joint computational times for all the values of the parameter.
%    paramOptimization: Logical variable to know if the parameter has been optimized.
%    methodOptimization: Logical variable to know if the method has been optimized.
%    coefficientsMVmethods: Possible coefficients of the multivariate method.
%    computationalTimesRunsMVmethod: Matrix with the computational times of the multivariate methods.
%    application: Logical variable that indicates if the call comes from the applications script.


function plotResults(errorRatesGroup1, errorRatesGroup2, errorRates,...
    minimizingValuesErrorRatesGroup1, minimizingValuesErrorRatesGroup2, minimizingValuesErrorRates,...
    minimizingValuesMethodsErrorRatesGroup1, minimizingValuesMethodsErrorRatesGroup2, minimizingValuesMethodsErrorRates,...
    paramValues, paramValuesMinpower, methodsMinpower, percentMargin,...
    namesDATAtoCLASSmethods, computationalTimesRuns, computationalTimesRunsWholeParamValues,...
    paramOptimization, methodOptimization, coefficientsMVmethods, computationalTimesRunsMVmethod, application)

B = size(computationalTimesRuns,1);  % Number of measurements in the hypermatrixes
addTitle = ['(based on ', num2str(B), ' runs) '];  % Real value of B

if ~exist('coefficientsMVmethods','var') || isempty(coefficientsMVmethods)
  coefficientsMVmethods = zeros(B,length(namesDATAtoCLASSmethods),length(paramValues));
end
if ~exist('computationalTimesRunsMVmethod','var') || isempty(computationalTimesRunsMVmethod)
  computationalTimesRunsMVmethod = zeros(B,length(namesDATAtoCLASSmethods));
end
if ~exist('application','var')
  application = false;
end

numberErrors = size(errorRates,3);

% Trick for the boxplots: duplicating the unique measurement
if isequal(B,1)
  if ~application
    ErrorRatesGroup1 = zeros(2,length(namesDATAtoCLASSmethods),numberErrors);
    ErrorRatesGroup2 = zeros(2,length(namesDATAtoCLASSmethods),numberErrors);
    ErrorRates = zeros(2,length(namesDATAtoCLASSmethods),numberErrors);
    MinimizingValuesErrorRatesGroup1 = zeros(2,length(namesDATAtoCLASSmethods));
    MinimizingValuesErrorRatesGroup2 = zeros(2,length(namesDATAtoCLASSmethods));
    MinimizingValuesErrorRates = zeros(2,length(namesDATAtoCLASSmethods));
    MinimizingValuesMethodsErrorRatesGroup1 = zeros(2,1);
    MinimizingValuesMethodsErrorRatesGroup2 = zeros(2,1);
    MinimizingValuesMethodsErrorRates = zeros(2,1);
  end
  ParamValuesMinpower = zeros(2,length(namesDATAtoCLASSmethods),length(paramValues));
  MethodsMinpower = zeros(2,length(namesDATAtoCLASSmethods));
  ComputationalTimesRuns = zeros(2,length(namesDATAtoCLASSmethods),length(paramValues));
  ComputationalTimesRunsWholeParamValues = zeros(2,length(namesDATAtoCLASSmethods));
  CoefficientsMVmethods = zeros(2,length(namesDATAtoCLASSmethods),length(paramValues));
  ComputationalTimesRunsMVmethod = zeros(2,length(namesDATAtoCLASSmethods));
  for i1 = 1:2
    if ~application
      ErrorRatesGroup1(i1,:,:) = errorRatesGroup1;
      ErrorRatesGroup2(i1,:,:) = errorRatesGroup2;
      ErrorRates(i1,:,:) = errorRates;
      MinimizingValuesErrorRatesGroup1(i1,:,:) = minimizingValuesErrorRatesGroup1;
      MinimizingValuesErrorRatesGroup2(i1,:,:) = minimizingValuesErrorRatesGroup2;
      MinimizingValuesErrorRates(i1,:,:) = minimizingValuesErrorRates;
      MinimizingValuesMethodsErrorRatesGroup1(i1) = minimizingValuesMethodsErrorRatesGroup1;
      MinimizingValuesMethodsErrorRatesGroup2(i1) = minimizingValuesMethodsErrorRatesGroup2;
      MinimizingValuesMethodsErrorRates(i1) = minimizingValuesMethodsErrorRates;
    end
    ParamValuesMinpower(i1,:,:) = paramValuesMinpower;
    MethodsMinpower(i1,:) = methodsMinpower;
    ComputationalTimesRuns(i1,:,:) = computationalTimesRuns;
    ComputationalTimesRunsWholeParamValues(i1,:) = computationalTimesRunsWholeParamValues;
    CoefficientsMVmethods(i1,:,:) = coefficientsMVmethods;
    ComputationalTimesRunsMVmethod(i1,:) = computationalTimesRunsMVmethod;
  end
  if ~application
    errorRatesGroup1 = ErrorRatesGroup1;
    errorRatesGroup2 = ErrorRatesGroup2;
    errorRates = ErrorRates;
    minimizingValuesErrorRatesGroup1 = MinimizingValuesErrorRatesGroup1;
    minimizingValuesErrorRatesGroup2 = MinimizingValuesErrorRatesGroup2;
    minimizingValuesErrorRates = MinimizingValuesErrorRates;
    minimizingValuesMethodsErrorRatesGroup1 = MinimizingValuesMethodsErrorRatesGroup1;
    minimizingValuesMethodsErrorRatesGroup2 = MinimizingValuesMethodsErrorRatesGroup2;
    minimizingValuesMethodsErrorRates = MinimizingValuesMethodsErrorRates;
  end
  paramValuesMinpower = ParamValuesMinpower;
  methodsMinpower = MethodsMinpower;
  computationalTimesRuns = ComputationalTimesRuns;
  computationalTimesRunsWholeParamValues = ComputationalTimesRunsWholeParamValues;
  coefficientsMVmethods = CoefficientsMVmethods;
  computationalTimesRunsMVmethod = ComputationalTimesRunsMVmethod;
  if ~application
    clear ErrorRatesGroup1 ErrorRatesGroup2 ErrorRates MinimizingValuesErrorRatesGroup1...
        MinimizingValuesErrorRatesGroup2 MinimizingValuesErrorRates MinimizingValuesMethodsErrorRatesGroup1...
        MinimizingValuesMethodsErrorRatesGroup2 MinimizingValuesMethodsErrorRates
  end
  clear ParamValuesMinpower MethodsMinpower ComputationalTimesRuns...
      ComputationalTimesRunsWholeParamValues CoefficientsMVmethods ComputationalTimesRunsMVmethod
  B = 2;  % Number of measurements in the hypermatrixes
end

newParamValues = {};
for i1 = 1:length(paramValues)
  newParamValues = [newParamValues num2str(paramValues(i1))];
end
if isequal(size(errorRates,3),length(paramValues)+1)
  newParamValues = [newParamValues 'MV'];
end

for i1 = 1:length(namesDATAtoCLASSmethods)
  % Error rates for each global method
  if ~application
    if paramOptimization
      if methodOptimization
        yMin = min([minimizingValuesMethodsErrorRatesGroup1; minimizingValuesMethodsErrorRatesGroup2; minimizingValuesMethodsErrorRates]);
        yMax = max([minimizingValuesMethodsErrorRatesGroup1; minimizingValuesMethodsErrorRatesGroup2; minimizingValuesMethodsErrorRates]);
        yRange = yMax - yMin;
        if isequal(yRange,0)
          yRange = 0.01;
        end
        figure
        subplot(3,1,1)
          boxplot(minimizingValuesMethodsErrorRatesGroup1, 'labels', ' ', 'notch', 'on')
          xlabel('Group 1')
          title(['Misclassification error rates ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])  
        subplot(3,1,2)
          boxplot(minimizingValuesMethodsErrorRatesGroup2, 'labels', ' ', 'notch', 'on')
          xlabel('Group 2')
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])  
        subplot(3,1,3)
          boxplot(minimizingValuesMethodsErrorRates, 'labels', ' ', 'notch', 'on')
          xlabel('Global')
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
      else
        yMin = min(min([minimizingValuesErrorRatesGroup1; minimizingValuesErrorRatesGroup2; minimizingValuesErrorRates]));
        yMax = max(max([minimizingValuesErrorRatesGroup1; minimizingValuesErrorRatesGroup2; minimizingValuesErrorRates]));
        yRange = yMax - yMin;
        if isequal(yRange,0)
          yRange = 0.01;
        end
        figure
        subplot(3,1,1)
          boxplot(minimizingValuesErrorRatesGroup1(:,i1), 'labels', ' ', 'notch', 'on')
          xlabel('Group 1')
          title(['Misclassification error rates ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])  
        subplot(3,1,2)
          boxplot(minimizingValuesErrorRatesGroup2(:,i1), 'labels', ' ', 'notch', 'on')
          xlabel('Group 2')
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])  
        subplot(3,1,3)
          boxplot(minimizingValuesErrorRates(:,i1), 'labels', ' ', 'notch', 'on')
          xlabel('Global')
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])  
      end
      yMin = min(min(min(paramValuesMinpower)));
      yMax = max(max(max(paramValuesMinpower)));
      yRange = yMax - yMin;
      if isequal(yRange,0)
        yRange = 0.01;
      end
      figure
      boxplot(reshape(paramValuesMinpower(:,i1,:), B, length(paramValues)), 'labels', paramValues, 'notch', 'on')
      title(['Minimizing-power of the parameter values ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
      set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
    else
      yMin = min([min(min(min(errorRatesGroup1))) min(min(min(errorRatesGroup2))) min(min(min(errorRates)))]);
      yMax = max([max(max(max(errorRatesGroup1))) max(max(max(errorRatesGroup2))) max(max(max(errorRates)))]);
      yRange = yMax - yMin;
      if isequal(yRange,0)
        yRange = 0.01;
      end
      figure
      subplot(3,1,1)
        boxplot(reshape(errorRatesGroup1(:,i1,:), B, numberErrors), 'labels', newParamValues, 'notch', 'on')
        xlabel('Group 1')
        title(['Misclassification error rates ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
        set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
      subplot(3,1,2)
        boxplot(reshape(errorRatesGroup2(:,i1,:), B, numberErrors), 'labels', newParamValues, 'notch', 'on')
        xlabel('Group 2')
        set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
      subplot(3,1,3)
        boxplot(reshape(errorRates(:,i1,:), B, numberErrors), 'labels', newParamValues, 'notch', 'on')
        xlabel('Global')
        set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
      if length(paramValues)>1
        xMin = paramValues(1);
        xMax = paramValues(end);
        % xRange = xMax - xMin;
        figure
        subplot(3,1,1)
          plot(paramValues,mean(reshape(errorRatesGroup1(:,i1,1:length(paramValues)), B, length(paramValues)),1),'k','LineWidth',2)
          xlabel('Group 1')
          title(['Evolution of the mean error rates ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
          set(gca,'XLim',[xMin-0.25 xMax+0.25])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        subplot(3,1,2)
          plot(paramValues,mean(reshape(errorRatesGroup2(:,i1,1:length(paramValues)), B, length(paramValues)),1),'k','LineWidth',2)
          xlabel('Group 2')
          set(gca,'XLim',[xMin-0.25 xMax+0.25])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        subplot(3,1,3)
          plot(paramValues,mean(reshape(errorRates(:,i1,1:length(paramValues)), B, length(paramValues)),1),'k','LineWidth',2)
          xlabel('Global')
          set(gca,'XLim',[xMin-0.25 xMax+0.25])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
      end
    end
  end

  % Figure for the coefficients of the multivariate methods
  if (~paramOptimization) && (~all(all(coefficientsMVmethods(:,i1,:)==0)))
    yMin = min(min(min(coefficientsMVmethods)));
    yMax = max(max(max(coefficientsMVmethods)));
    yRange = yMax - yMin;
    if isequal(yRange,0)
      yRange = 0.01;
    end
    figure
    boxplot(reshape(coefficientsMVmethods(:,i1,:), B, length(paramValues)), 'labels', paramValues, 'notch', 'on')
    title(['Coefficients of the multivariate method ', addTitle, 'of ', namesDATAtoCLASSmethods{i1}])
    set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
  end

  % Figure for the computation times
  yMin = min([min(min(min(computationalTimesRuns))) min(min(computationalTimesRunsWholeParamValues))...
      min(min(computationalTimesRunsMVmethod))]);
  yMax = max([max(max(max(computationalTimesRuns))) max(max(computationalTimesRunsWholeParamValues))...
      max(max(computationalTimesRunsMVmethod))]);
  yRange = yMax - yMin;
  if isequal(yRange,0)
    yRange = 0.01;
  end
  figure
  if length(paramValues)>1
    subplot(2,1,1)
      boxplot(reshape(computationalTimesRuns(:,i1,:), B, length(paramValues)), 'labels', paramValues, 'notch', 'on')
      title(['Time spent for each value of the parameter ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
      set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
    subplot(2,1,2)
      if paramOptimization
        boxplot(computationalTimesRunsWholeParamValues(:,i1), 'labels', 'Vector of parameter values', 'notch', 'on')
        title(['Joint time spent for all values ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
        set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
      else
        if isequal(size(errorRates,3),length(paramValues)+1)
          boxplot([computationalTimesRunsWholeParamValues(:,i1) computationalTimesRunsMVmethod(:,i1)],...
              'labels', {'Vector of parameter values' 'MV method'}, 'notch', 'on')
          title(['Joint time spent for all values ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        else
          boxplot(computationalTimesRunsWholeParamValues(:,i1),...
              'labels', 'Vector of parameter values', 'notch', 'on')
          title(['Joint time spent for all values ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
          set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
        end
      end
  else
    boxplot(reshape(computationalTimesRuns(:,i1,:), B, length(paramValues)), 'labels', paramValues, 'notch', 'on')
    title(['Time spent for each value of the parameter ', addTitle, 'with ', namesDATAtoCLASSmethods{i1}])
    set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
  end
  
end

% Minimizing-power of the methods
if methodOptimization
  yMin = min(min(methodsMinpower));
  yMax = max(max(methodsMinpower));
  yRange = yMax - yMin;
  if isequal(yRange,0)
    yRange = 0.01;
  end
  figure
  boxplot(methodsMinpower, 'labels', namesDATAtoCLASSmethods, 'notch', 'on')
  title(['Minimizing-power of the methods ' addTitle])
  set(gca,'YLim',[yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
end

end
