% This is the matlab code for classifying data using Linear Discriminant
% Analysis. It creates a two transformation matrices, one specific to
% each of the classes
% 
% file : /ftp/pub/resources/courses/ece_8993_speech/homework/1998/utilities/lda_class_spec.m
%

% the procedure used in the function parallels the description of LDA in 
% K. Fukunaga, "Introduction to Statistical Pattern Recognition", Second Ed.,
% Academic Press, 1990.
%
% Fukunaga however does not describe class specific LDA

% load the set-1 data created using xmgr
%
load case1_set_1.dat;
s1_new = case1_set_1;

% get the size of set-2
%
[m,n] = size(s1_new);

% load the set-2 data created using xmgr
%
load case1_set_2.dat;
s2_new = case1_set_2;

% load the merged data
%
load case1_both.dat;
set_all = case1_both;

% get the size of set-2
%
[x,y] = size(s2_new);

% compute the mean of each class
%
set1_mean = mean(s1_new);
set2_mean = mean(s2_new);

set_mean(1,:) = set1_mean;
set_mean(2,:) = set2_mean;

% define the test set
%
x1=[-1,-1];
x2=[0,0];
x3=[0.5, 0.5];
x4=[0.5, -0.5];

% compute the mean of entire data set (eqn. 10.3, pp. 446)
% Note that we assume equiprobable class in this case and hence the factor
% 0.5
%
set_all_mean = 0.5*set1_mean + 0.5*set2_mean;

% compute the between class scatter (eqn. 10.2 pp. 446)
%
between_class_scatter = cov(set_mean);

% compute the covariance for each class
%
set1_cov = cov(s1_new);
set2_cov = cov(s2_new);

% compute the covariance of all the data
%
set_all_cov = cov(set_all);

% compute the within class scatter (eqn. 10.1, pp. 446)
%
within_class_scatter = 0.5*(set1_cov) + 0.5*(set2_cov);

% compute the optimizing criterion, in this case the ratio of between class 
% scatter to the within class scatter (pp. 451).
% note that for class specific LDA, within class scatter is nothing but 
% the covariance of the class
%
new_feat_f1_spec = inv(set1_cov)*between_class_scatter;
new_feat_f2_spec = inv(set2_cov)*between_class_scatter;

% compute the eigen vectors and eigen values
%
[eigen_vect_1,eigen_val_1] = eig(new_feat_f1_spec);
[eigen_vect_2,eigen_val_2] = eig(new_feat_f2_spec);

% this step is to take the eigen vectors corresponding to non zero eigen 
% values and throwing the rest. A good check at this point is to see
% if there are L-1 non-zero eigen values for an L-class problem which should 
% be true if we are doing every thing right.
%

% make a temporary matrix which will hold all the eigen values
% in one column and will be used for comparing and getting the
% corresponding eigen vectors
%
for k=1:n
	temp_eigen_val_1(k) = eigen_val_1(k,k);
	temp_eigen_val_2(k) = eigen_val_2(k,k);
end;

