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

% PURPOSE: This function applies different classification methods for multivariate data.

% INPUT VARIABLES:
%    numberMVtoCLASSmethod: Number of the method to be applied.
%    infoMVtoCLASSmethod: Necessary information for the classification method.
%    groupMV1: Matrix with the training multivariate data of group 1 as rows.
%    groupMV1test: Matrix with the testing multivariate data of group 1 as rows.
%    groupMV2: Matrix with the training multivariate data of group 2 as rows.
%    groupMV2test: Matrix with the testing multivariate data of group 2 as rows.

% OUTPUT VARIABLES:
%    errorRatesGroup1: Number of misclassified testing elements of group 1 (out of m1).
%    errorRatesGroup2: Number of misclassified testing elements of group 2 (out of m2).
%    outputMVtoCLASSmethod: Cell array with possible additional output information of the method.

% COMMENTS:
%   (1) The KNN method, based on the function 'knnclassify()', needs the Bioinformatics Toolbox of MATLAB.
%   (2) This function is based on code of Andrés-M. Alonso.

function [errorRatesGroup1 errorRatesGroup2 outputMVtoCLASSmethod] =...
    fromMVtoCLASS(numberMVtoCLASSmethod, infoMVtoCLASSmethod, groupMV1, groupMV1test, groupMV2, groupMV2test)

n1 = size(groupMV1,1);
n2 = size(groupMV2,1);
m1 = size(groupMV1test,1);
m2 = size(groupMV2test,1);

