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

% PURPOSE: This function plots, for two samples, the values of one vector of variables.

% INPUT VARIABLES:
%    groupMV1: Matrix with the values of the discriminat variables for the
%              training elements of sample 1 as rows.
%    groupMV1test: Matrix with the values of the discriminat variables for the
%                  testing elements of sample 1 as rows.
%    groupMV2: Matrix with the values of the discriminat variables for the
%              training elements of sample 2 as rows.
%    groupMV2test: Matrix with the values of the discriminat variables for the
%                  testing elements of sample 2 as rows.
%    figureTitle: Optional variable with the title of the figure.
%    percentMargin: Optional variable to control the distance from the more
%                   external data to the axes.
%    coefficients: Optional coefficients of the variables (as a column).

% COMMENTS:
%   (1) This function uses training and testing information. If desired, the user can
%       call it with the same values for both of them.

function plotDiscrimVariables(groupMV1, groupMV1test, groupMV2, groupMV2test, figureTitle, percentMargin, coefficients)

if ~exist('figureTitle','var') || isempty(figureTitle)
  figureTitle = 'Discriminant Variables';
end
if ~exist('percentMargin','var') || isempty(percentMargin)
  percentMargin = 5;
end
if ~exist('letter','var') || isempty(letter)
  letter = '';
end
if exist('coefficients','var')
  if isequal(length(size(coefficients)),3)  % Three dimensions
    coefficients = reshape(coefficients,size(coefficients,3),1);  % To build a vector
  end
  if all(coefficients==0)  % All are 0's  (e.g., for KNN)
    clear coefficients  % For posterior calculations
  end
end

