%  
%    	This program will perform LDA (Linear Discriminant Analysis) on the set of data.
%       complete matrix of data 9x3 matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function out=lda(class1, class2, data, test_data, output_file)

%  This set of code is dated as correct on November 24, 1997 by Bob Brown

%       Number of features in each vector
%       size will return the number of rows and collumns for the argument
%       I am using x as a dumby variable to hold the number of rows
%       feature_num is the number of columns which is the number of features

%       load data files into matlab


%       Order to reduce to is defined by q

%       Setup files for storing the results.

FID = fopen('results','w');
OUTPUT_TO_FILE=1;
if isequal(nargin,4)
   OUTPUT_TO_FILE=0;
end

[x,feature_num] = size(data);    


%       Number of features vectors

[num_test_vectors,x] = size(test_data);
[num_class1_vect,x] = size(class1);
[num_class2_vect,x] = size(class2);

%       compute mean of total data set and each class

data_mean=mean(data);
class1_mean=mean(class1);
class2_mean=mean(class2);

%       compute within class covariance for each class

class1_cov=cov(class1);
class2_cov=cov(class2);

%       compute between class covariance

data_cov=cov(data);

%       compute new feature vector that is inverse of class_cov multiplied by between class covariance

class1_features=inv(class1_cov)*data_cov;
class2_features=inv(class2_cov)*data_cov;

%       find eigen-vectors and eigen-values for each class

[eigen_vect_class1, eigen_val_class1]=eig(class1_features);
[eigen_vect_class2, eigen_val_class2]=eig(class2_features);

%       decide which eigen vectors are kept and which ones are thrown away.

%       In this next set of code I make a temporary 1xfeature_num matrix that will hold
%       the diagonal of the eigen value matrix (the eigen values)
%       Next, I loop from 1 to feature_num to find the maximum eigen value for each class
%       which corresponds to the maximum eigen vector.  This eigen vector will then be move to the
%       top of the ordered eigen vector matrix until it is completely ordered.

        %make temporary eigen value matrix for each class

for k=1:feature_num                                             %make temporary 1xfeature_num matrix
    temp_eigen_val_class1(k) = eigen_val_class1(k,k);           %to hold the eigen values so that we
    temp_eigen_val_class2(k) = eigen_val_class2(k,k);           %can order them
end;

        %order eigen vector matrix

for k=1:feature_num                                             
    [row1,col1] = find(temp_eigen_val_class1 == max(temp_eigen_val_class1));
    ordered_eigen_vect_class1(k,:) = eigen_vect_class1(col1,:);
    temp_eigen_val_class1(col1) = -10000;

    [row2,col2] = find(temp_eigen_val_class2 == max(temp_eigen_val_class2));
    ordered_eigen_vect_class2(k,:) = eigen_vect_class2(col2,:);
    temp_eigen_val_class2(col2) = -10000;
end;

%       keep the q largest valued rows for eigen vector matrix and store them to the reduced matrix
%       This is the transform for each class.

for n=1:2
    transform_class1(n,:) = ordered_eigen_vect_class1(n,:);             %Class1 transform
    transform_class2(n,:) = ordered_eigen_vect_class2(n,:);             %Class2 transform       
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%       compute reduced features for each class of training data
%       For loop to find reduced features for each vector 

%       The following loop computes the mean reduced class1 features by loop for the number of vectors in the 
%       class and multiplying the class1 transform by the transpose of the of each vector in class1
%       then summing and dividing by the number of vectors in class1
%       Dumby variable are used to hold vector from class1 (may not be nescessary) and the temporary mean
%       as we and compute and add up all the reduced class features.
%       Finally we reinitialize the temporary vector of class1 (this may no be nescessary)

%       CLASS1

temp_mean1=0;                   %initial variable to make matlab happy

for n=1:num_class1_vect
    temp_class1 = class1(n,:);
    reduced_class1_features = transform_class1*transpose(temp_class1);
    reduced_class1_features = transpose(reduced_class1_features);
    temp_mean1 = temp_mean1 + reduced_class1_features;
    temp_class1 = 0;
end;

%       Now we compute the mean of the reduced class1 features by dividing the temp_mean by the number of 
%       class1 vectors

mean_reduced_class1_features = temp_mean1/num_class1_vect;

%       CLASS2

temp_mean2=0;                   %initial variable to make matlab happy

for n=1:num_class2_vect
    temp_class2 = class2(n,:);
    reduced_class2_features = transform_class2*transpose(temp_class2);
    reduced_class2_features = transpose(reduced_class2_features);
    temp_mean2 = temp_mean2 + reduced_class2_features;
    temp_class2 = 0;
end;

%       Now we compute the mean of the reduced class2 features by dividing the temp_mean by the number of 
%       class2 vectors

mean_reduced_class2_features = temp_mean2/num_class2_vect;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%       TESTING

%       test_feature is the feature vector to be tested to determine if it is HSBE, MSBE, or LSBE.



for j=1:num_test_vectors        %begin for loop for number of test feature vectors

temp_test_feature_vect = test_data(j,:);

%       CLASS 1

reduced_test_feature_vect_class1 = transform_class1*transpose(temp_test_feature_vect);
reduced_test_feature_vect_class1 = transpose(reduced_test_feature_vect_class1);

class1_distance = reduced_test_feature_vect_class1 - mean_reduced_class1_features;



%       CLASS 2

reduced_test_feature_vect_class2 = transform_class2*transpose(temp_test_feature_vect);
reduced_test_feature_vect_class2 = transpose(reduced_test_feature_vect_class2);

class2_distance = reduced_test_feature_vect_class2 - mean_reduced_class2_features;


%       find magnitude of distance for each class

mag_class1_distance =0;
mag_class2_distance =0;

for n=1:2
    squared1 = class1_distance(n) * class1_distance(n);
    mag_class1_distance = mag_class1_distance + squared1;            %finds magnitude of distance

    squared2 = class2_distance(n) * class2_distance(n);              
    mag_class2_distance = mag_class2_distance + squared2;            %finds magnitude of distance

end;



%       determine which class is the minimum distance and which class the feature vector belongs to 

distance_vect=[mag_class1_distance, mag_class2_distance];
distance_vect(1)=sqrt(distance_vect(1));
distance_vect(2)=sqrt(distance_vect(2));

if isequal(OUTPUT_TO_FILE,1)
if min(distance_vect)==distance_vect(1)
   fprintf(output_file,'\nx%d belongs to set 1\n',j);
else
   fprintf(output_file,'\nx%d belongs to set 2\n',j);
end
fprintf(output_file,'Euclidean distances:\n');
fprintf(output_file,'x%d & set1 mean: %2.4f\n',j, distance_vect(1));
fprintf(output_file,'x%d & set2 mean: %2.4f\n\n',j, distance_vect(2));
end
out = distance_vect;

end;    %end for loop for the number of feature vectors

