chmm.cpp

简单明了的HMM-GMM-KEAMS调用接口

/***************************************************************************

Module Name:

	Continuous Observation Hidden Markov Model with Gaussian Mixture

History: 

	2003/12/13	Fei Wang

***************************************************************************/
#include <stdafx.h>

#include <crtdbg.h>
#include <math.h>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <strstream>
#include "CHMM.h"
using namespace std;


CHMM::CHMM(int stateNum, int dimNum, int mixNum)
{
	m_stateNum = stateNum;
	m_maxIterNum = 100;
	m_endError = 0.001;

	Allocate(stateNum, dimNum, mixNum);

	for (int i = 0; i < m_stateNum; i++)
	{
		// The initial probabilities
		m_stateInit[i] = 1.0 / m_stateNum;

		// The transition probabilities
		for (int j = 0; j <= m_stateNum; j++)
		{
			m_stateTran[i][j] = 1.0 / (m_stateNum + 1);
		}
	}
}
CHMM::~CHMM()
{
	Dispose();
}

void CHMM::Allocate(int state, int dim, int mix)
{
	m_stateModel = new GMM*[state];
	m_stateInit = new double[state];
	m_stateTran = new double*[state];

	for (int i = 0; i < state; i++)
	{
		m_stateModel[i] = new GMM(dim, mix);
		m_stateTran[i] = new double[state + 1]; // Add a final state
	}
}

void CHMM::Dispose()
{
	for (int i = 0; i < m_stateNum; i++)
	{
		delete m_stateModel[i];
		delete[] m_stateTran[i];
	}
	delete[] m_stateModel;
	delete[] m_stateTran;
	delete[] m_stateInit;
}

void CHMM::Zero()
{
	for (int i = 0; i < m_stateNum; i++)
	{
		// The initial probabilities
		m_stateInit[i] = 0;

		// The transition probabilities
		for (int j = 0; j < m_stateNum + 1; j++)
		{
			m_stateTran[i][j] = 0;
		}
	}
}

void CHMM::Norm()
{
	double count = 0;
	int i,j;
	for ( j = 0; j < m_stateNum; j++)
	{
		count += m_stateInit[j];
	}
	for ( j = 0; j < m_stateNum; j++)
	{
		m_stateInit[j] /= count;
	}

	for (i = 0; i < m_stateNum; i++)
	{
		count = 0;
		for ( j = 0; j < m_stateNum; j++)
		{
			count += m_stateTran[i][j];
		}
		if (count > 0)
		{
			for ( j = 0; j < m_stateNum + 1; j++)
			{
				m_stateTran[i][j] /= count;
			}
		}
	}
}

double CHMM::GetStateInit(int i)
{
	_ASSERT(i >= 0 && i < m_stateNum);
	return m_stateInit[i];
}

double CHMM::GetStateFinal(int i)
{
	_ASSERT(i >= 0 && i < m_stateNum);
	return m_stateTran[i][m_stateNum];
}

double CHMM::GetStateTrans(int i, int j)
{
	_ASSERT(i >= 0 && i < m_stateNum && j >= 0 && j < m_stateNum);
	return m_stateTran[i][j];
}

GMM* CHMM::GetStateModel(int i)
{
	_ASSERT(i >= 0 && i < m_stateNum);
	return m_stateModel[i];
}

double CHMM::GetProbability(std::vector<double*>& seq)
{
	vector<int> state;
	return Decode_LHJ(seq, state);
}
double CHMM::GetProbability_LHJ(std::vector<double*>& seq,vector<int>& state)
{
	return Decode_LHJ(seq, state);
}
double CHMM::GetProbability2(std::vector<double*>& seq,vector<int> seq_state)
{
	vector<int> state;
	return Decode_LHJ2(seq, state,seq_state);
}

