function  [ w_fda threshold m1_lda m2_lda  Y_latent_train ]=pr_MLDA(featuredata,trainlabel)

% *************************************************************************
% This function is used to  calculate the maximum uncertainty linear
% discriminant analysis (MLDA) classifier based on training
% data which have passed the feature selection.This method employed a
% maximum entropy covariance selection method instead of the within-class scatter matrix
% in fisher linear discriminant analysis.The MLDA is also a linear
% discriminant method. y(x)=w.x+b.It maps the selected feature space into a
% one-dimensional feature space.

% Input:

%          featuredata:  the  training data which have passed the feature
%                        selection.A matrix n*s, n is the number of training
%                        subjects and s is the number of the selected features.

%           trainlabel:  the training label corresponding the training subject.
%                        A vector n*1, n is the number of training subjects.
%                        Note that the program requires the label of patient
%                        is -1 and the label of control is 1 to help us calculate
%                        the classification accuracy, sensitivity and
%                        specificity  uniformly.

% Output:

%              w_fda:  the weight vector of the selected features. A vector
%                      s*1, s is the number of the selected features.

%          threshold:  the separation of two classes data based on training
%                      data in one-dimensional feature space.

%             m1_lda:  the mean value of the patient data in one-dimensional
%                      feature space

%             m2_lda:  the mean value of the control data in one-dimensional
%                      feature space.

%     Y_latent_train:  the value of the trainging data one-dimensional
%                      feature space. A vector n*1, n is the number of
%                      training subject.

% References:

%Thomaz, C.E., Boardman, J.P., Hill, D.L.G., Hajnal, J.V., Edwards, A.D.,
%Rutherford, M.A., Gillies, D.F., Rueckert, D., 2004. Whole Brain
%Voxel-based Analysis Using Registration and Multivariate Statistics. 

%Thomaz, C., Boardman, J., Counsell, S., Hill, D., Hajnal, J., Edwards, A.,
%Rutherford, M., Gillies, D., Rueckert, D., 2007. A multivariate
%statistical analysis of the developing human brain in preterm infants. Image and Vision Computing 25, 981-994. 

% Written by Zhengjia Dai, State Key Laboratory of Cognitive Neuroscience
% and Learning, Beijing Normal University, Beijing, P.R.China, 100875. 
% Contact information:
% Zhengjia Dai: daizhengjia13@gmail.com
% Yong He : yong.h.he@gmail.com; yong.he@bnu.edu.cn
% September 24th,2011
% *************************************************************************

train_num=size(featuredata,1);

% finding the patient's data in training data
class1=featuredata(find(trainlabel==-1),:);

% finding the control's data in training data
class2=featuredata(find(trainlabel==1),:);

m1=mean(class1);

m2=mean(class2);

s1=zeros(size(class1,2));

s2=zeros(size(class2,2));

%calculating within-class scatter matrix of patient group
n1=size(class1,1);
for j=1:n1
    s1=s1+(class1(j,:)-m1)'*(class1(j,:)-m1);
end

%calculating within-class scatter matrix of control group
n2=size(class2,1);
for j=1:n2
    s2=s2+(class2(j,:)-m2)'*(class2(j,:)-m2);
end

% calculating the within-class scatter matrix
sw=s1+s2;

% maximum entropy covariance selection
sp=sw/(train_num-2);
m_lamda=trace(sp)/(train_num-1);
[TZXL TZZ]=eig(sp);
aaa=diag(TZZ);
aaa(find(aaa<=m_lamda))=m_lamda;
TZZ=diag(aaa);
sww=(TZXL*TZZ*TZXL')*(train_num-2);

% calculating the weight vector of feature
w_fda=inv(sww)*(m1-m2)';
w_fda=w_fda/norm(w_fda);


Y_latent_train=featuredata*w_fda;

m1_lda=mean(Y_latent_train(find(trainlabel==-1)));

m2_lda=mean(Y_latent_train(find(trainlabel==1)));

threshold=(n1*m1_lda+n2*m2_lda)/(n1+n2);
