function [stateSeq,logProb] = logvitlh(llhoods,transitions)
% [stateSeq,logProb] = logvitlh(llhoods,transitions)
%    Find Viterbi path through a set of local match likelihoods.
%    llhoods is an array with each row corresponding to a particular state 
%    and each column describing one time frame; the values are the 
%    log-likelihoods of that state given the features at that time (i.e. the 
%    log of the 'b' values in the classic Rabiner formulation).  
%    transitions is a transition matrix.  
%    First state (first row of llhoods, transitions) is assumed 
%    as entry state and last as (obligatory) exit.
%    Find the most likely path by dynamic programming, return it.
%    First and last states are not included in returned path.
%    logProb returns actual cost of best path.
% 2003-03-13 dpwe@ee.columbia.edu based on sacha's EPFL logvit.m

transitions(transitions<1e-100) = 1e-100;
logTrans = log(transitions);

[numStates, numPts] = size(llhoods);
nMinOne = numStates - 1;

% Initialize the delta and psi vectors for the emitting states
for i=2:nMinOne,
  delta(i) = logTrans(1,i) + llhoods(i,1);
  psi(1,i) = 1;
end;
delta = delta(:);

% Recursion
for t = 2:numPts,
  deltaBefore = delta;
  for i = 2:nMinOne,
    [maxDelta,index] = max( deltaBefore(2:nMinOne) + logTrans(2:nMinOne,i) );
%    if index == 0; disp('argh'); end;
    delta(i) = maxDelta + llhoods(i,t);
    psi(t,i) = index + 1;
  end;
end;

% Termination
stateSeq(numPts+2) = numStates; % Exit state
[maxDelta,index] = max( delta(2:nMinOne) + logTrans(2:nMinOne,numStates) );
logProb = maxDelta;
stateSeq(numPts+1) = index + 1;

% Backtracking
for t = numPts:-1:2,
  stateSeq(t) = psi(t,stateSeq(t+1));
end;
stateSeq(1) = 1; % Entry state

% No, cut out first and last points
stateSeq = stateSeq(2:(numPts+1));