//对每个观测，先判断分割的结果是否valid,若invalid，则对观测概率赋一个统一的较低的值
double CHMM::Decode_LHJ2(vector<double*>& seq, vector<int>& state,vector<int> seq_state)
{
	// Viterbi
	int size = (int)seq.size();
	int i,j,t;
	double* lastLogP = new double[m_stateNum];
	double* currLogP = new double[m_stateNum];
	int** path = new int*[size];
	

	FILE *fp;
	fp=fopen("d:\\work\\DVCoach\\HMM\\prob.txt","a");


	double *like=new double[size];
	for(i=0;i < size; i++)
		like[i]=-100;

	double FixProb=0.00001;
	double l;
	// Init
	path[0] = new int[m_stateNum];
	for (i = 0; i < m_stateNum; i++)
	{
		//currLogP[i] = LogProb(m_stateInit[i]) + LogProb(m_stateModel[i]->GetProbability(seq[0]));
		if(!seq_state[0])  //invalid
			l=FixProb;
		else
            l=m_stateModel[i]->GetProbability(seq[0]);
		if(like[0]<l)
			like[0]=l;

		fprintf(fp,"%f ",l);
		l=LogProb(l);
		currLogP[i] = LogProb(m_stateInit[i]) + l;
		path[0][i] = -1;
	}
	fprintf(fp,"\n");

	// Recursion
	for ( t = 1; t < size; t++)
	{
		path[t] = new int[m_stateNum];

		double* temp = lastLogP;
		lastLogP = currLogP;
		currLogP = temp;

		for ( i = 0; i < m_stateNum; i++)  
		{
			currLogP[i] = -1e308;

			// Searching the max for last state.//DP, 到 i 的最大概率
			for ( j = 0; j < m_stateNum; j++)
			{
				double l = lastLogP[j] + LogProb(m_stateTran[j][i]);
				if (l > currLogP[i])
				{
					currLogP[i] = l;
					path[t][i] = j;
				}
			}
			//currLogP[i] += LogProb(m_stateModel[i]->GetProbability(seq[t]));
			if(!seq_state[t])
				l=FixProb;
			else
				l=m_stateModel[i]->GetProbability(seq[t]);

			if(like[t]<l)
				like[t]=l;

			fprintf(fp,"%f ",l);
			l=LogProb(l);
			currLogP[i] += l;
		}
		fprintf(fp,"\n");
	}
	FILE *fp2;
	fp2=fopen("d:\\work\\DVCoach\\HMM\\prob2.txt","a");
	for ( i = 0; i < size; i++)
		fprintf(fp2,"%f \n",like[i]);
	delete like;

	fclose(fp);
	fclose(fp2);
	

	// Termination
	int finalState = 0;
	double prob = -1e308;
	for ( i = 0; i < m_stateNum; i++)
	{
		if (currLogP[i] > prob)
		{
			prob = currLogP[i];
			finalState = i;
		}
	}

	// Decode
	state.push_back(finalState);
	for ( t = size - 2; t >=0; t--)
	{
		int stateIndex = path[t+1][state.back()];
		state.push_back(stateIndex);
	}

	// Reverse the state list
	reverse(state.begin(), state.end());

	// Clean up
	delete[] lastLogP;
	delete[] currLogP;
	for ( i = 0; i < size; i++)
	{
		delete[] path[i];
	}
	delete[] path;

	prob = exp(prob / size);
	return prob;
}