% create the transform (comprising of all eigen vectors with a corresponding
% non zero eigen value
%
i = 1;
j = 1;
for k=1:n	
  if temp_eigen_val_1(k) ~= 0 
	transform_1_spec(:,i) = eigen_vect_1(:,k);
	i = i + 1;
 end;

 if temp_eigen_val_2(k) ~= 0
 	transform_2_spec(:,j) = eigen_vect_2(:,k);
	j = j + 1;
 end;
end;

% compute the transformed data using the transform specific to the class
%
new_transformed_set_1 = transform_1_spec(:,1)'*(s1_new)';
new_transformed_set_2 = transform_2_spec(:,1)'*(s2_new)';
mean_trans_set_1 = mean(new_transformed_set_1);
mean_trans_set_2 = mean(new_transformed_set_2);

% plot these new data points in the original space.
% here we project the transformed data points on to the axis defined
% by the eigen vector corresponding to the non-zero eigen value
%

% since we have only one non-zero eigen value for a 2 class problem, the 
% transform will be a 2 X 1 matrix which also represents a vector (or 
% equivalently a straight line. the transformed data points lie on this 
% straight line
%
slope_1 = transform_1_spec(2,1)/transform_1_spec(1,1);
slope_2 = transform_2_spec(2,1)/transform_2_spec(1,1);

for k=1:100

	new_plot_set_1(k,1) = -sqrt((new_transformed_set_1(:,k)*new_transformed_set_1(:,k))/(1+(slope_1*slope_1)));
	new_plot_set_1(k,2) = slope_1*new_plot_set_1(k,1);
	new_plot_set_2(k,1) = sqrt((new_transformed_set_2(:,k)*new_transformed_set_2(:,k))/(1+(slope_2*slope_2)));
	new_plot_set_2(k,2) = slope_2*new_plot_set_2(k,1);
end;

plot(s1_new(:,1), s1_new(:,2), 'ko');
hold on;
plot(s2_new(:,1), s2_new(:,2), 'k*');

%plot(new_plot_set_1(:,1), new_plot_set_1(:,2), 'ko');
%plot(new_plot_set_2(:,1), new_plot_set_2(:,2), 'k*');


% find the decision region 
%
% initialize the size of the sample data
%
index1 = 1;
index2 = 1;

% plot between two somewhat arbitrary x values that fit the data chosen for
% this experiment
%
for x = 6:-0.1:-6

        % plot between two somewhat arbitrary y values that fit the data chosen
        % for this experiment
        %
        for y = 6:-0.1:-6

                dist1=dist(mean_trans_set_1,(transform_1_spec(:,1)' * [x,y]'));
                dist2=dist(mean_trans_set_2,(transform_2_spec(:,1)' * [x,y]'));

                % see if the sign of the distance difference has changed. if so
                % then add the point to the plot data and increment the index
                % value
                %
                if ((dist1 < dist2)) 
                        plot_values_lda_spec1(index1,:) = [x,y];
                        index1 = index1 + 1;
                end;
                if ((dist2 < dist1)) 
			 plot_values_lda_spec2(index2,:) = [x,y];
                        index2 = index2 + 1;
                end;
        end
end

% plot the data points found on the current figure
%
plot(plot_values_lda_spec1(:,1), plot_values_lda_spec1(:,2), 'g.');
plot(plot_values_lda_spec2(:,1), plot_values_lda_spec2(:,2), 'b.');
hold off;

% transform the test vectors and compute the euclidean distance
% from the mean of each of the classes
%
dist_1 = [(transform_1_spec(:,1)'*x1') - (mean_trans_set_1)];
dist_2 = [(transform_2_spec(:,1)'*x1') - (mean_trans_set_2)];

dist_3 = [(transform_1_spec(:,1)'*x2') - (mean_trans_set_1)];
dist_4 = [(transform_2_spec(:,1)'*x2') - (mean_trans_set_2)];

dist_5 = [(transform_1_spec(:,1)'*x3') - (mean_trans_set_1)];
dist_6 = [(transform_2_spec(:,1)'*x3') - (mean_trans_set_2)];

dist_7 = [(transform_1_spec(:,1)'*x4') - (mean_trans_set_1)];
dist_8 = [(transform_2_spec(:,1)'*x4') - (mean_trans_set_2)];

% calculate the magnitude of euclidean distance
%
% for test vector 1
%
orig_dist_1 = dist_1 * transpose(dist_1);
mag_dist_1 = sqrt(orig_dist_1);

% for test vector 2
%
orig_dist_2 = dist_2 * transpose(dist_2);
mag_dist_2 = sqrt(orig_dist_2);

% for test vector 3
%
orig_dist_3 = dist_3 * transpose(dist_3);
mag_dist_3 = sqrt(orig_dist_3);

% for test vector 4
%
orig_dist_4 = dist_4 * transpose(dist_4);
mag_dist_4 = sqrt(orig_dist_4);

% for test vector 5
%
orig_dist_5 = dist_5 * transpose(dist_5);
mag_dist_5 = sqrt(orig_dist_5);

% for test vector 6
%
orig_dist_6 = dist_6 * transpose(dist_6);
mag_dist_6 = sqrt(orig_dist_6);

% for test vector 7
%
orig_dist_7 = dist_7 * transpose(dist_7);
mag_dist_7 = sqrt(orig_dist_7);

% for test vector 8
%
orig_dist_8 = dist_8 * transpose(dist_8);
mag_dist_8 = sqrt(orig_dist_8);

% open a output file for printing
% the results
%

% compare the magnitude of the distance
% and classify the set to which the test vector belong
%
% classify for the test vector 1
%
if ((mag_dist_1) < (mag_dist_2)) 
	fprintf('The test set x1 belongs to set 1\n');
else
   	fprintf('The test set x1 belongs to set 2\n');
end;

% classify for the test vector 2
%
if ((mag_dist_3) < (mag_dist_4)) 
   	fprintf( 'The test set x2 belongs to set 1\n');
else
   	fprintf('The test set x2 belongs to set 2\n');
end;

% classify for the test vector 3
%
if ((mag_dist_5) < (mag_dist_6))
   	fprintf('The test set x3 belongs to set 1\n');
else
   	fprintf('The test set x3 belongs to set 2\n');
end;

% classify for the test vector 4
%
if ((mag_dist_7) < (mag_dist_8)) 
   	fprintf('The test set x4 belongs to set 1\n');
else
   	fprintf('The test set x4 belongs to set 2\n');
end;

% end of file
%


