double get_cost();	/* user-defined cost function */

vector vit( vector _y, int _M, int _istate)
{
/*
 * generic viterbi routine
 *
 * R. Perry, June 1998
 *
 * Reference:
 *
 *   G. Forney, "The Viterbi Algorithm", Proceedings of the IEEE, Vol. 61,
 *   No. 3, March 1973, pp. 268-278.
 *
 * Arguments:
 *
 *   y - input data samples
 *   M - number of states
 *   istate - optional initial state
 *
 * Uses global variables:
 *
 *   HUGE_VAL - large constant representing infinite cost
 *
 * Returns:
 *
 *   path_index - indices of lowest-cost path
 *
 */

/*
 * leave V and w global in case we want to examine them
 *
 * matrix V;
 * vector w;
 */

vector w_new, path_index;
int i_min;
double w_min, p;

int i, j, k, N;

N = cols(_y);

V = -ones( _M, N+1);		/* viterbi path array and costs */

if( narg > 2 && _istate > 0) {	/* known initial state */
  V( _istate, 1) = _istate;
  w = HUGE_VAL*ones( 1, _M);
  w( _istate) = 0;
} else {			/* unknown initial state */
  V(:,1) = [1:_M]';
  w = zeros( 1, _M);
}


/*** main Viterbi loop ***/

  for( j = 1; j <= N; ++j) {	/* for each column, i.e. time index */
  
    w_new = HUGE_VAL*ones( 1, _M);	/* set big values to start with */
  
    for( i = 1; i <= _M; ++i) {		/* for each row, i.e. state index */
  
      if( w(i) == HUGE_VAL)		/* HUGE_VAL marks unreachable states */
        continue;
  
      for( k = 1; k <= _M; ++k) {	/* for each possible next state */
  
	  p = get_cost( _y(j), i, k);

          if( p == HUGE_VAL)		/* unreachable state */
	    continue;
  
        if( p + w(i) < w_new(k)) {	/* found better path */
          w_new(k) = p + w(i);
          V( k, j+1) = i;
        }
      }
    }
  
    w = w_new;
  }

/*** get optimal path ***/
  
  i_min = find( w == min(w));
  i_min = i_min(1);
  w_min = w(i_min);
  
  path_index = zeros(1,N);

  for( j = N; j >= 1; --j) {
    path_index(j) = i_min;
    i_min = V( i_min, j+1);
  }
      
  return path_index;
}