double CHMM::Decode_LHJ(vector<double*>& seq, vector<int>& state)
{
	// Viterbi
	int size = (int)seq.size();
	double* lastLogP = new double[m_stateNum];
	double* currLogP = new double[m_stateNum];
	int** path = new int*[size];
	int i,j;

	FILE *fp;
	fp=fopen("d:\\work\\DVCoach\\HMM\\prob.txt","a");


	double *like=new double[size];
	for(i=0;i < size; i++)
		like[i]=-100;

	// Init
	path[0] = new int[m_stateNum];
	for ( i = 0; i < m_stateNum; i++)
	{
		//currLogP[i] = LogProb(m_stateInit[i]) + LogProb(m_stateModel[i]->GetProbability(seq[0]));
		double l=m_stateModel[i]->GetProbability(seq[0]);
		if(like[0]<l)
			like[0]=l;

		fprintf(fp,"%f ",l);
		l=LogProb(l);
		currLogP[i] = LogProb(m_stateInit[i]) + l;
		path[0][i] = -1;
	}
	fprintf(fp,"\n");

	// Recursion
	for (int t = 1; t < size; t++)
	{
		path[t] = new int[m_stateNum];

		double* temp = lastLogP;
		lastLogP = currLogP;
		currLogP = temp;

		for ( i = 0; i < m_stateNum; i++)  
		{
			currLogP[i] = -1e308;

			// Searching the max for last state.//DP, 到 i 的最大概率
			for ( j = 0; j < m_stateNum; j++)
			{
				double l = lastLogP[j] + LogProb(m_stateTran[j][i]);
				if (l > currLogP[i])
				{
					currLogP[i] = l;
					path[t][i] = j;
				}
			}
			//currLogP[i] += LogProb(m_stateModel[i]->GetProbability(seq[t]));
			double l=m_stateModel[i]->GetProbability(seq[t]);
			if(like[t]<l)
				like[t]=l;

			fprintf(fp,"%f ",l);
			l=LogProb(l);
			currLogP[i] += l;
		}
		fprintf(fp,"\n");
	}
	FILE *fp2;
	fp2=fopen("d:\\work\\DVCoach\\HMM\\prob2.txt","a");
	for ( i = 0; i < size; i++)
		fprintf(fp2,"%f \n",like[i]);
	delete like;

	fclose(fp);
	fclose(fp2);
	

	// Termination
	int finalState = 0;
	double prob = -1e308;
	for ( i = 0; i < m_stateNum; i++)
	{
		if (currLogP[i] > prob)
		{
			prob = currLogP[i];
			finalState = i;
		}
	}

	// Decode
	state.push_back(finalState);
	for ( t = size - 2; t >=0; t--)
	{
		int stateIndex = path[t+1][state.back()];
		state.push_back(stateIndex);
	}

	// Reverse the state list
	reverse(state.begin(), state.end());

	// Clean up
	delete[] lastLogP;
	delete[] currLogP;
	for ( i = 0; i < size; i++)
	{
		delete[] path[i];
	}
	delete[] path;

	prob = exp(prob / size);
	return prob;
}
double CHMM::Decode(vector<double*>& seq, vector<int>& state)
{
	// Viterbi
	int size = (int)seq.size();
	double* lastLogP = new double[m_stateNum];
	double* currLogP = new double[m_stateNum];
	int** path = new int*[size];
	int i,j,t;

	// Init
	path[0] = new int[m_stateNum];
	for ( i = 0; i < m_stateNum; i++)
	{
		currLogP[i] = LogProb(m_stateInit[i]) + LogProb(m_stateModel[i]->GetProbability(seq[0]));
		path[0][i] = -1;
	}

	// Recursion
	for ( t = 1; t < size; t++)  //对每一个观测，求属于每个状态的当前最大累加概率
	{
		path[t] = new int[m_stateNum];

		double* temp = lastLogP;
		lastLogP = currLogP;
		currLogP = temp;

		for ( i = 0; i < m_stateNum; i++)
		{
			currLogP[i] = -1e308;

			// Searching the max for last state.
			for ( j = 0; j < m_stateNum; j++)
			{
				double l = lastLogP[j] + LogProb(m_stateTran[j][i]);
				if (l > currLogP[i])
				{
					currLogP[i] = l;
					path[t][i] = j;
				}
			}

			currLogP[i] += LogProb(m_stateModel[i]->GetProbability(seq[t]));
		}
	}

	// Termination
	int finalState = 0;
	double prob = -1e308;
	for ( i = 0; i < m_stateNum; i++)
	{
		if (currLogP[i] > prob)
		{
			prob = currLogP[i];
			finalState = i;
		}
	}

	// Decode
	state.push_back(finalState);
	for ( t = size - 2; t >=0; t--)
	{
		int stateIndex = path[t+1][state.back()];
		state.push_back(stateIndex);
	}

	// Reverse the state list
	reverse(state.begin(), state.end());

	//dump the state 
	FILE *fp=fopen("$state.txt","a");
	for( i=0;i<state.size();i++)
	{
	    int stateIndex=state[i];
		fprintf(fp,"%5d",stateIndex);
	}