switch size(groupMV1,2)
  case 1
    xMin = 1;
    xMax = max([length(groupMV1test) length(groupMV2test)]);
    xRange = xMax - xMin;
    if isequal(xRange,0)
      xRange = 0.01;
    end
    yMin = min([min(groupMV1test), min(groupMV2test), mean(groupMV1,1), mean(groupMV2,1)]);
    yMax = max([max(groupMV1test), max(groupMV2test), mean(groupMV1,1), mean(groupMV2,1)]);
    yRange = yMax - yMin;
    if isequal(yRange,0)
      yRange = 0.01;
    end
    plot(groupMV1test,'bo','MarkerSize',3,'LineWidth',1)  % Testing elements of group 1
    hold on
    plot(groupMV2test,'rs','MarkerSize',3,'LineWidth',1)  % Testing elements of group 2
    line([1 xMax],[mean(groupMV1test,1) mean(groupMV1test,1)],'Color','b','LineStyle',':','LineWidth',1)  % Mean of the testing elements of group 1
    line([1 xMax],[mean(groupMV2test,1) mean(groupMV2test,1)],'Color','r','LineStyle',':','LineWidth',1)  % Mean of the testing elements of group 2
    line([1 xMax],[mean(groupMV1,1) mean(groupMV1,1)],'Color','b','LineWidth',1)  % Mean of the training elements of group 1
    line([1 xMax],[mean(groupMV2,1) mean(groupMV2,1)],'Color','r','LineWidth',1)  % Mean of the training elements of group 2
    hold off
    xlabel('Testing elements')
    if exist('p','var') && ~isempty(p)
      ylabel(['Variable x' int2str(p) letter])
    else
      ylabel(['Variable x' letter])
    end
    title(figureTitle)
    set(gca, 'XLim', [xMin-xRange*percentMargin/100 xMax+xRange*percentMargin/100])
    set(gca, 'YLim', [yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
  case 2
    xMin = min([min(groupMV1test(:,1)), min(groupMV2test(:,1)), mean(groupMV1(:,1),1), mean(groupMV2(:,1),1)]);
    xMax = max([max(groupMV1test(:,1)), max(groupMV2test(:,1)), mean(groupMV1(:,1),1), mean(groupMV2(:,1),1)]);
    xRange = xMax - xMin;
    if isequal(xRange,0)
      xRange = 0.01;
    end
    yMin = min([min(groupMV1test(:,2)), min(groupMV2test(:,2)), mean(groupMV1(:,2),1), mean(groupMV2(:,2),1)]);
    yMax = max([max(groupMV1test(:,2)), max(groupMV2test(:,2)), mean(groupMV1(:,2),1), mean(groupMV2(:,2),1)]);
    yRange = yMax - yMin;
    if isequal(yRange,0)
      yRange = 0.01;
    end
    plot(groupMV1test(:,1),groupMV1test(:,2),'bo','MarkerSize',3,'LineWidth',1)  % Testing elements of group 1
    if exist('coefficients','var')
      globalCentroid = mean([groupMV1; groupMV2],1);
      paramLine = sqrt(xRange^2+yRange^2);
      lambda = -paramLine:(min([xRange yRange])/100):paramLine;
      % % Subspace
      coefficients = coefficients/sqrt(sum(coefficients.^2));
      subspace = nncopy(globalCentroid,length(lambda),1) +...
          nncopy(lambda',1,length(coefficients)).*nncopy(coefficients',length(lambda),1);
      subspace(subspace(:,1) < xMin-0.7*xRange*percentMargin/100,:) = [];
      subspace(subspace(:,1) > xMax+0.7*xRange*percentMargin/100,:) = [];
      subspace(subspace(:,2) < yMin-0.7*yRange*percentMargin/100,:) = [];
      subspace(subspace(:,2) > yMax+0.7*yRange*percentMargin/100,:) = [];
      % A slight graphical translation for the two lines not to overlap
      subspace(:,1) = subspace(:,1) - xRange/500;
      subspace(:,2) = subspace(:,2) - yRange/500;
    end
    hold on
    plot(groupMV2test(:,1),groupMV2test(:,2),'rs','MarkerSize',3,'LineWidth',1)  % Testing elements of group 2
    plot(mean(groupMV1test(:,1),1),mean(groupMV1test(:,2),1),'bo','MarkerSize',8,'LineWidth',1)  % Centroid of the testing elements of group 1
    plot(mean(groupMV2test(:,1),1),mean(groupMV2test(:,2),1),'rs','MarkerSize',8,'LineWidth',1)  % Centroid of the testing elements of group 2
    plot(mean(groupMV1(:,1),1),mean(groupMV1(:,2),1),'bo','MarkerSize',8,'LineWidth',3)  % Centroid of the training elements of group 1
    plot(mean(groupMV2(:,1),1),mean(groupMV2(:,2),1),'rs','MarkerSize',8,'LineWidth',3)  % Centroid of the training elements of group 2
    if exist('coefficients','var')  % Projecting directions
      plot(subspace(:,1),subspace(:,2),'Color','m','LineStyle','--','LineWidth',2)
    end
    hold off
    xlabel(['Variable x1' letter])
    ylabel(['Variable x2' letter])
    title(figureTitle)
    set(gca, 'XLim', [xMin-xRange*percentMargin/100 xMax+xRange*percentMargin/100])
    set(gca, 'YLim', [yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
  case 3
    xMin = min([min(groupMV1test(:,1)), min(groupMV2test(:,1)), mean(groupMV1(:,1),1), mean(groupMV2(:,1),1)]);
    xMax = max([max(groupMV1test(:,1)), max(groupMV2test(:,1)), mean(groupMV1(:,1),1), mean(groupMV2(:,1),1)]);
    xRange = xMax - xMin;
    if isequal(xRange,0)
      xRange = 0.01;
    end
    yMin = min([min(groupMV1test(:,2)), min(groupMV2test(:,2)), mean(groupMV1(:,2),1), mean(groupMV2(:,2),1)]);
    yMax = max([max(groupMV1test(:,2)), max(groupMV2test(:,2)), mean(groupMV1(:,2),1), mean(groupMV2(:,2),1)]);
    yRange = yMax - yMin;
    if isequal(yRange,0)
      yRange = 0.01;
    end
    zMin = min([min(groupMV1test(:,3)), min(groupMV2test(:,3)), mean(groupMV1(:,3),1), mean(groupMV2(:,3),1)]);
    zMax = max([max(groupMV1test(:,3)), max(groupMV2test(:,3)), mean(groupMV1(:,3),1), mean(groupMV2(:,3),1)]);
    zRange = zMax - zMin;
    if isequal(zRange,0)
      zRange = 0.01;
    end
    plot3(groupMV1test(:,1),groupMV1test(:,2),groupMV1test(:,3),'bo','MarkerSize',3,'LineWidth',1)  % Testing elements of group 1
    if exist('coefficients','var')
      globalCentroid = mean([groupMV1; groupMV2],1);
      paramLine = sqrt(xRange^2+yRange^2+zRange^2);
      lambda = -paramLine:(min([xRange yRange zRange])/100):paramLine;
      % % Line from the classical discriminant analysis
      coefficients = coefficients/sqrt(sum(coefficients.^2));
      subspace = nncopy(globalCentroid,length(lambda),1) +...
          nncopy(lambda',1,length(coefficients)).*nncopy(coefficients',length(lambda),1);
      subspace(subspace(:,1) < xMin-0.5*xRange*percentMargin/100,:) = [];
      subspace(subspace(:,1) > xMax+0.5*xRange*percentMargin/100,:) = [];
      subspace(subspace(:,2) < yMin-0.5*yRange*percentMargin/100,:) = [];
      subspace(subspace(:,2) > yMax+0.5*yRange*percentMargin/100,:) = [];
      subspace(subspace(:,3) < zMin-0.5*zRange*percentMargin/100,:) = [];
      subspace(subspace(:,3) > zMax+0.5*zRange*percentMargin/100,:) = [];
      % A slight graphical translation for the two lines not to overlap
      subspace(:,1) = subspace(:,1) - xRange/500;
      subspace(:,2) = subspace(:,2) - yRange/500;
      subspace(:,3) = subspace(:,3) - zRange/500;
    end
    hold on
    plot3(groupMV2test(:,1),groupMV2test(:,2),groupMV2test(:,3),'rs','MarkerSize',3,'LineWidth',1)  % Testing elements of group 2
    plot3(mean(groupMV1test(:,1),1),mean(groupMV1test(:,2),1),mean(groupMV1test(:,3),1),'bo','MarkerSize',8,'LineWidth',1)  % Centroid of the testing elements of group 1
    plot3(mean(groupMV2test(:,1),1),mean(groupMV2test(:,2),1),mean(groupMV2test(:,3),1),'rs','MarkerSize',8,'LineWidth',1)  % Centroid of the testing elements of group 2
    plot3(mean(groupMV1(:,1),1),mean(groupMV1(:,2),1),mean(groupMV1(:,3),1),'bo','MarkerSize',8,'LineWidth',3)  % Centroid of the training elements of group 1
    plot3(mean(groupMV2(:,1),1),mean(groupMV2(:,2),1),mean(groupMV2(:,3),1),'rs','MarkerSize',8,'LineWidth',3)  % Centroid of the training elements of group 2
    if exist('coefficients','var')  % Projecting directions
      plot3(subspace(:,1),subspace(:,2),subspace(:,3),'Color','m','LineStyle','--','LineWidth',2)
    end
    hold off
    xlabel(['Variable x1' letter])
    ylabel(['Variable x2' letter])
    zlabel(['Variable x3' letter])
    title(figureTitle)
    set(gca, 'XLim', [xMin-xRange*percentMargin/100 xMax+xRange*percentMargin/100])
    set(gca, 'YLim', [yMin-yRange*percentMargin/100 yMax+yRange*percentMargin/100])
    set(gca, 'ZLim', [zMin-zRange*percentMargin/100 zMax+zRange*percentMargin/100])
    set(gca, 'XGrid', 'on')
    set(gca, 'YGrid', 'on')
end

end
