% function [lik]=hmmcost(X, alphabet, K, E, P, Pi)
% 
% simple Hidden Markov Model - variable lengths and discrete symbols
%
% X - cell array of strings
% alphabet - string of characters
% K - number of states (default 2)
% E - observation emission probabilities
% P - state transition probabilities
% Pi - initial state prior probabilities

function [lik]=hmmcost(X, alphabet, K, E, P, Pi)

LA = length(alphabet); 

epsi=1e-6;

% number of sequences
N=length(X(1,:));

% length of each sequence
T=ones(1,N);
for n=1:N,
  T(n)=length(X{n});
end;

TMAX = max(T);

fprintf('\n********************************************************************\n');
fprintf('Testing %i sequences of maximum length %i from an alphabet of size %i\n',N,TMAX,LA);
fprintf('HMM with %i hidden states\n',K);
fprintf('********************************************************************\n');

B=zeros(TMAX,K);

LL=[];
lik=0;

  %%%% FORWARD-BACKWARD
  

  Gammainit=zeros(1,K);
  Gammasum=zeros(1,K);
  Gammaksum = zeros(LA,K);
  Scale=zeros(TMAX,1);
  sxi=zeros(K,K);

  for n=1:N
    alpha=zeros(T(n),K);
    beta=zeros(T(n),K);
    gamma=zeros(T(n),K);
    gammaksum = zeros(LA,K);  

    % Inital values of B = Prob(output|s_i), given data X

    Xcurrent=X{n};
    
    for i=1:T(n)
      m = findstr(alphabet,Xcurrent(i));
      if (m == 0)
	fprintf('symbol not found\n');  
	return;
      end
      B(i,:) = E(m,:);
    end;

    scale=zeros(T(n),1);
    alpha(1,:)=Pi(:)'.*B(1,:);
    scale(1)=sum(alpha(1,:));
    alpha(1,:)=alpha(1,:)/scale(1);
    for i=2:T(n)
      alpha(i,:)=(alpha(i-1,:)*P).*B(i,:);
      scale(i)=sum(alpha(i,:));
      alpha(i,:)=alpha(i,:)/scale(i);
    end;

    beta(T(n),:)=ones(1,K)/scale(T(n));
    for i=T(n)-1:-1:1
      beta(i,:)=(beta(i+1,:).*B(i+1,:))*(P')/scale(i); 
    end;

    gamma=(alpha.*beta)+epsi;
    gamma=rdiv(gamma,rsum(gamma));

    gammasum=sum(gamma);

    for i = 1:T(n)
      m = findstr(alphabet,Xcurrent(i));
      gammaksum(m,:) = gammaksum(m,:) + gamma(i,:);
    end;

    for i=1:T(n)-1
      t=P.*( alpha(i,:)' * (beta(i+1,:).*B(i+1,:)));
      sxi=sxi+t/sum(t(:));
    end;

    Gammainit=Gammainit+gamma(1,:);
    Gammasum=Gammasum+gammasum;
    Gammaksum = Gammaksum + gammaksum;

    for i=1:T(n)-1
      Scale(i,:) = Scale(i,:) + log(scale(i,:));
    end;
    Scale(T(n),:) = Scale(T(n),:) + log(scale(T(n),:));

  end;
  
  lik=sum(Scale);
  fprintf('test log likelihood = %f \n',lik);