% viterbi(mu,A,E,X) uses the Viterbi algorithm to find the most likely hidden state sequence associated with the data in (matrix) X

function [Pi] = viterbi(mu,A,E,X)

m = length(mu);                        % number of hidden states
n = length(X(1,:));                    % number of observed values

v = zeros(m,n+1);                      % place to store probabilities
                                       % row is hidden state
                                       % column is observation number
v(:,1) = log(mu');                     % log of initial probabilities
LA = log(A);                           % replace probabilities by logarithms
LE = log(E);                           % replace probabilities by logarithms

for i=2:n,                             % loop through all observations
  [ma,j] = max(v(:,i-1)*ones(1,m) + LA);
  v(:,i) = ma' + sum(LE(:,X(:,i)),2);  % add over all replications
  point(:,i) = j';                     % points back to most likely state
end

[ma,j] = max(v(:,n));                  % find the most likely last h. state

Pi(n) = j;

for i=n:-1:2,
  Pi(i-1) = point(Pi(i),i);            % look back to find most likely h.s.
end

break

% ---------------------------------------------------------------------------
% Same thing without using vectors over hidden states

w = zeros(m,n+1);                      % place to store probabilities
                                       % row is hidden state
                                       % column is observation number
w(:,1) = log(mu');                     % log of initial probabilities

for i=2:n,                             % observation number
  for L=1:m,                           % hidden states
    [ma,j] = max(w(:,i-1) + LA(:,L));  % largest prob of ending in h.s. L
    w(L,i) = ma;
    for k=1:length(X(:,1)),
      w(L,i) = w(L,i) + LE(L,X(k,i));  % prob of emission of observed state
    end
    point(L,i) = j;
  end
end

[ma,j] = max(w(:,n));

Pi2(n) = j;

for i=n:-1:2,
  Pi2(i-1) = point(Pi2(i),i);
end

[max(max(abs(w-v))) max(abs(Pi-Pi2))]

show(Pi,'FL');
show(Pi2,'FL');

Pi = Pi2;