% Initial number of states for training HMM
% You can set it to what you think is suatable

N=  ;


% Number of observation symbols
%This depends on how you model the inputs from the Scribble

M=  ;


% Number of symbols to recognize
% In general you can set it to 1 if you just 
% wanted to test individual trials 

N_symbols=   ;



%----------------------------------------
%----------------------------------------
%----------------------------------------
%When classifying a new character do the following:
%----------------------------------------
%----------------------------------------
%----------------------------------------

% read a new letter 
% represent it in the way you decided to
% classify it using  character = HMM_classify(A, B, p0, Z);


%function [character]=HMM_classify(A,B,p0,Z)

% A the transition matrix, (NxN symbols).
% B the observation symbol probability matrix, (NxM symbols).
% p0 the initial state distribution vector, (N symbols).
% Z the observerd sequence, vector, (1xT).

for i=1:N

   [alfa prob] = HMM_forward(A, B, p0, Z);
   %SEE NEXT FUNCTION 

      if (prob > max)
         max = prob;
         if (i == 1) character = '1';
         elseif (i == 2) character = '2';
         elseif (i == 3) character = '3';
         elseif (i == 4) character = '4';

% End of function


%----------------------------------------
%----------------------------------------
%----------------------------------------
%----------------------------------------
%    YOU NEED TO IMPLEMENT "FORWARD" ALGORITHM
%    YOU ALREADY KNOW A,B,p0 
%----------------------------------------
%----------------------------------------
%----------------------------------------
%----------------------------------------

function [alfa, prob]=hmm_forward(A,B,p0,z)

% A the transition matrix, (NxN).
% B the observation symbol probability matrix, (NxM).
% p0 the initial state distribution vector, (Nx1).
% z the observerd sequence, vector, (1xT).
% prob vector, (1xT).
% alfa, the scaled alfa variable, matrix, (NxT).

alfa_1 = p0 * B;

for t=2:length(z)
   alfa_t = B * sum(alfa_(t-1) * A)
   end
   
% End of function


%----------------------------------------
%----------------------------------------
%----------------------------------------
%----------------------------------------
%    YOU NEED TO IMPLEMENT "BACKWARD" ALGORITHM
%    YOU ALREADY KNOW A,B,p0 
%----------------------------------------
%----------------------------------------
%----------------------------------------
%----------------------------------------


function [beta]=HMM_backward(A,B,p0,z)

% A the transition matrix, (NxN).
% B the observation symbol probability matrix, (NxM).
% p0 the initial state distribution vector, (Nx1).
% z the observerd sequence, vector, (1xT).
% beta, the scaled beta variable, matrix, (NxT).

beta_T=ones(length(p0),length(z));

for t=T-1:1
   beta_t = sum (beta_(t+1) * A * B)
end



%----------------------------------------
%    AFTER YOU DECIDED HOW TO MODEL
%    THE SYSTEM and HAVE 
%    FORWARD AND BACKWARD ALGORITHMS
%    YOU HAVE TO LEARN/TRAIN HMMs
%----------------------------------------



% Set up initial HMM matrices values
% OBSERVE THAT THIS IS JUST AN EXAMPLE 
% AND THAT A, B, p0 SHOULD BE SET A BIT SMARTER!!!!!!

   A = zeros(N);
   p0 = zeros(N,1);
   B = ones(N,M);

      
% Train it TTT times

for i=1:TTT
   % SEE NEXT FUNCTION 
   [A, B, p0] = HMM_training(A,B,p0,Z,lZ);
end


% STORE FINAL A, B, p0 for each letter 



   
%----------------------------------------
%----------------------------------------
%----------------------------------------
%   NEED TO IMPLEMENT  [A, B, p0] = HMM_training(A,B,p0,Z,lZ);
%----------------------------------------
%----------------------------------------
%----------------------------------------


function [A_new, B_new, p0_new]=HMM_training(A,B,p0,z,lz)

% A the transition matrix, (NxN).
% B the observation symbol probability matrix, (NxM).
% p0 the initial state distribution vector, (Nx1).

% z the concatenated sequence of observation sequences, 
% E.g z=[z1 z2], (1xlength(z))
% lz the length of each observation sequence, E.g lz=[length(z1) length(z2)]
% YOU CAN DO THE ABOVE IN A DIFFERENT WAY!


  gamma = alfa* beta /(sum(alfa * beta))
  
  eta = alfa * A * B * beta / (sum(alfa * beta))


  p0 = gamma

  A_new = sum(eta) / sum (sum (eta))
  B_new (v_k)= sum_(v_k)(gamma)  / sum (gamma) 








