function  [stateseq,probs,maxprobstateseq] = revhmm (A,B,classprobs)

% REVHMM : Runs Viterbi Algorithm with Posterior Probabilities
%
% We have a sequence of T timesteps, at each timestep is a data point that 
%  has to be classified into one of N classes. For each data point x, you 
%  have used some other classifier and got P(y|x) for each of the N classes y.
%
% A is a N x N matrix of transition probabilities
%   A(i,j) is the probability that a point of class i is followed by one of class j
% B is a N x T matrix of posterior probabilities
%   B(n,t) = probability that t-th point is in class n
% classprobs is a N-element vector
%   classprobs(n) = a priori probability of class n
%   This is uniform if not supplied (and if you reweight your classes, you don't need to supply
%   it)
%
% stateseq is a T x 1 matrix of best labels from Viterbi decoding
% probs is a N x T matrix of probabilities computed along the way

[N,T]=size(B);
if (N == 0)
  error ('N should be positive');
end
  
if nargin<3
  classprobs = repmat(1/N,N,1);
end

if ((N~=size(A,1)) || (N~=size(A,2)))
  error (sprintf('A should be of size %d x %d\n',N,N));
end

probs = zeros(N,T);
probs(:,1) = B(:,1);
for t=2:T
  for j=1:N
    score=0;
    bestsnj=-inf;
    for n=1:N
      snj=A(n,j) * probs(n,t-1);
      score = score+snj;
      if snj>bestsnj
	prevstate(j,t)=n;
	bestsnj=snj;
      end
    end
    probs(j,t)=score*B(j,t) / classprobs(j) ;   % incorrect Bayesianally but works better sometimes
%    probs(j,t)=score*B(j,t);
%    probs(j,t)=score*bestsnj;                  % consider only path from previous best state, not summing over all previous states
  end
  probs(:,t) = probs(:,t) / sum(probs(:,t));
end
[m,mi]=max(probs(:,T));

stateseq = zeros(T,1);

stateseq(T)=mi;
for t=T-1:-1:1
  stateseq(t) = prevstate(stateseq(t+1),t+1);
end

[tmp,maxprobstateseq] = max(probs);