// Hmm.cpp: implementation of the CHmm class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "hmm_1.h"
#include "Hmm.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CHmm::CHmm()
{
	m_pA = NULL;
	m_pB = NULL;
	m_pPI = NULL;

	m_pInput = NULL;
	m_pRltMatrix = NULL;
	m_plamda = NULL;

	m_iDimA = 0;
	m_iDimB = 0;
	m_iInputLen = 0;
}

CHmm::~CHmm()
{
}

/*
	int iDimA		轉換矩陣的維數
	int iDimB		輸出符號數
	double **pA		狀態轉換矩陣iDimA*iDimA
	double **pB		輸出符號的概率 iDimA*iDimB
	double *pPI		初始狀態概率
*/
void CHmm::SetHmmParam(int iDimA, int iDimB, double **pA, double **pB, double *pPI)
{
	m_pA = pA;
	m_iDimA = iDimA;

	m_pB = pB;
	m_iDimB = iDimB;
	
	m_pPI = pPI;
}

/*
 *	向前算法
 *	int *pInput			檢測序列
 *	int iInputLen		序列長度
 *	double **pRltMatrix	中間結果
 */
double CHmm::Forward(int *pInput, int iInputLen, double **pRltMatrix)
{
	for(int i=0; i<iInputLen; i++)
	{
		if(pInput[i]>=m_iDimB)
		{
			AfxMessageBox("輸入的數據超界");
			return 0.0;
		}
	}


	m_pRltMatrix = pRltMatrix;

	double rlt = 0;
	m_pInput = pInput;
	m_iInputLen = iInputLen;

//	NewMatrix();

	double *alfa;
	alfa = new double[m_iDimA];

	for(i=0; i<m_iDimA; i++)
	{
		alfa[i] = forward_recursion(iInputLen-1, i);
		//m_pRltMatrix[i][0] = alfa[i];
		rlt+=alfa[i];
	}



	delete [] alfa;
	return rlt;
}

double CHmm::forward_recursion(int iIndex, int j)
{
	double rlt= 0;
	double *alfa_last;
	alfa_last = new double[m_iDimA];

	if(iIndex==1)
	{

		for(int i=0; i<m_iDimA; i++)
		{
			alfa_last[i] = m_pPI[i]*m_pB[i][m_pInput[0]];
			m_pRltMatrix[i][0] = alfa_last[i];
		}
	}
	else
	{
		for(int i=0; i<m_iDimA; i++)
		{
			alfa_last[i] = forward_recursion(iIndex-1, i);
		}
	}

	for(int i=0; i<m_iDimA; i++)
	{
		rlt+=alfa_last[i]*m_pA[i][j];
	}

	rlt *= m_pB[j][m_pInput[iIndex]];

	m_pRltMatrix[j][iIndex] = rlt;
	//m_pRltMatrix[j][iIndex+1] = rlt;
	delete [] alfa_last;
	return rlt;
}


/*
 *	向后算法
 *	int *pInput			檢測序列
 *	int iInputLen		序列長度
 *	double **pRltMatrix	中間結果
 */
double CHmm::Backward(int *pInput, int iInputLen, double **pRltMatrix)
{
	for(int i=0; i<iInputLen; i++)
	{
		if(pInput[i]>=m_iDimB)
		{
			AfxMessageBox("輸入的數據超界");
			return 0.0;
		}
	}


	m_pRltMatrix = pRltMatrix;

	double rlt = 0;
	m_pInput = pInput;
	m_iInputLen = iInputLen;

//	NewMatrix();

	double *beta;
	beta = new double[m_iDimA];

	for(i=0; i<m_iDimA; i++)
	{
		beta[i] = backward_recursion(0, i)*m_pPI[i]*m_pB[i][m_pInput[0]];
		m_pRltMatrix[i][0] = beta[i];
		rlt+=beta[i];
	}


	return rlt;
}
double CHmm::backward_recursion(int iIndex, int j)
{
	double rlt = 0.0;
	double *beta_last;
	beta_last = new double[m_iDimA];

	if(iIndex==m_iInputLen-1)
	{
		rlt = 1.0;
		//rlt = 2.0;

	}
	else
	{
		for(int i=0; i<m_iDimA; i++)
		{
			beta_last[i] = backward_recursion(iIndex+1, i);

			rlt+=beta_last[i]*m_pA[j][i]*m_pB[i][m_pInput[iIndex+1]];
		}
	}


	//保存中間結果
	m_pRltMatrix[j][iIndex+1] = rlt;


	return rlt;
}


/*
 *	韋特比算法
 *	int *pInput			檢測序列
 *	int iInputLen		序列長度
 *	int *pBestPath		結果最佳序列
 *	double **pRltMatrix	中間結果
 */
double CHmm::Viterbi(int *pInput, int iInputLen, int *pBestPath, double **pRltMatrix, int **plamda)
{

	for(int i=0; i<iInputLen; i++)
	{
		if(pInput[i]>=m_iDimB)
		{
			AfxMessageBox("輸入的數據超界");
			return 0.0;
		}
	}

	m_pRltMatrix = pRltMatrix;
	m_plamda = plamda;


	double rlt = 0.0;

	m_pInput = pInput;
	m_iInputLen = iInputLen;


	double *deta;
	deta = new double[m_iDimA];

	for(i=0; i<m_iDimA; i++)
	{
		deta[i] = viterbi_recursion(iInputLen-1, i);
		if(rlt<deta[i])
		{
			rlt=deta[i];
			pBestPath[iInputLen-1] = i+1;
		}
	}

	
	for(i=iInputLen-2; i>=0; i--)
	{
		pBestPath[i] = m_plamda[pBestPath[i+1]-1][i+1];
	}

	delete [] deta;
	return rlt;
}
//韋特比算法 的遞推函數
double CHmm::viterbi_recursion(int iIndex, int j)
{
	double rlt= 0;
	double *deta_last;
	deta_last = new double[m_iDimA];

	if(iIndex==1)
	{

		for(int i=0; i<m_iDimA; i++)
		{
			deta_last[i] = m_pPI[i]*m_pB[i][m_pInput[0]];

			//保存中間結果
			m_pRltMatrix[i][0] = deta_last[i];
			
		}
	}
	else
	{
		for(int i=0; i<m_iDimA; i++)
		{
			deta_last[i] = viterbi_recursion(iIndex-1, i);
		}
	}

	double temp;
	for(int i=0; i<m_iDimA; i++)
	{
		temp = deta_last[i]*m_pA[i][j];
		if(rlt<temp)
		{
			rlt=temp;

			//保存最佳路徑
			m_plamda[j][iIndex] = i+1;
		}
	}

	rlt *= m_pB[j][m_pInput[iIndex]];

	//保存中間結果
	//m_pRltMatrix[iIndex][j] = rlt;
	m_pRltMatrix[j][iIndex] = rlt;

	delete [] deta_last;
	return rlt;
}