switch numberMVtoCLASSmethod
  case 0  % Identity -- trick to call any classification method
    errorRatesGroup1 = groupMV1test;
    errorRatesGroup2 = groupMV2test;
    outputMVtoCLASSmethod = infoMVtoCLASSmethod;
  case 1  % Method WI: classification based on Fisher's method
    varsIncluded = infoMVtoCLASSmethod{1};
    stdvarsIncluded = infoMVtoCLASSmethod{2};
    numVars = size(groupMV1,2);
    [outclass, err, w] = classifyMod([groupMV1test(:,varsIncluded); groupMV2test(:,varsIncluded)],...
        [groupMV1(:,varsIncluded); groupMV2(:,varsIncluded)], [zeros(n1,1); ones(n2,1)]);
    % Error rates
    errorRatesGroup1 = sum(outclass(1:m1)); % Number of misclassified testing elements of group 1 (out of m1).
    errorRatesGroup2 = sum(1-outclass((m1+1):(m1+m2))); % % Number of misclassified testing elements of group 2 (out of m2).
    % Coefficients
    %   To obtain a vector with the first component being positive. The sign of
    %   all the components can be changed, and '-w' determines the same
    %   one-dimensional subspace. We change the sign to facilitate the
    %   interpretation of the results
    if isequal(sign(w(1)),-1)
      w = -w;
    end
    w1 = zeros(numVars,1);
    w1(varsIncluded) = w.*stdvarsIncluded;
    w2 = zeros(numVars,1);
    w2(varsIncluded) = w;
    outputMVtoCLASSmethod{1} = w1;  % Standardized version of the coefficients
    outputMVtoCLASSmethod{2} = w2;  % Nonstandardized version of the coefficients
    outputMVtoCLASSmethod{3} = outclass(1:m1);  % Misclassified elements of testing sample 1
    outputMVtoCLASSmethod{4} = 1-outclass((m1+1):(m1+m2));  % Misclassified elements of testing sample 2
  case 2  % Method WD: classification based on weighted distance
    varsIncluded = infoMVtoCLASSmethod{1};
    stdvarsIncluded = infoMVtoCLASSmethod{2};
    Sw = infoMVtoCLASSmethod{3};
    numVars = size(groupMV1,2);
    Dmeans = mean(groupMV1(:,varsIncluded))' - mean(groupMV2(:,varsIncluded))';
    w0 = ones(length(varsIncluded),1);
    w0 = w0/(w0'*Sw*w0);
    A = [];
    bb = [];
    Aeq = [];
    beq = [];
    % Our constraint
    lb = zeros(length(varsIncluded),1);
    ub = [];
    % The optimization problem
    [w,fval,exitflag1] = fmincon(@(w) objectiveFunction(w,Dmeans),...
        w0, A, bb, Aeq, beq, lb, ub,@(w) constraintsFunction(w,Sw),...
        optimset('Algorithm','active-set','Display','off')); % Notice that the coefficients of the excluded variables are 0
    % Error rates
    C = w'*(mean(groupMV1(:,varsIncluded))' + mean(groupMV2(:,varsIncluded))')/2;
    outclass = (w'*[groupMV1test(:,varsIncluded); groupMV2test(:,varsIncluded)]' > C);
    equalDist = (w'*[groupMV1test(:,varsIncluded); groupMV2test(:,varsIncluded)]' == C);
    if find(equalDist)
      % The elements such that w'*groupsMVtest(:,varsIncluded)' == C are
      % modified so as to assign them at random...
      outclass(equalDist) = randi(2,length(find(equalDist)),1)-1;
    end
    errorRatesGroup1 = sum(outclass(1:m1)); % Number of misclassified testing elements of group 1 (out of m1).
    errorRatesGroup2 = sum(1-outclass((m1+1):(m1+m2))); % Number of misclassified testing elements of group 2 (out of m2).
    % Coefficients
    w1 = zeros(numVars,1);
    w1(varsIncluded) = w.*stdvarsIncluded;
    w2 = zeros(numVars,1);
    w2(varsIncluded) = w;
    outputMVtoCLASSmethod{1} = w1;  % Standardized version of the coefficients
    outputMVtoCLASSmethod{2} = w2;  % Nonstandardized version of the coefficients
    outputMVtoCLASSmethod{3} = outclass(1:m1);  % Misclassified elements of testing sample 1
    outputMVtoCLASSmethod{4} = 1-outclass((m1+1):(m1+m2));  % Misclassified elements of testing sample 2
  case 3  % Method: KNN (multivariate K-nearest neighbours)
    numberNeighbours = infoMVtoCLASSmethod{1};
    outclass = knnclassify([groupMV1test; groupMV2test], [groupMV1; groupMV2], [zeros(n1,1); ones(n2,1)],numberNeighbours);
    errorRatesGroup1 = sum(outclass(1:m1)); % Number of misclassified testing elements of group 1 (out of m1).
    errorRatesGroup2 = sum(1-outclass((m1+1):(m1+m2))); % Number of misclassified testing elements of group 2 (out of m2).
    outputMVtoCLASSmethod{1} = outclass(1:m1);  % Misclassified elements of testing sample 1
    outputMVtoCLASSmethod{2} = 1-outclass((m1+1):(m1+m2));  % Misclassified elements of testing sample 2
end

end


% PURPOSE: This function implements the objective function of the Fisher's Lineal
%          Discriminant Analysis (LDA).
% INPUT VARIABLES:
%    w: Vector of variables (the parameters that determines the direction in the LDA).
%    Dmeans: Vector with the difference between the means of the two groups.
% OUTPUT VARIABLES:
%    R: Evaluation of the objective function.
function R = objectiveFunction(w, Dmeans)
R = -(w'*Dmeans)^2;
end

% PURPOSE: This function implements the usual constraint of the Fisher's Lineal Discriminant
%          Analysis (LDA) to determine a unique eigenvalue (from the infinite number of them).
% INPUT VARIABLES:
%    w: Vector of variables (the parameters that determines the direction in the LDA).
%    Sw: Matrix (within-group scatter matrix of variability in the LDA).
% OUTPUT VARIABLES:
%    ceq: (Empty variable.)
%    ceq: Evaluation of the constraining function that is equalized to zero.
function [c,ceq] = constraintsFunction(w, Sw)
c = [];
ceq = w'*Sw*w - 1;
end
