%Author: Ashish Kapoor MIT Media Lab, 10/20/02
%Please send corrections to: ash@media.mit.edu

clear; clc;

num_states = 7;
num_obs_symbols = 4;

%*******************************************************************************************
%Load Training Data
fprintf(1, 'Loading Training and Testing Data ....\n\n');

Train_class_1=cell(0);
for i=1:30
    fname=sprintf('hmmdata\\train\\class1\\sequence%2.2d',i);
    obs=load(fname);
    Train_class_1{i}=obs;
end
Train_class_2=cell(0);
for i=1:20
    fname=sprintf('hmmdata\\train\\class2\\sequence%2.2d',i);
    obs=load(fname);
    Train_class_2{i}=obs;
end

%Load Testing Data
Test_class_1=cell(0);
for i=1:10
    fname=sprintf('hmmdata\\test\\class1\\sequence%2.2d',i);
    obs=load(fname);
    Test_class_1{i}=obs;
end
Test_class_2=cell(0);
for i=1:10
    fname=sprintf('hmmdata\\test\\class2\\sequence%2.2d',i);
    obs=load(fname);
    Test_class_2{i}=obs;
end
%*********************************************************************************************

%Initializing Parameters for HMMs
fprintf(1, 'Initializing HMMs with %d states .... \n\n', num_states);

pi_class_1=zeros(num_states,1); pi_class_1(1)=1;
pi_class_2=zeros(num_states,1); pi_class_2(1)=1;

a_class_1=zeros(num_states, num_states);
a_class_2=zeros(num_states, num_states);
%Randomizing the initialization of transition matrix while
%keeping the HMM strictly left-to-right
for i=1:num_states-1
    a_class_1(i,i)=rand; a_class_1(i,i+1)=1-a_class_1(i,i);
    a_class_2(i,i)=rand; a_class_2(i,i+1)=1-a_class_2(i,i);
end
a_class_1(num_states, num_states) = 1;
a_class_2(num_states, num_states) = 1;

%Initializing observation matrices with equal probability of observing each symbol 
b_class_1=rand(num_states, num_obs_symbols);
b_class_2=rand(num_states, num_obs_symbols);
for i=1:size(b_class_1,1)
    b_class_1(i,:)=b_class_1(i,:)/sum(b_class_1(i,:));
end
for i=1:size(b_class_2,1)
    b_class_2(i,:)=b_class_2(i,:)/sum(b_class_2(i,:));
end

%**********************************************************************************************

%Learning the parameters using Baulm-Welch
fprintf(1, 'Learning Class 1 HMM parameters .... \n\n');
[pi_class_1, a_class_, b_class_1, llk_array_1] = learn_hmm_param(pi_class_1, a_class_1, b_class_1, num_states, Train_class_1, 0.1);

fprintf(1, '\nLearning Class 2 HMM parameters .... \n\n');
[pi_class_2, a_class_2, b_class_2, llk_array_2] = learn_hmm_param(pi_class_2, a_class_2, b_class_2, num_states, Train_class_2, 0.1);

fig=figure;
set(fig,'doublebuffer', 'on');

subplot(2, 1, 1); plot(llk_array_1); xlabel('Number of Iterations'); ylabel('Loglikelihood'); title('Class 1');
subplot(2, 1, 2); plot(llk_array_2); xlabel('Number of Iterations'); ylabel('Loglikelihood'); title('Class 2');

%**********************************************************************************************

%Testing on the dataset
fprintf(1, 'Testing ....\n');

correct_class_1=0;
for i=1:size(Test_class_1,2)
    p_hmm_1=log_prob_obs_hmm(pi_class_1, a_class_1, b_class_1, num_states, Test_class_1{i});
    p_hmm_2=log_prob_obs_hmm(pi_class_2, a_class_2, b_class_2, num_states, Test_class_1{i});
    if p_hmm_1 >= p_hmm_2
        correct_class_1=correct_class_1+1;
    end
end

correct_class_2=0;
for i=1:size(Test_class_2,2)
    p_hmm_1=log_prob_obs_hmm(pi_class_1, a_class_1, b_class_1, num_states, Test_class_2{i});
    p_hmm_2=log_prob_obs_hmm(pi_class_2, a_class_2, b_class_2, num_states, Test_class_2{i});
    if p_hmm_2 >= p_hmm_1
        correct_class_2=correct_class_2+1;
    end
end

fprintf(1,'\n                    Confusion Matrix for the Test Data\n');
fprintf(1,'              Recognized as Class 1          Recognized as Class 2\n');
fprintf(1,'Class 1                %d                              %d\n', correct_class_1, size(Test_class_1,2)-correct_class_1);
fprintf(1,'Class 2                %d                              %d\n', size(Test_class_2,2)-correct_class_2, correct_class_2);