function[] = mmda(input_train_file,input_test_file,output_file,num_dim);

% This matlab code does the Modified Multiple Discriminant Analysis (MMDA)
% as in theory
%
% open the training file 
% training file should be class_id followed by the features
%
fp = fopen(input_train_file,'r');
prev_class_id = 0;
index = 0;
while(~feof(fp)) 
  a = fscanf(fp, '%f',num_dim + 1);
     
  class_id = a(1);
  if (class_id == prev_class_id)
    class(class_id,index,:) = a(2:num_dim+1)';
  else
    index = 1;
    class(class_id,index,:) = a(2:num_dim+1)';
  end
  index = index + 1;
  prev_class_id = class_id;
end
fclose(fp);

% get the number of classes, number of training samples
% 
get_size = size(class)
num_classes = get_size(1);
num_train_samples = get_size(2);
vec_length = get_size(3);

for i = 1:1:vec_length
  overall_mean(i) = 0.0;
end

% calculate the within-class scatter matrix
%
for i = 1:1:num_classes
  % get the mean for each class
  %
  mean_class(i,:,:) = mean(class(i,:,:));

  % get the within-class covariance matrix for each class
  %
  sw(i,:,:) = cov(shiftdim(class(i,:,:)));

  overall_mean = overall_mean + num_train_samples * mean_class(i,:,:);
end

total_num_data = num_classes * num_train_samples;
overall_mean = overall_mean / total_num_data;

for i = 1:1:vec_length
  for j = 1:1:vec_length
    SW(i,j) = 0.0;
  end
end

% calculate the within-class scatter matrix
%
for i = 1:1:num_classes
  SW(:,:) = SW(:,:) + shiftdim(sw(i,:,:));
end

for i = 1:1:vec_length
  for j = 1:1:vec_length
    sb(i,j) = 0.0;
  end
end

% calculate the total between-class covariance matrix
%
for i = 1:1:num_classes
  val1 = (mean_class(i,:,:) - overall_mean)';
  val2 = val1';
  sb(:,:) = sb(:,:) + (num_train_samples * val1 * val2);
end

% get the transformation matrix which is obtained
% by solving for the eigen values of the between class
% scatter matrix
% W's column are the eigenvectors and correspond to the 
% i-th transformation
%
[W,Y] = eigs(sb,num_classes - 1);

vec_length = num_classes - 1;

% now recalculate the within-class variance of each dimesnion
% to be used later for classification
%
for i = 1:1: vec_length
  r(i) = W(:,i)' * SW(:,:) * W(:,i);
end

% do the classification
%
% compute the transformed features
%
for i = 1:1:num_classes
  for j = 1:1:num_train_samples
      y(i,j,:) = W(:,:)' * shiftdim(class(i,j,:));
  end
end

% compute the mean
%
for i = 1:1:num_classes
   trans_mean(i,:,:) = mean(y(i,:,:));
end

% calculate the distance between the input feature vector and
% any of the classes
%
% open the test file
%
fp = fopen(input_test_file,'r');
fp1 = fopen(output_file,'w');
k = 1;
while(~feof(fp))
  test_vec = fscanf(fp, '%f', num_dim);
  for i = 1:1:num_classes
    D(i) = 0;
    for j = 1:1:vec_length
      D(i) = D(i) + abs(((W(:,j)' * test_vec) - trans_mean(i,j))/sqrt(r(j)));
    end
  end
  [val(k),ind(k)] = min(D);
  fprintf(fp1,'%d\n',ind(k));
  k = k + 1;
end
fclose(fp);
fclose(fp1);
