/*
**   File：Viterbi.cpp
**   功能：給定HMM和觀察序列，求最可能的狀態
*/

#include "StdAfx.h"
#include <math.h>
#include "vector.h"
//#include "hmm.h"

using namespace std;
//#define VITHUGE  100000000000.0

/**************************************************************************
** 函數名稱：Viterbi
** 功能：Viterbi算法
** 參數：phmm：HMM結構指針
**       T：觀察值的個數
**       O：觀察序列
**       delta，psi為中間變量
**       q：求得的最佳狀態序列
**       pprob：概率
**/
//  A 
void main()
{   
	double **A,**B;
	int 	i, j;	   /* 狀態下標 */
	int   	t;	      /* 時間下標 */	
	int     T,N;      //時間和序列長度
	vector<double>  start;       //hmm中的pi 初始狀態
	int     *q;       //最后記錄的是最佳狀態序列
	vector<int>  object;  //狀態序列
	int  	maxvalind;    // 最佳狀態序號
	double  beststate;    //最后回溯用 記錄最佳狀態序列
	double	maxval, val;
	double  **delta;
	int     **psi;
	/* 1. 初始化  */
	
	for (i = 1; i <= N; i++) 
	{
		delta[1][i] = start[i] * B[i][object[1]];  //delta符號
		psi[1][i] = 0;  
	}	

	/* 2. 遞歸 */
	
	for (t = 2; t <= T; t++) 
	{
		for (j = 1; j <= N; j++) 
		{
			maxval = 0.0;
			maxvalind = 1;	
			for (i = 1; i <= N; i++) 
			{
				val = delta[t-1][i]*A[i][j];
				if (val > maxval) {
					maxval = val;	
					maxvalind = i;	
				}
			}	
			delta[t][j] = maxval*B[j][object[t]];
			psi[t][j] = maxvalind;  //  記錄argmax(delta*A)
		}
	}
	/* 3. 終止 */

    beststate = 0.0;
	q[T] = 1;
	for (i = 1; i <= N; i++)   //找到最佳的序列
	{
		if (delta[T][i] > beststate) 
		{
			beststate = delta[T][i];	
			q[T] = i;
		}
	}  	
	/* 4. Path (state sequence) backtracking */
	for (t = T - 1; t >= 1; t--)
		q[t] = psi[t+1][q[t+1]];
	//q[t]最后存儲的是最佳的觀測序列
}